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

背景信息

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

准备工作

1. 创建E-HPC集群。具体操作，请参见使用向导创建集群。

   配置集群时，软硬件参数配置如下：

   参数	说明
   硬件参数	部署方式为精简，包含1个管控节点和1个计算节点，规格如下： 管控节点：采用ecs.c7.large实例规格，该规格配置为2 vCPU，4 GiB内存。 计算节点：采用ecs.ebmc5s.24xlarge实例规格，该规格配置为96 vCPU、192 GiB内存。
   软件参数	镜像选择CentOS 7.6公共镜像，调度器选择pbs，打开VNC开关。

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

   集群用户用于登录集群，进行编译软件、提交作业等操作，配置用户权限时，权限组请选择sudo权限组。

操作步骤

1. 登录E-HPC集群。

   登录时，请使用具有sudo权限的用户。具体操作，请参见登录集群。

2. 安装BWA、GATK、Samtools软件。

   关于如何安装BWA、GATK、Samtools，请参见BWA、GATK、Samtools。

3. 下载并解压人类参考基因组、人类基因测序样本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/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

4. 构建索引。

   1. 构建bwa比对所需的参考基因组的index数据库。

      bwa index b37_human_g1k_v37.fasta

   2. 创建fasta序列格式索引。

      samtools faidx b37_human_g1k_v37.fasta 

5. 比对reads。

   将样本测序数据reads与人类参考基因组进行比对，并将输出文件转化为bam格式，可有效节省磁盘空间。

   time bwa mem -t 52 \ -R '@RG\tID:ehpc\tPL:illumina\tLB:library\tSM:b37' b37_human_g1k_v37.fasta \
   NA20318_ERR250971_1.filt.fastq NA20318_ERR250971_2.filt.fastq \
   | samtools view -S -b - > ERR250971.bam

6. 将比对记录按照顺序从小到大进行排序。

   time samtools sort -@ 52 -O bam -o ERR250971.sorted.bam ERR250971.bam

7. 去除重复序列。

   1. 去除DNA序列PCR扩增产生的重复reads序列。

      time gatk MarkDuplicates \
      -I ERR250971.sorted.bam -O ERR250971.markdup.bam \
      -M ERR250971.sorted.markdup_metrics.txt

   2. 构建markdup测序bam文件索引。

      samtools index ERR250971.markdup.bam

8. 重新校正碱基质量值。

   1. 利用已有的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

   2. 对原始碱基质量值进行调整。

      time gatk ApplyBQSR \
      --bqsr-recal-file ERR250971.BQSR.table \
      -R b37_human_g1k_v37.fasta \
      -I ERR250971.markdup.bam \
      -O ERR250971.BQSR.bam

   3. 构建BQSR测序bam文件索引。

      time samtools index ERR250971.BQSR.bam

9. 变异检查。

   1. 生成参考基因组dict文件。

      time gatk CreateSequenceDictionary  \
      -R b37_human_g1k_v37.fasta \
      -O b37_human_g1k_v37.dict

   2. 输出变异检测结果vcf文件。

      time gatk HaplotypeCaller \
      -R b37_human_g1k_v37.fasta  \
      -I ERR250971.BQSR.bam \
      -O ERR250971.HC.vcf.gz

测序耗时说明