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

背景信息

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

操作步骤

1. 创建并登录集群。

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

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

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

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

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

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

      具体操作，请参见创建用户。

   4. 在集群页面，找到Genome集群，单击远程连接。

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

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

测序耗时