en>Malcolma: added Category:Category theory; removed {{uncategorized}} using HotCat

2013-12-08T11:08:54Z

added Category:Category theory; removed {{uncategorized}} using HotCat

New page

{{Orphan|date=October 2013}}

'''SNV calling from NGS data''' refers to a range of methods for identifying the existence of [[single nucleotide polymorphism|single nucleotide variants]] (SNVs) from the results of [[Next generation sequencing|NGS]] experiments. These are computational techniques, and are in contrast to specialist experimental methods based on known population-wide single nucleotide polymorphisms (see [[SNP genotyping]]). Due to the increasing abundance of NGS data, these techniques are becoming increasingly popular for performing SNP genotyping, with a wide variety of algorithms designed for specific experimental designs and applications.<ref name = "Nielsen 2011">{{cite journal
|title=Genotype and SNP calling from next-generation sequencing data
|author=Nielsen, Rasmus and Paul, Joshua S and Albrechtsen, Anders and Song, Yun S
|journal=Nature Reviews Genetics
|volume=12
|number=6
|pages=443–451
|year=2011
|publisher=Nature Publishing Group
}}</ref> In addition to the usual application domain of SNP genotyping, these techniques have been successfully adapted to identify rare SNPs within a population,<ref name = "Bansal 2010">{{cite journal
|title=A statistical method for the detection of variants from next-generation resequencing of DNA pools
|author=Bansal, Vikas
|journal=Bioinformatics
|volume=26
|number=12
|pages=i318-i324
|year=2010
|publisher=Oxford Univ Press
}}</ref> as well as detecting [[somatic]] SNVs within an individual using multiple tissue samples.<ref name = "Roth 2012">{{cite journal
|title=JointSNVMix: a probabilistic model for accurate detection of[somatic mutations in normal/tumour paired next-generation sequencing data
|author=Roth, Andrew and Ding, Jiarui and Morin, Ryan and Crisan, Anamaria and Ha, Gavin and Giuliany, Ryan and Bashashati, Ali and Hirst, Martin and Turashvili, Gulisa and Oloumi, Arusha and others
|journal=Bioinformatics
|volume=28
|number=7
|pages=907–913
|year=2012
|publisher=Oxford Univ Press
}}</ref>

== Methods for detecting germline variants ==

Most NGS based methods for SNV detection are designed to detect [[germline]] variations in the individual's genome. These are the mutations that an individual biologically inherits from their parents, and are the usual type of variants searched for when performing such analysis (except for certain specific applications where [[#Methods for detecting somatic variants|somatic mutations]] are sought). Very often, the searched for variants occur with some (possibly rare) frequency, throughout the population, in which case they may be referred to as [[single nucleotide polymorphisms]] (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:<ref name="Nielsen 2011"/>
# Filtering the set of NGS reads to remove sources of error/bias
# Aligning the reads to a reference genome
# 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
# Filtering the predicted results, often based on metrics relevant to the application
The usual output of these procedures is a [[Variant Call Format|VCF]] file.

=== Probabilistic methods ===

[[File:Heterozygous SNV call, from aligned NGS reads.png|thumb|right|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 valled (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 (genetics)|locus]] (position on the genome) the number of occurrences of each distinct nucleotide among the reads alinged 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.<ref name="Martin 2010">{{cite journal
|title=SeqEM: an adaptive genotype-calling approach for next-generation sequencing studies
|author=Martin, Eden R and Kinnamon, DD and Schmidt, Michael A and Powell, EH and Zuchner, S and Morris, RW
|journal=Bioinformatics
|volume=26
|number=22
|pages=2803–2810
|year=2010
|publisher=Oxford Univ Press
}}</ref> 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.<ref name="Nielsen 2011"/>

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 [[Maximum a posteriori estimation|MAP]] estimate.

Probabilistic methods for variant calling are based on [[Bayes Theorem|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:
:<math>
\begin{align}
P(G|D) &= \frac{P(D|G) P(G)}{P(D)}\\[8pt]
&= \frac{P(D|G)\,P(G)}{\sum\limits_{i=1}^{n} P(D|G_i)\,P(G_i)}\\[8pt]
\end{align}</math>
In the above equation:
*<math>D</math> refers to the observed data; that is, the aligned reads
*<math>G</math> is the genotype whose probability is being calculated
*<math>G_i</math> refers to the ''i''th 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 <math>P(G)</math>, the error model used to model the probabilities <math>P(D|G)</math>, and the partitioning of the overall genotypes into separate sub-genotypes, whose probabilities can be individually estimated in this framework.<ref name="You 2012">{{cite journal
|title=SNP calling using genotype model selection on high-throughput sequencing data
|author=You, Na and Murillo, Gabriel and Su, Xiaoquan and Zeng, Xiaowei and Xu, Jian and Ning, Kang and Zhang, Shoudong and Zhu, Jiankang and Cui, Xinping
|journal=Bioinformatics
|volume=28
|number=5
|pages=643–650
|year=2012
|publisher=Oxford Univ Press
}}</ref>

==== 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 [[Hardy Weinberg Equilibrium]].<ref name ="Li 2009">{{cite journal
|title=SNP detection for massively parallel whole-genome resequencing
|author=Li, Ruiqiang and Li, Yingrui and Fang, Xiaodong and Yang, Huanming and Wang, Jian and Kristiansen, Karsten and Wang, Jun
|journal=Genome research
|volume=19
|number=6
|pages=1124–1132
|year=2009
|publisher=Cold Spring Harbor Lab
}}</ref> 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.<ref name="Martin 2010">{{cite jounral
|title=SeqEM: an adaptive genotype-calling approach for next-generation sequencing studies
|author=Martin, Eden R and Kinnamon, DD and Schmidt, Michael A and Powell, EH and Zuchner, S and Morris, RW
|journal=Bioinformatics
|volume=26
|number=22
|pages=2803-2810
|year=2010
|publisher=Oxford Univ Press
}}</ref>

==== Error models for data observations ====

The error model used in creating a probabilistic method for variant calling is the basis for calculating the <math>P(D|G)</math> 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 <math>P(D|G)</math> 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 score|Phred]] quality scores) have been incorporated in these calculations, taking into account the expected error rate in each individual read at a locus.<ref name="Li 2008">{{cite journal
|title=Mapping short DNA sequencing reads and calling variants using mapping quality scores
|author=Li, Heng and Ruan, Jue and Durbin, Richard
|journal=Genome research
|volume=18
|number=11
|pages=1851–1858
|year=2008
|publisher=Cold Spring Harbor Lab
}}</ref> 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.<ref name="Li 2009"/>

==== 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 [[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]].<ref name="Abecasis 2010">{{cite journal
|title=A map of human genome variation from population-scale sequencing.
|author=Abecasis, GR and Altshuler, David and Auton, A and Brooks, LD and Durbin, RM and Gibbs, Richard A and Hurles, Matt E and McVean, Gil A and Bentley, DR and Chakravarti, A and others
|journal=Nature
|volume=467
|number=7319
|pages=1061–1073
|year=2010
}}</ref>

=== Heuristic based algorithms ===

As an alternative to probabilistic methods, [[Heuristic (computer science)|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.<ref name="Koboldt 2012">{{cite journal
|title=VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing
|author=Koboldt, Daniel C and Zhang, Qunyuan and Larson, David E and Shen, Dong and McLellan, Michael D and Lin, Ling and Miller, Christopher A and Mardis, Elaine R and Ding, Li and Wilson, Richard K
|journal=Genome research
|volume=22
|number=3
|pages=568–576
|year=2012
|publisher=Cold Spring Harbor Lab
}}</ref>

=== 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 align the NGS reads to. In human genetics studies, high quality references are available, from sources such as the HapMap project,<ref name="Gibbs 2003">{{cite journal
|title=The international HapMap project
|author=Gibbs, Richard A and Belmont, John W and Hardenbol, Paul and Willis, Thomas D and Yu, Fuli and Yang, Huanming and Ch'ang, Lan-Yang and Huang, Wei and Liu, Bin and Shen, Yan and others
|journal=Nature
|volume=426
|number=6968
|pages=789–796
|year=2003
|publisher=Nature Publishing Group
}}</ref> 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 [[sequence assembly|assembled]] in order to create a reference sequence for alignment.<ref name="Nielsen 2011"/>

=== 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 [[Polymerase chain reaction|PCR]]). Such biases can result in dubious variant calls - for instance if a read containing a sequencing error at some locus is duplicated due to a 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.<ref name="Nielsen 2011"/>

== Methods for detecting somatic variants ==

In addition to methods that align reads from individual sample(s) to a reference genome in order to detect [[germline]] genetic variants, reads from multiple tissue samples within a single individual can be aligned and compared in order to detect somatic variants. These variants correspond to [[genetic mutation|mutations]] that have occurred [[de novo]] within groups of [[somatic cell]]s within an individual (that is, they are not present within the individual's germline cells). This form of analysis has been frequently applied to the study of [[cancer]], where many studies are designed around investigating the profile of somatic mutations within cancerous tissues. Such investigations have resulted in diagnostic tools that have seen clinical application, and are used to improve scientific understanding of the disease, for instance by the discovery of new cancer-related genes, identification of involved [[gene regulatory networks]] and [[metabolic pathways]], and by informing models of how tumors grow and evolve.<ref name = "Shyr 2013">{{cite journal
|title=Next generation sequencing in cancer research and clinical application
|author=Shyr, Derek and Liu, Qi and others
|journal=Biological procedures online
|volume=15
|number=4
|year=2013
}}</ref>

=== 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.<ref name="Roth 2012"/>

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.<ref name="Roth 2012"/><ref name="Larson 2012">{{cite journal
|title=SomaticSniper: identification of somatic point mutations in whole genome sequencing data
|author=Larson, David E and Harris, Christopher C and Chen, Ken and Koboldt, Daniel C and Abbott, Travis E and Dooling, David J and Ley, Timothy J and Mardis, Elaine R and Wilson, Richard K and Ding, Li
|journal=Bioinformatics
|volume=28
|number=3
|pages=311–317
|year=2012
|publisher=Oxford Univ Press
}}</ref> In addition there has recently been some investigation in [[machine learning]] based techniques for performing this analysis.<ref name="Ding 2012">{{cite journal
|title=Feature-based classifiers for somatic mutation detection in tumour--normal paired sequencing data
|author=Ding, Jiarui and Bashashati, Ali and Roth, Andrew and Oloumi, Arusha and Tse, Kane and Zeng, Thomas and Haffari, Gholamreza and Hirst, Martin and Marra, Marco A and Condon, Anne and others
|journal=Bioinformatics
|volume=28
|number=2
|pages=167–175
|year=2012
|publisher=Oxford Univ Press
}}</ref>

== List of available software ==
*[http://soap.genomics.org.cn/index.html SOAPsnp]
*[http://128.32.118.212/thorfinn/realSFS realSFS]
*[http://samtools.sourceforge.net Samtools]
*[http://www.broadinstitute.org/gsa/wiki/index.php/The_Genome_Analysis_Toolkit GATK]
*[http://faculty.washington.edu/browning/beagle/beagle.html Beagle]
*[http://mathgen.stats.ox.ac.uk/impute/impute_v2.html IMPUTE2]
*[http://genome.sph.umich.edu/wiki/Thunder MaCH]
*[http://compbio.bccrc.ca/software/snvmix SNVmix]
*[http://varscan.sourceforge.net VarScan]
*[http://gmt.genome.wustl.edu/somatic-sniper Somaticsniper]
*[http://compbio.bccrc.ca/software/jointsnvmix JointSNVMix]

== References ==

{{reflist|2}}

[[Category:DNA sequencing]]

Glossary of stack theory - Revision history

en>Malcolma: added Category:Category theory; removed {{uncategorized}} using HotCat