# Comparing mutation frequency between a case and a pool of controls

I"m working in genomics and trying to come up with the appropriate statistical test for my question.

To call mutations in a tumor's DNA, we use sequencing that samples from the total population of cells in a tumor. For each genomic position, we then get two numbers: The number of reads that support a mutation and the number of reads that support the wildtype sequence (and summing those gives the total number of reads).

If I want to identify mutations that are unique to the tumor sample (as compared to the normal control), I can set up a contingency table and use Fisher's exact test. (e.g. Is 21/100 mutation-supporting reads at this site significantly different from 2/98 in the control?)

What comparable test is appropriate when I have a pool of several normal samples? I want to test whether the frequency of mutation at a particular site in my tumor sample is significantly higher than it is across 10 normal samples. So I'm comparing one case (21/100) to many controls: (1/94, 3/85, 0/100 ... ). The null hypothesis is that there is no difference between case and control.

• Addendum: if I'm assuming that the normal samples are quantifying the 'noise' rate in the method, is it reasonable to combine their data? (and thereby 'hide' the information that they're independent estimates of the rate) What do I lose by doing so? Commented Jul 29, 2015 at 21:58
• This problem has been gone over extensively for the TCGA project. Do a web search for the MuTect and MutSig programs developed by the Broad Institute, for example. These explicitly model error patterns in sequence reads and so forth. Unless all of your tumor samples have same-patient controls, you will get into trouble if you start comparing tumor variants against separate normal cases, as many "normal" variants have yet to be identified.
– EdM
Commented Jul 29, 2015 at 22:52
• I'm familiar with both of those programs (and am an author on many of the TCGA papers). Neither addresses the question I'm asking, which is about quantifying error rates in ultra-deep sequencing, not about simply identifying somatic mutations). Mutect is roughly equivalent to the 2x2 Fisher test I describe - I need the statistical expansion to a pool. Commented Jul 30, 2015 at 0:53
• Because you didn't apparently sample normal tissue from the persons with cancer, you do not get to use paired designs. Only 10 normal samples seems rather underpowered for an unpaired design, so I'm not sure you can say what a "wild type" really is.
– DWin
Commented May 4, 2016 at 17:50

Logistic regression with samples as the independent variable would be a useful way to proceed. Reshape the data so that there is one row per read, with a 0 value for wild type and 1 for mutant in the first column, and a second column identifying the sample. Your null hypothesis then is a particular pre-specified contrast of the tumor sample against the mean of the controls, avoiding multiple-comparison issues.

This has some advantages: it reduces to chi-square when there is only 1 tumor and 1 control, it takes differences in total counts among samples into account, it generalizes to multinomial regression if you want to consider all 4 bases at the position instead of just mutant/wild-type, and it allows for testing differences among the controls. See this page for further discussion of this approach.

That said, I'm still concerned about biological interpretation of the results if you don't have a normal-DNA match for the tumor.

• Thanks for the suggestion - reading up on it now. FWIW, I understand your concern, but I know that this mutation is somatic, as it's called in a primary tumor at high VAF. The patient has gone into morphologic remission now, and the question is whether the mutation is still present at a detectable level (using multiple samples instead of just a single matched normal is intended to give a more accurate assessment of the error rate is at that position, which includes sequence context, etc) Commented Jul 30, 2015 at 15:21
• If you are interested in whether the mutation is gone, don't forget that therapy can change the VAF by relative selection of subclones from a heterogeneous tumor. So the mutation might still be present even if not "detectable" at a given read depth. Consider also using the "Sensitivity Calculation" at the end of the MuTect paper "Online Methods" section to determine how low a VAF you can reasonably detect. Use the controls to estimate the context-specific base-error rate (and thus the Phred score) for that calculation.
– EdM
Commented Jul 30, 2015 at 16:01

My answer seems like more than a comment but not quite an answer. Sorry if I'm posting in the wrong format.

I'll take the other side of "how do you know it's somatic?" If you want to use normal samples to estimate context-specific technical error rates, do you risk (and/or do you care about) low-frequency residual clones confounding that estimate? I suppose it depends on whether or not you can assume non-zero residual disease can exist in these normal samples at a frequency in the ballpark of your error rate (~0.1-1%). Are your normal samples from the same tissue type? If you have on the order of 10 I'm assuming these are all hematological. Ideally, I'd consider estimating these error rates from completely independent tissues, or, more practically, from different patients (samples processed in the same way, but from individuals who never had evidence of harboring the variant in question). Just an opinion.

Regarding the actual statistical approach, I wonder if a likelihood framework would be useful here. Your data are generated by a binomial process, where the number of non-reference reads, X, out of n sequencing reads is given by X~Bin(n,p). Here p is a reflection the the true variant allele frequency/fraction (VAF) and your technical error rate:

p=P(true non-reference)(1-error rate) + P(true reference)(error rate)

So the assumption that VAF=0 in control samples reduces to p=(1)(error rate).

Your tumor sample and N control samples are N+1 independent experiments:

X[i]~Bin(n[i],p[i])

Your null hypothesis is that p[i] is the same for all samples (estimate for p comes from pooling all tumor and control samples). Your alternative is that p[i] is different for the tumor sample versus your control samples (maximum likelihood estimate for tumor p[i] is derived from the tumor sample, estimate for control p[i] is derived from pooled control samples). Your likelihood statistic is then:

D=-2ln(p(observed non-reference counts X[1], X[2], ..., X[N+1] | jointly estimated p)) + 2ln(p(observed non-reference counts X[1], X[2], ..., X[N+1] | p estimated independently for tumor vs. control samples))

You can then calculate a p-value for this test statistic either empirically or using a chi-square distribution with 1 degree of freedom.