bernatgel
commented
2 years ago Hi Diya,
TL;DR: You should add ymax=50000 to your call to kpPlotBAMDensity. Full explanation below.
I see what your problem is. You have two problems here, actually, a biology/data problem and a karyoploteR one.
The biology/data problem is that you have a veeery high peak in your Pag3770_10 chromosome, about 45X higher than the surrounding regions. If you are using a custom genome (you have assembled or it's an imperfect genome for a non-model organism) and your data is whole-genome DNA-sequencing, your peak most probably arises from an erroneous genome assembly. It could be possible that a bunch of repeats located in different parts of the genome were converted into a single element here in your Pag3770_10 chromosome and you are seeing the sum of the coverage of all those repeats in this single spot. You can sometimes see similar things in the model organisms too near centromeres and telomeres. I don't think there's anything you can do here except make an even better genome draft, if possible, but repetitive regions are a tough nut to crack.
The other problem you have is a karyoploteR problem, a problem with the axis and actually a problem I've seen several times, so it's clear it my fault in designing a and documenting kpAxis and not yours in using it! Anyway, the problem here is that your axis scales are wrong and arbitrary. The thing is that karyoploteR executes the commands one by one, strictly one after the other and independently! So issuing a kpAxis([...], ymax=50000) will not affect in any way kpPlotBAMDensity. If you do not specify ymax in kpPlotBAMDensityit will autoadjust the ymax parameter to the maximum value in the plotted region (31578 for chromosome 4 and 946938 in chromosome 10, in your example). To correctly set your axis according to your plotting, you have two options:
Manually set your axis and data plotting ymax:
kp<-kpPlotBAMDensity(kp, data="../gatk/3770_merged_coordsort.bam", window.size=1e4, ymax=50000, r0=at$r0, r1=at$r1, col="darkorange")
kpAxis(kp, ymax=50000, r0=at$r0, r1=at$r1, cex=0.5)or set your axis depending on your actual data
kp<-kpPlotBAMDensity(kp, data="../gatk/3770_merged_coordsort.bam", window.size=1e4, r0=at$r0, r1=at$r1, col="darkorange")
kpAxis(kp, ymax=max(kp$latest.plot$computed.values$density), r0=at$r0, r1=at$r1, cex=0.5)These two approaches have their own problems. For the first one, you need to know the range of your data and you might end up with data out of range (more on that later). In your example, the peak with a density of 900K will be well over the maximum 50K limit.
The second approach will create correct representations also, but with different limits for each plot, making it more difficult to interpret and sometimes leading to incorrect interpretations as in your question.
In your case, since the high peak seems an error in your data, I would probably go for the first option (setting ymax=50000 in both commands). Since clipping is on by default in zoomed plots, the peak should get cropped at 50000 and you'd get what you need. The only problem, is that you'll not have any indication of the cropping, but you could add one with a kpSegments or something similar, plotting a red line over the regions with a density higher than 50000. I think it would make sense to add this directly to kpPlotBAMDensity (and kpPlotBigWig and kpPlotBAMCoverage and...) but I can not provide a timeframe for it to be done.
So... after this long message, I hope this all makes sense and hope it helps you get the images you need. If you need any further clarification, please ask! And if this is not much of a trouble, could you comment later telling me which route you took? It will help me understand how everybody is using karyoploteR and how to improve the package and its documentation!
Bernat
diyasen2021
Hello,
I am trying to plot read coverage using kpPlotBAMDensity but I find large differences between chromosomes using this method. For instance, I run the code below.
and I obtain the following image for chromosome 4.
where the mean read density is ~10,000 summary(kp$latest.plot$computed.values$density) Min. 1st Qu. Median Mean 3rd Qu. Max. 38 9770 10064 10128 10358 31758
and by changing to chromosome 10, I obtain the following image
where the mean read density is ~18,000 summary(kp$latest.plot$computed.values$density) Min. 1st Qu. Median Mean 3rd Qu. Max. 11412 17548 18337 29892 19102 946938
Although chromosome 4 seems okay, the plot for chromosome 10 is dominated by the single peak near the beginning of the chromosome and this seems to be suppressing the read density along the rest of the chromosome, thereby leading to an inaccurate representation. What could possibly be going wrong here? Any help is much appreciated.
Thanks for your time.
Diya