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
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

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

GATK4重测序数据怎么分析？

GATK4 × WGS

什么是WGS重测序？

分析流程

参数设置

原始数据质控过滤

参考基因组序列比对

bam排序转sam文件

picard去除重复

samtools建立索引

GATK4变异检测

HaplotypeCaller

CombineGVCFs

GenotypeGVCFs

Merge合并VCF

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

热门文章

最新文章

相关课程

相关电子书

相关实验场景

热门

活动广场

任务中心

开发者评测

高校计划

乘风者计划

训练营

阿里云MVP

话题

直播

下载

镜像站

技术资料

插件

GATK4重测序数据怎么分析？

GATK4 × WGS

什么是WGS重测序？

分析流程

参数设置

原始数据质控过滤

参考基因组序列比对

bam排序转sam文件

picard去除重复

samtools建立索引

GATK4变异检测

HaplotypeCaller

CombineGVCFs

GenotypeGVCFs

Merge合并VCF

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

热门文章

最新文章

相关课程

相关电子书

相关实验场景