[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