limit the coverage bar height

[Rohit-Satyam]

Hi Developers!!

I am trying to limit the yaxis range of coverage because some of the regions have extraordinarily high coverage which make other regions look deprived of any reads. This is giving us a false sense that the coverage of genome is poor. However, when I try using ymax and ymin to limit the bar height I am unable to control the bar height. Isn't there a way to scale the coverage??

Here is the code that I used:

module load libxml2/2.9.6/gnu-6.4.0
module load R/4.1.1/gnu-6.4.0
library(Biostrings)
library(karyoploteR)
p <- readDNAStringSet("PlasmoDB-56_PvivaxP01_Genome.fasta")
## Creating custom Granges###################
t <- data.frame(p@ranges)
t$names <- stringr::str_split(t$names,pattern = " | ", simplify = TRUE)[,1]
t$start <- 1
df <- data.frame(chr=t$names,start=t$start,end=t$width)[1:16,]
gr <- makeGRangesFromDataFrame(df)
#####################################
jpeg(file="saving_plot1.jpeg",width=1800, height=2000,res=200)
kp <- plotKaryotype(genome = gr,cex =0.8)
kp <- kpAddBaseNumbers(kp, tick.dist = 1e6, cex=0.5, units="Mb", tick.len = 3, add.units = "TRUE")
tt <- kpPlotBAMCoverage(kp, data="M.sorted.bam",max.valid.region.size = 3e7,col="red",border="red",lwd=2,ymin=10,ymax=1000,clipping=TRUE)
kpAxis(kp, ymin=10,ymax=1000,cex=0.5)
dev.off()

Kindly help!!

[airichli]

I also met the same problem

[Rohit-Satyam]

@bernatgel Can you please help us here!!

[Rohit-Satyam]

I tried using Bigwig file obtained from deeptools command:

bamCoverage -b file.bam -o file.bw -p 20 --outFileFormat bigwig --normalizeUsing RPKM

And then used kpPlotBigWig as given below. However, the code below assign same Y-axis limit (0-160) to all the chromosomes (technically the maximum coverage per chromosome should be different)

## Reading reference genome
p <- readDNAStringSet("curtisi.fa")
## Creating custom Granges###################
t <- data.frame(p@ranges)
t$names <- stringr::str_split(t$names,pattern = " | ", simplify = TRUE)[,1]
t$start <- 1
# Removing contigs
df <- t[!grepl("contig", t$names),]
df <- data.frame(chr=df$names,start=df$start,end=df$width)

gr <- makeGRangesFromDataFrame(df)

#I just have one bigwig file
big.wig.files <- dir(path = "~/Documents/try",
                    pattern = ".bw",
                    all.files = T,
                    full.names = T)

kp <- plotKaryotype(genome = gr,cex =0.8)
kp <- kpAddBaseNumbers(kp, tick.dist = 1e6, cex=0.5, units="Mb", tick.len = 5, add.units = "TRUE")


out.at <- autotrack(1:length(big.wig.files), 
                    length(big.wig.files), 
                    margin = 0.3, 
                    r0 = 0.3,
                    r1 = 1)
i=1
bigwig.file <- big.wig.files[i]
  
  at <- autotrack(i, length(big.wig.files), r0 = out.at$r0, r1 = out.at$r1, margin = 0.2)
  
  # Plot bigwig
  kp <- kpPlotBigWig(kp, 
                     data = bigwig.file, 
                     ymax = "per.chr", ymin = 0,
                     r0 = at$r0, 
                     col="red",
                     r1 = at$r1)
  computed.ymax <- ceiling(kp$latest.plot$computed.values$ymax)
  
  # Add track axis
  kpAxis(kp, 
         ymin = 0 
         ymax = computed.ymax, 
         numticks = 2,
         r0 = at$r0, 
         r1 = at$r1,
         cex = 0.5)

[bernatgel]

Hi @Rohit-Satyam and @airichli

This is a missing option in karyoploteR, and one that should be included. One way to achieve this should be via clipping, but due to limitations in the way the R processes clippings and the internal structure of karyoploteR, this optimal option is the most difficult. Another option would be to include a new parameter, something along the lines of "max.y.plot" or "clip.data.at.ymax" or something a bit more descriptive, and simply modify the coverage so anything over a certain threshold is changed to the maximum allowed coverage. This would allow you to create the image in the first post without the big towers, and simply adjust the ymax as you wish. Implementing this should be fairly easy, but I cannot do it right now and the next few weeks (I'm away from the lab). But the idea would be to add something like:

bam.cov <- apply(bam.cov, function(x) {
              x[x>ymax] <- ymax
              return(x) 
          })

to line 109 in https://github.com/bernatgel/karyoploteR/blob/master/R/kpPlotBAMCoverage.R

If you can do this and use the modified version (you do not need to rebuild the package, even! simply run the modified function definition after loading karyoploteR) and it should work.

Hope this helps in the meantime!

Bernat