目录

1. 一个流程

   1. 工作流程
   2. GVCF文件示例
   3. Phred-Scaled Likelyhoods

2. 主要算法

   1. Roadmap
   2. Algorithms
   3. Mark Duplicates
   4. Base Quality Score Recalibration (BQSR)
3. 背景知识
   1. VCF格式
      1. 组织结构
      2. 归一化
      3. 一些Metrics的含义
   2. Reference Genome
   3. PF
   4. Inbreeding Coefficient
   5. Spanning or Overlapping Deletions
   6. Lane-Library-Sample-Cohort
   7. Heterozygosity
   8. Phred-Scaled Quality Scores
   9. Read Groups

---

GATK Germline Variant Discovery

GATK提供了一个可运行的示例流程，用于在多个样本的队列中联合调用种系 SNP 和插入缺失。该工作流程适用于全基因组或外显子组数据。使用三个全基因组示例片段来演示 HaplotypeCaller 的 GVCF 工作流程，用于联合变异分析。使用 GenomicsDB 数据库结构，根据家族谱系进行基因型细化（genotype refinement），并评估了细化的效果。

Run on Colab

工作步骤

1. 开始于处理后的BAM文件；

2. 对单个样本使用HaplotypeCaller，设置为GVCF模式（-ERC GVCF）；获得单个样本的GVCF（推荐后缀.g.vcf），其中包含变体等位基因的每个可能的基因型的可能性，包括一个symbolic allele;

   • This will produce a GVCF file that contains likelihoods for each possible genotype for the variant alleles, including a symbolic allele.

   • VCF详解

3. 合并GVCF，官方推荐GenomicsDB；

4. 联合基因型分型；

   • 用CalculateGenotypePosteriors细化基因型调用。“ped”文件格式是指广泛使用的连锁谱系数据格式。我们可以使用 trio.ped 文件中提供的系谱信息。其次，我们可以使用人口先验；我们使用来自gnomAD的群体等位基因频率资源。

   • 更多关于Genotype Refinement

5. 比较一致性和差别；

   • 有一些不同的GATK/Picard工具来比较位点和基因型的一致性。CollectVariantCallingMetrics产生摘要和详细指标。摘要指标提供队列级的变异指标，而细节指标则对调用集中每个样本的变异指标进行分割，并增加一些字段。

   • 更多关于Metrics

6. 输出校正的VCF文件。

GVCF文件示例

##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele not already represented at this location by REF and ALT">
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##contig=<ID=20,length=63025520,assembly=GRCh37>
##source=HaplotypeCaller
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NA12878
20  10000000    .   T   <NON_REF>   .   .   END=10000008    GT:DP:GQ:MIN_DP:PL  0/0:16:48:16:0,48,497
20  10000009    .   A   <NON_REF>   .   .   END=10000009    GT:DP:GQ:MIN_DP:PL  0/0:16:32:16:0,32,485
20  10000010    .   C   <NON_REF>   .   .   END=10000017    GT:DP:GQ:MIN_DP:PL  0/0:16:48:16:0,48,472

Phred-scaled Likelyhoods

PL values are phred-scaled likelihood scores, normalized such that the most likely genotype will have a score of 0. So the approximate likelihoods are 10^(-PL/10). In this case, for PL values of 103,0,26 the likelihoods would be:

• 10^(-10.3) approximately 5.0E-11
• 10^0 approximately 1
• 10^(-2.6) approximately 0.0025

So the heterozygous case is the most likely and is indicated as such by the PL values.

Algorithms in GATK Germline Pipeline

Roadmap

(FASTQ)Preprocess(BAM) ➡ Variant Calling(via GVCF to VCF) ➡ Filteration & Recalibration(VCF for Follow-up analysis)

Algorithms

• Mark Duplicates
• Base Quality Score Recalibration (BQSR)
• HaplotypeCaller
• Joint-Call Cohort
• Refine Genotypes
• Filter and Annotate Variants

Mark Duplicates

• A single DNA fragment may be sequenced several times
• Example: duplicates can arise during library construction using PCR
• These duplicate reads are not informative and cannot be considered as evidence for or against a putative variant
• Without this step, you risk having over-representation in your sequence of areas preferentially amplified during PCR
• Duplicate reads can also result from a single amplification cluster, incorrectly detected as multiple clusters by the optical sensor of the sequencing instrument (Optical Duplicate)

Base Quality Score Recalibration (BQSR)

Quality Scores:

• Per-base estimates of error emitted by the sequencer
• Expresses the level of confidence for each base called
• Use standard Pred scores: Q20 is a general cutoff for high quality and represents 99% certainty that a base was called correctly
• 99% certainty means 1 out of 100 expected to be wrong.
• In a small dataset of 1M reads with a read length of 50, this means 50M bases total. With 99% confidence, this means 50,000 possible erroneous bases.

How Does It Work:

1. You provide GATK Base Recalibrator with a set of known variants.
2. GATK Base Recalibrator analyzes all reads looking for mismatches between the read and reference, skipping those positions which are included in the set of known variants (from step 1).
3. GATK Base Recalibrator computes statistics on the mismatches (identified in step 2) based on the reported quality score, the position in the read, the sequencing context (ex: preceding and current nucleotide).
4. Based on the statistics computed in step 3, an empirical quality score is assigned to each mismatch, overwriting the original reported quality score.

背景知识

VCF格式

VCF是一种文本文件格式（很可能以压缩的方式存储）。它包含元信息（meta-information lines）、标题行（header line）和数据（data lines）3个部分。当前版本是4.x。

例：

##fileformat=VCFv4.1
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##GATKCommandLine.HaplotypeCaller=<ID=HaplotypeCaller,Version=3.4-3-gd1ac142,Date="Mon May 18 17:36:4
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##contig=<ID=chr1,length=249250621,assembly=b37>
##reference=file:human_genome_b37.fasta
#CHROM      POS ID  REF ALF QUAL    FILTER  INFO    FORMAT  align.bam
AF086833    60  .   A   T   54  .   DP=43   GT:PL   0/1:51,0,48 
AF086833    55  .   TGAGGA  TGA 182 .   INDEL;IDV=42    GT:PL   0/1:215,0,255
AF086833    60  .   A   C,T 39.8    .   DP=126  GT:PL   1/2:99,87,47,86,0,49

组织结构

一个VCF文件包含至少8列以tab分隔的常规信息（主要是位置及质量等信息），第9列及以后表示额外信息（这部分很重要）。

• CHROM：变异发生的染色体或者contig
• POS：变异发生的基因组坐标（对缺失来讲，显示的是缺失开始的位置）
• ID：一般是dbSNP的ID（可有可无）
• REF：Forward strand（正链）上的参考等位基因
• ALT：正链上对应的改变的等位基因（可以有多个）
• QUAL：REF/ALT 变异位点存在的概率（和FASTQ的质量值、SAM的MAPQ一样，都是Phred值 -10 * log(1-p)）
• FILTER：一般都要经过适当的过滤后才能成为后续可用的variant callset（变异位点数据集）。此处可能有三种情况：
  • 没有通过过滤的变异位点
  • "PASS": 表示全部通过了过滤
  • ".": 表示对这个位点没有进行任何过滤
• INFO：以tag=value的形式给出位点注释信息；多对信息以分号分隔
• FORMAT：针对样本的注释。以冒号分隔，并且对应后面各列（每列代表一个样本）的信息

归一化

以第二条记录为例，这是一个缺失（deletion），它有多种不同的表示形式，

• TGAGGA->TGA
• GAGGA->GA
• GGA->

为了避免后续麻烦，需进行归一化：

一个归一化后的变体一定是简约的（parsimonious）和左对齐的（left aligned）。

VCF是最常用的数据格式之一，一些软件对归一化方法有额外要求，使用时需要视情况而定。

Parsimony

Left alignment

一些Metrics的含义

GT(Genotype)

GT表示样本中这个位点的基因型，其中0表示参考REF，1表示变异ALT的第一个entry，2表示ALT的第二个entry（以此类推）。

以二倍体为例，GT可能的值及其解释：

• 0/0 表示样本是纯合子，并且和参考的等位基因一样
• 0/1表示样本是杂合子，有一个参考的等位基因，一个变异的等位基因
• 1/2 样本是杂合子，两个都是变异的等位基因
• 1/1 样本是纯合子，且两个都是变异的第一个等位基因
• 2/2样本是纯合子，且两个都是变异的第二个等位基因（以此类推）

PL(Phred-Scaled Genotype Likelihood)

PL用来衡量不同基因型可能发生的概率。假设某基因型的可能性为P=1，则其phred值为-10log(P)，即为0。因此，0对应最可能的基因型。 例如，PL为51,0,48，则依次代表

• phred(P(0/0)) = 51
• phred(P(0/1)) = 0
• phred(P(1/1)) = 48

GQ(Genotype Quality)

GQ用于衡量GT的可信度，它和PL相似，也是Phred值，但不相同。GQ一般取PL中第二个小的值（除非第二小的PL大于99）。

在GATK分析中，GQ最大就限制在99，因为超过99其实是没有什么意义的，并且还多占字符。因此，如果GATK中发现PL值中第二个小的值比99还要大，软件就将GQ标为99。用GQ值就可以得到第一位和第二位之间相差多大。

其他

• AD（DepthPerAlleleBySample): unfiltered allele depth 就是有多少reads出现了某个等位基因（其中也包含了没有经过variant caller过滤的reads），但是不包括没意义的reads（就是那些统计结果不过关的，没法说服软件相信这个等位基因）。在二倍体样本中表示为：逗号分隔的两个值，前者为对应的ref基因型，后者为对应alt的基因型
• DP（Coverage，即reads覆盖度，指一些reads被过滤后的覆盖度）： filtered depth, at the sample level 只有通过variant caller软件过滤后的reads才能计算入内，但是DP也纳入了那些经过过滤但没有意义的reads（uninformative reads），这一点又和AD不同
• AC（Allele Count）：该Allele数目
• AF(Allele Frequency)：该Allele频率
• AN(Allel Number)：Allele总数目

Reference Genome

• Analysis set：分析集参考基因组具有特殊的功能，以适应序列读数对齐。这种类型的基因组参考可能与你用来浏览基因组的参考不同。见关于人类基因组参考构建的文件，以了解一个例子。
• Contig：一个contig是一个没有 "物理 "空白的连续序列（"N "个碱基的延伸在这里不被认为是空白），比如一个染色体。也可以是不完全组装中的支架，细菌基因组中的质粒，诸如此类。
• Alternate contigs、alternate scaffolds或alternate loci允许在过于复杂而不能用单一表示的区域中表示不同的单倍型。关于ALT等位基因组的目的和用途的更多讨论，请参见人类基因组参考构建的文件。
• Primary assembly：初级组装指的是(i)组装的染色体，(ii)未定位的（已知属于特定染色体，但顺序或方向未知）和(iii)未放置的（染色体未知）序列的集合。它代表一个非冗余的单倍体基因组。
• PAR：代表假性自体区。哺乳动物X和Y染色体上的PAR区允许性染色体之间的重组。因为PAR序列共同创造了一个二倍体或假自体序列区，所以在基因组组装中，X和Y染色体的序列被有意的相同。分析组基因组进一步对Y染色体PAR区域中的两个进行硬屏蔽，以便使读数只映射到X染色体PAR区域。
• Assemblies：不同的assemblies会改变基因座的坐标，并不经常发布。在人类基因组方面，Hg19和GRCh38/hg38代表两个不同的主要组合。比较来自不同基因组的数据，需要调整基因组坐标以匹配基因座，有时是不完美的。
• Patches：补丁是为一个特定的集合定期发布的区域性修复。它们的目的是在不破坏染色体坐标的情况下，改善汇编的代表性或增加信息。有两种类型的补丁，固定的和新的，代表不同类型的序列变化。
  • Fix patches：固定补丁代表的序列将在下一个主要汇编版本中取代主要汇编序列。在解释数据时，固定斑块应优先于染色体。
  • Novel patches：新的补丁代表交替的位点。在解释数据时，将新补丁视为群体序列变异。

PF

Illumina测序仪执行一个内部质量过滤程序，称为贞洁度过滤器(chastity filter)，通过该过滤器的reads被称为PF，即pass-filter。根据Illumina的说法，贞洁度被定义为最亮的碱基强度除以最亮和第二亮的碱基强度之和的比率。如果在前25个周期中，没有超过1个碱基调用的贞洁度值低于0.6，那么reads簇(cluster)就会通过过滤。这个过滤过程将最不可靠的簇从图像分析结果中去除。

Inbreeding Coefficient

近交系数是一个变异位点的过剩杂合度。它可以作为不良映射（poor mapping）的代理（近交系数高的位点通常是基因组中映射不良的位置，处于该区域的读数与该区域不匹配，因为它们属于其他地方）。为了正确计算这个注释，至少需要10个样本（最好是更多）。

近交系数依赖于 Hardy-Weinberg 原理背后的数学。

完整原文

Spanning or Overlapping Deletions

我们使用术语——跨越删除（spanning deletion）或重叠删除（overlapping deletion）——来指代跨越感兴趣位置的缺失。

跨越缺失的存在影响我们如何在它跨越的任何位点代表携带缺失的那些样本的基因型，无论是杂合还是纯合变异形式。 VCF v4.3 规范的第 8 页第 5 项保留了 * 等位基因以引用重叠删除。这不要与用于表示象征性替代等位基因(symbolic alternate alleles)的括号星号 <*> 混淆。

Lane-Library-Sample-Cohort

四个主要的测序数据组织单位：

• Lane: 测序的基本机器单位。该通道反映了高通量测序机的基本独立运行。对于Illumina机器，这就是物理测序通道。
• Library: 一个DNA制备单位，在某些时候被物理地汇集在一起。多个通道可以从同一文库的等分试样中运行。DNA库及其制备是被测序的自然单位。例如，如果文库的复杂性有限，那么许多序列会被重复，并会导致各泳道的高重复率。如果使用RNAseq，文库的准备过程包括反转录成cDNA。
• Sample: 来自单一个体的生物样本。多个具有不同性质的文库可以从原始样本的DNA来源中构建。在我们的文件中，我们将样本视为独立的个体，我们试图确定其基因组序列。请注意，从这个角度来看，尽管肿瘤/正常样本来自同一个人，但它们是不同的。
• Cohort: 正在一起分析的样本的集合。这个组织单位是最主观的，非常具体地取决于测序项目的设计目标。对于像1000基因组这样的群体发现项目，分析队列是每个群体中约100个个体。对于有许多深度测序样本的外显子组项目（如有800个EOMI样本的ESP），我们将整套样本分成~50个个体的队列，进行多样本分析。

Heterozygosity

杂合度。在GATK基因分型中，我们使用 "预期杂合度 "值来计算一个基因座为非参考的先验概率。鉴于预期杂合度hets，我们计算N个样本在一个位点上是同源参考的概率为1-sum_i_2N（hets / i）。为人类提供的默认值是hets=1e-3；0.001的值意味着从生物群体中随机选择的两条染色体将以1000bp的比率相互差异。在这种情况下，hets类似于群体遗传学中的参数theta。如果需要，hets参数值可以被修改。

请注意，这个数值与任何给定样本具有杂合基因型的可能性无关，在GATK中，这纯粹是由观察数据P(D | AB)在可能存在AB杂合基因型的模式下的概率决定的。这种AB基因型的后验概率将使用hets先验，但是GATK在确定一个位点是多态的概率时只使用这个后验概率。因此，改变hets参数只是增加了一个位点在所有样本中被称为非参考的机会，但实际上根本没有改变输出基因型的可能性，因为这些不是后验概率。改变GATK是否考虑杂合基因型的可能性的一个数量是倍性，它描述了物种中每个个体携带多少个染色体副本。

Phred-Scaled Quality Scores

计算公式： Q = -10 log E

Phred Quality Score	Error	Accuracy (1 - Error)
10	1/10 = 10%	90%
20	1/100 = 1%	99%
30	1/1000 = 0.1%	99.9%
40	1/10000 = 0.01%	99.99%
50	1/100000 = 0.001%	99.999%
60	1/1000000 = 0.0001%	99.9999%

Read Groups

在实践中，一个read group指由测序仪一次运行产生的一组reads。

读取组在SAM/BAM/CRAM文件中是由官方SAM规范中定义的一些标签来识别的。这些标签，如果分配得当，不仅可以让我们区分样本，还可以区分与文物有关的各种技术特征。有了这些信息在手，我们就可以在重复标记和基础重新校准的步骤中减轻这些人工制品的影响。

要查看 BAM 文件的读取组信息，请使用以下命令。

samtools view -H sample.bam | grep '^@RG'

输出：

@RG ID:H0164.2  PL:illumina PU:H0164ALXX140820.2    LB:Solexa-272222    PI:0    DT:2014-08-20T00:00:00-0400 SM:NA12878  CN:BI

GATK所要求的read group的字段的含义：

• ID： Read group identifier
• PU： Platform Unit
• SM： Sample
• PL： Platform/technology used to produce the read
• LB： DNA preparation library identifier