[BioC] summarizeOverlaps with mapq filter?

Sang Chul Choi scchoi at alaska.edu
Thu Jul 4 03:24:14 CEST 2013


Hi,

I wonder if I could use "summarizeOverlaps" to count reads with mapping quality (or MAPQ) larger than some value, say, 30.  I took the example script in the function's help to specify what function signature I want to use.   I used to use "readBamGappedAlignments" to read a BAM file, and filter out reads with MAPQ less than a specified value, and use the filtered GappedAlignments object to count reads mapped on genomic features.  But, I like the feature of summarizeOverlaps that I can use to read a list of BAM files.  But, I do not know use the summarizeOverlaps function to use only reads with MAPQ larger than some specified value.

If this is not possible currently, then I might have to create another BAM file of reads with MAPQ larger than some value and use the BAM file to count overlaps using summarizeOverlaps function.  I do not like this two-step procedure.

Thank you for your help!

SangChul

---------------------------------
   library(Rsamtools)

     fls <- list.files(system.file("extdata",package="GenomicRanges"),
                       recursive=TRUE, pattern="*bam$", full=TRUE)
     names(fls) <- basename(fls)
     bf <- BamFileList(fls, index=character())
     features <- GRanges(
         seqnames = c(rep("chr2L", 4), rep("chr2R", 5), rep("chr3L", 2)),
         ranges = IRanges(c(1000, 3000, 4000, 7000, 2000, 3000, 3600,
                            4000, 7500, 5000, 5400),
                          width=c(rep(500, 3), 600, 900, 500, 300, 900,
                                  300, 500, 500)), "-",
         group_id=c(rep("A", 4), rep("B", 5), rep("C", 2)))

     se <- summarizeOverlaps(features, bf)
---------------------------------

===============================
Sang Chul Choi Ph.D.
Research Associate

Sang Chul Choi (222 WRRB IAB)
902 Koyukuk Dr (PO Box 757000)
311 Irving I, University of Alaska Fairbanks
Fairbanks, AK 99775-7000

Office +1-907-474-5071
Fax +1-907-474-6967
Google Phone +1-607-542-9362
===============================



More information about the Bioconductor mailing list