View Single Post
Old 11-05-2015, 11:42 PM   #1
MAPK
Junior Member
 
Location: Here (at lab!)

Join Date: Nov 2015
Posts: 2
Default SNP coverage using Pile up From RSamtools

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

Last edited by MAPK; 11-05-2015 at 11:45 PM.
MAPK is offline   Reply With Quote