Announcement

Collapse

Welcome to the New Seqanswers!

Welcome to the new Seqanswers! We'd love your feedback, please post any you have to this topic: New Seqanswers Feedback.
See more
See less

SNP coverage using Pile up From RSamtools

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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, 11:45 PM.
Working...
X