使用BWA、GATK、Samtools软件进行基因测序

背景信息

生命科学领域内基因测序技术的飞速发展，使得人类发现的基因序列以指数级增长，对于如此数量庞大的基因进行同源性搜寻、比对、变异检查等，往往伴随着巨大的数据处理和并行计算。高性能计算可以提供强大的算力支持，使用多种调度器提高并发效率，使用GPU进行计算加速等。

本文以经典及普及的二代全基因组测序WGS（Whole Genome Sequencing）流程为例，结合二代测序软件GATK，介绍人类全基因组测序的通用流程。在实际生信分析中，需要结合不同的业务需求进行相应的调整，如增加变异检测质控及过滤等，您可以结合实际场景进行调整。

在进行基因测序时，您可以使用BWA构建索引及比对记录，再使用Samtools对比对记录进行排序，然后使用GATK去除重复序列、重新校正碱基质量值、变异检查。

• BWA（Burrows-Wheeler-Alignment Tool）是一款将DNA序列映射到参考基因组上的软件，例如比对人类基因组。

• GATK（The Genome Analysis Toolkit）是一款二代重测序数据分析软件，是基因分析的工具集。主要用于去除重复序列、重新校正碱基质量值、变异检查等。

• Samtools是用于处理sam和bam格式的工具软件，能够查看二进制文件、转换文件格式、对文件排序及合并，可以结合sam格式文件中的flag、tag等信息，对比对结果进行统计汇总。

准备工作

1. 创建一个E-HPC集群。具体操作，请参见创建标准版集群。

   本文使用的集群配置示例如下：

   配置项	配置
   系列	标准版
   部署模式	公共云集群
   集群类型	SLURM
   节点配置	包含1个管理节点、1个登录节点和1个计算节点，规格如下： 管理节点：采用ecs.c7.large实例规格，该规格配置为2 vCPU，4 GiB内存。 登录节点：采用ecs.c7.large实例规格，该规格配置为2 vCPU，4 GiB内存。 计算节点：采用ecs.g7.16xlarge实例规格，该规格配置为64 vCPU、256 GiB内存。
   集群镜像	CentOS 7.6公共镜像

2. 创建集群用户。具体操作，请参见用户管理。

   集群用户用于登录集群，进行编译软件、提交作业等操作，本文创建的用户示例如下：

   • 用户名：testuser

   • 用户组：sudo权限组

3. 安装软件。

   待安装软件详情如下：

   软件名称	版本	下载地址
   BWA	4.2.0.0	BWA
   GATK	0.7.17	GATK
   Samtools	1.14	Samtools

步骤一：连接E-HPC集群

1. 选择任意一种方式连接E-HPC集群。具体操作，请参见连接集群。

2. 执行以下命令，创建作业执行目录，本文以/home/data/human为例。

   说明

   为确保有足够的存储空间用于下载和处理数据，建议将作业执行目录配置在集群挂载的共享存储目录/home下。

   sudo mkdir -p /home/data/human

步骤二：下载数据并进行预处理

1. 执行以下命令，创建并打开作业脚本文件，脚本文件命名为data_pre.slurm。

   sudo vim data_pre.slurm

   脚本内容示例如下：

   #!/bin/bash -l
   #SBATCH -J data_pre
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH --output=data_pre.log
   
   cd /home/data/human/
   
   # 下载并解压人类参考基因组、人类基因测序样本1、人类基因测序样本2
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_human_g1k_v37.fasta
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_1000G_phase1.indels.b37.vcf
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_Mills_and_1000G_gold_standard.indels.b37.vcf
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/b37_dbsnp_138.b37.vcf.zip
   unzip b37_dbsnp_138.b37.vcf.zip
   
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz
   gunzip gatk-examples_example1_NA20318_ERR250971_1.filt.fastq.gz
   
   wget https://public-ehpc-package.oss-cn-hangzhou.aliyuncs.com/lifescience/gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz
   gunzip gatk-examples_example1_NA20318_ERR250971_2.filt.fastq.gz
   
   # 使用bgzip对文件进行压缩
   bgzip -f b37_1000G_phase1.indels.b37.vcf
   bgzip -f b37_Mills_and_1000G_gold_standard.indels.b37.vcf
   bgzip -f b37_dbsnp_138.b37.vcf
   
   # 使用tabix建立索引
   tabix -p vcf  b37_1000G_phase1.indels.b37.vcf.gz
   tabix -p vcf  b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz
   tabix -p vcf  b37_dbsnp_138.b37.vcf.gz

2. 执行以下命令，提交作业。

   sbatch data_pre.slurm

   返回示例如下，表示已成功提交作业。

   

3. 执行以下命令，查看作业状态。

   squeue

4. 待作业执行完成后，执行以下命令，查看data_pre.log文件结果。

   cat data_pre.log

步骤三：构建参考基因组索引

1. 执行以下命令，创建并打开作业脚本文件，脚本文件命名为bwa_index.slurm。

   sudo vim bwa_index.slurm

   脚本内容示例如下：

   #!/bin/bash -l
   #SBATCH -J bwa_index
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH --output=bwa_index.log
   
   cd /home/data/human/
   
   # 构建BWA比对所需的参考基因组的index数据库
   time bwa index b37_human_g1k_v37.fasta
   
   # 创建fasta序列格式索引
   samtools faidx b37_human_g1k_v37.fasta

2. 执行以下命令，提交作业。

   sbatch --dependency=afterok:jobid1 bwa_index.slurm

3. 执行以下命令，查看作业状态。

   squeue

4. 待作业执行完成后，执行以下命令，查看bwa_index.log文件结果。

   cat bwa_index.log

步骤四：执行基因测序分析

1. 执行以下命令，创建并打开作业脚本文件，脚本文件命名为wgs_main_process.slurm。

   sudo vim wgs_main_process.slurm

   脚本内容示例如下：

   重要

   请将第六行#SBATCH -w compute000字段中compute000内容替换为集群实际运行节点。

   #!/bin/bash -l
   #SBATCH -J wgs_main_process
   #SBATCH --partition comp
   #SBATCH -N 1
   #SBATCH -n 64
   #SBATCH -w compute000
   #SBATCH --output=wgs_main_process.log
   
   cd /home/data/human/
   
   # 将样本测序数据reads与人类参考基因组进行比对，并将输出文件转化为bam格式，以节省磁盘空间
   time bwa mem -t 64 \
   -R '@RG\tID:ehpc\tPL:illumina\tLB:library\tSM:b37' \
   b37_human_g1k_v37.fasta \
   gatk-examples_example1_NA20318_ERR250971_1.filt.fastq \
   gatk-examples_example1_NA20318_ERR250971_2.filt.fastq | samtools view -S -b - > ERR250971.bam
   
   # 将比对记录按照顺序从小到大进行排序
   time samtools sort -@ 64 -O bam -o ERR250971.sorted.bam ERR250971.bam
   
   # 去除DNA序列PCR扩增产生的重复reads序列
   time gatk MarkDuplicates \
   -I ERR250971.sorted.bam -O ERR250971.markdup.bam \
   -M ERR250971.sorted.markdup_metrics.txt
   
   # 构建markdup测序bam文件索引
   samtools index ERR250971.markdup.bam
   
   # 生成参考基因组dict文件
   time gatk CreateSequenceDictionary  \
   -R b37_human_g1k_v37.fasta \
   -O b37_human_g1k_v37.dict
   
   # 利用已有的snp及indels数据库，建立相关模型，构建重校准表
   time gatk BaseRecalibrator \
   -R b37_human_g1k_v37.fasta \
   -I ERR250971.markdup.bam  \
   --known-sites b37_1000G_phase1.indels.b37.vcf.gz \
   --known-sites b37_Mills_and_1000G_gold_standard.indels.b37.vcf.gz \
   --known-sites b37_dbsnp_138.b37.vcf.gz -O ERR250971.BQSR.table
   
   # 对原始碱基质量值进行调整
   time gatk ApplyBQSR \
   --bqsr-recal-file ERR250971.BQSR.table \
   -R b37_human_g1k_v37.fasta \
   -I ERR250971.markdup.bam \
   -O ERR250971.BQSR.bam
   
   # 构建BQSR测序bam文件索引
   time samtools index ERR250971.BQSR.bam
   
   # 输出变异检测结果vcf文件
   time gatk HaplotypeCaller \
   -R b37_human_g1k_v37.fasta  \
   -I ERR250971.BQSR.bam \
   -O ERR250971.HC.vcf.gz

2. 执行以下命令，提交作业。

   sbatch --dependency=afterok:jobid2 wgs_main_process.slurm

3. 执行以下命令，查看作业状态。

   squeue

4. 待作业执行完成后，执行以下命令，查看wgs_main_process.log文件结果。

   cat wgs_main_process.log

测序耗时说明

本文使用的计算节点为ecs.g7.16xlarge，仅为演示全基因组测序流程，整体大约耗时6~7小时，下表展示了基因测序各流程大约使用的时间供您参考，并非最优性能。