Genomic density and Rainfall plot

Rainfall plots are a powerful tool to visualize the distribution and clustering of genomic regions across the genome. Each point in a rainfall plot represents a genomic region, with:

The x-axis showing its genomic coordinate.
The y-axis showing the log-transformed minimal distance to its two neighboring regions (log10(dist)).

Clusters of nearby regions appear as dense “rainfalls” of points, making this approach ideal for identifying mutation hotspots, differentially methylated regions (DMRs), or other non-random patterns of localization.

However, when region density is high, overplotting can obscure meaningful interpretation. To address this, we overlay a genomic density track that summarizes the fraction of a genomic window covered by regions—providing a smoother representation of local enrichment.

Code

library(ggalign)
#> Loading required package: ggplot2
#> 
#> Attaching package: 'ggalign'
#> The following object is masked from 'package:ggplot2':
#> 
#>     element_polygon

The data comes from the circlize package and includes both hyper- and hypo-methylated DMRs in the human genome.

Code

load(system.file(package = "circlize", "extdata", "DMR.RData", mustWork = TRUE))
bed_data <- dplyr::bind_rows(hyper = DMR_hyper, hypo = DMR_hypo, .id = "DMR") |>
    dplyr::relocate(DMR, .after = dplyr::last_col())

if the amount of regions in a cluster is high, dots will overlap, and direct assessment of the number and density of regions in the cluster will be impossible. To overcome this limitation, additional tracks are added which visualize the genomic density of regions (defined as the fraction of a genomic window that is covered by genomic regions).

circle_genomic(
    "hg19",
    radial = coord_radial(inner.radius = 0.2, rotate.angle = TRUE),
    direction = "inward",
    theme = theme(
        plot.margin = margin(t = 10, b = 10, l = 2, unit = "mm")
    )
) -
    # Remove radial axes for a cleaner circular look
    scheme_theme(
        axis.text.r = element_blank(),
        axis.ticks.r = element_blank()
    ) +

    # Cytoband (Ideogram) Layer ------------------
    ggalign(size = 0.1) +
    geom_rect(
        aes(
            xmin = start, xmax = end, ymin = 0, ymax = 1,
            fill = gieStain
        )
    ) +
    scale_fill_manual(
        values = c(
            gpos100 = rgb(0, 0, 0, maxColorValue = 255),
            gpos = rgb(0, 0, 0, maxColorValue = 255),
            gpos75 = rgb(130, 130, 130, maxColorValue = 255),
            gpos66 = rgb(160, 160, 160, maxColorValue = 255),
            gpos50 = rgb(200, 200, 200, maxColorValue = 255),
            gpos33 = rgb(210, 210, 210, maxColorValue = 255),
            gpos25 = rgb(200, 200, 200, maxColorValue = 255),
            gvar = rgb(220, 220, 220, maxColorValue = 255),
            gneg = rgb(255, 255, 255, maxColorValue = 255),
            acen = rgb(217, 47, 39, maxColorValue = 255),
            stalk = rgb(100, 127, 164, maxColorValue = 255)
        ),
        guide = "none"
    ) +
    scale_x_continuous(
        breaks = scales::breaks_pretty(2),
        labels = scales::label_bytes()
    ) +
    guides(
        r = "none", r.sec = "axis",
        theta = guide_axis_theta(angle = 0)
    ) +
    theme(axis.text.theta = element_text(size = 6)) +
    scale_y_continuous(limits = c(0, 1), oob = scales::oob_keep) +
    # Add chromosome labels
    geom_text(aes(middle, label = seqnames, y = 2.2),
        vjust = 0,
        data = function(d) {
            data <- ggalign_attr(d, "ranges")
            data$middle <- (data$start + data$end) / 2
            data
        }
    ) +

    # add rainfall plot -----------------------
    ggalign(bed_data) +
    geom_point(
        aes(middle, log10(dist), color = DMR),
        data = function(d) {
            d <- dplyr::bind_rows(!!!lapply(
                split(d, ~DMR), function(dd) genomic_dist(dd)
            ))
            d$middle <- (d$start + d$end) / 2
            dplyr::rename(d, seqnames = chr)
        }
    ) +
    scale_color_brewer(
        palette = "Dark2",
        guide = guide_legend(
            position = "inside",
            theme = theme(legend.position.inside = c(0.5, 0.5))
        )
    ) +

    # add density plot -----------------------
    ggalign(DMR_hyper, size = 0.5) +
    geom_density(
        aes(middle, density, color = after_scale(fill)),
        fill = RColorBrewer::brewer.pal(3, "Dark2")[1L],
        stat = "identity",
        data = function(d) {
            d <- genomic_density(d)
            d$middle <- (d$start + d$end) / 2
            d
        }
    ) +
    ggalign(DMR_hypo, size = 0.5) +
    geom_density(
        aes(middle, density, color = after_scale(fill)),
        fill = RColorBrewer::brewer.pal(3, "Dark2")[2L],
        stat = "identity",
        data = function(d) {
            d <- genomic_density(d)
            d$middle <- (d$start + d$end) / 2
            d
        }
    ) &
    facet_sector(vars(seqnames), drop = FALSE) &
    theme(panel.background = element_rect(color = "black", fill = NA))