全基因组测序(WGS)流程及实践

全基因组测序(WGS)流程及实践

在开始正文之前，声明一下，此篇教程参考黄树嘉老师在知乎上发的专栏《从零开始完整学习全基因组测序数据分析》。本文的目的是记录笔记以便随时浏览查阅。还有一点就是，流程的具体形式是次要的，重要的是怎么利用这个工具来解决生物学问题。用哲学的话来说，如何对数据进行分析和解读是根本之“道”，而具体的某个软件、某个参数只是“术”，不要依赖或者深陷于某个工具，而是要认清自己要解决的问题是什么，所需要的结果是什么，然后再巧妙地选择合适的技术去实现。

⚠️：本文代码皆在MacOS系统终端执行，并以E.coli K12为例讲解WGS流程。

目录：
0.准备阶段
1.原始数据质控
2.read比对，排序和去除重复序列
3.Indel区域重（“重新”的“重”）比对
4.碱基质量值重校正
5.变异检测
6.变异结果质控和过滤

0.准备阶段

下载并安装以下软件

BWA, Samtools, GATK4.0, sratoolkit, tabix, fastp, Trimmomatic。

下载E.coli K12的参考基因组序列

# 下载E.coli K12的参考基因组序列
curl -O ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz

# 解压并重命名为E.coli_K12_MG1655.fa
bgzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa

# 用samtools创建索引，生成E.coli_K12_MG1655.fa.fai文件
# 通过索引，可以方便地获取fasta文件中任意位置的序列。
# e.g. samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200
samtools faidx E.coli_K12_MG1655.fa

下载E.coli K12的测序数据

需要用到NCBI的官方工具包sratoolkit，直接下载.fastq格式文件，当然也可以在NCBI上下载好.sra文件并用sratoolkit转换成我们所需的.fastq格式文件。这个数据来自Illumina MiSeq测序平台，read长度是300bp，测序类型是双末端测序(Pair-End)。

「注」: 从被测DNA序列两端各测序一次的模式，这就被称为双末端测序(Pair-End Sequencing，简称PE测序)。

# 方式一
# 下载.sra格式
curl -O ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra
# 转换为.fastq格式
fastq-dump --split-files SRR1770413.sra

# 方式二：直接下载.fastq格式
fastq-dump --split-files SRR1770413

# 用bgzip（不推荐gzip）将其压缩为.gz文件，节省空间
bgzip -f SRR1770413_1.fastq
bgzip -f SRR1770413_2.fastq

来看一下当前E.coli K12相关数据目录结构：

./wgs_practice/  
├── bin  
├── input  
│   ├── E.coli  
│       ├── reference  
│       │   ├── E.coli_K12_MG1655.fa  
│       │   ├── E.coli_K12_MG1655.fa.fai  
│       │   └── work.log.sh  
│       └── rawdata  
│           ├── SRR1770413_1.fastq.gz  
│           ├── SRR1770413_2.fastq.gz  
│           └── work.log.sh     
└── output

✳ 小扩展：FASTA和FASTQ
准备好了这几个文件之后，你也许心里会有困惑，fasta文件和fastq文件有什么区别呢？这两个都是存储核苷酸序列信息(DNA/RNA)或者蛋白质的氨基酸序列信息(Amino Acid sequence, aa)最常使用的两种文本文件。

FASTA
1) 来源：由William. R. Pearson 和 David. J. Lipman在1988年所编写的“FASTA”的比对软1) 来源：由William. R. Pearson 和 David. J. Lipman在1988年所编写的“FASTA”的比对软件，目的是用于生物序列数据的处理。自那之后，生物科学家们把FASTA作为这种存储有！顺！序！的序列数据的文件后缀。这包括我们常用的参考基因组序列、蛋白质序列、编码DNA序列(coding DNA sequence，简称CDS）、转录本序列等文件都是如此，文件后缀除了.fasta之外，也常用.fa或者.fa.gz(gz压缩)。

2) 构成：序列头信息(有时包括一些其它的描述信息)和具体的序列数据。头信息独占一行，以大于号(>)开头作为识别标记，其中除了记录该条序列的名字之外，有时候还会接上其它的信息。注意的点：①头信息格式没有严格地限制。所以有业内不成文的规则：用一个空格把头信息分为两个部分：第一部分是序列名字，它和大于号(>)紧接在一起；第二部分是注释信息，这个可以没有，就看具体需要。e.g. >gene_00284728 length=231;type=dna ②FASTA由于是文本文件，它里面的内容是否有重复是无法自检的，在使用之前需要我们进行额外的检查。只需检查序列名字是否有重复即可。紧接的下一行是具体的序列内容，直到另一行碰到另一个大于号(>)开头的新序列或者文件末尾。

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC
AGATAAAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACCATTACCACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGG

FASTQ
1) 来源：产生自测序仪的原始测序数据，文件大小依照不同的测序量(或测序深度)而有很大差异，文件后缀通常都是.fastq，.fq或者.fq.gz(gz压缩)。

2) 构成：每四行成为一个独立的单元，我们称之为read。

第一行：以‘@’开头，是这一条read的名字，这个字符串是根据测序时的状态信息转换过来的，中间不会有空格，它是每一条read的唯一标识符。

第二行：测序read的序列，N代表的是测序时那些无法被识别出来的碱基。

第三行：以“+”开始，可以储存一些附加信息，一般是空的。

第四行：测序read的质量值，它描述的是每个测序碱基的可靠程度，用ASCII码表示，越大说明测序的质量越好。公式：Q = -10log(p_error)。即，质量值是测序错误率的对数(10为底数)乘以-10(并取整)。比如，如果该碱基的测序错误率是0.01，那么质量值就是20(俗称Q20)。为什么要用ASCII码来代表，直接用数字不行吗？行！但很难看。并非所有的ASCII码都是可见的字符，为了能够让碱基的质量值表达出来，必须避开所有这些不可见字符。最简单的做法就是加上一个固定的整数(33)！我们把加33的的质量值体系称之为Phred33。

@SRR1770413.1 1 length=301
CACCCGGCATCAGGTGCGGTACTTTTGCGCCTCCCAGCCGGACCGGCCCTGCGGCGTAATACGCGGCACTTTCCCCCCTCCAGCCGGTCCGCCTTTCCCGCTCCTCCTGCCTGTTCGCGCCCAGCCTCACATCCCTCGCTGCCTGCGTATCCAGCTCACTCTCCCTGGTTGCCGCCTACATCCCCTTCCCACTGGCGCTCTGCCGCTGCGTGTTCTTCCGCCTGTCCAGCATCTTCATCTGCTCCCTCCCGCTGTTCCACCCCTTTGCACCCCCCCTCTGCCCCTCCTGCTCGCCAGCCCC
+
CCCCCECFCEFC@8F8C77B7BFEFD,C+@@@BCB##########################################################################################################################################################################################################################################################################
@SRR1770413.2 2 length=301
TGTCTCAATCAATAAAACATCGTCAATTTCCGTGACACCGTCACTTTTCTTTTCTACTTTCATTGGGCTAAAACTGTCGCTTACCGCTTCCAGGCCAACAACTGCCTCACGTTCCCCCCGCCACTCGTACCCGCCCAGCCCACTTTTCTTTCCCAACCGTCCCCCAATCCCCCGTTCCTTTTTCCCCCGCCCAGGTAACAAACCCCCCTCCTCCCACCAAGCCTTAAACACCCAACCGGTTACCGCCACATTGCCTTCCTGAAGCTCGCAACCCCCCCCCTCTGCACCCCCGTCCCCCCTC
+
A@CCBF9FGGD9FF9C9E<EF7CE,,CEEFE,C,,:,;C+;B,8C@CCFEEEECC,<C@@@,<C,,,6;,,,,;@,<C,C@C,CC,678B@,,,:@,,9,,::,:B,?,A,:9AA##########################################################################################################################################################################################

1.原始测序数据质控

质控的意义：以illumina为首的二代测序基本都是运用边合成边测序的技术。碱基的合成依靠的是化学反应，这使得碱基链可以不断地从5'端一直往3'端合成并延伸下去。但在这个合成的过程中随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，这就会带来一个问题——越到后面碱基合成的错误率就会越高。而测序数据的质量好坏会影响我们的下游分析。

在分析测序之前先搞清楚如下两个地方：①原始数据是通过哪种测序平台产生的，它们的错误率分布是怎么样的，是否有一定的偏向性和局限性，是否会显著受GC含量的影响等。②评估它们有可能影响哪些方面的分析。

一般我们可以从如下几个方面来分析：

read各个位置的碱基质量值分布
好的结果是碱基质量值都基本都在大于30，而且波动很小。
碱基的总体质量值分布
对于二代测序，最好是达到Q20的碱基要在95%以上(最差不低于90%)，Q30要求大于85%(最差也不要低于80%)。
read各个位置上碱基分布比例，目的是为了分析碱基的分离程度
最好平均在1%以内。
GC含量分布
对于人类来说，我们基因组的GC含量一般在40%左右。
read各位置的N含量
N在测序数据中一般是不应该出现的。
read是否还包含测序的接头序列
在测序之前需要构建测序文库，测序接头就是在这个时候加上的，其目的一方面是为了能够结合到flowcell上，另一方面是当有多个样本同时测序的时候能够利用接头信息进行区分。一般的WGS测序是不会测到这些接头序列的。
read重复率，这个是实验的扩增过程所引入的

关于如何认识fq数据的事情就说到这，接下来谈谈如何去除测序接头和低质量序列。fastp软件可以很轻松地实现上述功能。除了质控报告，fastq还兼具过滤，接头处理，滑窗质量剪裁，PE数据的碱基校正，全局剪裁，polyG剪裁，分子标签UMI处理，输出文件切分等功能。

使用fastp对原始数据进行质控

fastp \
-i SRR1770413_1.fastq.gz \
-o SRR1770413_1.clean.fastq.gz \
-I SRR1770413_2.fastq.gz \
-O SRR1770413_2.clean.fastq.gz \
-z 4 \
-f 5 -t 5 -F 5 -T 5 \
-5 -W 5 -M 20 \
-Q \
-l 50 \
-c \
-w 4

# 参数说明
# 1. I/O options   即输入输出文件设置
## -i, --in1    输入read1文件名  
## -o, --out1   输出read1文件名  
## -I, --in2    输入read2文件名  
## -O, --out2   输出read2文件名  
## -z, --compression    输出压缩格式。给定一个数字1-9，调整压缩比率和效率的平衡  

# 2. adapter trimming options   过滤序列接头参数设置  
## -A, --disable_adapter_trimming    默认已开启，设置该参数则关闭  

# 3. global trimming options   剪除序列起始和末端的低质量碱基数量参数  
## -f,--trim_front1  
## -t, --trim_tail1  
## -F, --trim_front2  
## -T, --trim_tail2  

# 4. per read cutting by quality options   划窗裁剪  
## 具体的原理就是通过滑动一定长度的窗口，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除。  
### 不是挖掉read中的这部分低质量序列，像切菜一样，直接从低质量区域开始把这条read后面的所有其它碱基全！部！剁！掉！否则就是在人为改变实际的基因组序列情况。  
## -5, --cut_by_quality5  
## -3, --cut_by_quality3  
## -W, --cut_window_size    设置窗口大小  
## -M, --cut_mean_quality    设置窗口碱基平均质量阈值  

# 5. quality filtering options   根据碱基质量来过滤序列  
## -Q, --disable_quality_filtering    控制是否去除低质量，默认自动去除，设置-Q关闭  
## -q, --qualified_quality_phred    设置低质量的标准，默认Q15  
## -u, --unqualified_percent_limit    低质量碱基所占百分比，默认40代表40%  

# 6. length filtering options   根据序列长度来过滤序列  
## -l, --length_required    设规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉  

# 7. base correction by overlap analysis options   通过overlap来校正碱基  
## -c, --correction    是对overlap的区域进行纠错，只适用于pairend reads  

# 8. threading options   设置线程数  
## -w, --thread 设置线程数

或使用Trimmomatic对原始数据进行质控

java -jar /path/trimmomatic-0.39.jar PE \
-threads 4 \
-phred33 \
SRR1770413_1.fastq.gz SRR1770413_2.fastq.gz \
SRR1770413_1.clean.fastq.gz SRR1770413_1.unpaired.fastq.gz SRR1770413_2.clean.fastq.gz SRR1770413_2.unpaired.fastq.gz \
ILLUMINACLIP:/path/TruSeq3-PE.fa:2:30:10:8:True \
SLIDINGWINDOW:5:20 \
LEADING:5 \
TRAILING:5 \
MINLEN:50

# 参数说明
# 1.java -jar trimmomatic-0.39.jar PE \
## PE:双端测序，单端数据用"SE"  

# 2. [-threads <threads>] \
## 设置线程数  

# 3. [-phred33|-phred64] \
## 设置碱基的质量格式(默认-phred64，自v0.32版本之后可自动识别是phred33还是phred64)  

# 4. [-trimlog <trimLogFile>] \
## 设置trimmommatic工具处理的日志文件为'trim.log'，每两行为一对reads信息  
## 生成的日志文件'trim.log'包含5列信息：
###  1) read名称
###  2) 切除后剩余的read长度
###  3) 切除后第一个碱基所在位置(0-base)，也就是起始位置切除了几个碱基
###  4) 切除后最后一个碱基所在位置
###  5) read末尾切除了几个碱基
## 由于该生成文件较大，如后续步骤不需该文件信息，可考虑不设置

# 5. [<inputFile1> <inputFile2>] \
## 输入文件  

# 6. [<outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] \
## 输出文件  

# 7. [ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>] \
## 作用:过滤reads中的测序接头和引物序列  
## fastaWithAdaptersEtc：指定包含接头和引物序列(所有被视为污染的序列)的fasta文件路径  
## seedMismatches：比对时接头序列时所允许的最大错配数  
## palindrome clip threshold：PE的两条read同时和PE的adapter序列比对，匹配度加起来超过这个阈值才进行切除  
## simple clip threshold：任意adapter序列与某read比对最低分阈值，大于这个阈值，就进行切除  
## minAdapterLength：指定palindrome模式下可以切除的接头序列最短长度，默认为8  
## keepBothReads：只作用于palindrome模式，设置":True"保留反向链，/
### 这里是说去除接头序列后，由于正反链包含相同的序列信息(尽管序列是反向互补的)，默认情况(":false")下会去除反向链。

# 8. [SLIDINGWINDOW:<windowSize>:<requiredQuality>] \
## 作用:从reads的5'端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗
## widowSize:设置窗口大小  
## requiredQuality：设置窗口碱基平均质量阈值  

# 9. [LEADING:<quality>] \
## quality:规定read开头的碱基是否要被切除的质量阈值  

# 10. [TRAILING:<quality>] \
## quality:规定read末尾的碱基是否要被切除的质量阈值  

# 11. [MINLEN:<length>]
## length:设规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉
time trimm PE -phred33 \
SRR1770413_1.fastq.gz SRR1770413_2.fastq.gz \
../cleandata/SRR1770413_1.clean.fastq.gz ../cleandata/SRR1770413_1.unpaired.fastq.gz \
../cleandata/SRR1770413_2.clean.fastq.gz ../cleandata/SRR1770413_2.unpaired.fastq.gz \
ILLUMINACLIP:/Users/xtingmao/Applications/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10:8:True \
SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50 \
&& echo "** fq QC done **"

2. read比对，排序和去除重复序列

序列比对

NGS测序下来的短序列(read)存储于FASTQ文件里面，FASTQ文件中紧挨着的两条read之间没有任何位置关系，它们都是随机来自于原本基因组中某个位置的短序列而已。因此，我们需要先把这一大堆的短序列捋顺，一个个去跟该物种的参考基因组比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。

首先，我们需要为参考基因组的构建索引——这其实是在为参考序列进行Burrows Wheeler变换（wiki: 块排序压缩），以便能够在序列比对的时候进行快速的搜索和定位。

# 为参考序列构建BWA比对所需的FM-index(比对索引)
bwa index E.coli_K12_MG1655.fa

完成之后，你会看到类似如下几个以E.coli_K12_MG1655.fasta为前缀的文件：

E.coli_K12_MG1655.fa.amb
E.coli_K12_MG1655.fa.ann
E.coli_K12_MG1655.fa.bwt
E.coli_K12_MG1655.fa.pac
E.coli_K12_MG1655.fa.sa

这些就是在比对时真正需要被用到的文件。这一步完成之后，我们就可以将read比对至参考基因组了。
利用bwa mem比对模块将E.coli K12质控后的测序数据定位到其参考基因组上，同时通过管道('|' 操作符)将比对数据流引到samtools转换为BAM格式(SAM的二进制压缩格式)，然后重定向('>'操作符)输出到文件中保存下来。命令如下：

bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' \
../input/E.coli/reference/E.coli_K12_MG1655.fa \
../input/E.coli/cleandata/SRR1770413_1.clean.fastq.gz \
../input/E.coli/cleandata/SRR1770413_2.clean.fastq.gz \
| samtools view -Sb - > ../output/E.coli/map/E_coli_K12.bam

# Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
# -t，线程数，我们在这里使用4个线程
# -R 接的是Read Group的字符串信息，以@RG开头，制表符(\t)分离各选项
## 1) ID，这是Read Group的分组ID，一般设置为测序的lane ID
## 2) PL，指的是所用的测序平台
## 3) SM，样本ID
# 用samtools将它转化为BAM文件，为了有效节省磁盘空间，而且BAM会更加方便于后续的分析
# -b后面的 '-'，它代表是上面管道引流过来的数据

「Tips」BWA MEM比对模块是有一定适用范围的：它是专门为长read比对设计的，目的是为了解决，第三代测序技术这种能够产生长达几十kb甚至几Mbp的read情况。一般只有当read长度≥70bp的时候，才推荐使用，如果比这个要小，建议使用BWA ALN模块。

排序

FASTQ文件里面这些被测序下来的read是随机分布于基因组上面的，第一步的比对是按照FASTQ文件的顺序把read逐一定位到参考基因组上之后，随即就输出了，它不会也不可能在这一步里面能够自动识别比对位置的先后位置重排比对结果。因此，比对后得到的结果文件中，每一条记录之间位置的先后顺序是乱的，我们后续去重复等步骤都需要在比对记录按照顺序从小到大排序下来才能进行，所以这才是需要进行排序的原因。命令如下：

samtools sort \
-@ 4 \
-m 4G \
-O bam \
-o ../output/E.coli/sorted/E_coli_K12.sorted.bam \
../output/E.coli/map/E_coli_K12.bam

# Usage: samtools sort [options...] [in.bam]
# -@，用于设定排序时的线程数，我们设为4
# -m，限制排序时最大的内存消耗，这里设为4GB
# -O 指定输出为bam格式
# -o 是输出文件的名字，这里叫E_coli_K12.sorted.bam
# 最后就是输入文件——E_coli_K12.bam

# 删除不必要文件(可选)
## rm -f ../output/E.coli/map/E_coli_K12.bam

标记PCR重复

构建测序文库时细胞量的不足和打断过程中引起的部分降解，会使整体或者局部的DNA浓度过低。如果直接取样很可能会漏掉原本基因组上的一些DNA片段，导致测序不全。PCR扩增原本的目的就是为了增大微弱DNA序列片段的密度，但由于整个反应都在一个试管中进行，因此其他一些密度并不低的DNA片段也会被同步放大，在取样上机时，这些DNA片段就很可能被重复取到相同的几条去进行测序。最直接的后果就是同时增大了变异检测结果的假阴和假阳率。主要原因有：①DNA打断的那一步会发生一些损失，主要表现是会引发一些碱基发生颠换变换(嘌呤变嘧啶或者反之)，PCR扩大了这个信号。②PCR过程中也会带来新的碱基错误。③对于真实的变异，PCR反应可能会对包含某一个碱基的DNA模版扩增更加剧烈(这个现象称为PCR Bias)。因此，如果反应体系是对含有reference allele的模板扩增偏向强烈，那么变异碱基的信息会变小，从而会导致假阴。

PCR对真实的变异检测和个体的基因型判断都有不好的影响。GATK、Samtools、Platpus等这种利用贝叶斯原理的变异检测算法都是认为所用的序列数据都不是重复序列(即将它们和其他序列一视同仁地进行变异的判断，所以带来误导)，因此必须要进行标记(去除)或者使用PCR-Free的测序方案。

这里我们采用gatk的MarkDuplicates模块来实现这个功能，即只把重复序列在输出的新结果中标记出来，但不删除。命令如下：

gatk MarkDuplicates \
-I ../output/E.coli/sorted/E_coli_K12.sorted.bam \
-O ../output/E.coli/marked/E_coli_K12.sorted.markdup.bam \
-M ../output/E.coli/marked/E_coli_K12.sorted.markdup_metrics.txt

# 删除不必要文件(可选)
# rm -f ../output/E.coli/sorted/E_coli_K12.sorted.bam

创建比对索引文件

上面一步完成之后，我们需要为E_coli_K12.sorted.markdup.bam创建索引文件，它的作用能够让我们可以随机访问这个文件中的任意位置，而且后面的“局部重比对”步骤也要求这个BAM文件一定要有索引，命令如下：

samtools index ../output/E.coli/marked/E_coli_K12.sorted.markdup.bam

3. Indel区域重（“重新”的“重”）比对

merge

有时在进行"局部区域重比对"这一步骤之前还有一个merge的操作，将同个样本的所有比对结果合并成唯一一个大的BAM文件。之所以会有这种情况，是因为有些样本测得非常深，其测序结果需要经过多次测序(或者分布在多个不同的测序lane中)才全部获得，这个时候我们一般会先分别进行比对并去除重复序列后再使用samtools进行合并。

merge的例子如下：

samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

局部重比对

局部重比对的目的是将BWA比对过程中所发现有潜在序列插入或者序列删除(insertion和deletion，简称Indel)的区域进行重新校正。为什么需要进行这个校正呢？其根本原因来自于参考基因组的序列特点和BWA这类比对算法本身，注意这里不是针对BWA，而是针对所有的这类比对算法，包括bowtie等。这类在全局搜索最优匹配的算法在存在Indel的区域及其附近的比对情况往往不是很准确。另一个重要的原因是在这些比对算法中，对碱基错配和开gap的容忍度是不同的。具体体现在罚分矩阵的偏向上，例如，在read比对时，如果发现碱基错配和开gap都可以的话，它们会更偏向于错配。但是这种偏向错配的方式，有时候却还会反过来引起错误的开gap！这就会导致基因组上原本应该是一个长度比较大的Indel的地方，被错误地切割成多个错配和短indel的混合集，这必然会让我们检测到很多错误的变异。而且，这种情况还会随着所比对的read长度的增长(比如三代测序的Read，通常都有几十kbp)而变得越加严重。

GATK的局部重比对模块，应用Smith-Waterman算法，可以极其有效地实现对全局比对结果的校正和调整，最大程度低地降低由全局比对算法的不足而带来的错误。除此之外，还会对这个区域中的read进行一次局部组装，把它们连接成为长度更大的序列，这样能够更进一步提高局部重比对的准确性。

具体的做法是：

第一步，RealignerTargetCreator ，目的是定位出所有需要进行序列重比对的目标区域。候选的重比对区除了要在样本自身的比对结果中寻找之外，还应该把已知Indel区域也包含进来。

第二步，IndelRealigner，对所有在第一步中找到的目标区域运用算法进行序列重比对，最后得到捋顺了的新结果。

那么既然Indel局部重比对这么好，这么重要，似乎看起来在任何情况下都应该是必须的。然鹅，回答是否定的！但否定是有前提的！那就是我们后面的变异检测必须是使用GATK，而且必须使用GATK的HaplotypeCaller模块，仅当这个时候才可以减少这个Indel局部重比对的步骤。原因是GATK的HaplotypeCaller中，会对潜在的变异区域进行相同的局部重比对！但是其它的变异检测工具或者GATK的其它模块就没有这么干了！所以切记！

因为在GATK 4.0中没发现IndelRealigner这个功能，所以针对本文中的例子，这一步省略掉了。

4. 碱基质量值重校正(BQSR)

在WGS分析中，变异检测是一个极度依赖测序碱基质量值的步骤。因为这个质量值是衡量我们测序出来的这个碱基到底有多正确的重要(甚至是唯一)指标。BQSR(Base Quality Score Recalibration)主要是通过机器学习的方法构建测序碱基的错误率模型，然后对这些碱基的质量值进行相应的调整。

图中，横轴(Reported quality score)是测序结果在Base calling之后报告出来的质量值，也就是我们在FASTQ文件中看到的那些；纵轴(Empirical quality score)代表的是“真实情况的质量值”。

原理：
试想一下，如果我们看到某一个碱基报告的质量值是20时，那么它的预期错误率是1%，反过来想，就等于是说如果有100个质量值都是20的碱基，那么从统计上讲它们中将只有1个是错的！做了这个等效变换之后，我们的问题就可以转变成为寻找错误碱基的数量了。我们知道个体与个体之间的差异其实是很小的，那么在一个群体中发现的已知变异，在某个个体身上也很可能是同样存在的。因此，这个时候我们可以对比对结果进行直接分析，首先排除掉所有的已知变异位点，然后计算每个(报告出来的)质量值下面有多少个碱基在比对之后与参考基因组上的碱基是不同的，这些不同碱基就被我们认为是错误的碱基，它们的数目比例反映的就是真实的碱基错误率，换算成Phred score之后，就是纵轴的Empirical quality score了。

具体的做法是：
第一步，BaseRecalibrator，计算出所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件(sample_name.recal_data.table)。

第二步，ApplyBQSR，这一步利用第一步得到的校准表文件(sample_name.recal_data.table)重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新BAM文件。

❕注意，因为BQSR实际上是为了(尽可能)校正测序过程中的系统性错误，因此，在执行的时候是按照不同的测序lane或者测序文库来进行的，这个时候@RG信息(BWA比对时所设置的)就显得很重要了，算法就是通过@RG中的ID来识别各个独立的测序过程。

针对本文中的例子，没有执行BQSR是因为E.coli K12没有那些必须的known变异集(或者有但没找到)，所以无法进行。

放个human的例子在这：

gatk BaseRecalibrator \
-R /path/to/reference/Homo_sapiens_assembly38.fasta \
-I /path/to/sample_name.sorted.markdup.bam \
--known-sites /path/to/gatk/bundle/1000G_phase1.indels.hg38.vcf \
--known-sites /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf \
--known-sites /path/to/gatk/bundle/dbsnp_146.hg38.vcf \
-O /path/to/sample_name.sorted.markdup.recal_data.table

gatk ApplyBQSR \
--bqsr-recal-file /path/to/sample_name.sorted.markdup.recal_data.table \
-R /path/to/reference/Homo_sapiens_assembly38.fasta \
-I /path/to/sample_name.sorted.markdup.bam \
-O /path/to/sample_name.sorted.markdup.BQSR.bam

5. 变异检测

这是目前所有WGS数据分析流程的一个目标——获得样本准确的变异集合。这里变异检测的内容一般会包括：SNP、Indel，CNV和SV等，这个流程中我们只做其中最主要的两个：SNP和Indel。

我们这里使用GATK HaplotypeCaller模块对样本中的变异进行检测。它会先推断群体的单倍体组合情况，计算各个组合的几率，然后根据这些信息再反推每个样本的基因型组合。因此它不但特别适合应用到群体的变异检测中，而且还能够依据群体的信息更好地计算每个个体的变异数据和它们的基因型组合。

一般来说，在实际的WGS流程中对HaplotypeCaller的应用有两种做法，差别只在于要不要在中间生成一个gVCF：
(1) 直接进行HaplotypeCaller，这适合于单样本，或者那种固定样本数量的情况，也就是执行一次HaplotypeCaller之后就老死不相往来了。否则你会碰到仅仅只是增加一个样本就得重新运行这个HaplotypeCaller的坑爹情况(即，N+1难题)，而这个时候算法需要重新去读取所有人的BAM文件，这将会是一个很费时间的痛苦过程；
(2) 每个样本先各自生成gVCF，然后再进行群体joint-genotype。这其实就是GATK团队为了解决(1)中的N+1难题而设计出来的模式。gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了。

推荐选择第二种方式，变异检测不是一个样本的事情，有越多的同类样本放在一起joint calling结果将会越准确，而如果样本足够多的话，在低测序深度的情况下也同样可以获得完整并且准确的结果，而这样的分步方式是应对多样本的好方法。命令如下：

# 0. 为E.coli K12的参考序列生成一个.dict文件
gatk CreateSequenceDictionary \
-R ../input/E.coli/reference/E.coli_K12_MG1655.fa \
-O ../input/E.coli/reference/E.coli_K12_MG1655.dict

# 1. 生成中间文件gvcf
gatk HaplotypeCaller \
-R ../input/E.coli/reference/E.coli_K12_MG1655.fa \
--emit-ref-confidence GVCF \
-I ../output/E.coli/marked/E_coli_K12.sorted.markdup.bam \
-O ../output/E.coli/gvcf/E_coli_K12.g.vcf

# 2. 通过gvcf检测变异
gatk GenotypeGVCFs \
-R ../input/E.coli/reference/E.coli_K12_MG1655.fa \
-V ../output/E.coli/gvcf/E_coli_K12.g.vcf \
-O ../output/E.coli/gvcf/E_coli_K12.vcf

#3 压缩 
bgzip -f ../output/E.coli/gvcf/E_coli_K12.vcf

#4 构建tabix索引
tabix -p vcf ../output/E.coli/gvcf/E_coli_K12.vcf.gz

对于单个样本来说，除了上面这种方式，还有三种完成变异检测的方式，这里也以人为例，放上其他方式的命令供大家选择。

## 方式一，直接调用HaplotypeCaller输出样本VCF，面对较大的输入文件时，速度较慢
gatk HaplotypeCaller \
-R /path/to/reference/Homo_sapiens_assembly38.fasta \
-I /path/to/sample_name.sorted.markdup.BQSR.bam \
-O /path/to/sample_name.HC.vcf.gz

## 方式二，输出这个样本每个染色体的vcf，然后在合并所有的染色体结果，目的是提高速度，这不是必须的，仅是通过分染色体获得速度的提升
chrom=( chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY chrM )
for i in ${chrom[@]}; do
    gatk HaplotypeCaller \
    -R /path/to/Homo_sapiens_assembly38.fasta \
    -I /path/to/sample_name.sorted.markdup.BQSR.bam \
    -L $i \
    -O /path/to/sample_name.HC.${i}.vcf.gz
done && wait
merge_vcfs=""
for i in ${chrom[@]}; do
    merge_vcfs=${merge_vcfs}" -I /path/to/sample_name.HC.${i}.vcf.gz \\"\n
done && gatk MergeVcfs ${merge_vcfs} -O /path/to/sample_name.HC.vcf.gz

## 方式三，输出每个染色体的gvcf，然后对每个染色体单独进行GenotypeGVCFs，目的是提高速度，这不是必须的，仅是通过分染色体获得速度的提升
chrom=( chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY chrM )
for i in ${chrom[@]}; do
    gatk HaplotypeCaller \
    --emit-ref-confidence GVCF \
    -R /path/to/reference/Homo_sapiens_assembly38.fasta \
    -I /path/to/sample_name.sorted.markdup.BQSR.bam \
    -L $i \
    -O /path/to/sample_name.HC.${i}.g.vcf.gz && \
    gatk GenotypeGVCFs \
    -R /path/to/reference/Homo_sapiens_assembly38.fasta \
    -V /path/to/sample_name.HC.${i}.g.vcf.gz \
    -O /path/to/sample_name.HC.${i}.vcf.gz
done && wait
merge_vcfs=""
for i in ${chrom[@]}; do
    merge_vcfs=${merge_vcfs}" -I /path/to/sample_name.HC.${i}.vcf.gz \\"\n
done && gatk MergeVcfs ${merge_vcfs} -O /path/to/sample_name.HC.vcf.gz

为了保持一致，来看一下现在的目录结构，供大家对照。

./wgs_practice/
├── bin
├── input
│   └── E.coli
│       ├── reference
│       │   ├── E.coli_K12_MG1655.dict
│       │   ├── E.coli_K12_MG1655.fa
│       │   ├── E.coli_K12_MG1655.fa.amb
│       │   ├── E.coli_K12_MG1655.fa.ann
│       │   ├── E.coli_K12_MG1655.fa.bwt
│       │   ├── E.coli_K12_MG1655.fa.fai
│       │   ├── E.coli_K12_MG1655.fa.pac
│       │   ├── E.coli_K12_MG1655.fa.sa
│       │   └── work.log.sh
│       ├─── rawdata
│       │   ├── SRR1770413_1.fastq.gz
│       │   ├── SRR1770413_2.fastq.gz
│       │   └── work.log.sh
│       └──── cleandata
│           ├── SRR1770413_1.clean.fastq.gz
│           ├── SRR1770413_2.clean.fastq.gz
│           └── work.log.sh
├── output
│   └── E.coli
│       ├── E_coli_K12.g.vcf
│       ├── E_coli_K12.g.vcf.idx
│       ├── E_coli_K12.sorted.markdup.bam
│       ├── E_coli_K12.sorted.markdup.bam.bai
│       ├── E_coli_K12.sorted.markdup_metrics.txt
│       ├── E_coli_K12.vcf
│       └── E_coli_K12.vcf.idx
└── work.log.sh

6.变异结果质控和过滤

这是我们这个流程中的最后一步了。什么是质控？质控的含义和目的是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。首选的质控方案是GATK VQSR(Variant Quality Score Recalibration)，它通过机器学习的方法利用多个不同的数据特征训练一个模型(高斯混合模型)对变异数据进行质控。

下面是大家经常都会用到的VQSR基本命令，以human为例:

## VariantRecalibrator用来进行模型计算，获得数据的情况
## ApplyVQSR则是根据我们设定的ts_filter_level参数，最终过滤得到数据

## 首先是SNP mode
time gatk VariantRecalibrator \
   -R $reference/Homo_sapiens_assembly38.fasta \
   -V $outdir/poplation/${outname}.HC.vcf.gz \
   -resource:hapmap,known=false,training=true,truth=true,prior=15.0 $GATK_bundle/hapmap_3.3.hg38.vcf \
   -resource:omini,known=false,training=true,truth=false,prior=12.0 $GATK_bundle/1000G_omni2.5.hg38.vcf \
   -resource:1000G,known=false,training=true,truth=false,prior=10.0 $GATK_bundle/1000G_phase1.snps.high_confidence.hg38.vcf \
   -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $GATK_bundle/dbsnp_146.hg38.vcf \
   -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
   -mode SNP \
   -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
   -rscriptFile $outdir/poplation/${outname}.HC.snps.plots.R \
   --tranches-file $outdir/poplation/${outname}.HC.snps.tranches \
   -O $outdir/poplation/${outname}.HC.snps.recal && \
time $gatk ApplyVQSR \
   -R $reference/Homo_sapiens_assembly38.fasta \
   -V $outdir/poplation/${outname}.HC.vcf.gz \
   --ts_filter_level 99.0 \
   --tranches-file $outdir/poplation/${outname}.HC.snps.tranches \
   -recalFile $outdir/poplation/${outname}.HC.snps.recal \
   -mode SNP \
   -O $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz && echo "** SNPs VQSR done **"

## 然后是Indel mode
time $gatk VariantRecalibrator \
   -R $reference/Homo_sapiens_assembly38.fasta \
   -input $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz \
   -resource:mills,known=true,training=true,truth=true,prior=12.0 $GATK_bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf \
   -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $GATK_bundle/dbsnp_146.hg38.vcf \
   -an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
   -mode INDEL \
   --max-gaussians 6 \
   -rscriptFile $outdir/poplation/${outname}.HC.snps.indels.plots.R \
   --tranches-file $outdir/poplation/${outname}.HC.snps.indels.tranches \
   -O $outdir/poplation/${outname}.HC.snps.indels.recal && \
time $gatk ApplyVQSR \
   -R $reference/Homo_sapiens_assembly38.fasta \
   -input $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz \
   --ts_filter_level 99.0 \
   --tranches-file $outdir/poplation/${outname}.HC.snps.indels.tranches \
   -recalFile $outdir/poplation/${outname}.HC.snps.indels.recal \
   -mode INDEL \
   -O $outdir/poplation/${outname}.HC.VQSR.vcf.gz && echo "** SNPs and Indels VQSR (${outname}.HC.VQSR.vcf.gz finish) done **"

小拓展：
-resource后跟参数含义：
训练集名字，这个名字是可以随便改动的，但是为了便于交流，一般还是默认按照数据集的名字来设置(如上面的例子)；
known：该数据是否作为已知变异数据集，用于对变异数据的标注；
training：该数据是否作为模型训练的数据集，用于训练VQSR模型；
truth：该数据是否作为验证模型训练的真集数据，这个数据同时还是VQSR训练bad model时自动进行参数选择的重要数据；
prior：该数据集在VQSR模型训练中的权重，或者叫Prior likelihood(这里转化为Phred-scale，比如20代表的是0.99)。

对于SNP来说，VQSR的训练集数据目前主要有四个：HapMap，Omni，1000G和dbSNP。
第一个是HapMap，它来自国际人类单倍体型图计划，这个数据集包含了大量家系数据，并且有非常严格的质控和严密的实验验证，因此它的准确性是目前公认最高的。
第二个是Omni，这个数据源自Illumina的Omni基因型芯片，大概2.5百万个位点，它也是一个高可信的变异结果。
第三个是1000G，它是千人基因组计划(1000 genomes project)质控后的变异数据，一般不作为模型验证的真集数据。
第四个是dbSNP，这是一个绝对不可以作为训练集位点的数据——收集的数据实际都是研究者们发表了相关文章提交上来的变异，这些变异很多是没做过严格验证的，很多甚至还是假的，在没被反复验证之前，是不可信的，dbSNP的唯一作用就是用于标注我们的变异集中哪些是已经在其它研究中出现过的。

对于Indel来说，VQSR模型的训练参数只有两个：Mill和dbSNP。
第一个是Mills，对于Indel来说能正在算得上真集的并不多，Mills_and_1000G_gold_standard.indels.hg38.vcf算是其中一个，并被专门做过验证。
第二个还是dbSNP，这个就和上面SNP模式下的作用是一样的。不过假如这一步对Indel进行VQSR的VCF数据是顺着上面SNP VQSR后下来的话，那么这个dbSNP的参数可以省略。

不幸的是，使用VQSR需要具备以下两个条件：
第一，需要一个精心准备的已知变异集，它将作为训练质控模型的真集。人类的基因组数据，就有一套适合用来训练过滤模型的已知变异集(dbSNP，1000G，Hapmap和omini等)。GATK的bundle主要就是对这些高质量的数据集做了精心的处理和选择，然后把它们作为VQSR时的真集位点。强调一点，是真集的『位点』而不是真集的『数据』！因为，VQSR并不是用这些变异集里的数据来训练的，而是用我们自己的变异数据。实际上，这些已知变异集的意义是告诉我们群体中哪些位点存在着变异，如果在其他人的数据里能观察到落入这个集合中的变异位点，那么这些被已知集包括的变异就有很大的可能是正确的。也就是说，我们可以从数据中筛选出那些和真集『位点』相同的变异，把它们当作是真实的变异结果。接着，进行VQSR的时候，程序就可以用这个筛选出来的数据作为真集数据来训练，并构造模型了。第二，要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。

很遗憾我们这个E.coli K12的变异结果并不适合通过VQSR来进行过滤，那么碰到这种情况的时候该怎么办？这个时候，我们就不得不选择硬过滤的方式来质控了。所谓硬过滤其实就是通过人为设定一个或者若干个指标阈值（也可以叫数据特征值），然后把所有不满足阈值的变异位点采用一刀切掉的方法。我们可以直接使用GATK VQSR所用的指标来执行硬过滤——毕竟这些指标都是经过精挑细选的，既然VQSR就是用这些指标来训练质控模型的，那么它们就可以在一定程度上描述每个变异的质量，我们用这些指标设置对应的阈值来进行硬过滤也将是合理的。

VQSR使用的数据指标有6个：
QualByDepth（QD）
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)

- QualByDepth（QD）

QD是变异质量值(Quality)除以覆盖深度(Depth)得到的比值。这里的变异质量值就是VCF中QUAL的值——用来衡量变异的可靠程度，这里的覆盖深度是这个位点上所有含有变异碱基的样本的覆盖深度之和，通俗一点说，就是这个值可以通过累加每个含变异的样本（GT为非0/0的样本）的覆盖深度（VCF中每个样本里面的DP）而得到。

- FisherStrand（FS）

FS是一个通过Fisher检验的p-value转换而来的值，它要描述的是测序或者比对时对于只含有变异的read以及只含有参考序列碱基的read是否存在着明显的正负链特异性（Strand bias，或者说是差异性）。这个差异反应了测序过程不够随机，或者是比对算法在基因组的某些区域存在一定的选择偏向。如果测序过程是随机的，比对是没问题的，那么不管read是否含有变异，以及是否来自基因组的正链或者负链，只要是真实的它们就都应该是比较均匀的，也就是说，不会出现链特异的比对结果，FS应该接近于零。

- StrandOddsRatio (SOR)

SOR原理和FS类似，但是用了不同的统计检验方法计算差异程度，用的是symmetric odds ratio test，它更适合于高覆盖度的数据。

- RMSMappingQuality (MQ)

MQ这个值是所有比对至该位点上的read的比对质量值的均方根(先平方、再平均、然后开方)。它和平均值相比更能够准确地描述比对质量值的离散程度。

- MappingQualityRankSumTest (MQRankSum)

对杂合位点进行的不同碱基之间比对质量的曼惠特尼秩和检验结果，通过ref和alt碱基的比对质量差异来评估位点的可信度。

- ReadPosRankSumTest (ReadPosRankSum)

对杂合位点进行的秩和检验，看不同的碱基是否倾向于出现在reads上的特定位置(例如接近reads的起始或者终止)。

对于本文的例子，执行硬过滤，命令如下：

# 使用SelectVariants，选出SNP
gatk SelectVariants \
-select-type SNP \
-V ../output/E.coli/gvcf/E_coli_K12.vcf.gz \
-O ../output/E.coli/sel/E_coli_K12.snp.vcf.gz

# 为SNP作硬过滤
gatk VariantFiltration \
-V ../output/E.coli/sel/E_coli_K12.snp.vcf.gz \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "Filter" \
-O ../output/E.coli/sel/E_coli_K12.snp.filter.vcf.gz

# 使用SelectVariants，选出Indel
gatk SelectVariants \
-select-type INDEL \
-V ../output/E.coli/gvcf/E_coli_K12.vcf.gz \
-O ../output/E.coli/sel/E_coli_K12.indel.vcf.gz

# 为Indel作过滤
gatk VariantFiltration \
-V ../output/E.coli/sel/E_coli_K12.indel.vcf.gz \
--filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "Filter" \
-O ../output/E.coli/sel/E_coli_K12.indel.filter.vcf.gz

# 重新合并过滤后的SNP和Indel
gatk MergeVcfs \
-I ../output/E.coli/sel/E_coli_K12.snp.filter.vcf.gz \
-I ../output/E.coli/sel/E_coli_K12.indel.filter.vcf.gz \
-O ../output/E.coli/sel/E_coli_K12.filter.vcf.gz

# 删除无用中间文件
rm -f ../output/E.coli/sel/E_coli_K12.snp.vcf.gz* \
../output/E.coli/sel/E_coli_K12.snp.filter.vcf.gz* \
../output/E.coli/sel/E_coli_K12.indel.vcf.gz* \
../output/E.coli/sel/E_coli_K12.indel.filter.vcf.gz*

「Tips」不管我们的指标和阈值设置的多么完美，硬过滤的硬伤都是一刀切，它并不会根据多个维度的数据自动设计出更加合理的过滤模式。硬过滤作为一个简单粗暴的方法，不到不得已的时候不推荐使用，即使使用了，也一定要明白每一个过滤指标和阈值都意味着什么。

7.小结

到这里，这篇文章就结束了。文章基本包含了WGS实践的一个基本流程，但也仅仅只是一个流程的搭建。我们的个人能力不能只是会跑一个流程，我想，更深层次地，是要去掌握如何分析数据的能力，如何发现问题和解决数据问题等的能力，理解问题的本质，选择合适的工具，而不是反过来被工具所束缚。毕竟，工具是死的，人是活的，用聪明的脑袋去灵活地使用工具。

参考：

[1] 使用fastp进行数据质控

[2] NGS数据过滤之Trimmomatic详细说明

[3] fastp

[4] 从零开始完整学习全基因组测序数据分析：第2节 FASTA和FASTQ

[5] 从零开始完整学习全基因组测序数据分析：第3节数据质控

[6] 从零开始完整学习全基因组测序数据分析：第4节构建WGS主流程

[7] 从零开始完整学习全基因组测序数据分析：第5节理解并操作BAM文件

[8] GATK4.0和全基因组数据分析实践（上）

[9] GATK4.0和全基因组数据分析实践（下）

[10] GATK4全基因组数据分析最佳实践，我以这篇文章为标志，终结当前WGS系列数据分析的流程主体问题 | 完全代码

[11] 我应该如何正确设置GATK VQSR的模型训练参数？

该文档首发于微信公众号：浮生小记x。

编辑于 2021-12-10 · 著作权归作者所有

赞同 66