SNV calling from NGS data

A set of hypothetical NGS reads are shown, aligned against a reference sequence. At the annotated locus, the reads contain a mixture of A/G nucleotides, against the A reference allele. Depending on the prior genotype probabilities, and the chosen error model, this may be called as a heterozygous SNV (genotype AG predicted), the G nucleotides may be classified as errors and no variant called (genotype AA predicted), or alternatively the A nucleotides may be classified as errors and a homozygous SNV called (genotype GG predicted).

In an ideal error free world with high read coverage, the task of variant calling from the results of a NGS data alignment would be simple; at each locus (position on the genome) the number of occurrences of each distinct nucleotide among the reads aligned at that position can be counted, and the true genotype would be obvious; either AA if all nucleotides match allele A, BB if they match allele B, or AB if there is a mixture. However, when working with real NGS data this sort of naive approach is not used, as it cannot account for the noise in the input data.[4] The nucleotide counts used for base calling contain errors and bias, both due do the sequenced reads themselves, and the alignment process. This issue can be mitigated to some extent by sequencing to a greater depth of read coverage, however this is often expensive, and many practical studies require making inferences on low coverage data.[1]

Probabilistic methods aim to overcome the above issue, by producing robust estimates of the probabilities of each of the possible genotypes, taking into account noise, as well as other available prior information that can be used to improve estimates. A genotype can then be predicted based on these probabilities, often according to the MAP estimate.

Probabilistic methods for variant calling are based on Bayes' theorem. In the context of variant calling, Bayes' Theorem defines the probability of each genotype being the true genotype given the observed data, in terms of the prior probabilities of each possible genotype, and the probability distribution of the data given each possible genotype. The formula is:

${\displaystyle {\begin{aligned}P(G\mid D)&={\frac {P(D\mid G)P(G)}{P(D)}}\\[8pt]&={\frac {P(D\mid G)\,P(G)}{\sum \limits _{i=1}^{n}P(D\mid G_{i})\,P(G_{i})}}\end{aligned}}}$

In the above equation:

• ${\displaystyle D}$ refers to the observed data; that is, the aligned reads
• ${\displaystyle G}$ is the genotype whose probability is being calculated
• ${\displaystyle G_{i}}$ refers to the ith possible genotype, out of n possibilities

Given the above framework, different software solutions for detecting SNVs vary based on how they calculate the prior probabilities ${\displaystyle P(G)}$, the error model used to model the probabilities ${\displaystyle P(D\mid G)}$, and the partitioning of the overall genotypes into separate sub-genotypes, whose probabilities can be individually estimated in this framework.[5]

The calculation of prior probabilities depends on available data from the genome being studied, and the type of analysis being performed. For studies where good reference data containing frequencies of known mutations is available (for example, in studying human genome data), these known frequencies of genotypes in the population can be used to estimate priors. Given population wide allele frequencies, prior genotype probabilities can be calculated at each locus according to the Hardy–Weinberg equilibrium.[6] In the absence of such data, constant priors can be used, independent of the locus. These can be set using heuristically chosen values, possibly informed by the kind of variations being sought by the study. Alternatively, supervised machine-learning procedures have been investigated that seek to learn optimal prior values for individuals in a sample, using supplied NGS data from these individuals.[4]

The error model used in creating a probabilistic method for variant calling is the basis for calculating the ${\displaystyle P(D\mid G)}$ term used in Bayes' theorem. If the data was assumed to be error free, then the distribution of observed nucleotide counts at each locus would follow a Binomial Distribution, with 100% of nucleotides matching the A or B allele respectively in the AA and BB cases, and a 50% chance of each nucleotide matching either A or B in the AB case. However, in presence of noise in the read data this assumption is violated, and the ${\displaystyle P(D\mid G)}$ values need to account for the possibility that erroneous nucleotides are present in the aligned reads at each locus.

A simple error model is to introduce a small error to the data probability term in the homozygous cases, allowing a small constant probability that nucleotides which don't match the A allele are observed in the AA case, and respectively a small constant probability that nucleotides not matching the B allele are observed in the BB case. However more sophisticated procedures are available which attempt to more realistically replicate the actual error patterns observed in real data in calculating the conditional data probabilities. For instance, estimations of read quality (measured as Phred quality scores) have been incorporated in these calculations, taking into account the expected error rate in each individual read at a locus.[7] Another technique that has successfully been incorporated into error models is base quality recalibration, where separate error rates are calculated – based on prior known information about error patterns – for each possible nucleotide substitution. Research shows that each possible nucleotide substitution is not equally likely to show up as an error in sequencing data, and so base quality recalibration has been applied to improve error probability estimates.[6]

In the above discussion, it has been assumed that the genotype probabilities at each locus are calculated independently; that is, the entire genotype is partitioned into independent genotypes at each locus, whose probabilities are calculated independently. However, due to linkage disequilibrium the genotypes of nearby loci are in general not independent. As a result, partitioning the overall genotype instead into a sequence of overlapping haplotypes allows these correlations to be modelled, resulting in more precise probability estimates through the incorporation of population-wide haplotype frequencies in the prior. The use of haplotypes to improve variant detection accuracy has been applied successfully, for instance in the 1000 Genomes Project.[8]

As an alternative to probabilistic methods, heuristic methods exist for performing variant calling on NGS data. Instead of modelling the distribution of the observed data and using Bayesian statistics to calculate genotype probabilities, variant calls are made based on a variety of heuristic factors, such as minimum allele counts, read quality cut-offs, bounds on read depth, etc. Although they have been relatively unpopular in practice in comparison to probabilistic methods, in practice due to their use of bounds and cut-offs they can be robust to outlying data that violate the assumptions of probabilistic models.[9]

An important part of the design of variant calling methods using NGS data is the DNA sequence used as a reference to which the NGS reads are aligned. In human genetics studies, high quality references are available, from sources such as the HapMap project,[10] which can substantially improve the accuracy of the variant calls made by variant calling algorithms. As a bonus, such references can be a source of prior genotype probabilities for Bayesian-based analysis. However, in the absence of such a high quality reference, experimentally obtained reads can first be assembled in order to create a reference sequence for alignment.[1]

Various methods exist for filtering data in variant calling experiments, in order to remove sources of error/bias. This can involve the removal of suspicious reads before performing alignment and/or filtering of the list of variants returned by the variant calling algorithm.

Depending on the sequencing platform used, various biases may exist within the set of sequenced reads. For instance, strand bias can occur, where there is a highly unequal distribution of forward vs reverse directions in the reads aligned in some neighborhood. Additionally, there may occur an unusually high duplication of some reads (for instance due to bias in PCR). Such biases can result in dubious variant calls – for instance if a fragment containing a PCR error at some locus is over amplified due to PCR bias, that locus will have a high count of the false allele, and may be called as a SNV – and so analysis pipelines frequently filter calls based on these biases.[1]

Until recently, software tools for carrying out this form of analysis have been heavily underdeveloped, and were based on the same algorithms used to detect germline variations. Such procedures are not optimized for this task, because they do not adequately model the statistical correlation between the genotypes present in multiple tissue samples from the same individual.[3]

More recent investigations have resulted in the development of software tools especially optimized for the detection of somatic mutations from multiple tissue samples. Probabilistic techniques have been developed that pool allele counts from all tissue samples at each locus, and using statistical models for the likelihoods of joint-genotypes for all the tissues, and the distribution of allele counts given the genotype, are able to calculate relatively robust probabilities of somatic mutations at each locus using all available data.[3][12] In addition there has recently been some investigation in machine learning based techniques for performing this analysis.[13][14][15][16]

In 2021, the Sequencing Quality Control Phase 2 Consortium[17] has published a number of studies that investigated the effects of sample preparations, sequencing library kits, sequencing platforms, and bioinformatics workflows on the accuracy of somatic SNV detection[18] based on a pair of tumor-normal cell lines that the Consortium has established as the reference samples, data, and call sets.[19]