have many error likelihoods, how to combine to get a confidence or p value?

Question

I'm working in bioinformatics and its been a long time since I dusted on my statistics. Basically I'm working on variant calling which amounts to sequencing a large number of sequence reads and analyzing them to see how they differ from a known sample. The same position in a genome will be covered many times like this:

AAAAAAAAAAAAAAAA - Reference sample

AAAAAAAAAAAAAAAA - Sequenced read1
AAAAAAAAAAAAAAAA - Sequenced read2
AAAAAAATAAAAAAAA - Sequenced read3
 AAAAAATAAAAAAAAA - Sequenced read4

The whole game is to figure out when the reads represent a sample that differs from the reference. In the example above there is a single position where half of the reads suggest there is a variant of A to T near the middle.

One way to associate a p-value with the likelihood that the middle position is in fact different from the reference sample is to compute, or collect a sequence error rate, E, and use the total depth and depth of just the T's in a binomial test to see if we have a significant non-reference signal that is higher than expected given the error rate.

I'm working on a new method for this where I am trying to combine multiple types of sequence error. Basically, for any position in my sample I have a list like:

The likelihood of seeing an A to T error in position 8 in a read (E1)
The likelihood of seeing an A to T error in position 7 in a read (E2)
...

and so on where I am calculating many types of error metrics, and they're not completely independent of each other.

Now in a real world example I have maybe 100 sequence reads covering a position and for each read that covers it I have an error probability of seeing that base substitution due to sequence error as a function of the within-read position.

Can you point in the right direction to figure out how to combine a list of error likelihoods into a p-value?

One basic thing I can think of is to just multiply my error likelihoods to get the overall likelihood that these sequence reads reported a sequence at a position by chance after which I can apply a threshold to the final likelihoods to float my putative true-positives to the top of the heap, but I feel as though there must be a better way to think about this.

Answer 1

If errors along a single sequence read or across sequence reads are not independent, then as @jsk said in a comment you can't just multiply probabilities. If you have a detailed model of how the sequence-read errors are related than you could in principle devise a way to combine those error probabilities to give a mathematically correct error estimate.

In practice, however, I wonder whether you are trying to re-invent something for which reliable methods already exist. My concern has to do both with the statistical content and with other practical considerations in variant calling, with the latter perhaps more of a problem.

For example, sequence data typically include per-base per-read Phred scores, related to the probability of a technical base-call error. These error probabilities are calculated from several properties of the raw data as calibrated against well-characterized sequences, with algorithms updated as technology changes.

There certainly can be sequence-specific differences in read-error probabilities, an issue examined for example in this paper. If you have evidence that Phred scores aren't working well in particular sequence contexts, that could be of some interest. Otherwise, the Phred scores provide a well accepted measure of per-base per-read error probabilities.

The MuTect method for variant calling compares a model having no variant at a site against a model having a variant allele at the site, using individual Phred scores to account for sequencing noise at the site of interest and at other sites along the reads. It does assume independence of sequencing errors across reads. It should be possible to build on the MuTect variant-detection approach, explained in detail on the paper's Online Methods, if you have a particular model of non-independence in mind.

There are practical issues in sequencing that can be even more important than technical read errors. For example, as the online methods of the MuTect paper say:

A common source of false positive mutation calls is contamination of the tumor DNA with DNA from other individuals [e.g., from those handling the sample]. Germ-line SNPs in the contaminating DNA appear as somatic mutations...even 2% contamination can give rise to 166 false positive calls per megabase and 10 false positive calls per megabase when excluding known SNP sites.

Compare that against the technical read error probabilities judged by Phred scores, with a reasonably good quality Phred score of 30 having a read-error probability of only 0.1%. The MuTect method thus incorporates a correction for potential DNA contamination in the sample. A potential oxoG artifact introduced via sample processing also has to be dealt with.

Then there are particular DNA sequences that tend to lead to false positives. These may come from sequence-specific errors either in preparing samples for sequencing or in technical read errors. Applying a mutation-calling algorithm to many known normal samples identifies such loci and determines a "panel of normals" against which to compare putative mutations in test samples. A mutation tentatively identified in a test sample is not called if it is within the panel of normals.

This panel-of-normals approach thus protects against systematic sequence-specific false-positive calls even if we don't know the mechanisms underlying the errors. This empirical approach seems likely to work better in a practice than a theoretical approach necessarily based on the handful of error mechanisms that we do understand.