Hi, I am using pile up from R Package, ‘Rsamtools’. I could not get the correct coverage for the SNP positions I am looking for in my BAM files. The issue is that I do not have the strand information (+ or -) and keep missing some of the reads in my output. Is there a way to resolve this issue?
Here are the parameters I am using for Pile up:
Code:
summary[1:5,]
[1,] "chr1" "115258748" "115258748"
[2,] "chr1" "115258748" "115258748"
[3,] "chr1" "115258747" "115258747"
[4,] "chr1" "115258747" "115258747"
[5,] "chr1" "115258747" "115258747"
data <- summary ### summary could be a subset of a.indel - . a.indel[wanted,core.ann.gr)
data.gr <- GRanges(seqnames =data[,"chr"],ranges = IRanges(start=as.numeric(data[,"start"]),end=as.numeric(data[,"end"])))
which <- data.gr
which
# GRanges object with 103 ranges and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 [115258748, 115258748]
# [2] chr1 [115258748, 115258748]
# [3] chr1 [115258747, 115258747]
# [4] chr1 [115258747, 115258747]
params <-ScanBamParam(which=which,flag=scanBamFlag(isUnmappedQuery=FALSE,isDuplicate=FALSE,isNotPassingQualityControls=FALSE),simpleCigar = FALSE,reverseComplement = FALSE,what=c("qname","flag","rname","seq","strand","pos","qwidth","cigar","qual","mapq") ) ### NOTE isValidVendorRead=FALSE shoudl be TRUE
param.pile <- PileupParam(max_depth=2500, min_base_quality=0, min_mapq=0,min_nucleotide_depth=1, min_minor_allele_depth=0,distinguish_strands=TRUE, distinguish_nucleotides=TRUE,ignore_query_Ns=TRUE, include_deletions=TRUE,cycle_bins=numeric() )
#minimum base quality 20