NGS 데이터 분석 - GATK 변이 발굴 (Germline)

NGS 데이터 분석 - GATK 전체 파이프라인 및 설치

NGS 데이터 분석 - GATK 데이터 전처리

 

[참고]

- Germline short variant discovery (SNPs + Indels) (link)

Cohort Data

Germline Short Variant Discovery (1) Cohort data

Call variants per-sample

HaplotypeCaller은 haplotypes의 local de-novo assembly를 통해 SNP과 indel의 calling을 동시에 실시할 수 있다. 즉, 변이가 있다고 생각되는 순간 mapping 정보를 버리고 해당 부위의 reads를 다시 assembly 하기 때문에 기존의 알고리즘(UnifiedGenotyper)보다 더욱 정확하다.

GVCF 모드는 대규모 variant calling에서 사용되며, HaplotypeCaller은 샘플별로 분석하여 중간 GVCF 파일을 생성한다. 이는 추후 여러 샘플을 활용한 joint genotyping에 활용되며, 대규모 코호트에서 매우 효율적으로 샘플수를 점진적으로 증가시킬 수 있다. 실제로는 이 과정은 샘플별로 데이터 전처리와 함께 진행되어 초기의 uBAM 파일에서 GVCF파일까지 생성하게 된다.

- HaplotypeCaller (in GVCF mode)

 

Consolidate GVCFs

여러 샘플에 걸친 GVCF 파일들을 병합하는 과정으로 이를 통해 joint genotyping의 규모와 속도를 늘릴 수 있다. 주의해야할 점은 이 과정에서는 GVCF를 단순히 합치는 것이지 이를 joint genotyping이라고 혼동해서는 안 된다. 버전 4 이전에서는 CombineGVCFs를 통해 합쳐졌고 버전 4에도 포함되었지만 이는 단순히 유산으로 남기기 위함이며, GenomicsDBImport의 성능이 훨씬 우수하며 GVCF 파일 대신 datastore을 생성한다.

- GenomicsDBImport

Joint-Call Cohort

GenotypeGVCFs는 샘플별 GVCFs을 모두 모아서 joint-called된 SNP과 indel을 생성한다. 코호트 단위의 분석을 통해서는 어려운 부위에서도 세밀하게 변이를 감지를 할 수 있고, squared-off 유전형 행렬을 생성하여 모든 샘플에 걸친 관심 있는 부위의 정보를 제공한다. 이 과정은 매우 빠르게 실행할 수 있으며(N+1 문제), 샘플이 추가될 때마다 언제든 다시 실행이 가능하다.

- GenotypeGVCFs

Filter Variants by Variant (Quality Score) Recalibration

GATK는 기본적으로 실제 변이를 놓치는 일이 없도록 하는 것을 중요시하지만, raw callset을 필터링하여 false positives를 줄이는 것이 필요하다. 이를 variant quality score recalibration (VQSR)이라고 하며, 머신러닝 기법을 활용하여 변이의 annotation 정보를 통해 실제 변이 여부를 확인하고, 각 변이에 QUAL보다 훨씬 정확한 VQSLOD score를 부여한다. 먼저 학습 변이를 통해 모델을 생성하고, 데이터에 모델을 적용하여 보정이 잘 된 확률을 각 변이에 부여한다. 그 다음 이 점수를 통해 필터링을 실시하여 원하는 수준에 도달한 callset을 찾는다.

이때, 알고리즘은 학습과정에서 매우 퀄리티가 높은 이미 알려진 변이들을 필요로 하며, 꽤 많은 데이터를 통해 변이의 좋고 나쁨을 학습해야 한다. 따라서 규모가 작은 데이터셋이나 많은 경우 학습이 어렵거나 불가능한 경우가 있다 (targeted sequencing data나 RNAseq, non-model organisms). 만약 어떤 이유든 VQSR을 실행할 수 없다면, 특정 annotation에 대해 일괄적인 thresholds를 설정해서 모든 변이에 동일하게 적용하는 hard-filtering 방법을 실시해야 한다. 현재 신경망 기반의 방법을 개발 중에 있으며 궁극적으로는 VQSR을 대체할 것이라고 한다.

- VariantRecalibrator

- ApplyVQSR

Single Sample Data

Germline Short Variant Discovery (2) Single Sample Data

단일 샘플의 경우 HaplotypeCaller의 기본 단일샘플 모드를 통해 BAM에서 VCF파일을 생성한다. 변이 필터링의 경우, CNNScoreVariants를 사용하여 각 변이에 모델의 quality 예측값을 부여하고, FIlterVariantTranches을 통해 목적에 맞게 SNP과 indel 민감도를 설정하여 제거한다.

- HaplotypeCaller (in single-sample mode)

- CNNScoreVariants

- FIlterVariantTranches

- Funcotator

 

HaplotypeCaller

Haplotype local re-assembly를 통한
Germline short variant discovery

- HaplotypeCaller (link)

- Input : BAM

- Output : VCF/GVCF

 

[Overview]

The previous step produced a table of per-read allele likelihoods for each candidate variant site under consideration. Now, all that remains to do is to evaluate those likelihoods in aggregate to determine what is the most likely genotype of the sample at each site. This is done by applying Bayes' theorem to calculate the likelihoods of each possible genotype, and selecting the most likely. This produces a genotype call as well as the calculation of various metrics that will be annotated in the output VCF if a variant call is emitted.

 

[Preliminary assumptions / limitations]

- Quality : Keep in mind that we are trying to infer the genotype of each sample given the observed sequence data, so the degree of confidence we can have in a genotype depends on both the quality and the quantity of the available data. By definition, low coverage and low quality will both lead to lower confidence calls. The GATK only uses reads that satisfy certain mapping quality thresholds, and only uses “good” bases that satisfy certain base quality thresholds (see documentation for default values).

- Ploidy : Both the HaplotypeCaller and GenotypeGVCFs assume that the organism of study is diploid by default, but the desired ploidy can be set using the -ploidy argument. The ploidy is taken into account in the mathematical development of the Bayesian calculation using a generalized form of the genotyping algorithm that can handle ploidies other than 2. Note that using ploidy for pooled experiments is subject to some practical limitations due to the number of possible combinations resulting from the interaction between ploidy and the number of alternate alleles that are considered. There are some arguments that aim to mitigate those limitations but they are not fully documented yet.

- Paired end reads : Reads that are mates in the same pair are not handled together in the reassembly, but if they overlap, there is some special handling to ensure they are not counted as independent observations.

- Single-sample vs multi-sample : We apply different genotyping models when genotyping a single sample as opposed to multiple samples together (as done by HaplotypeCaller on multiple inputs or GenotypeGVCFs on multiple GVCFs). The multi-sample case is not currently documented for the public but is an extension of previous work by Heng Li and others.

 

[Phase]

1. Define active regions. The program determines which regions of the genome it needs to operate on, based on the presence of significant evidence for variation.

2. Determine haplotypes by re-assembly of the active region. For each ActiveRegion, the program builds a De Bruijn-like graph to reassemble the ActiveRegion and identifies what are the possible haplotypes present in the data. The program then realigns each haplotype against the reference haplotype using the Smith-Waterman algorithm in order to identify potentially variant sites.

3. Determine likelihoods of the haplotypes given the read data. For each ActiveRegion, the program performs a pairwise alignment of each read against each haplotype using the PairHMM algorithm. This produces a matrix of likelihoods of haplotypes given the read data. These likelihoods are then marginalized to obtain the likelihoods of alleles per read for each potentially variant site.

4. Assign sample genotypes. For each potentially variant site, the program applies Bayes’ rule, using the likelihoods of alleles given the read data to calculate the posterior likelihoods of each genotype per sample given the read data observed for that sample. The most likely genotype is then assigned to the sample.

1. Define active regions

In this first step, the program traverses the sequencing data to identify regions of the genomes in which the samples being analyzed show substantial evidence of variation relative to the reference. The resulting areas are defined as “active regions”, and will be passed on to the next step. Areas that do not show any variation beyond the expected levels of background noise will be skipped in the next step. This aims to accelerate the analysis by not wasting time performing reassembly on regions that are identical to the reference anyway.

To define these active regions, the program operates in three phases. First, it computes an activity score for each individual genome position, yielding the raw activity profile, which is a wave function of activity per position. Then, it applies a smoothing algorithm to the raw profile, which is essentially a sort of averaging process, to yield the actual activity profile. Finally, it identifies local maxima where the activity profile curve rises above the preset activity threshold, and defines appropriate intervals to encompass the active profile within the preset size constraints.

Once this process is complete, the program applies a few post-processing steps to finalize the the active regions (see detailed doc above). The final output of this process is a list of intervals corresponding to the active regions which will be processed in the next step.

2. Determine haplotypes by local assembly of the active region.

The goal of this step is to reconstruct the possible sequences of the real physical segments of DNA present in the original sample organism. To do this, the program goes through each active region and uses the input reads that mapped to that region to construct complete sequences covering its entire length, which are called haplotypes. This process will typically generate several different possible haplotypes for each active region due to:

• real diversity on polyploid (including CNV) or multi-sample data
• possible allele combinations between variant sites that are not totally linked within the active region
• sequencing and mapping errors

In order to generate a list of possible haplotypes, the program first builds an assembly graph for the active region using the reference sequence as a template. Then, it takes each read in turn and attempts to match it to a segment of the graph. Whenever portions of a read do not match the local graph, the program adds new nodes to the graph to account for the mismatches. After this process has been repeated with many reads, it typically yields a complex graph with many possible paths. However, because the program keeps track of how many reads support each path segment, we can select only the most likely (well-supported) paths. These likely paths are then used to build the haplotype sequences which will be used for scoring and genotyping in the next step.

Once the haplotypes have been determined, each one is realigned against the original reference sequence in order to identify potentially variant sites. This produces the set of sites that will be processed in the next step. A subset of these sites will eventually be emitted as variant calls to the output VCF.

3. Evaluating the evidence for haplotypes and variant alleles

Now that we have all these candidate haplotypes, we need to evaluate how much evidence there is in the data to support each one of them. So the program takes each individual read and aligns it against each haplotype in turn (including the reference haplotype) using the PairHMM algorithm, which takes into account the information we have about the quality of the data (i.e. the base quality scores and indel quality scores). This outputs a score for each read-haplotype pairing, expressing the likelihood of observing that read given that haplotype.

Those scores are then used to calculate out how much evidence there is for individual alleles at the candidate sites that were identified in the previous step. The process is called marginalization over alleles and produces the actual numbers that will finally be used to assign a genotype to the sample in the next step.

[Calculating genotype likelihoods using Bayes' Theorem]

We use the approach described in Li 2011 to calculate the posterior probabilities of non-reference alleles (Methods 2.3.5 and 2.3.6) extended to handle multi-allelic variation.

The basic formula we use for all types of variation under consideration (SNPs, insertions and deletions) is:

$$ P(G|D) = \frac{ P(G) P(D|G) }{ \sum_{i} P(G_i) P(D|G_i) } $$

If that is meaningless to you, please don't freak out -- we're going to break it down and go through all the components one by one. First of all, the term on the left:

$$ P(G|D) $$

is the quantity we are trying to calculate for each possible genotype: the conditional probability of the genotype G given the observed data D.

Now let's break down the term on the right:

$$ \frac{ P(G) P(D|G) }{ \sum_{i} P(G_i) P(D|G_i) } $$

We can ignore the denominator (bottom of the fraction) because it ends up being the same for all the genotypes, and the point of calculating this likelihood is to determine the most likely genotype. The important part is the numerator (top of the fraction):

$$ P(G) P(D|G) $$

which is composed of two things: the prior probability of the genotype and the conditional probability of the data given the genotype.

The first one is the easiest to understand. The prior probability of the genotype G:

$$ P(G) $$

represents how probably we expect to see this genotype based on previous observations, studies of the population, and so on. By default, the GATK tools use a flat prior (always the same value) but you can input your own set of priors if you have information about the frequency of certain genotypes in the population you're studying.

The second one is a little trickier to understand if you're not familiar with Bayesian statistics. It is called the conditional probability of the data given the genotype, but what does that mean? Assuming that the genotype G is the true genotype,

$$ P(D|G) $$

is the probability of observing the sequence data that we have in hand. That is, how likely would we be to pull out a read with a particular sequence from an individual that has this particular genotype? We don't have that number yet, so this requires a little more calculation, using the following formula:

$$ P(D|G) = \prod{j} \left( \frac{P(D_j | H_1)}{2} + \frac{P(D_j | H_2)}{2} \right) $$

You'll notice that this is where the diploid assumption comes into play, since here we decomposed the genotype G into:

$$ G = H_1H_2 $$

which allows for exactly two possible haplotypes. In future versions we'll have a generalized form of this that will allow for any number of haplotypes.

Now, back to our calculation, what's left to figure out is this:

$$ P(D_j|H_n) $$

which as it turns out is the conditional probability of the data given a particular haplotype (or specifically, a particular allele), aggregated over all supporting reads. Conveniently, that is exactly what we calculated in Step 3 of the HaplotypeCaller process, when we used the PairHMM to produce the likelihoods of each read against each haplotype, and then marginalized them to find the likelihoods of each read for each allele under consideration. So all we have to do at this point is plug the values from that table into the equation above, and we can work our way back up to obtain:

$$ P(G|D) $$

for the genotype G.

4. Assigning per-sample genotypes

The previous step produced a table of per-read allele likelihoods for each candidate variant site under consideration. Now, all that remains to do is to evaluate those likelihoods in aggregate to determine what is the most likely genotype of the sample at each site. This is done by applying Bayes' theorem to calculate the likelihoods of each possible genotype, and selecting the most likely. This produces a genotype call as well as the calculation of various metrics that will be annotated in the output VCF if a variant call is emitted.

[Selecting a genotype and emitting the call record]

We go through the process of calculating a likelihood for each possible genotype based on the alleles that were observed at the site, considering every possible combination of alleles. For example, if we see an A and a T at a site, the possible genotypes are AA, AT and TT, and we end up with 3 corresponding probabilities. We pick the largest one, which corresponds to the most likely genotype, and assign that to the sample.

Note that depending on the variant calling options specified in the command-line, we may only emit records for actual variant sites (where at least one sample has a genotype other than homozygous-reference) or we may also emit records for reference sites. The latter is discussed in the reference confidence model documentation.

Assuming that we have a non-ref genotype, all that remains is to calculate the various site-level and genotype-level metrics that will be emitted as annotations in the variant record, including QUAL as well as PL and GQ. For more information on how the other variant context metrics are calculated, please see the corresponding variant annotations documentation.

 

[참고]

- HaplotypeCaller in a nutshell (link) 

- Assigning per-sample genotypes (HaplotypeCaller) (link)

• We recommend using a list of intervals to speed up analysis. See this document for details.

GVCF mode

This document describes the reference confidence model applied by HaplotypeCaller to generate a per-sample GVCF, invoked by -ERC GVCF or -ERC BP_RESOLUTION.

As explained here, HaplotypeCaller works by assembling the reads to create potential haplotypes, realigning the reads to their most likely haplotypes, and then projecting these reads back onto the reference sequence via their haplotypes to compute alignments of the reads to the reference. At that point, we can calculate the likelihoods of each possible genotype and emit variant calls.

What that article does not explain is how HaplotypeCaller additionally estimates the chance that some (unknown) non-reference allele is segregating at this position by examining the realigned reads that span the reference base. At this base we perform two calculations:

• Estimate the confidence that no SNP exists at the site by contrasting all reads with the REF base vs. all reads with any non-reference base.
• Estimate the confidence that no indel of size X (determined by command line parameter) could exist at this site by calculating the number of reads that provide evidence against such an indel, and from this value estimate the chance that we would not have seen the allele confidently.

Based on this, we emit the genotype likelihoods (PL) and compute the GQ (from the PLs) for the least confidence of these two models. We use a symbolic ALT allele, NON_REF, to hold the likelihood that the site is not homozygous reference, as well as allele-specific AD and PL field values.

We do this at all sites in the territory covered by the analysis, including homozygous-reference sites, both inside and outside the ActiveRegions determined by HaplotypeCaller.

1. Variant calling

Run the HaplotypeCaller on each sample's BAM file(s) (if a sample's data is spread over more than one BAM, then pass them all in together) to create single-sample gVCFs, with the option --emitRefConfidence GVCF, and using the .g.vcf extension for the output file.

Note that versions older than 3.4 require passing the options --variant_index_type LINEAR --variant_index_parameter 128000 to set the correct index strategy for the output gVCF.

2. Data aggregation step

A new tool called GenomicsDBImport is necessary to aggregate the GVCF files and feed in one GVCF to GenotypeGVCFs. You can read more about it here. You can also run CombineGVCFs if you are not able to use GenomicsDBImport.

3. Joint genotyping

Take the outputs from step 2 (or step 1 if dealing with fewer samples) and run GenotypeGVCFs on all of them together to create the raw SNP and indel VCFs that are usually emitted by the callers.

4. Variant recalibration

Finally, resume the classic GATK Best Practices workflow by running VQSR on these "regular" VCFs according to our usual recommendations.

That's it! Fairly simple in practice, but we predict this is going to have a huge impact in how people perform variant discovery in large cohorts. We certainly hope it helps people deal with the challenges posed by ever-growing datasets.

 

[참고]

- Calling variants on cohorts of samples using the HaplotypeCaller in GVCF mode (link)

- HaplotypeCaller Reference Confidence Model (GVCF mode) (link)