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.
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))) )