Enrichment between genomic regions

Simulated data

In a very simple scenario of having only one chromosome of size 250 Mbp.

First let’s create a function that draw random regions (ranges) in this chromosome.

library(IRanges)
randRegions <- function(sizes, max.pos = 2.5e+08, max.iter = 10) {
    gr = IRanges(runif(length(sizes), 0, max.pos - 
        sizes), width = sizes)
    dup = which(countOverlaps(gr, gr) > 1)
    iter = 1
    while (iter <= max.iter & length(dup) > 0) {
        gr[dup] = IRanges(runif(length(dup), 0, max.pos - 
            sizes[dup]), width = sizes[dup])
        dup = which(countOverlaps(gr, gr) > 1)
    }
    return(gr)
}

Now some regions will be our “repeats”: 10,000 regions from size 10 bp to 6 Kbp.

rep.r = randRegions(runif(10000, 10, 6000))
sum(width(rep.r))/2.5e+08

## [1] 0.1196666

They cover 11.97% of the chromosome.

Now if we have another set of regions and we want to know how much they overlap with the repeats we could use the hypergeometric test. With this test we assume that we are sampling bases in the genome and testing if it’s covered by a repeat. In that sense, we expect 11.97% of our regions to overlap a repeat. If we compare random regions there shouldn’t be a significant overlap and the distribution of the P-value should be flat.

testHG <- function(feat.r, nb = 1000, size = 1, nb.test = 3000, 
    total.b = 2.5e+08) {
    exp.b = sum(width(feat.r))
    sapply(1:nb.test, function(ii) {
        reg.r = randRegions(rep(size, nb))
        obs.ol = sum(overlapsAny(reg.r, feat.r))
        phyper(obs.ol, exp.b, total.b - exp.b, length(reg.r), 
            lower.tail = FALSE)
    })
}

ht.sim = rbind(data.frame(nb = 1000, size = 1, pv = testHG(rep.r, 
    1000, 1)), data.frame(nb = 1000, size = 1000, pv = testHG(rep.r, 
    1000, 1000)), data.frame(nb = 100, size = 1000, 
    pv = testHG(rep.r, 100, 1000)), data.frame(nb = 1000, 
    size = 100, pv = testHG(rep.r, 1000, 100)))

ht.sim %>% mutate(nbsize = paste0(nb, " x ", size, 
    "bp")) %>% group_by(nbsize) %>% arrange(pv) %>% 
    mutate(cumprop = (1:n())/n()) %>% ggplot(aes(x = pv, 
    y = cumprop, color = nbsize)) + geom_line() + theme_bw() + 
    geom_abline(linetype = 2) + ylab("cumulative proportion") + 
    xlab("P-value") + scale_color_brewer(palette = "Set1", 
    name = "regions") + theme(legend.justification = c(1, 
    0), legend.position = c(0.99, 0.01))

As expected, the hypergeometric test works well for region of 1 bp. Otherwise the distribution of the P-values is biased. The larger the regions the stronger the bias. To a lower extent, more regions also means more bias.

Using control regions with similar features

We want to control for the size distribution and the total number of regions tested. Instead of the hypergeometric test, we can get control regions and compare their overlap with the actual regions, using a logistic regression for example. The control regions must be randomly distributed in the genome but have the same size distribution as our original regions. In the logistic regression we compare the two binary variables: overlapping a repeat or not, being an original region or a control region.

testLR <- function(feat.r, nb = 1000, size = 1, nb.test = 3000) {
    sapply(1:nb.test, function(ii) {
        reg.r = randRegions(rep(size, nb))
        cont.r = randRegions(width(reg.r))
        df = rbind(data.frame(region = TRUE, ol = overlapsAny(reg.r, 
            feat.r)), data.frame(region = FALSE, ol = overlapsAny(cont.r, 
            feat.r)))
        pvs = glm(ol ~ region, data = df, family = binomial()) %>% 
            tidy %>% .$p.value
        pvs[2]
    })
}

lr.sim = rbind(data.frame(nb = 1000, size = 1, pv = testLR(rep.r, 
    1000, 1)), data.frame(nb = 1000, size = 1000, pv = testLR(rep.r, 
    1000, 1000)), data.frame(nb = 100, size = 1000, 
    pv = testLR(rep.r, 100, 1000)), data.frame(nb = 1000, 
    size = 100, pv = testLR(rep.r, 1000, 100)))

lr.sim %>% mutate(nbsize = paste0(nb, " x ", size, 
    "bp")) %>% group_by(nbsize) %>% arrange(pv) %>% 
    mutate(cumprop = (1:n())/n()) %>% ggplot(aes(x = pv, 
    y = cumprop, color = nbsize)) + geom_line() + theme_bw() + 
    geom_abline(linetype = 2) + ylab("cumulative proportion") + 
    xlab("P-value") + scale_color_brewer(palette = "Set1", 
    name = "regions") + theme(legend.justification = c(1, 
    0), legend.position = c(0.99, 0.01))

The distribution of the P-values is much better.

Controlling for correlated features

In the genome, the distribution of genes, repeats, functional regions, and others is not random. Different types of elements tend to be found together while others don’t. For example some repeats are located in GC-rich regions and others in AT-rich regions. Transposable elements don’t overlap exonic regions much. Their are hotspots of segmental duplications.

Sometimes we want to control for the overlap with one (or more) genomic features to test the independent association of another. For example, we known copy number variants (CNVs) are enriched in segmental duplications and transposable elements are also enriched in segmental duplications. We might want to test if CNVs are independently enriched in regions with transposable elements, controlling for the overlap with segmental duplications.

I tried to simulate a first set of regions that significantly overlaps our repeats and another one that significantly overlaps the first set. That way we should see a significant overlap with repeat when we test them separately, but the second one should be significant when we control for the first one.

corRegions <- function(sizes, feat.r, or = 2, max.iter = 10, 
    max.pos = 2.5e+08) {
    reg.r = randRegions(sizes)
    for (ii in 1:or) {
        reg.r = c(reg.r[overlapsAny(reg.r, feat.r)], 
            randRegions(sizes))
    }
    dup = which(countOverlaps(reg.r, reg.r) > 1)
    sizes = width(reg.r)
    iter = 1
    while (iter <= max.iter & length(dup) > 0) {
        reg.r[dup] = IRanges(runif(length(dup), 0, 
            max.pos - sizes[dup]), width = sizes[dup])
        dup = which(countOverlaps(reg.r, reg.r) > 1)
    }
    reg.r
}

## First set of regions
repcor.r = corRegions(rep(10000, 1000), rep.r)
cont.r = randRegions(width(repcor.r))
df = rbind(data.frame(region = TRUE, ol = overlapsAny(repcor.r, 
    rep.r)), data.frame(region = FALSE, ol = overlapsAny(cont.r, 
    rep.r)))
glm(ol ~ region, data = df, family = binomial()) %>% 
    tidy %>% kable

## Second set of regions
seed.r = randRegions(rep(10000, 6000))
seed.ol = overlapsAny(seed.r, repcor.r)
repcorcor.r = c(seed.r[seed.ol], seed.r[head(which(!seed.ol), 
    sum(!seed.ol) * 0.05)])
cont.r = randRegions(width(repcorcor.r))
df = rbind(data.frame(region = TRUE, ol = overlapsAny(repcorcor.r, 
    rep.r)), data.frame(region = FALSE, ol = overlapsAny(cont.r, 
    rep.r)))
glm(ol ~ region, data = df, family = binomial()) %>% 
    tidy %>% kable

df = rbind(data.frame(region = TRUE, ol = overlapsAny(repcorcor.r, 
    rep.r), repcor = overlapsAny(repcorcor.r, repcor.r)), 
    data.frame(region = FALSE, ol = overlapsAny(cont.r, 
        rep.r), repcor = overlapsAny(cont.r, repcor.r)))
glm(ol ~ region + repcor, data = df, family = binomial()) %>% 
    tidy %>% kable

randRegionsCons <- function(reg.r, feat.r, nb.seed = 1e+06) {
    seed.r = randRegions(rep(1, nb.seed))
    dist.df = distanceToNearest(seed.r, feat.r) %>% 
        as.data.frame
    reg.ol = overlapsAny(reg.r, feat.r)
    res.r = lapply(unique(width(reg.r)), function(size) {
        size.ii = which(width(reg.r) == size)
        res.r = IRanges()
        if (sum(reg.ol[size.ii]) > 0) {
            seed.ii = dist.df %>% filter(distance < 
                size/2) %>% .$queryHits %>% sample(sum(reg.ol[size.ii]))
            res.r = c(res.r, resize(seed.r[seed.ii], 
                size, fix = "center"))
        }
        if (sum(!reg.ol[size.ii]) > 0) {
            seed.ii = dist.df %>% filter(distance > 
                size/2) %>% .$queryHits %>% sample(sum(!reg.ol[size.ii]))
            res.r = c(res.r, resize(seed.r[seed.ii], 
                size, fix = "center"))
        }
        res.r
    })
    do.call(c, res.r)
}

contSize.r = randRegionsCons(repcorcor.r, repcor.r)
df = rbind(data.frame(region = TRUE, ol = overlapsAny(repcorcor.r, 
    rep.r), repcor = overlapsAny(repcorcor.r, repcor.r)), 
    data.frame(region = FALSE, ol = overlapsAny(contSize.r, 
        rep.r), repcor = overlapsAny(contSize.r, repcor.r)))
glm(ol ~ region, data = df, family = binomial()) %>% 
    tidy %>% kable

mean(overlapsAny(repcorcor.r, repcor.r))

## [1] 0.7812231

mean(overlapsAny(cont.r, repcor.r))

## [1] 0.1395349

mean(overlapsAny(contSize.r, repcor.r))

## [1] 0.7812231

Comparing different sets

What if we need to compare sets of regions A and B with a third one C. If the A and B are comparable in term of size and total number we could directly compare the overlap or an enrichment estimate (e.g. model estimate). If A and B have different size distribution or just total number of regions, these estimates may not be directly comparable. If they both overlap significantly with C, the previous test (control regions + logistic regression) should test them significant. Even the P-value might be affected by the difference in size/number between the two sets. But how should we compared them ? Which interpretable metric could we use to compare enrichment of two different sets or regions ?

A practical example would be two catalogs of CNVs, say from two different methods, that we want to compare to a functional annotation. If one catalogs has more CNVs, or has larger CNVs, how can we say which one overlaps better with the functional annotation ?

I simulate this scenario and compare a few metrics:

• The fold-change in overlap proportion: proportion overlapping / proportion overlapping in control regions.
• The diff-change in overlap proportion: proportion overlapping - proportion overlapping in control regions.
• The logistic regression estimate which are log odds ratio.

Well, the main question was how should I simulate this. I ended up simulating two sets with similar odds ratio so we already know which metric will work better… One of the value of simulation is to force us to define the question. Or at least think about it. In this example, forcing two different sets to have similar odds ratio seemed more natural than trying to double the proportion for example. The odds ratio seems more fair to me and might avoid the situation where we are more likely to observe a large effect size just because the regions are rarer/smaller.

Using a set of functional regions, I will try to compare a set of small CNVs and large CNVs. We expect more of the large CNVs to overlap the functional regions by chance.

fun.r = randRegions(rep(10, 30000))
cnv.sm = randRegions(rep(1000, 1000))
cnv.lg = randRegions(rep(10000, 1000))
mean(overlapsAny(cnv.sm, fun.r))

## [1] 0.112

mean(overlapsAny(cnv.lg, fun.r))

## [1] 0.688

testLR <- function(reg.r, feat.r) {
    cont.r = randRegions(width(reg.r))
    df.sm = rbind(data.frame(region = TRUE, ol = overlapsAny(reg.r, 
        feat.r)), data.frame(region = FALSE, ol = overlapsAny(cont.r, 
        feat.r)))
    rbind(data.frame(term = "fold-change", estimate = mean(overlapsAny(reg.r, 
        feat.r))/mean(overlapsAny(cont.r, feat.r)), 
        p.value = NA), data.frame(term = "diff-change", 
        estimate = mean(overlapsAny(reg.r, feat.r)) - 
            mean(overlapsAny(cont.r, feat.r)), p.value = NA), 
        glm(ol ~ region, data = df.sm, family = binomial()) %>% 
            tidy %>% select(term, estimate, p.value))
}

metrics.df = rbind(testLR(cnv.sm, fun.r) %>% mutate(region = "cnv.sm"), 
    testLR(cnv.lg, fun.r) %>% mutate(region = "cnv.lg"))
metrics.df %>% filter(term != "(Intercept)") %>% select(region, 
    term, estimate) %>% kable

null.df = lapply(1:1000, function(ii) {
    cnv.sm = randRegions(rep(1000, 1000))
    cnv.lg = randRegions(rep(10000, 1000))
    rbind(testLR(cnv.sm, fun.r) %>% mutate(region = "cnv.sm", 
        rep = ii), testLR(cnv.lg, fun.r) %>% mutate(region = "cnv.lg", 
        rep = ii))
})
null.df = do.call(rbind, null.df)

null.df %>% filter(term == "regionTRUE") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 0) + ggtitle("log odds ratio")

null.df %>% filter(term == "fold-change") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 1) + ggtitle("fold-change")

null.df %>% filter(term == "diff-change") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 0) + ggtitle("diff-change")

asso.df = lapply(1:1000, function(ii) {
    cnv.sm = randRegions(rep(1000, 1000))
    cnv.lg = randRegions(rep(10000, 1000))
    cnv.sm = c(cnv.sm[overlapsAny(cnv.sm, fun.r)], 
        randRegions(rep(1000, 1000)))
    cnv.lg = c(cnv.lg[overlapsAny(cnv.lg, fun.r)], 
        randRegions(rep(10000, 1000)))
    rbind(testLR(cnv.sm, fun.r) %>% mutate(region = "cnv.sm", 
        rep = ii), testLR(cnv.lg, fun.r) %>% mutate(region = "cnv.lg", 
        rep = ii))
})
asso.df = do.call(rbind, asso.df)

asso.df %>% filter(term == "regionTRUE") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 0) + ggtitle("log odds ratio")

asso.df %>% filter(term == "fold-change") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 1) + ggtitle("fold-change")

asso.df %>% filter(term == "diff-change") %>% ggplot(aes(x = estimate, 
    colour = region)) + geom_density() + theme_bw() + 
    geom_vline(xintercept = 0) + ggtitle("diff-change")