Hello! I'm trying to figure out how gviz produces the coverage plot for AlignmentTracks. I am finding a mismatch between the gviz coverage plot and one produced using the Rsamtools pileup function. When plotted (using the plottracks function for the AlignmentTrack and a simple scatter plot for Rsamtools pileup), the coverage looks the same, but the numerical values are different.

Thank you!

Code for Rsamtools:

which <- GRanges(chr,IRanges(z_start,end + zoom_f))
what <- c("rname","qname", "pos","qwidth")

scanparam <- ScanBamParam(which=which, what=what)

pileupparam <- PileupParam(max_depth=250, min_base_quality=13, min_mapq=0, min_nucleotide_depth=1, min_minor_allele_depth=0, distinguish_strands=FALSE, distinguish_nucleotides=FALSE, ignore_query_Ns=FALSE, include_deletions=TRUE, include_insertions=FALSE)

bamfile <- BamFile(x.bam)
pileup <- pileup(bamfile, scanBamParam = scanparam, PileupParam=pileupparam )
scatter.smooth(pileup$count)

Code for gviz:

gtrack <- GenomeAxisTrack()
alTracks <- AlignmentsTrack(x.bam, name = 'x' , isPaired = paired, col = "#808080")
plotTracks(c(list(gtrack,unlist(alTracks)), 
               from = 38053189, to = 38069893 ,chromosome = 'chr14', main = 'ESR1', cex.main = 4, transcriptAnnotation = "symbol", 
               type= "coverage", fill.coverage = "dodgerblue", cex = 1, sizes = c(1,1,1,rep(3,length(alTracks)))
    )