bernatgel / karyoploteR

karyoploteR - An R/Bioconductor package to plot arbitrary data along the genome
https://bernatgel.github.io/karyoploter_tutorial/
299 stars 42 forks source link

kpPlotRainfall error when given a path to VCF file that contains other variants than SNVs #153

Open fawaz-dabbaghieh opened 5 months ago

fawaz-dabbaghieh commented 5 months ago

Hi,

Thanks for the great library!

I am trying to generate some rainfall plots from somatic variants over chromosomes. I tested my code on a table generated from somatic SNVs and that works fine. However, I later noticed that kpPlotRainfall also takes a path to VCF as data, I tried that with VCF which contains more complex variants like indels, duplications, and BNDs, and I got the following error:

Error in `[[<-`(`*tmp*`, name, value = c(`NA` = "#888888")) : 
  1 elements in value to replace 0 elements
Calls: kpPlotRainfall -> $<- -> $<- -> [[<- -> [[<-
Execution halted

I am using custom genome and tracks (chm13), and generating the plots for each chromosome on a PDF page. This works fine when the data given to kpPlotRainfall is a GRanges object, but not when giving a VCF

# genome size and bands for chm13 t2t
sample = "sample1"
custom.genome <- toGRanges(genome_size)
custom.cytobands <- toGRanges(genome_bands)
chromosomes = c("chr1")
i <- 1

pdf("test.pdf", height=10, width=16)
# this is supposed to be part of a loop that plots for each chromosome
kp <- plotKaryotype(plot.type=1, genome = custom.genome[seqnames(custom.genome)==chromosomes[i]]
                    ,cytobands = custom.cytobands[seqnames(custom.cytobands)==chromosomes[i]],
                    chromosome = chromosomes[i])
kpAddCytobandLabels(kp, force.all = TRUE, srt = 90, col = 'orange')
title = paste0("[vcf] Somatic Mutations of ", sample, " and chromosome ", chromosomes[i])
kpAddMainTitle(kp, main=title, cex=1.5)
kpAddBaseNumbers(kp)
# only chromosome one VCF
kpPlotRainfall(kp, data = "sniffles_svs_chr1.vcf")
kpAxis(kp, ymax = 8, r0=0.01, r1=0.9)
kpAddLabels(kp, labels = c("Distance between mutations (log10)"), srt=90, pos=1, label.margin = 0.04,
            r0=0.01, r1=0.9)
dev.off()

To debug this, I kept removing different SVs (e.g. DEL, INS, INV, DUP, BND), until it worked with only SNVs left.

Is this the intended usage? Or am I doing something wrong?

Thanks for the help, Fawaz