Gviz: Create DataTrack object with GenomicAlignments::coverage() data?
1
0
Entering edit mode
stianlagstad ▴ 90
@stianlagstad-9723
Last seen 4.1 years ago

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.

gviz • 2.2k views
ADD COMMENT
1
Entering edit mode
@florianhahnenovartiscom-3784
Last seen 5.6 years ago
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 COMMENT
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 478 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6