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:
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