SNV calling from NGS data

SNV calling from NGS data is any of a range of methods for identifying the existence of

In addition to the usual application domain of SNP genotyping, these techniques have been successfully adapted to identify rare SNPs within a population,[2] as well as detecting somatic SNVs within an individual using multiple tissue samples.^[3]

Methods for detecting germline variants

Most NGS based methods for SNV detection are designed to detect

(SNPs). Technically the term SNP only refers to these kinds of variations, however in practice they are often used synonymously with SNV in the literature on variant calling. In addition, since the detection of germline SNVs requires determining the individual's genotype at each locus, the phrase "SNP genotyping" may also be used to refer to this process. However this phrase may also refer to wet-lab experimental procedures for classifying genotypes at a set of known SNP locations.

The usual process of such techniques are based around:^[1]

1. Filtering the set of NGS reads to remove sources of error/bias
2. Aligning the reads to a reference genome
3. Using an algorithm, either based on a statistical model or some heuristics, to predict the likelihood of variation at each locus, based on the quality scores and allele counts of the aligned reads at that locus
4. Filtering the predicted results, often based on metrics relevant to the application
5. SNP annotation to predict the functional effect of each variation.

The usual output of these procedures is a VCF file.

Probabilistic methods

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

. 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]

Prior genotype probability estimation

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

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]

Error models for data observations

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

, 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]

Partitioning of the genotype

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

Heuristic based algorithms

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]

Reference genome used for alignment

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]

Pre-processing and filtering of results

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]

Methods for detecting somatic variants

In addition to methods that align reads from individual sample(s) to a reference genome in order to detect

Recent developments

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]

List of available software

• VarNet
• Big Data Genomics: Avocado
• Beagle
• DeepVariant
• Freebayes
• GATK (including MuTect)
• IMPUTE2
• JointSNVMix
• MaCH
• Magnolia DCNN
• NeuSomatic
• NGSEP
• Pisces
• Platypus
• realSFS
• Reveel
• SAMtools
• SNVmix
• SOAPsnp
• SomaticSeq
• SomaticSniper
• Strelka
• VarDict
• VarScan

References