Question: R&Bioconductor: Extract coverage from bamfile given GRanges object?
1
gravatar for stianlagstad
3.9 years ago by
stianlagstad1.0k
Oslo, Norway
stianlagstad1.0k wrote:

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.

bioconductor R • 1.4k views
ADD COMMENTlink modified 3.9 years ago • written 3.9 years ago by stianlagstad1.0k
1
gravatar for stianlagstad
3.9 years ago by
stianlagstad1.0k
Oslo, Norway
stianlagstad1.0k wrote:

I think I figured it out:

library(Rsamtools)
library(RNAseqData.HNRNPC.bam.chr14)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

# Get reference to example bamfile
bamfile <-  RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
# Exmple transcript name
transcript <- "uc001yhh.1"

# Get exons for transcript
allExons <- GenomicFeatures::exons(
  TxDb.Hsapiens.UCSC.hg19.knownGene,
  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]
# Unlist
cov <- unlist(cov)
# Turn into numeric vector
cov <- as.numeric(cov)
# Plot
plot(cov)

coverage

ADD COMMENTlink modified 3.9 years ago • written 3.9 years ago by stianlagstad1.0k
1

Careful about what 'coverage' means when exons overlap; maybe you want to ScanBamParam(which=reduce(allExons)) to avoid double counting.

ADD REPLYlink written 3.9 years ago by Martin Morgan1.6k

For my real script I do that in an earlier step, but thanks!

ADD REPLYlink written 3.9 years ago by stianlagstad1.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1691 users visited in the last hour