Archive for August 10th, 2009

Multiple Testing in Genome-Wide Association Studies via Hidden Markov Models

Monday, August 10th, 2009

I just read a new paper out in “advanced access” in Bioinformatics:

Multiple Testing in Genome-Wide Association Studies via Hidden Markov Models

Wei et al.

Abstract

Motivation: Genome wide association studies (GWAS) interrogate common genetic variation across the entire human genome in an unbiased manner and hold promise in identifying genetic variants with moderate or weak effect sizes. However, conventional testing procedures, which are mostly p-value based, ignore the dependency and therefore suffer from loss of efficiency. The goal of this article is to exploit the dependency information among adjacent SNPs to improve the screening efficiency in GWAS.

Results: We propose to model the linear block dependency in the SNP data using hidden Markov Models. A compound decision-theoretic framework for testing HMM-dependent hypotheses is developed. We propose a powerful data-driven procedure (PLIS) that controls the false discovery rate (FDR) at the nominal level. PLIS is shown to be optimal in the sense that it has the smallest false negative rate (FNR) among all valid FDR procedures. By re-ranking significance for all SNPs with dependency considered, PLIS gains higher power than conventional p-value based methods. Simulation results demonstrate that PLIS dominates conventional FDR procedures in detecting disease associated SNPs. Our method is applied to analysis of the SNP data from a GWAS of type 1 diabetes. Compared to the BH procedure, PLIS yields more accurate results and has better reproducibility of findings.

Conclusion: The genomic rankings based on the our procedure are substantially different from the rankings based on the p-values. By integrating information from adjacent locations, the PLIS rankings benefit from the increased signal to noise ratio, hence our procedure often has higher statistical power and better reproducibility. This provide a promising direction in large-scale GWAS.

Summary

The topic is multiple testing correction in genome wide association studies (GWAS), which is probably one of the most important issues in such studies.  With the very large number of tests – typically hundreds of thousands to a million – you need to correct your significance value to avoid drowning in false positives.

The false discovery rate (FDR) method is a way of doing this, that essentially ranks the p-values and then picks the smallest while keeping the cumulative sum below the desired significance level.  Doing this in a GWAS ignores the dependency between tests caused by linkage disequilibrium, however, and this paper improves on this by taking the dependency into account.

They do this by fitting the data to a hidden Markov model where the hidden states are associated/not-associated and the emissions are z-values (the null for not associated and a mixture of normals for associated).  From this they can get a posterior probability of association/non-association for each marker, conditional on the test statistics for all markers in the genome: P(\theta_i \,|\, \mathbf{z}) where \theta_i is the state at marker i (0 if the marker is not associated with the disease and 1 if it is) and \mathbf{z} is the test statistics for all markers.

Now they consider all the P(\theta_i=0 \,|\, \mathbf{z}), that is all the posterior probabilities of not being associated, order them, and pick markers as long as the cumulative posterior probability is less than the significance threshold.

They say that this approach 1) guarantees that the false discovery rate is below the threshold and 2) that it is optimal in the sense that it is the method with that guarantee that has fewest false negatives, but they refer to an appendix for the proof of that, and that appendix is not in the paper, so I cannot really check that.  I would have loved to, though, since I want to know which assumptions about the data underlies this proof, but no matter…

Anyway, on with the summary.

They now validate the method with two simulation setups; one based on data simulated by a hidden Markov model, so matching the inference method and one with more realistic data.  For the first simulation study they show that the FDR guarantees are met and that the new method is more sensitive than those they compare it with.  For the second simulation study they essentially only show that it ranks true associations better than plain p-values.

They apply the method to a real data set and show that they are more successful in ranking markers that can be replicated in a replication cohort, again compared to plain p-values.

The good

First of all, I think it is an important problem to attack.  Doing so while taking the correlation between markers – and through that the correlation between their test statistics – is definitely the way to go.

Using hidden Markov models is also very sensible.  They are computationally efficient, usually easy to extend in various ways, and well founded in statistics so results are (relatively) easy to relate to.

The bad

I do have some problems with the method, though.  First some minor issues.

If the method really does compute the posterior probability of a marker being associated with the disease, then you would expect to consider all markers i where P(\theta_i=1\,|\,\mathrm{z})>P(\theta_i=0\,|\,\mathrm{z}) since those are the markers that are more likely to be assocated than not associated!  Picking only some of them means that for the rest of them you are essentially betting for a hypothesis less likely than the one you reject.

The issue here is, of course, that the value computed is in fact not the posterior probability of being associated.  The prior probability of association versus non-associated is probably not taken into account.  If it is, I couldn’t find any mentioning of it, at least.

If you include this prior belief you could just include the prior odds in your test and you would have a different approach to judging significance.  In practice it probably doesn’t matter much, so it is more an objection of aesthetics.

The ugly

Ok, now we come to the part about the paper I really didn’t like.  The simulation studies used to validate the method.

The first simulation study, I just can’t put much trust in. Not that I think there is anything dodgy in the results reported, but the simulations are from hidden Markov models matching the inference method, and there is no way that real data is generated that way.  There is nothing wrong with modelling the data as hidden Markov models – even if it is not generated by a process that even remotely resembles it – since it is just an analysis strategy anyway, but the simulation validation based on an unrealistic assumption is not particularly convincing…

The second simulation study is a more convincing setup, since here it is real LD data and a more realistic disease model.  However, here the FDR is not reported, only the sensitivity of seeing a marker in LD with a causal marker in top K of the ranking.  There is a better ranking with this approach than there is with just the p-values of the individual tests, but that says nothing about what the false discovery rate is.  So based on the results presented, I have no way of knowing what FDR to expect on realistic data (or how often I get a real hit below the FDR threshold for that matter).

I think this is a major problem with the validation of the method.  It is really only validated on data that is unlikely to resemble real GWAS data.  At least, the part of the method that has to do with FDR – the ranking results are okay.

Summary

Ok, I don’t want to end up sounding all negative.  It just looks that way since I ordered the criticism good -> bad -> ugly.

I stand by the criticism – I do think there are some problems with the validation – but all in all I like the method and I will definitely keep it in mind for my own future work.  I will look at that, as soon as they put up the source on CRAN (the paper just says that it will be made available there, but not when).

The main problem I have with the validation is really only the claims about false discovery rate.  Ok, since that is what the method is supposed to handle, that is a major problem, but as a method for ranking markers it looks pretty good.  Taking neighbouring markers into account in the analysis is what does this, I think, and is what we have also observed with our methods.

…and, you know, if you have a good ranking, maybe the false discovery rate isn’t all that important!  We only trust markers we validate in a replication data set anyway, and we will probably try to validate “top k” rather than all markers below a certain false discovery rate.  So in practical terms, the ranking is probably much more important than the multiple test correction.

Not so in the replication, of course, there you have to be strict about significance, but I’m not sure we need to be that strict in the initial discovery data set, as long as we don’t try to replicate thousands and thousands of markers, and we are probably not going to do that anyway.

I started my review by saying that correcting for multiple testing is very important in GWAS, and it is and it is very worthwhile to develop methods for it, but improving the ranking of markers is one of the problems that is even more important.
–
Wei, Z., Sun, W., Wang, K., & Hakonarson, H. (2009). MultipleTesting in Genome-Wide Association Studies via Hidden Markov Models Bioinformatics DOI: 10.1093/bioinformatics/btp476
222-225=-3

Posted in Paper reviews | No Comments »

Last week in the blogs

Monday, August 10th, 2009

Genetics and genomics

Physics

Research life

–

222-224=-2

Posted in Fun | No Comments »

Uncertainty and randomness, part II

Monday, August 10th, 2009

Ok, so in my previous post I considered tossing a (biased) coin and trying to figure out the probability of heads over tails.  The probability of tossing heads, p is a property of the coin and fixed, it is just unknown, so there will be some randomness in estimating it due to the randomness in coin tosses we use to estimate.

In the Bayesian analysis I used, I treated the parameter p as stochastic, because that is how that kind of analysis works, but it never really was.  In this post I will consider a random p.

I am going to consider an experiment where I am picking coins out of a bag, tossing them and recording whether they are heads or tails, and then throwing away the coin to pick another one.  I think that the coins are biased, but I don’t know how.  Can I figure that out?

Formalizing this mathematically, I am assuming that I am drawing coins each with a p probability of coming up head, where the ps are beta distributed with (unknown) parameters \alpha and \beta.  I get to toss each coin once to see if it turns up heads or tails, and then I have to throw the coin away, so for coin i I have a Bernoulli distribution with parameter p_i.

a <- 3 ; b <- 4
n <- 100
outcomes <- rbinom(n,1,rbeta(1,a,b))
s <- sum(outcomes)

The task is to figure out the distribution of ps, which means estimating hyper parameters \alpha and \beta.

Any ideas?

Doing the same as last time?

In the setup in my last post, we modelled the situation very similarly.  We drew a p at random and then tossed coins and got a distribution for p that way.  Could that work again?

We could pick a prior for p, say \mathrm{beta}(2,2) as last time, calculate the posterior based on the number of heads out of n tosses and see if that tends toward the real distribution of ps.

plot.analysis <- function(prior.a,prior.b ,real.a,real.b, n) {
  x <- seq(0,1,length.out=1000)

  # real p distribution
  real <- dbeta(x,real.a,real.b)

  # get prior
  prior <- dbeta(x,prior.a,prior.b)

  # perform experiment
  s <- sum(rbinom(n,1,rbeta(1,real.a,real.b)))
  # get posterior
  post <- dbeta(x,prior.a+s,prior.b+n-s)

  # plot results...
  plot(x, post, main='Prior, likelihood and posterior',
       type='l', xlab='', ylab='')
  lines(x, real, col='green')
  lines(x, prior, col='red')
}

opar <- par(mfrow=c(2,2), mar=c(2,2,2,2)+0.1)
plot.analysis(2,2, 3,4, 5)
plot.analysis(2,2, 3,4, 10)
plot.analysis(2,2, 3,4, 50)
plot.analysis(2,2, 3,4, 100)
par(opar)

(Ups, the title for these plots are wrong … there is the real distribution of ps but no likelihood … but I really cannot be bothered to re-do the plot right now…)

Ok, clearly that didn’t work.  Why?

For one thing, even if p is stochastic in this setup, there is only one p, but in the setup we are modelling there is one for each coin toss.

By the very nature of the beta distribution, the more heads and tails we connect, the “pointy’er” it will become, so if the real distribution of ps is relatively flat, it will look less and less like it.

A simple maximum likelihood solution

Since all we have to do to estimate the distribution of ps is to estimate the hyper parameters \alpha and \beta, we could just take the maximum likelihood approach.  We just need to work out the likelihood of our model and then maximize it with respect to the two.

If \pi(p;\alpha,\beta) is the distribution of ps, then the likelihood of ps and heads and tails is

\mathrm{lhd}(\alpha,\beta,p_1,\ldots,p_n) = \prod_i p_i^{h_i}(1-p_i)^{1-h_i} \pi(p_i)

where h_i is 1 if the ith coin toss is head and 0 otherwise.

We don’t want to maximize with respect to all the p_is – and with one coin toss for each we don’t have a lot of information to do that anyway – but we can just integrate them out, and since they are independent that gives us

\mathrm{lhd}(\alpha,\beta) = \prod_i \int p_i^{h_i}(1-p_i)^{1-h_i} \pi(p_i;\alpha,\beta) \,\mathrm{d}p_i

or

\log\mathrm{lhd}(\alpha,\beta) = s \log \int p\pi(p;\alpha,\beta) \,\mathrm{d}p + (n-s)\log \int (1-p)\pi(p;\alpha,\beta) \,\mathrm{d}p

where n is the number of coin tosses and s the number of heads.

With a beta prior for p

\pi(p;\alpha,\beta) = \frac{1}{B(\alpha,\beta)}p^{\alpha-1}(1-p)^{\beta-1}

where B(\alpha,\beta)=\int p^{\alpha-1}(1-p)^{\beta-1} \,\mathrm{d}p we get

\int p\pi(p;\alpha,\beta) \,\mathrm{d}p = B(\alpha+1,\beta)/B(\alpha,\beta)

and

\int (1-p)\pi(p;\alpha,\beta) \,\mathrm{d}p = B(\alpha,\beta+1)/B(\alpha,\beta)

so the log-likelihood is

\log\mathrm{lhd}(\alpha,\beta) = s\log B(\alpha+1,\beta) + (n-s)\log B(\alpha,\beta+1) - n\log B(\alpha,\beta).

You can probably optimize this with respect to \alpha and \beta analytically, but I’m just going to use a numerical optimisation to try it out:

a <- 3 ; b <- 4
n <- 100
s <- sum(rbinom(n,1,rbeta(1,a,b)))
logL <- function(a,b, n,s) {
  s*lbeta(a+1,b) + (n-s)*lbeta(a,b+1) - n*lbeta(a,b)
}
search <- nlm(function (x) -logL(x[1],x[2],n,s), c(3,4))
mle.a <- search$estimate[1]
mle.b <- search$estimate[2]
opar <- par(mfrow=c(1,2))
xs <- seq(1,10,length.out=500)
ys <- sapply(xs, function(x) logL(x,4,n,s))
plot(xs,ys,type='l',xlab=expression(alpha),main='Likelihood',ylab='logL')
abline(v=mle.a)
mtext(mle.a,at=mle.a)
ys <- sapply(xs, function(x) logL(3,x,n,s))
plot(xs,ys,type='l',xlab=expression(beta),main='Likelihood',ylab='logL')
abline(v=mle.b)
mtext(mle.b,at=mle.b)
par(opar)

and we can plot the true p density to compare it with the inferred:

x <- seq(0,1,length.out=100)
true <- dbeta(x,3,4)
inferred <- dbeta(x,mle.a,mle.b)
plot(x,true,type='l',col='red',
     main='True and inferred distribution',
     xlab='', ylab='',
     ylim=range(true,inferred))
lines(x,inferred, col='blue')

You might get different results here; when I experimented with it I didn’t always get results this close to the true values, but I picked a nice result for the plots here…

–

222-223=-1

Posted in Fun, Teaching, Work | Tags: ,  | No Comments »