GATK4 × WGS

今天分享一篇笔记，主要内容是WGS重测序下机数据fq.gz文件到变异位点信息vcf文件的操作过程，包括序列质控过滤、基因组比对、检测变异、合并VCF等步骤。

本篇笔记只涉及代码步骤，均已测试过，不提供示例数据，有需要可以用公开发表数据进行学习。

关键词： WGS、GATK4、samtools、vcftools、fastp、picard。

什么是WGS重测序？

重测序是指在基因组学研究中，通过对DNA或RNA序列的测定，确定基因组或转录组的顺序。

通过重测序，我们可以揭示基因组中的突变、变异和重排等变化，从而深入研究生物体的遗传特征和功能。重测序技术的应用广泛，包括基因组测序、转录组测序、外显子测序等，为生物学研究、医学诊断和个体化治疗等领域提供了重要的工具和数据基础。

分析流程

首先建立一个项目文件夹，将原始的测序数据放在数据目录下，通常采用二代高通量测序，每个样品下机数据为两个fq文件，成对存在，以下代码流程需要根据自己的实际需要进行修改，此处只作为示例，部分方法借鉴了网上公开信息和软件手册。

参数设置

● sample_name_1:样品一名称或者代号，例如A
● sample_name_2:样品二名称或者代号，例如B
● INPUT_S1_1:样品一的第一条测序fq.gz文件名称，例如Rawdata1_1.fq.gz
● INPUT_S1_2:样品一的第二条测序fq.gz文件名称，例如Rawdata1_2.fq.gz
● INPUT_S2_1:样品二的第一条测序fq.gz文件名称，例如Rawdata2_1.fq.gz
● INPUT_S2_2:样品二的第二条测序fq.gz文件名称，例如Rawdata2_2.fq.gz
● OUTPUT:输出结果文件名称，例如xxxx.vcf.gz

设置样品名称和具体的文件路径信息，后续分析将会读取这些信息，避免重复输入修改，提高效率。

原始数据质控过滤

dir="/Data/"
fastp -i ${dir}/01_rawdata/Unknown_1.fq.gz \
      -I ${dir}/01_rawdata/Unknown_2.fq.gz \
      -o ${dir}/02_cleardata/good_1.fq.gz \
      -O ${dir}/02_cleardata/good_2.fq.gz \
      -h R1.html -j R1.json --thread 23

这段代码使用 fastp 程序，用于对测序数据进行质量控制（QC）和预处理，以下是代码中每个部分的简要解释：

-i ${dir}/01_rawdata/Unknown_1.fq.gz：这是输入文件参数，指定了第一对（R1）原始测序数据文件的路径。${dir} 是一个变量，表示某个目录的路径。
-I ${dir}/01_rawdata/Unknown_2.fq.gz：这是第二对（R2）原始测序数据文件的路径。
-o ${dir}/02_cleardata/good_1.fq.gz：这是输出文件参数，指定了处理后的第一对（R1）数据的输出路径。数据经过质量控制和预处理后，将被保存在这个文件中。
-O ${dir}/02_cleardata/good_2.fq.gz：这是处理后的第二对（R2）数据的输出路径。
-h R1.html：生成质量报告的 HTML 文件的名称，这里命名为 R1.html，其中记录了质量控制的各种统计信息和图形。
-j R1.json：生成包含质量控制信息的 JSON 格式文件的名称，这里命名为 R1.json。
--thread 23：指定线程数，表示程序将使用 23 个线程来加速处理。

主要原理是在指定的输入路径中找到原始测序数据文件，使用 fastp 对这些数据进行质量控制和预处理，然后将处理后的数据和质量统计信息保存到指定的输出路径中。生成的质量报告以 HTML 和 JSON 格式保存，程序会使用 23 个线程来加速处理过程。

参考基因组序列比对

bwa mem -t 23 -R '@RG\tID:R1\tSM:R1\tPL:illumina' \
   /NGS/Ref.fa \
   /Data/02_cleardata/good_1.fq.gz \
   /Data/02_cleardata/good_2.fq.gz \
   > /03_bwa_sam/R1.mem.sam

这段代码使用了 bwa mem 工具来执行基因组比对，其主要功能是将测序数据映射到参考基因组上。以下是代码的简要解释：

bwa mem：这是一个用于执行基因组比对的命令行工具。它通过将测序数据映射到参考基因组上来确定测序数据的起源位置。
-t 23：这个参数指定线程数，即使用 23 个线程来加速比对过程。
-R '@RG\tID:R1\tSM:R1\tPL:illumina'：这是一个描述测序数据的元数据标签。其中，@RG 表示"Read Group"，ID:R1 表示读组的唯一标识符，SM:R1 表示样本名称，PL:illumina 表示测序平台为 Illumina。
/NGS/Ref.fa：这是参考基因组的路径。程序将测序数据与这个基因组进行比对。
/Data/02_cleardata/good_1.fq.gz 和 /Data/02_cleardata/good_2.fq.gz：这是经过质量控制和预处理的测序数据的路径。这些数据将与参考基因组进行比对。
> /03_bwa_sam/R1.mem.sam：这个符号（>）将比对的结果输出到指定路径下的文件中，这里是 /03_bwa_sam/R1.mem.sam。生成的文件将是 SAM 格式，其中包含了测序数据与参考基因组的比对信息。

使用 bwa mem 工具将经过质控和预处理的测序数据与参考基因组进行比对。比对的结果将以 SAM 格式保存在指定路径下的文件中，同时使用 23 个线程来加速比对过程，元数据标签描述了测序数据的一些信息，如样本名称、测序平台等。

bam排序转sam文件

#!/usr/bin/bash
dir="/work/" # 这行定义了一个名为 dir 的变量，存储路径为 "/work/"，表示工作目录的路径。
samtools sort -@20 -m 2G -o R1.sort.bam R1.mem.sam
samtools sort -@20 -m 2G -o R2.sort.bam R2.mem.sam
# 使用 samtools 工具对名为 "R1.mem.sam" 的 SAM 格式文件进行排序，并将排序后的结果保存为 "R1.sort.bam" 的 BAM 文件。选项 -@20 指定了线程数为 20，-m 2G 指定了单个线程的最大内存为 2GB。
# @为线程数，m为单个线程的内存，因此需要注意20×2G=40G不能超过机器最大限制

这段代码旨在将经过比对的测序数据文件进行排序并生成 BAM 文件，利用 samtools 工具对两个经过比对的 SAM 文件进行排序，并生成对应的 BAM 文件。通过设置线程数和单线程内存限制，脚本确保了排序过程在给定的资源限制下进行，避免了资源超限问题，排序结果将分别保存为 "R1.sort.bam" 和 "R2.sort.bam" 文件。

picard去除重复

picard -Xmx100g  MarkDuplicates \ # 这个参数设置了工具的最大可用内存为 100GB，用于处理大规模数据
   I=$input1 O=$out1  \
   CREATE_INDEX=false  REMOVE_DUPLICATES=true \
   M=${out1}.marked_dup_metrics.txt

使用 Picard 工具的 MarkDuplicates 功能来处理输入的 BAM 文件。它标记并移除其中的重复测序数据，并生成一个输出 BAM 文件，同时还会生成一个包含标记重复度度量信息的文本文件，处理过程中分配最大 100GB 的内存。

samtools建立索引

#!/usr/bin/bash
dir="/work/"
input1="R1.sort.rmdup.bam"
input2="R2.sort.rmdup.bam"
# -c 这个选项表示要为 BAM 文件创建一个 CRAM 格式的索引。CRAM 是一种高效的压缩格式，用于存储比对数据。
samtools index -c ${dir}file_sam_rmsort/${input1} &
samtools index -c ${dir}file_sam_rmsort/${input2} &

使用 samtools 工具来为一个已经排序的 SAM 或 BAM 文件创建索引。

GATK4变异检测

#!/usr/bin/bash
dir="/work/"
gatk="/software/GATK/gatk-4.4.0.0/gatk"
ref="/Ref/ref.fa"
input_1=${dir}file_sam_rmsort/R1.sort.rmdup.bam
input_2=${dir}file_sam_rmsort/R2.sort.rmdup.bam
out_1=${dir}file_gvcf/R1.g.vcf.gz
out_2=${dir}file_gvcf/R2.g.vcf.gz
REF=${ref}
SM1="R1"
SM2="R2"
cd $dir
mkdir tmpdir
chroms=($(grep '>' $REF |sed 's/>//' | tr '\n' ' '))

这段脚本首先定义了工作目录、GATK路径、参考基因组文件等各项参数，然后针对两个输入的已去重且排序的BAM文件，使用GATK工具，对每个样本分别生成带有变异信息的gVCF格式文件。

生成的gVCF文件将在工作目录下的“file_gvcf”文件夹中，分别命名为“R1.g.vcf.gz”和“R2.g.vcf.gz”，包含了每个样本在每个染色体上的变异信息。

此外，脚本还创建了一个名为“tmpdir”的临时目录用于中间文件，并获取了参考基因组中的染色体信息。

最终，脚本将生成的gVCF文件存储在指定的输出路径中，然后依次进行如下子步骤处理：

HaplotypeCaller

for chr in ${chroms[@]}
do
        if [ ! -f ./file_gvcf2/${SM1}.${chr}.vcf.gz ]; then
        ${gatk} --java-options "-Xmx50g -XX:ParallelGCThreads=4 -Djava.io.tmpdir=./tmpdir" \
HaplotypeCaller -R ${ref} -I ${input_1} -OVI False \
--intervals ${chr}  --native-pair-hmm-threads 4 -ERC GVCF -O ./file_gvcf2/${SM1}.${chr}.vcf.gz \
1>./file_gvcf2/${SM1}.${chr}.vcf.gz.log 2>&1 &
        fi
done

这段循环代码遍历给定的染色体列表，对每个染色体使用GATK工具的HaplotypeCaller功能进行变异检测，生成每个样本的gVCF格式变异文件。

CombineGVCFs

for chr in ${chroms[@]}; do
        ${gatk} --java-options \
        "-Xmx50g -XX:ParallelGCThreads=4" \
        CombineGVCFs -R ${ref} \
        -V ${dir}file_gvcf2/${SM1}.${chr}.vcf \
        -V ${dir}file_gvcf2/${SM2}.${chr}.vcf \
        -O ${dir}file_chr_vcf/${chr}.g.vcf > /dev/null &
done

此行代码通过${gatk}来调用GATK工具CombineGVCFs功能，使用--java-options参数设置Java虚拟机的选项。其中-Xmx50g指定Java最大可用内存为50GB，-XX:ParallelGCThreads=4设置并行垃圾回收线程数为4。

GenotypeGVCFs

for chr in ${chroms[@]}; do
        ${gatk} --java-options 、\
        "-Xmx50g -XX:ParallelGCThreads=4" \
        GenotypeGVCFs -R ${ref} \
        -V ${dir}file_chr_vcf/${chr}.g.vcf \
        -O ${dir}file_allvcf/${chr}.vcf > /dev/null &
done

接着，调用GenotypeGVCFs功能，以给定的参考基因组${ref}和已处理的gVCF文件${dir}file_chr_vcf/${chr}.g.vcf为输入，生成对应染色体${chr}的VCF格式变异文件，结果存储在${dir}file_allvcf/${chr}.vcf路径下。/dev/null将标准输出重定向为空，&将任务放入后台运行。

Merge合并VCF

${gatk} MergeVcfs \
$(for i in ${chroms[@]}; \
do \
   echo "-I ${dir}file_allvcf/${i}.vcf "; \
done) \
-O ${dir}file_allvcf/final_Yr45.vcf

最后调用MergeVcfs合并VCF文件，得到最终的结果，到此流程结束！感谢您的阅读，以下补充一点儿重测序数据分析软件安装和环境部署的教程。

如何安装分析所需的软件？

在Linux系统中安装和部署miniconda3，并创建WGS重测序分析环境，并安装GATK4、samtools、vcftools、fastp、picard软件包，可以按照以下步骤进行操作：

1. 下载miniconda3安装包：
   wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
2. 安装miniconda3：
   bash Miniconda3-latest-Linux-x86_64.sh
3. 创建WGS重测序分析环境：
   conda create -n wgs_analysis
4. 激活环境：
   conda activate wgs_analysis
5. 安装GATK4：
   conda install -c bioconda gatk4
6. 安装samtools：
   conda install -c bioconda samtools
7. 安装vcftools：
   conda install -c bioconda vcftools
8. 安装fastp：
   conda install -c bioconda fastp
9. 安装picard：
   conda install -c bioconda picard

行文至此，暂时收笔，如果喜欢本篇内容欢迎点赞转发，如果您也对重测序分析感兴趣，欢迎扫描下方二维码加入生信分析交流群，免费交流广告禁入。