On 06/07/2012 11:34 PM, Yu Chuan Tai wrote: > Hi Martin, > > Is it correct that the below code calculate the number of reads overlapping > with an amplicon, and overlapping means at least 1 base overlap, and it > doesn't have to be fully within the amplicon? In the case that a read > overlaps > with 2 amplicons, will it be counted twice? When I used this approach to > calculate > amplicon-level read counts, I found the number of read counts > overlapping with > all the amplicons is larger than the total number of read counts in the > BAM file, > and wonder if that's b/c a read could be counted more than once? yes > I found that the code below gives more read counts than using > summarizeOverlaps(). > I think it's b/c the latter counts a read at most once. see ?summarizeOverlaps for how counting occurs with different modes. > > If I want to calculate the coverage of SNVs/INDELs outputted from samtools, > is it correct that using summarizeOverlaps() will under-estimate the > coverage, since > a read may overlap with several SNVs/INDELs? It is 'easy' using findOverlaps() or countOverlaps() to create counting schemes that are different from those implemented in samtools / scanBam, summarizeOverlaps, etc. Martin > > Thanks! > Yu Chuan > > > param = ScanBamParam(what="seq", which=GRanges("seq1",IRanges(100, > 500))) > > dna = scanBam(fl, param=param)[[1]][["seq"]] > > length(dna) # 365 reads overlap region > > alphabetFrequency(dna, collapse=TRUE, baseOnly=TRUE) # 2838 + 3003 GC > > > > On Thu, 7 Jun 2012, Martin Morgan wrote: > >> On 06/07/2012 07:54 AM, Yu Chuan Tai wrote: >>> Hi Martin, >>> >>> Thanks! I will look into the links below. By 'better support for >>> paired-end reads in the 'devel' version', which package are you >>> referring to? >> >> Mostly GenomicRanges, e.g., readGappedAlignmentPairs, building on >> additional facilities in Rsamtools. Herve is responsible for this. >> >> Martin >> >>> >>> Best, >>> Yu Chuan >>> >>> On Thu, 7 Jun 2012, Martin Morgan wrote: >>> >>>> On 06/06/2012 09:53 PM, Yu Chuan Tai wrote: >>>>> Hi Martin, >>>>> >>>>> More questions on your approaches below. If my BAM files are >>>>> generated by Bowtie2 on pair-end fastq files, for scanBamFlag(), >>>>> should I set isPaired=TRUE? Do I need to worry about other input >>>>> arguments for scanBamFlag() or ScanBamParam(), if I want to >>>>> calculate coverage properly? >>>> >>>> It really depends on what you're interested in doing; see for instance >>>> the post by Herve the other day >>>> >>>> https://stat.ethz.ch/pipermail/bioconductor/2012-June/046052.html >>>> >>>>> >>>>> Also, summarizeOverlaps() doesn't seem to handle paired-end reads. >>>>> How to get around this, or it won't affect coverage calculation? >>>> >>>> There is better support for paired-end reads in the 'devel' version of >>>> Biocondcutor; see >>>> >>>> http://bioconductor.org/developers/useDevel/ >>>> >>>> whether and what aspects of paired-endedness are important depends on >>>> how you are using your coverage. >>>> >>>>> >>>>> Finally, is there any way to calculate base-specific coverage at any >>>>> genomic locus or interval in Rsamtools? Thanks! >>>> >>>> I tried to answer this in your other post. >>>> >>>> Martin >>>> >>>>> >>>>> Best, Yu Chuan >>>>> >>>>>> More specifically, after >>>>>> >>>>>> library(Rsamtools) example(scanBam) # defines 'fl', a path to a >>>>>> bam file >>>>>> >>>>>> for a _single_ genomic range >>>>>> >>>>>> param = ScanBamParam(what="seq", which=GRanges("seq1", >>>>>> IRanges(100, 500))) dna = scanBam(fl, param=param)[[1]][["seq"]] >>>>>> length(dna) # 365 reads overlap region alphabetFrequency(dna, >>>>>> collapse=TRUE, baseOnly=TRUE) # 2838 + 3003 GC >>>>>> >>>>>> though you'd likely want to specify several regions (vector >>>>>> arguments to GRanges) and think about flags (scanBamFlag() and the >>>>>> flag argument to ScanBamParam), read mapping quality, reads >>>>>> overlapping more than one region, etc. (summarizeOverlaps >>>>>> implements several counting strategies, but it is 'easy' to >>>>>> implement arbitrary approaches). >>>>>> >>>>>>> >>>>>>> Martin >>>>>>> >>>>>>>> >>>>>>>> Thanks for any input! >>>>>>>> >>>>>>>> Best, Yu Chuan >>>>>>>> >>>>>>>> _______________________________________________ Bioconductor >>>>>>>> mailing list Bioconductor at r-project.org >>>>>>>> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the >>>>>>>> archives: >>>>>>>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>>>>>> >>>>>>> >>>>>> >>>>>> >>>>>> >>>>>>>> >>>>>>>> >>>> -- >>>>>> Computational Biology Fred Hutchinson Cancer Research Center 1100 >>>>>> Fairview Ave. N. PO Box 19024 Seattle, WA 98109 >>>>>> >>>>>> Location: M1-B861 Telephone: 206 667-2793 >>>>>> >>>> >>>> >>>> -- >>>> Computational Biology / Fred Hutchinson Cancer Research Center >>>> 1100 Fairview Ave. N. >>>> PO Box 19024 Seattle, WA 98109 >>>> >>>> Location: Arnold Building M1 B861 >>>> Phone: (206) 667-2793 >>>> >> >> >> -- >> Computational Biology / Fred Hutchinson Cancer Research Center >> 1100 Fairview Ave. N. >> PO Box 19024 Seattle, WA 98109 >> >> Location: Arnold Building M1 B861 >> Phone: (206) 667-2793 >> -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793