I have RNASeq data and wanted to plot coverage on both positive and negative strands of the species I am studying. I wanted to plot it like this (https://davetang.org/muse/2013/09/07/creating-a-coverage-plot-in-r/), but I was only able to plot for 1 chromosome. My organism has 16 chromosomes and I wanted to plot read coverage for all chromosomes on one plot. It tried several ways in R but since the X axis is really big (1:40 million bases), I am not getting the correct plot. Is there a better option to do this? Thanks!
You can use karyoploteR as Bastien suggested, but if you can a whole genome coverage plot I think the best approach would be to perform a bit of preprocessing instead of going directly to kpPlotCoverage.
I would recommend using something such as deepTools bamCoverage to create a BigWig file with the (possibly stranded) coverage data. And then use karyoploteR's kpPlotBigWig to plot them. If you want to work only in R, you could use kpPlotBAMDensity directly on the BAM files, but it does not (yet) support read strand filtering. These methods are documented in the R documentation (?kpPlotBigWig and ?kpPlotBamDensity) but they have not been yet added to the tutorial.
Important note: kpPlotBigWig does not work on windows machines, so if you are on Windows you should go for the second option.
Take a look at KaryoploteR, especially PlotCoverage
Does anyone know if I can do exactly the same but by using pybedtools or other things?
While not quite the same thing take a look at Qualimap and deepTools as alternates.