5
$\begingroup$

I have genomic .bed file data of 4 different types; type A,B,C,D. These are some examples-

Type A:

1   101380000   101710000   A
1   110085000   110320000   A

Type B:

1   100930000   101335000   B
1   10430000    10560000    B

Type C:

1   101775000   101825000   C
1   103570000   104070000   C

Type D:

1   153325000   153475000   D
1   156355000   156415000   D

This data is for chr1, but I such data for all the chromosomes. I want to draw a chromosome map (possibly using the Gviz package?) which can assign different colors to the 4 types of genomic ranges and make me an idiogram like plot for the chromosomes I wish to see, so I know where my ranges physically are, on the chromosome.

Can this be done using the Gviz package in R? If not, what is the best way and code to do this?

$\endgroup$
  • $\begingroup$ On my (retina) screen some of these ranges would be less than a pixel in width. Are you sure this is the best approach? $\endgroup$ – heathobrien Feb 21 '18 at 11:20
  • $\begingroup$ @heathobrien Aah then I guess perhaps not :/ Can you show me how to do it anyway, then I can think about it. My bed files are ~1000-4000 lines in number so maybe I can see something. $\endgroup$ – rishi Feb 21 '18 at 12:22
  • $\begingroup$ @heathobrien on smaller chromosomes perhaps $\endgroup$ – rishi Feb 21 '18 at 12:28
  • 1
    $\begingroup$ I haven't done that specifically, but it should be easy with ggbio. If you want to see the distribution of your ranges across a chromosome, I suggest plotting the centres of the ranges as points. There are a bunch of examples here $\endgroup$ – heathobrien Feb 21 '18 at 12:48
  • $\begingroup$ You could also use karyoploteR too instead of Gviz, see the vignette in Bioconductor to see the resulting plots $\endgroup$ – llrs Feb 21 '18 at 13:40
10
$\begingroup$

You can use karyoploteR to plot your regions on an ideogram quite easily.

Disclosure: I'm one of the authors of the tool.

For this example I'll start creating 4 GRanges with random regions, but with your data you should load your bed into R and split it into four GRanges.

library(karyoploteR)


typeA <- createRandomRegions(nregions = 100, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 20, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 1000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 100, length.mean = 2e6, length.sd = 5e6)

colors <- c(A="#FFD700", B="#00BFFF", C="#FF7F50", D="#90EE90")

And we can plot them. We have to start creating an empty ideogram with plotKaryotype (we are working with human here, but you can use any genome you want) and then plot the regions in the GRanges with kpPlotRegions. As a bonus, if two regions overlap it will plot them one above the other.

kp <- plotKaryotype(genome = "hg19")

kpPlotRegions(kp, data=typeA, r0=0, r1=0.2, col=colors["A"])
kpPlotRegions(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"])
kpPlotRegions(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"])
kpPlotRegions(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"])
legend(x = "right", legend = names(colors), fill = colors, border=NA)

And you'll get something similar to this

enter image description here

In this example we have fairly large regions and not a lot of them. As @heathobrien pointed out, plotting the actual regions in the whole genome plot might not be the best visualization option. You can either use zoom to create different plots of smaller parts of the genome or change your visualization to show the density of regions along the genome with kpPlotDensity.

Plotting the density and tweaking a bit the visualization (and changing the example to have 10 times more regions) with this codei

typeA <- createRandomRegions(nregions = 1000, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 200, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 10000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 1000, length.mean = 2e6, length.sd = 5e6)

kp <- plotKaryotype(genome = "hg19", plot.type = 4, 
                    ideogram.plotter = NULL, labels.plotter = NULL, main="Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)

kpPlotDensity(kp, data=typeA, r0=0, r1=0.2, col=colors["A"], window.size     = 10e6)
kpAddLabels(kp, labels = "Type A", r0=0, r1=0.2)
kpPlotDensity(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"], window.size = 10e6)
kpAddLabels(kp, labels = "Type B", r0=0.25, r1=0.45)
kpPlotDensity(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"], window.size = 10e6)
kpAddLabels(kp, labels = "Type C", r0=0.5, r1=0.7)
kpPlotDensity(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"], window.size = 10e6)
kpAddLabels(kp, labels = "Type D", r0=0.75, r1=0.95)

it would look something like this

Whole Genome Plot of density of regions using karyoploteR

The plots could be tweaked a lot changing various parameters and you could add other data on top of it if needed. You can find more information on the Bioconductor landing page or at the karyoploteR's Tutorial and Examples page.

$\endgroup$
  • 3
    $\begingroup$ @bernatgel thank you so much for the elaborate answer! I used karyoploteR as the comments suggested, but your code gave me a better answer, which I can now tweak about. Thanks! $\endgroup$ – rishi Mar 1 '18 at 11:46
  • $\begingroup$ I moved the discussion to chat so as not keep this excellent answer clean and not to bother bernatgel with the discussion. $\endgroup$ – terdon Mar 1 '18 at 15:56
  • $\begingroup$ The downloaded source packages are in ‘/tmp/RtmpKhPald/downloaded_packages’ Warning message: package ‘karyoploteR’ is not available (for R version 3.4.3) Can you give me a solution for this? Do I upgrade my R? $\endgroup$ – rishi Mar 6 '18 at 13:58
  • 1
    $\begingroup$ @rishi karyoploteR is available from R3.4 onward, so it should be available for 3.4.3. Are you installing it with biocLite? source("bioconductor.org/biocLite.R") biocLite("karyoploteR") $\endgroup$ – bernatgel Mar 6 '18 at 15:35
  • 1
    $\begingroup$ @rishi Can you send me the exact code you are using to install karyoploteR and the output of sessionInfo()? maybe by email would be best: bgel at igtp.cat and I'll try to understand what's going on $\endgroup$ – bernatgel Mar 7 '18 at 9:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.