自己在研究生的第一个项目就是处理重测序数据，在做的过程中参考了很多资料，现在把整个重测序的上游分析流程分享给大家，一些关键环节已经帮大家踩过坑了，相信跟着学习实践一定能够收获很多。制作不易，希望大家多多关注点赞支持！！

0.准备阶段

在开始之前，我们需要做一些准备工作，主要是部署好相关的软件和工具。我们在这个重测序数据分析过程中用到的所有软件都是开源的，它们的代码全部都能够在github上找到，具体如下：

BWA: 这是最权威，使用最广的NGS数据比对软件；conda install bowtie2

Samtools: 是一个专门用于处理比对数据的工具；conda install -c bioconda bwa #-c 即channel，下载生信软件都要加上-c bioconda software

Picard: 它是目前最著名的组学研究中心-Broad研究所开发的一款强大的NGS数据处理工具，功能方面和Samtools有些重叠，但更多的是互补，它是由java编写的，我们直接下载最新的.jar包就行了。conda install -c bioconda picard

GATK :目前GATK有3.x和4.x两个不同的版本，代码在github上也是分开的。4.x是今年新推出的，在核心算法层面并没太多的修改，但使用了新的设计模式，做了很多功能的整合，是更适合于大规模集群和云运算的版本，后续GATK团队也将主要维护4.x的版本，而且它的代码是100%开源的，这和3.x只有部分开源的情况不同。看得出GATK今年的这次升级是为了应对接下来越来越多的大规模人群测序数据而做出的改变，但现阶段4.x版本还不稳定，真正在使用的人和机构其实也还不多。短期来看，3.x版本还将在业内继续使用一段时间；其次，3.x对于绝大分部的分析需求来说是完全足够的。我们在这里以GATK4作为流程的重要工具进行分析流程的构建。conda install -c bioconda gatk4

事实上，对于构造重测序分析流程来说，以上这个四个工具就完全足够了。它们的安装都非常简单，除了BWA和Samtools由C编写的，安装时需要进行编译之外，另外两个只要保证系统中的java是1.8.x版本及以上的，那么直接下载jar包就可以使用了。操作系统方面推荐linux（集群）或者Mac OS。

1.原始数据质控

这一部分考虑到现在公司给出的测序数据已经进行了质控，到我们手里基本都是clean data,网上关于质控的资源也比较多，基本上通过fastqc和fastp就可以进行质控和数据清洗。所以这里不再赘述。

2.数据预处理

序列比对

先问一个问题：为什么需要比对？

我们已经知道NGS测序下来的短序列（read）存储于FASTQ文件里面。虽然它们原本都来自于有序的基因组，但在经过DNA建库和测序之后，文件中不同read之间的前后顺序关系就已经全部丢失了。因此，FASTQ文件中紧挨着的两条read之间没有任何位置关系，它们都是随机来自于原本基因组中某个位置的短序列而已。

因此，我们需要先把这一大堆的短序列捋顺，一个个去跟该物种的 参考基因组【注】比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。这 也是核心流程真正意义上的第一步，只有完成了这个序列比对我们才有下一步的数据分析。

【注】参考基因组：指该物种的基因组序列，是已经组装成的完整基因组序列，常作为该物种的标准参照物，比如人类基因组参考序列，fasta格式。

序列比对本质上是一个寻找最大公共子字符串的过程。大家如果有学过生物信息学的话，应该或多或少知道BLAST，它使用的是动态规划的算法来寻找这样的子串，但在面对巨量的短序列数据时，类似BLAST这样的软件实在太慢了！因此，需要更加有效的数据结构和相应的算法来完成这个搜索定位的任务。

我们这里将用于流程构建的BWA就是其中最优秀的一个，它将BW(Burrows-Wheeler)压缩算法和后缀树相结合，能够让我们以较小的时间和空间代价，获得准确的序列比对结果。

以下我们就开始流程的搭建。

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

$ bwa index genome.fasta #对基因组文件进行bwa索引

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

这些就是在比对时真正需要被用到的文件。这一步完成之后，我们就可以将read比对至参考基因组了：

.├── genome.fasta.amb
 ├── genome.fasta.ann
 ├── genome.fasta.bwt
 ├── genome.fasta.pac
 └── genome.fasta.sa

这些就是在比对时真正需要被用到的文件。这一步完成之后，我们就可以将read比对至参考基因组了：

$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/genome.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b > sample.bam

大伙如果以前没使用过这个比对工具的话，那么可能不明白上面参数的含义。我们这里调用的是bwa的mem比对模块，在解释这样做之前，我们不妨先看一下bwa mem的官方用法说明，它就一句话：

Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]

其中，[options]是一系列可选的参数，暂时不多说。这里的 < idxbase>要输入的是参考基因组的BW索引文件，我们上面通过bwa index构建好的那几个以human.fasta为前缀的文件便是；< in1.fq>和 [in2.fq]输入的是质控后的fastq文件。但这里输入的时候为什么会需要两个fq（in1.fq和in2.fq）呢？我们上面的例子也是有两个：read_1.fq.gz和read_2.fq.gz。这是因为这是双末端测序（也称Pair-End）的情况，那什么是“双末端测序”呢？这两个fq之间的关系又是什么？这个我需要简单解析一下。

我们已经知道NGS是短读长的测序技术，一次测出来的read的长度都不会太长，那为了尽可能把一个DNA序列片段尽可能多地测出来，既然测一边不够，那就测两边，于是就有了一种从被测DNA序列两端各测序一次的模式，这就被称为双末端测序（Pair-End Sequencing，简称PE测序）。如下图是Pair-End测序的示意图，中间灰色的是被测序的DNA序列片段，左边黄色带箭头和右边蓝色带箭头的分别是测序出来的read1和read2序列，这里假定它们的长度都是100bp。虽然很多时候Pair-End测序还是无法将整个被测的DNA片段完全测通，但是它依然提供了极其有用的信息，比如，我们知道每一对的read1和read2都来自于同一个DNA片段，read1和read2之间的距离是这个DNA片段的长度，而且read1和read2的方向刚好是相反的（这里排除mate-pair的情形）等，这些信息对于后面的变异检测等分析来说都是非常有用的。

Pair-End 测序

另外，在read1在fq1文件中位置和read2在fq2文件中的文件中的位置是相同的，而且read ID之间只在末尾有一个’/1’或者’/2’的差别。

read1 ID和read2 ID的差别

既然有双末端测序，那么与之对应的就有单末端测序（Single End Sequecing，简称SE测序），即只测序其中一端。因此，我们在使用bwa比对的时候，实际上，in2.fq是非强制性的（所以用方括号括起来），只有是双末端测序的数据时才需要添加。

回到上面我们的例子，大伙可以看到我这里除了用法中提到的参数之外，还多了2个额外的参数，分别是：-t，线程数，我们在这里使用4个线程；-R 接的是 Read Group的字符串信息，这是一个非常重要的信息，以@RG开头，它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。在Read Group中，有如下几个信息非常重要：

(1) ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；

(2) PL，指的是所用的测序平台，这个信息不要随便写！特别是当我们需要使用GATK进行后续分析的时候，更是如此！这是一个很多新手都容易忽视的一个地方，在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台，当然，如果实在不知道，那么必须设置为UNKNOWN，名字方面不区分大小写。如果你在分析的时候这里没设置正确，那么在后续使用GATK过程中可能会碰到类似如下的错误：

ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.

这个时候你需要对比对文件的header信息进行重写，就会稍微比较麻烦。

我们上面的例子用的是PL:illumina。如果你的数据是CG测序的那么记得不要写成CG！而要写COMPLETE。

(3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本。

(4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；

除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后，在RG字符串中要用制表符（\t）将它们分开。

最后在我们的例子中，我们将比对的输出结果直接重定向到一份sample_name.sam文件中，这类文件是BWA比对的标准输出文件，它的具体格式我会在下一篇文章中进行详细说明。但SAM文件是文本文件，一般整个文件都非常巨大，因此，为了有效节省磁盘空间，一般都会用samtools将它转化为BAM文件（SAM的特殊二进制格式），而且BAM会更加方便于后续的分析。所以我们上面比对的命令可以和samtools结合并改进为：

$ bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam

我们通过管道(“|”)把比对的输出如同引导水流一样导流给samtools去处理，上面samtools view的-b参数指的就是输出为BAM文件，这里需要注意的地方是-b后面的’-‘，它代表就是上面管道引流过来的数据，经过samtools转换之后我们再重定向为sample_name.bam。

关于BWA的其他参数，我这里不打算对其进行一一解释，在绝大多数情况下，采用默认是合适的做法。

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

对bam文件进行进行重新排序（reorder）

由BWA生成的sam文件时按字典式排序法进行的排序（lexicographically）进行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在进行callsnp的时候是按照染色体组型（karyotypic）进行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要对原始sam文件进行reorder。可以使用picard-tools中的ReorderSam完成。

Picard的SortSam需指定一个tmp目录，用于存放中间文件，中间文件会很大，above 10G.注意指定目录的空间大小。

排序

以上，我们就完成了read比对的步骤。接下来是排序：

排序这一步我们也是通过使用samtools来完成的，命令很简单：

Usage: samtools sort [options...] [in.bam]

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

$ time samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam

其中，-@，用于设定排序时的线程数，我们设为4；-m，限制排序时最大的内存消耗，这里设为4GB；-O 指定输出为bam格式；-o 是输出文件的名字，这里叫sample_name.sorted.bam。我会比较建议大伙在做类似分析的时候在文件名字将所做的关键操作包含进去，因为这样即使过了很长时间，当你再去看这个文件的时候也能够立刻知道当时对它做了什么；最后就是输入文件——sample_name.bam。

【注意】排序后如果发现新的BAM文件比原来的BAM文件稍微小一些，不用觉得惊讶，这是压缩算法导致的结果，文件内容是没有损失的。

去除重复序列（或者标记重复序列）

在排序完成之后我们就可以开始执行去除重复（准确来说是 去除PCR重复序列）的步骤了。

首先，我们需要先理解什么是重复序列，它是如何产生的，以及为什么需要去除掉？要回答这几个问题，我们需要再次理解在建库和测序时到底发生了什么。

我们知道在NGS测序之前都需要先构建测序文库：通过物理（超声）打断或者化学试剂（酶切）切断原始的DNA序列，然后选择特定长度范围的序列去进行PCR扩增并上机测序。

因此，这里重复序列的来源实际上就是由PCR过程中所引入的。因为所谓的PCR扩增就是把原来的一段DNA序列复制多次。可是为什么需要PCR扩增呢？如果没有扩增不就可以省略这一步了吗？

情况确实如此，但是很多时候我们构建测序文库时能用的细胞量并不会非常充足，而且在打断的步骤中也会引起部分DNA的降解，这两点会使整体或者局部的DNA浓度过低，这时如果直接从这个溶液中取样去测序就很可能漏掉原本基因组上的一些DNA片段，导致测序不全。而PCR扩增的作用就是为了把这些微弱的DNA多复制几倍乃至几十倍，以便增大它们在溶液中分布的密度，使得能够在取样时被获取到。所以这里大家需要记住一个重点，PCR扩增原本的目的是为了增大微弱DNA序列片段的密度，但由于整个反应都在一个试管中进行，因此其他一些密度并不低的DNA片段也会被同步放大，那么这时在取样去上机测序的时候，这些DNA片段就很可能会被重复取到相同的几条去进行测序（下图为PCR扩增示意图）。

PCR扩增示意图：PCR扩增是一个指数扩增的过程，图中原本只有一段双链DNA序列，在经过3轮PCR后就被扩增成了8段

看到这里，你或许会觉得，那没必要去除不也应该可以吗？因为即便扩增了多几次，不也同样还是原来的那一段DNA吗？直接用来分析对结果也不会有影响啊！难道不是吗？

会有影响，而且有时影响会很大！最直接的后果就是同时增大了变异检测结果的假阴和假阳率。主要有几个原因：

• DNA在打断的那一步会发生一些损失，主要表现是会引发一些碱基发生颠换变换（嘌呤-变嘧啶或者嘧啶变嘌呤），带来假的变异。PCR过程会扩大这个信号，导致最后的检测结果中混入了假的结果；
• PCR反应过程中也会带来新的碱基错误。发生在前几轮的PCR扩增发生的错误会在后续的PCR过程中扩大，同样带来假的变异；
• 对于真实的变异，PCR反应可能会对包含某一个碱基的DNA模版扩增更加剧烈（这个现象称为PCR Bias）。如果反应体系是对含有reference allele的模板扩增偏向强烈，那么变异碱基的信息会变小，从而会导致假阴。

PCR对真实的变异检测和个体的基因型判断都有不好的影响。GATK、Samtools、Platpus等这种利用贝叶斯原理的变异检测算法都是认为所用的序列数据都不是重复序列（即将它们和其他序列一视同仁地进行变异的判断，所以带来误导），因此必须要进行标记（去除）或者使用PCR-Free的测序方案（这个方案目前正变得越来越流行，特别是对于RNA-Seq来说尤为重要，现在著名的基因组学研究所——Broad Institute，基本都是使用PCR-Free的测序方案）。

那么具体是如何做到去除这些PCR重复序列的呢？我们可以抛开任何工具，仔细想想，既然PCR扩增是把同一段DNA序列复制出很多份，那么这些序列在经过比对之后它们一定会定位到基因组上相同的位置，比对的信息看起来也将是一样的！于是，我们就可以根据这个特点找到这些重复序列了！

事实上，现有的工具包括Samtools和Picard中去除重复序列的算法也的确是这么做的。不同的地方在于，samtools的rmdup是直接将这些重复序列从比对BAM文件中删除掉，而Picard的MarkDuplicates默认情况则只是在BAM的FLAG信息中标记出来，而不是删除，因此这些重复序列依然会被留在文件中，只是我们可以在变异检测的时候识别到它们，并进行忽略。

考虑到尽可能和现在主流的做法一致（但我并不是说主流的做法就一定是对的，要分情况看待，只是主流的做法容易被做成生产流程而已），我们这里也用Picard来完成这个事情：

picard MarkDuplicates I=sample_name.sorted.bam O=sample_name.sorted.markdup.bam M=sample_name.markdup_metrics.txt

这里只把重复序列在输出的新结果中标记出来，但不删除。如果我们非要把这些序列完全删除的话可以这样做：java -jar 把参数REMOVE_DUPLICATES设置为ture，那么重复序列就被删除掉，不会在结果文件中留存。我比较建议使用第一种做法，只是标记出来，并留存这些序列，以便在你需要的时候还可以对其做分析。

picard MarkDuplicates REMOVE_DUPLICATES=true I=sample_name.sorted.bam O=sample_name.sorted.markdup.bam M=sample_name.markdup_metrics.txt

创建索引

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

samtools index sample_name.sorted.markdup.bam

完成之后，会生成一份sample_name.sorted.markdup.bam.bai文件，这就是上面这份BAM的index。

同时，还要对参考基因组进行GATK的索引,也就是准备参考基因组.fai和.dict文件

gatk CreateSequenceDictionary -R genome.fa -O genome.dict && samtools faidx genome.fa

生成gvcf文件并合并

有多个个体的情况下，首先我们用HaplotypeCaller命令给每一个个体生成gvcf文件，然后用CombineGVCFs按染色体合并gvcf文件。

我是用shell生成的批量脚本然后放到服务器上一起跑的，就会比较快。批量写脚本什么语言都可以，就按照你的需求循环改变名称、个体号、染色体号等各种字符就行

#1 生成gvcf文件
gatk  HaplotypeCaller -R ref.fa -I sample_name.sorted.markdup.bam --emit-ref-confidence GVCF --min-base-quality-score 10 -O chrx.g.vcf.gz
#2 gvcf文件按染色体合并
ls chrx.g.vcf.gz > chrx_gvcf.list 
gatk CombineGVCFs -R ref.fa -V chrx_gvcf.list   -L X(染色体号) -O chrx.merged.g.vcf.gz

生成基因型文件(这一步变成vcf了）

gatk GenotypeGVCFs -R ref.fa -V chrx.merged.g.vcf.gz -O chrx.genotype.vcf.gz

过滤

上一步得到的是rawdata，我们还需要对原始数据做变异质控。质控是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。

SNP过滤有两种情况，一种是仅根据位点质量信息（测序深度，回帖质量等）对SNP进行粗过滤。另一种是考虑除了测序质量以外的信息进行的过滤。 接下来我会分别介绍到这两种过滤。

从测序位点质量上看，在GATK HaplotypeCaller之后，首选的质控方案是GATK VQSR， 原理是利用自身数据和已知变异位点集的overlap，通过GMM模型构建一个分类器来对变异数据进行打分，从而评估每个位点的可行度。

虽然经常看到VQSR的教程，但是很可惜目前应该只有人类上可以做（还是需要高质量的已知变异集），所以我们其他物种只能做hard-filtering硬过滤了。

挑选snp

gatk SelectVariants -R genome.fa -variant sample.vcf -O sample_snp.vcf -select-type SNP

过滤snp

gatk  VariantFiltration -variant sample_snp.vcf -O sample_snp_filter.vcf -R genome.fa