This is a test version of Biostars. For the public version, visit https://www.biostars.org.
fold enrichment between genomic regions/loci

Hi there,

I have the following problem. I called variants in two genomic regions (genome-wide and centromere) following the same approach; now, the two have on average three order of magnitude of size difference and, therefore, I wish to understand whether any of them has an enrichment for overall variants identified.

To do so I started with a very preliminary table — see below:

|   id  |  chr_len  | size_no-cent | size_cent | variants_gw | genome-wide | variants_ce | centromere |
|:-----:|:---------:|:------------:|:---------:|:-----------:|:-----------:|:-----------:|:----------:|
|  chr1 | 245130294 |    239630887 |   5499407 |      179329 |        1336 |        9017 |        610 |
|  chr2 | 242677997 |    240267333 |   2410664 |      190309 |        1263 |        1632 |       1477 |
|  chr3 | 198760610 |    192905570 |   5855040 |      161069 |        1198 |        2568 |       2280 |
|  chr4 | 192622708 |    188490246 |   4132462 |      161336 |        1168 |         145 |      28500 |
|  chr5 | 182259134 |    177862654 |   4396480 |      148078 |        1201 |        5360 |        820 |
|  chr6 | 172045469 |    167796642 |   4248827 |      152031 |        1104 |        2592 |       1639 |
|  chr7 | 163647842 |    154845695 |   8802147 |      135038 |        1147 |        5582 |       1577 |
|  chr8 | 145680301 |    142875109 |   2805192 |      120813 |        1183 |         964 |       2910 |
|  chr9 | 137521496 |    128036976 |   9484520 |      102612 |        1248 |        6120 |       1550 |
| chr10 | 134988382 |    130691620 |   4296762 |       60229 |        2170 |        5415 |        793 |
| chr11 | 134265560 |    128811401 |   5454159 |      110559 |        1165 |        5910 |        923 |
| chr12 | 133236650 |    130160056 |   3076594 |      108783 |        1197 |        2672 |       1151 |
| chr13 | 101016140 |     96437082 |   4579058 |       84539 |        1141 |        2099 |       2182 |
| chr14 |  97561019 |     90289179 |   7271840 |       73720 |        1225 |        6723 |       1082 |
| chr15 |  95272790 |     91581871 |   3690919 |       68285 |        1341 |        2603 |       1418 |
| chr16 |  89874169 |     86069415 |   3804754 |       71579 |        1202 |         660 |       5765 |
| chr17 |  82991830 |     79862750 |   3129080 |       70596 |        1131 |        2198 |       1424 |
| chr18 |  76653175 |     74405679 |   2247496 |       62193 |        1196 |        1658 |       1356 |
| chr19 |  60904190 |     56011623 |   4892567 |       57681 |         971 |        2592 |       1888 |
| chr20 |  65674320 |     60115194 |   5559126 |       51930 |        1158 |        2716 |       2047 |
| chr21 |  38797553 |     34439617 |   4357936 |       34352 |        1003 |        2170 |       2008 |
| chr22 |  45250002 |     38279301 |   6970701 |       34998 |        1094 |        5273 |       1322 |
|  chrX | 227216509 |    224241275 |   2975234 |       82401 |        2721 |        1058 |       2812 |

where I do have the chromosome number (id), its length, its size (excluding centromere), its centromere size, the total number of variants identified genome-wide, the ratio between the length and the former (genome_wide), the total number of variants identified at the centromere, and the ratio between the centromere size and the former (centromere).

I really like the idea of plotting this in an enrichment analysis fashion (like a volcano plot); however, I noticed that the input for such plots requires a p-value, a p-adjusted value as well as any form of fold change (I have seen log2 and log10 so far) for all entries — chromosomes in my case.
The best I can do is to compute a single p-value comparing variants genome-wide and at the centromere by doing a parametric two-tail t-test between groups with unequal variance (which resulted in a non-significant difference between the two: 0,17748). I have seen a common tool to compute such values for genes is DESeq2, so is there a way to process my input to get such information and/or any other alternative down this line if relevant?

Any advice/suggestion is much appreciated; thanks in advance!

statistical-test fold-enrichment variant-calling
• 182 views • link
updated by SeqBench • SeqBench • written by Matteo Ungaro • Matteo Ungaro

1 answer

You don't need DESeq2 for this, and the t-test is the wrong shape for the data.

The t-test treats your 23 per-chromosome ratios as 23 exchangeable observations and throws away the counts behind them, so a centromere with 9017 variants and one with 145 carry equal weight. What you actually have is counts observed over intervals of different length, which is a rate comparison. A Poisson rate-ratio test per chromosome gives you exactly the two numbers you want:

res <- t(apply(df, 1, function(r) {
  pt <- poisson.test(c(r[["variants_ce"]], r[["variants_gw"]]),
                     c(r[["size_cent"]],   r[["size_no_cent"]]))
  c(rate_ratio = unname(pt$estimate), p = pt$p.value)
}))
padj <- p.adjust(res[, "p"], method = "BH")

rate_ratio is your fold enrichment (variants per bp in centromere / variants per bp elsewhere), so the volcano is just log2(rate_ratio) on x and -log10(padj) on y. With counts this large almost everything will be "significant", so read the effect size, not the p-value.

If you want to account for the fact that variant density genuinely varies between chromosomes beyond Poisson noise, fit a quasi-Poisson or negative binomial GLM with offset(log(size)) and region as the predictor, which handles the overdispersion without pretending each chromosome is an independent replicate of the same rate.

One caveat that matters more than the statistics: centromeres are megabase-scale satellite repeat, and short-read mappability there is poor. Whatever difference you find is at least partly technical -- inflated by mismapping in some regions, deflated by dropout in others -- so it's worth reporting callable/covered bases per region instead of raw interval size, and using that as the denominator.

Related, your chr4 row looks off: 145 variants over 4.13 Mb is ~1 per 28.5 kb, while every other centromere is in the 800-2900 bp range. That's a 10-30x outlier in the opposite direction from the rest. I'd check coverage and whether that region was masked or filtered before reading anything biological into the table.

• 0 views • link
written by SeqBench • SeqBench

@Leo thank you so much for the explanation, I see how that makes sense!

About this

One caveat that matters more than the statistics: centromeres are megabase-scale satellite repeat, and short-read mappability there is poor. Whatever difference you find is at least partly technical -- inflated by mismapping in some regions, deflated by dropout in others -- so it's worth reporting callable/covered bases per region instead of raw interval size, and using that as the denominator.

I'm aware of it, and I'm using HiFi long reads — not that it fully accounts for the problem but at least it can be mitigated; additionally, I'm not quite aligning to centromeres from a single FASTA genome but rather to centromeres' pangenomes containing both haplotypes to consider haplotype-aware variation for this sample.
In graph space reads will be aligned to the more similar centromeric sequence between the two haplotypes which, should be the most accurate way determine from which haplotype and how many variants it bears.

I hear you and, still, myself consider this highly experimental but I wanted to give it a try since, otherwise, any linear alternative would generates not only artefacts but also a huge amount of false positives/negatives coming from the mismapping and dropouts.

About this

Related, your chr4 row looks off: 145 variants over 4.13 Mb is ~1 per 28.5 kb, while every other centromere is in the 800-2900 bp range. That's a 10-30x outlier in the opposite direction from the rest. I'd check coverage and whether that region was masked or filtered before reading anything biological into the table.

The region was not masked nor filtered, I also noticed for chr4 that huge discrepancy but that is what the approach returns... the coverage is fine — actually on the high end when compared to other chromosomes after extracting centromeric-specific reads.
Off the top of my head, my first intuition is that those reads don't align quite well to the centromere graph of chr4 for many possible reasons e. g. that region isn't well resolved in the assemblies, bigger discrepancy in size between the two haplotypes, etc.

On a different note, when running your code on my df I had to restrict it only to numeric columns and then parse back the centromeres identifiers after calculating the p.adjust values; I hope this is fine. Also, more interestingly, at what level can I fit the offset(log(size)), I assume during plotting? Apologies, stupid question, I run it on the df fields for size (size_no-cent and size_cent) before using the code block you shared; however, I didn't manage to integrate the region as a predictor. Should it be something as follow:

glm(id ~ offset(log(size)), family=poisson(), data=df)

Thanks again!

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

The offset goes in the model fit, not the plot. It's what turns counts into rates, so what you're actually fitting is log(count/size).

Your formula won't run because id is the chromosome label rather than a count, and there's no region term in it. You need the data long rather than wide, two rows per chromosome:

long <- data.frame(
  chrom  = rep(df$id, 2),
  count  = c(df$variants_ce, df$variants_gw),
  size   = c(df$size_cent,   df$`size_no-cent`),
  region = rep(c("centromere", "rest"), each = nrow(df))
)

fit <- glm(count ~ region + offset(log(size)), family = quasipoisson, data = long)
exp(coef(fit))["regionrest"]

That last line is your rate ratio for rest vs centromere, with the overdispersion accounted for. Restricting to numeric columns and rejoining the IDs afterwards is fine.

On chr4, one thing specific to graph alignment worth considering: if your chr4 centromere graph already represents this sample's haplotypes well, those bases simply stop being called as variants. A better-matched graph produces fewer variants, not more, so low count with good coverage might be saying the chr4 graph is the good one rather than the broken one. Your other idea, a large size discrepancy between haplotypes splitting reads across two paths, would leave a different signature -- reduced read support per path rather than per region -- so looking at per-path support instead of total coverage should separate the two.

• 0 views • link
written by SeqBench • SeqBench

@Leo fantastic that clarifies everything.

I agree with your statements about the graph; in fact, as a first analysis I run a benchmark between linear vs. graph approach to ascertain which one was a better representation of the genomic landscape of the sample in question. Upon inspection, I also figured out what happen to the chr4 graph: it has been circularized for some reason...

The algorithm to build centromeres' pangenomes is very new, and I'll report this back to the developers; with this in mind, one issue might now be the exact definition of start/end positions during alignment in turns resulting in that low variants count, what do you think?

One last thing, you mentioned this in your first message:

With counts this large almost everything will be "significant", so read the effect size, not the p-value.

which I understand and, indeed, witness; however, with your last code I get a nested list (fit) – containing such information under the effects variable – but for one I wish to know how to interpret the 29.1173 and, secondly, why you apply the exp(coef()) functions to it?

Should I read the resulting 1.112012 in terms of significance? Thanks a lot for following-up and for the nice explanation!

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

Circularisation would do it, and I'd guess it's most of the explanation rather than a side detail. If the two haplotypes got collapsed into one circular structure, the differences between them are now baked into the graph topology instead of existing as alternate paths a read can disagree with, so they stop being callable as variants. That's the same mechanism I described, and circularisation is a concrete cause for it. The start/end ambiguity you're asking about is real but second-order: reads spanning the artificial origin get their alignments broken there, so you lose a window's worth of calls around one arbitrary point. That's tens of variants at most, not thousands. Worth reporting the circularisation itself to the developers though.

On the numbers, don't read fit$effects -- that's an internal by-product of the QR decomposition rather than anything interpretable, which is why 29.1173 doesn't correspond to anything meaningful. Use summary(fit), or coef(fit) if you just want the coefficients.

The exp() is because Poisson and quasipoisson use a log link, so coefficients come out on the log scale. coef(fit)["regionrest"] is log(rate ratio), and exponentiating puts it back on the ratio scale.

So 1.112 means variant density outside centromeres is about 11% higher than inside, averaged over all your chromosomes. That's an effect size, not a significance statement -- the p-value for it is in summary(fit). And it's a nice illustration of the earlier point: with counts in the hundreds of thousands that 11% will almost certainly come back with a tiny p-value, but 11% is not obviously biologically interesting. Worth deciding up front what size of difference would actually mean something to you.

One caveat: since that model only has region in it, 1.112 is one average across all chromosomes. It isn't comparable to the per-chromosome rate ratios from poisson.test, which is where the chr4 weirdness shows up.

• 0 views • link
written by SeqBench • SeqBench

Log in to answer this question.

More posts like this