使用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等信息,对比对结果进行统计汇总。
准备工作
操作步骤
登录E-HPC集群。
登录时,请使用具有sudo权限的用户。具体操作,请参见登录集群。
安装BWA、GATK、Samtools软件。
下载并解压人类参考基因组、人类基因测序样本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
构建索引。
构建bwa比对所需的参考基因组的index数据库。
bwa index b37_human_g1k_v37.fasta
创建fasta序列格式索引。
samtools faidx b37_human_g1k_v37.fasta
比对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
将比对记录按照顺序从小到大进行排序。
time samtools sort -@ 52 -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
重新校正碱基质量值。
利用已有的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
变异检查。
生成参考基因组dict文件。
time gatk CreateSequenceDictionary \ -R b37_human_g1k_v37.fasta \ -O b37_human_g1k_v37.dict
输出变异检测结果vcf文件。
time gatk HaplotypeCaller \ -R b37_human_g1k_v37.fasta \ -I ERR250971.BQSR.bam \ -O ERR250971.HC.vcf.gz
测序耗时说明
该测序流程可结合实际场景作进一步的优化,以提升性能。例如,在染色体切割时结合调度器提高并发率,通过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 |