bernatgel / karyoploteR

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

kpPlotRegions with bed file error? #92

Open TrentPrall opened 3 years ago

TrentPrall commented 3 years ago

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:

Screen Shot 2020-11-17 at 4 53 14 PM

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

| [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.

TrentPrall commented 3 years ago

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

bernatgel commented 3 years ago

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

bernatgel commented 3 years ago

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 (bgel@igtp.cat) and I'll take a look.

Bernat

TrentPrall commented 3 years ago

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