第4节构建WGS主流程

这是WGS数据分析的流程图。流程的目的是准确检测出每个样本（这里特指人）基因组中的变异集合，也就是人与人之间存在差异的那些DNA序列。

一、准备阶段

这个WGS数据分析过程中用到的所有软件都是开源的，它们的代码全部都能够在 github 上找到，具体如下：

BWA（Burrow-Wheeler Aligner）： 这是最权威，使用最广的NGS数据比对软件，目前已经更新到0.7.16版本；

Samtools：是一个专门用于处理比对数据的工具，由BWA的作者（lh3）所编写；

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

GATK：同样是 Broad 研究所开发的，是目前业内最权威、使用最广的基因数据变异检测工具。值得注意的是，目前 GATK 有 3.x 和 4.x 两个不同的版本，代码在 github 上也是分开的。

4.x 是今年新推出的，在核心算法层面并没太多的修改，但使用了新的设计模式，做了很多功能的整合，是更适合于大规模集群和云运算的版本，后续GATK 团队也将主要维护 4.x 的版本，而且它的代码是100%开源的。

 3.x 只有部分是开源的，其次，3.x 对于绝大分部的分析需求来说是完全足够的。我们在这里也以 GATK3.8（最新版本）作为流程的重要工具进行分析流程的构建。

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

二、原始数据质控

数据的质控，之前文章已经介绍过，这里不再赘述。

三、数据预处理

序列比对

（1）为什么需要比对？

① NGS测序数据比对的必要性：

NGS（Next-Generation Sequencing）测序得到的短序列（reads）存储在FASTQ文件中。尽管这些reads来自有序的基因组，但由于测序和建库过程，文件中的read顺序完全随机，丢失了原始基因组中的位置信息。为了解析数据，我们需要通过比对找到每条read在参考基因组中的确切位置，并按顺序排列。这就是“比对”的核心任务，也是后续分析的基础。

② 比对的原理：

比对的目标是寻找短序列（reads）与参考基因组之间的最大公共子字符串。传统工具如BLAST使用动态规划完成比对，但面对NGS海量数据，其速度无法满足需求。为此，诸如BWA等工具采用更高效的算法和数据结构，例如Burrows-Wheeler变换（BWT）和后缀树，实现快速且高效的比对。

③ 比对的目的：

为了重新恢复这些短序列在基因组中的位置，我们需要将它们逐一与参考基因组进行比较，确定每条read在基因组中的位置并排序。这一过程称为序列比对，是NGS数据分析的第一步。只有完成比对，后续的分析（如变异检测、基因表达分析等）才能进行。

④ 核心概念

1）参考基因组：

指该物种的完整基因组序列，通常作为标准参照。

通常以FASTA格式存储。

2）序列比对的本质：

比对的核心是找到reads与参考基因组之间的最大公共子串。

BLAST等早期工具基于动态规划算法实现这一任务，但面对NGS产生的海量数据，速度和效率已无法满足需求。

(2)现代比对工具：BWA

① 高效比对的实现：

BWA（Burrows-Wheeler Aligner）结合了Burrows-Wheeler压缩算法和后缀树，极大地提高了比对速度和效率，在短时间内完成精确的序列比对。

第一步：构建参考基因组的索引

在进行比对前，需要为参考基因组创建索引，这是对序列进行Burrows-Wheeler变换（BWT）的过程，以便快速搜索和定位。

代码示例：bwa index human.fasta

以人类基因组为例（长度约3Gb），索引构建通常需要3小时左右。完成后，生成的索引文件如下：

---

实操：比对步骤

1.参考基因组索引构建

在比对前，需为参考基因组创建索引，以优化搜索效率。

代码示例：

$ bwa index human.fasta

例如，人类基因组约3Gb，索引构建可能耗时数小时，生成的文件包括human.fasta的多个前缀文件。

2.双末端测序（PE）比对

以BWA的mem模块为例，比对双末端测序数据：

$ 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 > sample_name.sam

注释：

参数解析：

-t 4：使用4线程，提高运行速度。

-R：指定Read Group（RG）信息，用于分组和后续分析：

ID：分组ID（如测序lane ID）。

PL：测序平台（如illumina）。

SM：样本ID。

LB：文库名称，区分不同组。

/path/to/human.fasta：参考基因组的索引文件。

read_1.fq.gz和read_2.fq.gz：双末端测序的两个FASTQ文件。

关于双末端测序：

Pair-End测序（PE）：对DNA片段两端分别测序，生成read1和read2。其方向相反且来自同一片段，提供更多信息（如片段长度、方向等），便于后续变异分析。

单末端测序（SE）：仅测序一端，仅需提供一个FASTQ文件。

---

输出优化：SAM到BAM转换

由于SAM文件是文本格式且占用空间较大，推荐将其转换为更高效的BAM格式：$ 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 -S -b：将SAM转换为BAM。

>：保存输出为sample_name.bam。

注意事项

模块适用性：

BWA MEM模块适用于长read（≥70bp），尤其针对第三代测序数据。

短read（<70bp）推荐使用BWA ALN模块。

PL参数设置：

确保测序平台设置正确（如ILLUMINA）。错误设置可能导致后续分析工具（如GATK）报错。

通过以上步骤，完成了NGS数据的核心环节：序列比对，为后续分析（如变异检测）奠定基础。

---

接下来是排序：

排序是对比对后生成的 BAM 文件进行位置排序的必要步骤，主要目的是确保比对结果按基因组位置从小到大排列，以便后续分析工具（如去重复标记、变异检测等）能正常运行。下面是排序命令的逐步解析及补充说明：

排序命令详解：$ time samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam

注释：

参数说明：

time

用于显示命令运行所花费的时间，方便评估效率。非必需，但建议使用。

samtools sort

samtools的排序子命令，专门用于对 BAM 文件进行排序。

-@ 4

指定排序使用的线程数，这里设置为 4 个线程。根据硬件性能调整以优化速度。

-m 4G

限制每个线程的内存消耗，设定为 4GB。总内存消耗约为4G × 线程数，需确保机器内存足够。

-O bam

指定输出格式为 BAM（二进制格式）。其他选项包括 SAM（文本格式）和 CRAM（更高效的压缩格式）。

-o sample_name.sorted.bam

输出文件的名称，建议采用有意义的命名格式（如添加sorted表明已排序）。

sample_name.bam

输入的 BAM 文件，即比对生成的未排序文件。

为什么需要排序？

位置无序：BWA 比对生成的 BAM 文件是基于 FASTQ 文件的读取顺序逐条比对的，导致记录的基因组位置是随机分布的。

后续分析需求：许多分析工具（例如samtools markdup、bcftools mpileup）都要求 BAM 文件按照基因组位置排序，才能高效处理比对结果。

数据处理优化：排序后的 BAM 文件在数据检索时更加高效，特别是结合索引文件（如 BAI）使用。

结果说明

排序生成的sample_name.sorted.bam文件通常比原始文件稍小，这是因为排序过程中压缩算法的优化。

使用samtools index为排序后的 BAM 文件生成索引文件（.bai格式）可以进一步提升数据检索和分析速度。

例如，生成索引的命令如下：

$ samtools index sample_name.sorted.bam

---

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

在完成排序后，去除重复（准确地说，是去除PCR 复制引入的重复序列）是下一步重要的处理环节。以下是关于重复序列的来源及其去除的详细说明。

什么是重复序列？

重复序列是指测序过程中产生的、定位到参考基因组上同一位置并具有相同起点的比对片段。这类重复通常不是因为生物学上的重复片段，而是因技术原因（如PCR 扩增）人为引入的。

重复序列的来源

DNA 片段化

测序前需要将 DNA 随机打断成较短的片段（通过超声波或酶切），形成测序文库。

打断后的片段浓度可能过低，无法满足测序需求。

PCR 扩增

为了增加片段浓度，从而确保每个片段都有足够的数量能够被测序仪读取，PCR 扩增被广泛应用。

PCR 的缺点是无法区分片段的初始浓度，因此所有片段都会被放大，导致高初始浓度片段可能被扩增过多。

取样测序

扩增后的 DNA 片段进入测序仪时，会随机取样。高密度的片段更容易被多次取样，从而生成大量冗余的重复序列。

为什么要去除重复序列？

减少偏倚

PCR 扩增过程可能会放大特定片段，导致这些片段在测序数据中占比过高，从而影响下游分析的准确性（例如突变检测和拷贝数变异分析）。

提高数据质量

去除重复能使数据更好地代表样本的真实情况，减少因技术原因导致的错误。

节省存储空间

重复片段会占用不必要的存储空间，增加计算资源的负担。

如何去除重复？

在排序完成的 BAM 文件基础上，常用工具如samtools markdup或Picard可以标记并去除重复。以下是一个示例命令：$ samtools markdup -r -s sample_name.sorted.bam sample_name.rmdup.bam

注释：

参数详解：

-r

删除重复片段。如果省略，则只标记重复片段而不删除。

-s

在标准输出中记录统计信息，便于分析重复比例。

sample_name.sorted.bam

排序后的输入 BAM 文件。

sample_name.rmdup.bam

去除重复后的输出 BAM 文件。

注意事项

非 PCR 引入的重复

某些生物学现象（例如基因组高重复区）也可能导致测序中出现重复片段。这种情况下，需根据分析需求决定是否保留。

多重比对和去重复

比对到多个位置的片段（如低复杂度区域）可能会干扰重复标记，需综合考虑是否去除。

重复率评估

去除重复后，建议用工具（如 FastQC 或 samtools stats）评估重复率，通常重复率在 10%-30% 是合理范围。

为什么需要去除 PCR 重复序列？

假变异的放大

在建库和打断步骤中，DNA 片段可能发生碱基改变（如嘌呤-嘧啶转换）。

PCR 扩增会对这些改变的片段放大，使它们在数据中显得比真实比例更多。最终，变异检测算法可能会误判这些为真实变异。

PCR 引入的新错误

PCR 扩增并非完美复制过程，会随机引入新的碱基错误。

早期扩增的错误会被不断放大，导致测序数据中伪变异比例升高。

PCR 偏好性（PCR Bias）

在扩增过程中，某些 DNA 模板可能被优先扩增。例如：

参考等位基因（reference allele）常因扩增效率更高而占据主导，稀有的变异等位基因信号被淹没，导致假阴性结果。

误导变异检测算法

主流变异检测工具（如GATK、Samtools、Platypus）假设输入的测序数据中，所有片段都是独立的，并具有随机覆盖性。

如果未去除 PCR 重复，这些工具可能会将高覆盖区域错误解释为生物学意义上的重复或变异。

结论：不去除 PCR 重复会显著提高假阳性和假阴性率，尤其是在变异检测、基因型判别和 RNA-Seq 等应用中。

如何去除 PCR 重复？

核心思想：

PCR 重复片段因来源于同一模板，其比对起始位置（start position）和比对方向（strand orientation）是相同的。通过对 BAM 文件的比对结果进行筛选，可以标记或删除这些重复片段。

方法实现

1.基于比对信息的去重

通常，使用工具如samtools markdup或Picard自动检测和去除重复。

示例命令（samtools markdup）：$ samtools markdup -r sample_name.sorted.bam sample_name.rmdup.bam

示例命令（Picard MarkDuplicates）：

$ java -jar picard.jar MarkDuplicates /

      I=sample_name.sorted.bam /

      O=sample_name.rmdup.bam /

      M=sample_name.metrics.txt /

      REMOVE_DUPLICATES=true

参数解析：

-r / REMOVE_DUPLICATES=true：删除标记为重复的片段。

M：生成重复标记的统计文件（如重复率）。

2.标记重复而不删除

如果需要保留重复信息供后续分析使用，可仅标记重复：

bash

复制代码

$ samtools markdup sample_name.sorted.bam sample_name.markdup.bam

确实，从直觉上看，重复序列似乎只是原始 DNA 的复制，保留也不会显著影响结果。然而，正如你所提到的，重复序列带来的问题是深层次的，会直接影响变异检测结果的准确性。下面，我们详细解答为什么必须去除重复，以及如何实现这一过程。

为什么需要去除 PCR 重复序列？

假变异的放大

在建库和打断步骤中，DNA 片段可能发生碱基改变（如嘌呤-嘧啶转换）。

PCR 扩增会对这些改变的片段放大，使它们在数据中显得比真实比例更多。最终，变异检测算法可能会误判这些为真实变异。

PCR 引入的新错误

PCR 扩增并非完美复制过程，会随机引入新的碱基错误。

早期扩增的错误会被不断放大，导致测序数据中伪变异比例升高。

PCR 偏好性（PCR Bias）

在扩增过程中，某些 DNA 模板可能被优先扩增。例如：

参考等位基因（reference allele）常因扩增效率更高而占据主导，稀有的变异等位基因信号被淹没，导致假阴性结果。

误导变异检测算法

主流变异检测工具（如GATK、Samtools、Platypus）假设输入的测序数据中，所有片段都是独立的，并具有随机覆盖性。

如果未去除 PCR 重复，这些工具可能会将高覆盖区域错误解释为生物学意义上的重复或变异。

结论：不去除 PCR 重复会显著提高假阳性和假阴性率，尤其是在变异检测、基因型判别和 RNA-Seq 等应用中。

如何去除 PCR 重复？

核心思想：

PCR 重复片段因来源于同一模板，其比对起始位置（start position）和比对方向（strand orientation）是相同的。通过对 BAM 文件的比对结果进行筛选，可以标记或删除这些重复片段。

方法实现

1.基于比对信息的去重

通常，使用工具如samtools markdup或Picard自动检测和去除重复。

示例命令（samtools markdup）：

bash

复制代码

$ samtools markdup -r sample_name.sorted.bam sample_name.rmdup.bam

示例命令（Picard MarkDuplicates）：

bash

复制代码

$ java -jar picard.jar MarkDuplicates /      I=sample_name.sorted.bam /      O=sample_name.rmdup.bam /      M=sample_name.metrics.txt /      REMOVE_DUPLICATES=true

参数解析：

-r / REMOVE_DUPLICATES=true：删除标记为重复的片段。

M：生成重复标记的统计文件（如重复率）。

2.标记重复而不删除

如果需要保留重复信息供后续分析使用，可仅标记重复：

bash

复制代码

$ samtools markdup sample_name.sorted.bam sample_name.markdup.bam

PCR-Free 测序方案

在高质量数据需求的场景下，可以采用PCR-Free 测序，直接跳过 PCR 扩增步骤。这种方法：

避免扩增引入的假信号。

减少偏倚，获得更真实的等位基因信息。

目前，如 Broad Institute 等机构的 RNA-Seq 或基因组测序工作中，已大量使用 PCR-Free 测序方案。

事实上，现有的工具，如Samtools和Picard，确实采用了相似的方法来去除重复序列。不同之处在于，Samtools的`rmdup`命令会直接从比对后的BAM文件中删除这些重复序列，而Picard的`MarkDuplicates`工具在默认情况下仅在BAM文件的FLAG信息中对重复序列进行标记，而不进行删除。因此，这些重复序列依然保留在文件中，但我们可以在进行变异检测时识别并忽略它们。