<p>The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.</p>
<p>My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.</p>
<p>This is what I think might work</p>
<ol>
<li>Calculate your true OR, qtl snps versus all others for say TFBS</li>
<li>Divide your snps according to MAF into 10 bins / stratas</li>
<li>Within each bin, we permute the label whether the snp is a qtl</li>
<li>Calculate OR of selection / others for TFBS</li>
<li>Repeat 1-3 N times to collect the permutated OR</li>
<li>p value would be number of (permutated OR &gt; observed OR+1)/(N+1)</li>
</ol>
<p>Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:</p>
<pre><code>snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))
</code></pre>
<p>We have the MAF for all snps, with those in qtl being different from others and put them in bins.</p>
<pre><code>snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))
</code></pre>
<p>We sort the snps according to MAF to make it easier:</p>
<pre><code>o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]
</code></pre>
<p>And we create the bin:</p>
<pre><code>bin_MAF = cut(snp_MAF,10)
</code></pre>
<p>Lastly a function to calculate OR:</p>
<pre><code>calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}
</code></pre>
<p>So our observed OR (in log) is:</p>
<pre><code>obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859
</code></pre>
<p>Now we collect the OR under permutation:</p>
<pre><code>N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)
</code></pre>
<p><a href="https://i.sstatic.net/sK6zH.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/sK6zH.png" alt="enter image description here" /></a></p>
<p>Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation</p>

CC BY-SA 4.0
2020-09-11T15:05:47.680Z
2020-09-11T15:05:47.680Z
StupidWolf
<p>I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).</p>
<p>A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:</p>
<ol>
<li>Take a single <strong>random sample of all SNPs</strong> with size 10<em>n</em> (<em>n</em> is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).</li>
<li>Do a <strong>bootstrapping</strong>, taking multiple random samples (excluding qtl-SNPs), each of size <em>n</em>, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a <em>p</em>-value.</li>
</ol>
<p>Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?</p>
<p>Thanks.</p>

CC BY-SA 4.0
2020-09-04T14:43:19.000Z
also posted in https://stats.stackexchange.com/questions/485997/odds-ratio-and-enrichment-of-snps-in-gene-regions. Avoid posting the same question twice. I would suggest thinking about the nature of  your question which forum is more suitable
StupidWolf
CC BY-SA 4.0
2020-09-04T15:02:00.543Z
@StupidWolf the MAF distribution of the full 7M and the qtl-SNPs is different. That's why I don't want to consider the full panel and instead use a more comparable 'background' of SNPs to have a meaningful estimation of the enrichment.
sergiovm
<r><snp><gwas><p-values><qtls>
Odds ratio and enrichment of SNPs in gene regions?
# Odds ratio and enrichment of SNPs in gene regions?

(Asked by: username_0 on 2020-09-04)

I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).

A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:

1. Take a single **random sample of all SNPs** with size 10*n* (*n* is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).
2. Do a **bootstrapping**, taking multiple random samples (excluding qtl-SNPs), each of size *n*, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a *p*-value.

Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?

Thanks.

**Tags:** r, snp, gwas, p-values, qtls
> **Comments:**
> **username_1**: also posted in https://stats.stackexchange.com/questions/485997/odds-ratio-and-enrichment-of-snps-in-gene-regions. Avoid posting the same question twice. I would suggest thinking about the nature of  your question which forum is more suitable  
> **username_0**: @StupidWolf the MAF distribution of the full 7M and the qtl-SNPs is different. That's why I don't want to consider the full panel and instead use a more comparable 'background' of SNPs to have a meaningful estimation of the enrichment.  

---

## 1 answers

### Answer by username_1 on 2020-09-11. Score: 3 | Accepted answer

The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.
My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.
This is what I think might work

Calculate your true OR, qtl snps versus all others for say TFBS
Divide your snps according to MAF into 10 bins / stratas
Within each bin, we permute the label whether the snp is a qtl
Calculate OR of selection / others for TFBS
Repeat 1-3 N times to collect the permutated OR
p value would be number of (permutated OR > observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

We have the MAF for all snps, with those in qtl being different from others and put them in bins.
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

We sort the snps according to MAF to make it easier:
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

And we create the bin:
bin_MAF = cut(snp_MAF,10)

Lastly a function to calculate OR:
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

So our observed OR (in log) is:
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

Now we collect the OR under permutation:
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)


Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation


---

