Odds ratio and enrichment of SNPs in gene regions?

Question

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 10n (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.

Answer 1

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

1. Calculate your true OR, qtl snps versus all others for say TFBS
2. Divide your snps according to MAF into 10 bins / stratas
3. Within each bin, we permute the label whether the snp is a qtl
4. Calculate OR of selection / others for TFBS
5. Repeat 1-3 N times to collect the permutated OR
6. 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