阿里云首页 弹性高性能计算E-HPC

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

本文以BWA、GATK、Samtools软件为例介绍如何使用E-HPC集群进行基因测序计算。

背景信息

人类对生命真相的探索与追问从未停歇,而求索过程,必然离不开科学技术的支持。生命科学领域内基因测序技术的飞速发展,人类发现的基因序列以指数级增长,对于如此数量庞大的基因进行同源性搜寻、比对、变异检查等,往往伴随着巨大的数据处理和并行计算。高性能计算可以提供强大的算力支持,使用多种调度器提高并发效率,使用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. 创建并登录集群。

    1. 登录弹性高性能计算控制台

    2. 创建一个名为Genome的集群。

      具体操作,请参见创建集群。请注意以下配置参数:

      • 计算节点:选择ecs.ebmc5s.24xlarge

      • VNC:打开VNC开关,打开后可以自动部署远程可视化窗口。

    3. 创建一个名为Genome的普通用户。

      具体操作,请参见创建用户

    4. 集群页面,找到Genome集群,单击远程连接

    5. 远程连接页面,输入Genome的用户名、密码和节点端口,单击ssh连接

  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 最佳实践 使用BWA、GATK、Samtools软件进行基因测序