Reflecting coverage differences when ploting several AlignmentsTracks using plotTrack

[JesseRop]

Dear developer,

Many thanks for the great package. It produces beautiful plots! I'm plotting several alignment tracks and would be grateful for advice on how to make the y-axis reflect the true depth and not be standardized to similar heights. As seen in the example below, the category with the lowest coverage has 300 reads while that with the highest has 2000 and yet they have the same height as the y-axis has been blown up for the lower coverage bam.

Below is my attempt at reproducing the issue using datasets provided in the package but I can't find 2 different bam files with coverage in the same region.

library(Gviz)
library(BSgenome.Hsapiens.UCSC.hg19)
sTrack <- SequenceTrack(Hsapiens)

afrom <- 2960000
ato <- 3160000
alTrack <- AlignmentsTrack(
  system.file(package = "Gviz", "extdata", "snps.bam"),
  isPaired = TRUE,
  stacking = "squish",
  name = 'Al1',
  cex.title=1.22,
  cex.axis = 1.2,
  fontcolor.title = 'black',
  background.title="white",
  col.axis = 'black',
  col.border.title="lightgray")

alTrack2 <- AlignmentsTrack(
  system.file(package = "Gviz", "extdata", "indels.bam"),
  isPaired = TRUE,
  stacking = "squish",
  name = 'Al2',
  cex.title=1.22,
  cex.axis = 1.2,
  fontcolor.title = 'black',
  background.title="white",
  col.axis = 'black',
  col.border.title="lightgray")

plotTracks(c(alTrack, alTrack2, sTrack), chromosome = "chr21", from = 44946590, to = 44946660, cex = 0.5, min.height = 8,
           sizes=c(5, 5,2),
           # sizes=c(rep(5, length(altracks)), c(1,3)),
           type = "coverage")

Kindly advice on how to scale the data tracks so that they all have the same y-axis limits to reflect the true coverage differences. Thanks!

[ivanek]

Hi @JesseRop, I am not sure if I got your question. Would setting of ylim=c(0,2700) in the plotTracks call do the trick for you? Best, Robert

[JesseRop]

Hi Robert, thank you very much for your prompt reply. Adding ylim=c(0,2700) in the plotTracks call should solve it for this particular example. However, I'm making hundreds of plots across the genome and was looking for a way to automate this without knowing the coverage depth beforehand.

[ivanek]

Hi @JesseRop,

Not in an easy way. Gviz reads data for each track separately and plots them directly.

You would have to read-in data from all files and create a list of Track objects. In that list you can iterate over all tracks to get the coverages and calculate the optimal ylim coordinates. (You can use Gviz:::.import.bam.alignments to import the data.)