In jmonlong/PopSV: Population-based detection of structural variants from Read-Depth signal

Constructing control regions for enrichment analysis

Introduction

When investigating enrichment in specific genomic features, it is common to compare the regions of interest with some controls regions. A simple approach to construct control regions could be to randomly select regions across the genome. Additionally, it is important for the control regions to have the same size distribution.

However, a random distribution across the genome is usually not realistic. Likely, regions in the genome were not tested because inaccessible or not included in the analysis. Moreover, you might want to control for some patterns and look for more.

For example, we first observed enrichment of CNVs in low-mappability regions. We then wanted to test additional enrichment in different repeat classes. Because repeats are enriched in low-mappability regions, repeats will likely be seen enriched in CNVs. We want to avoid spurious correlation and control for the low-mappability enrichment. By constructing control regions with the same low-mappability enrichment we can can test additional enrichment in the different repeat classes without being biased by the relation between low-mappability regions and repeats.

Constructing control regions with PopSV package

First we load the package and retrieve some annotations to play with.

library(PopSV)
library(AnnotationHub)
ah = AnnotationHub()
genes = ah[["AH49010"]] ## Genes
dgv = ah[["AH5120"]] ## SVs from DGV
dgv = dgv[sample.int(length(dgv), 1e4)] ## Reduce to 10K random SVs

We imported a gene annotation and 10 thousands SVs from DGV. If we want to construct control regions that fit the SV size and overlap with genes, we run:

dgv.cont = draw.controls(dgv, list(gene=genes), chr.prefix="chr")

Now let's verify that the size distribution is the same. By construction it should be exactly the same.

library(ggplot2)
size.df = rbind(data.frame(reg="dgv", size=width(dgv)),
                 data.frame(reg="control", size=width(dgv.cont)))
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge")
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge") + scale_x_log10()

And that the input and output regions overlap genes similarly. By default an approximation is used and might lead to small differences. To switch off the approximation see further down.

mean(overlapsAny(dgv, genes))
mean(overlapsAny(dgv.cont, genes))

draw.controls functions can accept any number of genomic features to control. Let's import two additional genomic annotation that we would like to control for our enrichment analysis: assembly gaps and segmental duplications.

gap = ah[["AH6444"]]
segdups = ah[["AH5121"]]
dgv.cont2 = draw.controls(dgv, list(gene=genes, gap=gap, sd=segdups), chr.prefix="chr")

Again, the size distribution must be the same:

## Same size distribution ?
size.df = rbind(data.frame(reg="dgv", size=width(dgv)),
                 data.frame(reg="control", size=width(dgv.cont2)))
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge")
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge") + scale_x_log10()

And the overlap with the three different genomic annotations similar.

## Same overlap with features ?
mean(overlapsAny(dgv, genes))
mean(overlapsAny(dgv.cont2, genes))
mean(overlapsAny(dgv, gap))
mean(overlapsAny(dgv.cont2, gap))
mean(overlapsAny(dgv, segdups))
mean(overlapsAny(dgv.cont2, segdups))

If we had used the first set of control regions (only genes overlap control) the gap and segmental duplication overlap proportions wouldn't match.

mean(overlapsAny(dgv.cont, segdups))
mean(overlapsAny(dgv.cont, gap))

Controlling for the distance to a feature

In addition to controlling for the overlap to a set of feature, draw.controls can also control for the distance to one feature. Although we can control for overlap to several feature the distance control is more complex to multiplex and for now we can only control for distance to one feature.

dgv.cont3 = draw.controls(dgv, list(gene=genes, sd=segdups), chr.prefix="chr", dist.gr=gap)

Again, the size distribution must be the same:

## Same size distribution ?
size.df = rbind(data.frame(reg="dgv", size=width(dgv)),
                 data.frame(reg="control", size=width(dgv.cont3)))
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge")
ggplot(size.df, aes(x=size, fill=reg)) + geom_histogram(position="dodge") + scale_x_log10()

And the overlap with the three different genomic annotations similar.

## Same overlap with features ?
mean(overlapsAny(dgv, genes))
mean(overlapsAny(dgv.cont3, genes))
mean(overlapsAny(dgv, segdups))
mean(overlapsAny(dgv.cont3, segdups))

Finally let's check that we now control for the distance to a gap.

## Same distance to gap ?
dist.df = rbind(data.frame(reg="dgv", dist=as.data.frame(distanceToNearest(dgv, gap))$distance),
    data.frame(reg="simple control", dist=as.data.frame(distanceToNearest(dgv.cont2, gap))$distance),
    data.frame(reg="distance control", dist=as.data.frame(distanceToNearest(dgv.cont3, gap))$distance))
ggplot(dist.df, aes(x=dist, fill=reg)) + geom_histogram(position="dodge") + xlab("distance to a gap")
ggplot(dist.df, aes(x=dist, colour=reg)) + stat_ecdf() + ylab("cumulative proportion of regions") + xlab("distance to a gap")

Approximation or exact overlap ?

In order to speed up the process, the input regions are grouped by size. Control regions are then built for each group using the average size. This approximation leads to small differences between feature overlap with input and output regions. It's not a problem usually, especially when the input regions have specific size patterns (e.g. CNVs from a binned genome).

If you want to force the control regions to have exactly the same overlap with the feature, use nb.class=Inf. However it will take more time to compute, especially if the size distribution of the input regions is very diverse. Of note the computation time also increases with additional features to control.

system.time((dgv.cont.approx = draw.controls(dgv, list(gene=genes), chr.prefix="chr")))
system.time((dgv.cont.exact = draw.controls(dgv, list(gene=genes), chr.prefix="chr", nb.class=Inf)))
mean(overlapsAny(dgv, genes))
mean(overlapsAny(dgv.cont.approx, genes))
mean(overlapsAny(dgv.cont.exact, genes))

R session

sessionInfo()

jmonlong/PopSV documentation built on Sept. 15, 2019, 9:29 p.m.