首页 >弹性高性能计算E-HPC >最佳实践 >使用BWA、GATK、Samtools软件进行基因测序

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

本文介绍如何使用E-HPC集群运行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集群。具体操作,请参见使用向导创建集群

    配置集群时,软硬件参数配置如下:

    参数

    说明

    硬件参数

    部署方式为精简,包含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,请参见BWAGATKSamtools

  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

测序耗时说明

本文使用的计算节点为ecs.ebmc5s.24xlarge,仅为演示全基因组测序流程,下表展示了基因测序各流程使用的时间,并非最优性能。
说明

该测序流程可结合实际场景作进一步的优化,以提升性能。例如,在染色体切割时结合调度器提高并发率,通过GATK分布式并发测序提高测序效率,结合FPGA加速测序等。如您有业务需求,请提交工单

流程

功能

运行时间(min)

bwa index

索引

50.60

bwa mem

比对

21.98

sort

排序

1.62

MarkDuplicates

去重

21.60

BaseRecalibrator

构建校准表

38.47

ApplyBQSR

重校准

17.22

CreateSequenceDictionary

创建字典

0.18

HaplotypeCaller

变异检查

175.93
阿里云首页 弹性高性能计算E-HPC 相关技术圈