Search
Question: Gviz: Create DataTrack object with GenomicAlignments::coverage() data?
0
gravatar for stianlagstad
20 months ago by
stianlagstad90
stianlagstad90 wrote:

I want to plot a transcript with it's coverage data as a histogram using DataTrack, but I'm a bit stuck.. Consider this code:

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]
# Set $transcript mcols
mcols(allExons)$transcript <- mcols(allExons)$TXNAME
# Create track
trTrack <- Gviz::GeneRegionTrack(allExons)
# Plot
Gviz::plotTracks(trTrack)
# Transcript plot works well

# Get coverage
cov <- GenomicAlignments::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)

# How do I create a DataTrack object of with these coverage values? I wish to plot the transcript with the coverage data I now have in cov.

Any ideas? I don't want to use an AlignmentsTrack. For this task I have to extract coverage data using coverage() with the bamfile.

ADD COMMENTlink modified 20 months ago by florian.hahne@novartis.com1.5k • written 20 months ago by stianlagstad90
1
gravatar for florian.hahne@novartis.com
20 months ago by
Switzerland

Please check 

? DataTrack

for details on how to initialise objects of the class. I suspect that cov, before casting to numeric,  is an RLE object, so something along these lines should work:

dt <- DataTrack(start=start(cov), end=end(cov), data=runValue(cov))

Florian

 
ADD COMMENTlink written 20 months ago by florian.hahne@novartis.com1.5k

Thank you, that brings me closer!

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]
# Set $transcript mcols
mcols(allExons)$transcript <- mcols(allExons)$TXNAME
# Create track
trTrack <- Gviz::GeneRegionTrack(
  allExons,
  chromosome = "chr14")
# Plot
Gviz::plotTracks(trTrack)
# Transcript plot works well

# Get coverage
cov <- GenomicAlignments::coverage(bamfile, param = Rsamtools::ScanBamParam(which = allExons))
# Coverage only for my transcript
cov <- cov[allExons]
# Unlist
cov <- unlist(cov)

# Create data track
dt <- DataTrack(
  start = start(cov),
  end = end(cov),
  data = runValue(cov),
  genome = "hg19",
  type = "h",
  chromosome = "chr14")
# Plot
Gviz::plotTracks(dt)
# Works!

# Try plotting both
Gviz::plotTracks(
  list(trTrack, dt),
  chromosome = "chr14",
  from = trTrack@start,
  to = trTrack@end)
# Coverage not showing

I'm guessing this is because I unlisted the cov earlier, so that it does't contain the coordinates I need in order to plot the data at the right place. Does anyone have any pointers?

ADD REPLYlink written 20 months ago by stianlagstad90
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.2.0
Traffic: 174 users visited in the last hour