R&Bioconductor: Extract coverage from bamfile given GRanges object?
• 2.4k views
Given a .bamfile and a GRanges object with locations of exons in a specific transcript, how can I extract coverage for only the regions in the GRanges object?
My end goal is to plot coverage for exons in a specific transcript.
I think I figured it out:
# Get reference to example bamfile
bamfile <- RNAseqData.HNRNPC.bam.chr14_BAMFILES
# Exmple transcript name
transcript <- "uc001yhh.1"
# Get exons for transcript
allExons <- GenomicFeatures::exons(
vals = list(tx_name = transcript),
columns = list("EXONNAME", "TXNAME", "GENEID"))
# Expand GRanges
allExons <- expand(allExons)
# Remove unwanted rows
allExons <- allExons[mcols(allExons)$TXNAME == transcript]
# Get coverage
cov <- coverage(bamfile, param = Rsamtools::ScanBamParam(which = allExons))
# Coverage only for my transcript
cov <- cov[allExons]
cov <- unlist(cov)
# Turn into numeric vector
cov <- as.numeric(cov)
Login before adding your answer.
Traffic: 1408 users visited in the last hour
Careful about what 'coverage' means when exons overlap; maybe you want to ScanBamParam(which=reduce(allExons)) to avoid double counting.
For my real script I do that in an earlier step, but thanks!