bernatgel
kpPlotRegions with bed file error?
Hi, First of all, amazing package. I'm having trouble creating a regions plot using a GRanges created from a bed file (I'd like something along these lines: https://bernatgel.github.io/karyoploter_tutorial/Tutorial/PlotCoverage/images/Figure2-1.png). I'm using a custom genome to plot against. I know it's loaded properly because i'm able to plot gene names from a .gff file mapped against the custom genome and run kpPlotBAMDensity using BAM files against the genome without issues. Here is the troublesome code:
load genome
mafa6.genome <- toGRanges(data.frame(chr=c("CM021939.1", "CM021940.1", "CM021941.1", "CM021942.1", "CM021943.1", "CM021944.1", "CM021945.1", "CM021946.1", "CM021947.1", "CM021948.1", "CM021949.1", "CM021950.1", "CM021951.1", "CM021952.1", "CM021953.1", "CM021954.1", "CM021955.1", "CM021956.1", "CM021957.1", "CM021958.1", "CM021959.1"), start=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), end=c(223606306, 194592313, 186444865, 171057148, 186553353, 179102756, 171798370, 144116383, 131032084, 96731059, 132457180, 130596009, 108762655, 125104124, 111111006, 79120393, 95081867, 73713002, 58824109, 75859114, 150377965)))
set plot parameters
pp <- getDefaultPlotParams(plot.type=2)
zoom to region
mhc.region <- toGRanges(data.frame("CM021942.1", 136475200, 141859400)) kp <- plotKaryotype(genome=mafa6.genome, plot.type=2, zoom=mhc.region, plot.params = pp)
set background colors
kpDataBackground(kp, data.panel=1) kpDataBackground(kp, data.panel=2)
plot gene intervals
gff.file <- "CM021942.1.gff" features <- import(gff.file) genes <- features[features$type=="gene"] kpPlotRegions(kp, data=genes[strand(genes)=="+"], avoid.overlapping = FALSE, data.panel="ideogram") kpPlotRegions(kp, data=genes[strand(genes)=="-"], avoid.overlapping = FALSE, data.panel="ideogram", col="lightcyan4")
plot read regions
bed <- "24740-MHCfull.filte.mafa6.sorted.bed" regions2 <- toGRanges(bed) kpPlotRegions(kp, data=regions2, data.panel=2, r0=1, r1=0, col="cadetblue")
However when I plot the bed regions they do not stack when overlapping and instead show up as a thin line on the bottom of the data panel:
Here is my genomic ranges for the bed. I am admittedly new to R so perhaps this is an issue?
head(regions2) GRanges object with 6 ranges and 2 metadata columns: seqnames ranges strand | name score <Rle> <IRanges> <Rle> |
[1] CM021939.1 15786-25427 + | 0c81a6d4-1186-49af-8ed6-919cbdbc051b 29 [2] CM021939.1 16320-17649 - | 29a19c7c-d41d-445a-b5f6-c32bcd3319fe 0 [3] CM021939.1 17913-19292 - | 29a19c7c-d41d-445a-b5f6-c32bcd3319fe 0 [4] CM021939.1 18333-19220 - | 74cdb513-34da-4cb6-8ade-8483d4c007ee 0 [5] CM021939.1 19038-19599 - | 1e2277f6-6550-4761-bc64-cdf4959a33a3 0 [6] CM021939.1 21497-22056 - | 1e2277f6-6550-4761-bc64-cdf4959a33a3 0
seqinfo: 544 sequences from an unspecified genome; no seqlengths
Any advice would be greatly appreciated and thank you in advance.
To clarify a bit more, if you look closely it does appear as if the ranges across the chromosome are being read in properly, but they aren't getting stacked properly
Hi TrentPrall,
I'll take a look at that. At a first glance your code seems correct, but I might be missing something. I'll get back to you when I find out what's happening
Bernat
Hi Trent, I tried your code with a bunch of random regions and it seems to work as expected, so I think this might have something to do with the distribution of the regions. I'd try with a wider region plotted (or ideally the whole chromosome), to see if there's a really tall stack of regions somewhere. Another thing you could try is setting num.layers=4 (or any other small number) to force each layer to accupy more vertical space. Or if you are comfortable with sharing this data you can send me the bed file ([email protected]) and I'll take a look.
Bernat
Hi Bernat,
Thank you for the quick response. I have tried your suggestions but have had no luck. I have sent you an email with data.
Thanks again, Trent