使用Spark进行WGS分析

发表于 更新于 分类于 阅读次数： 本文字数： 4.2k 阅读时长 ≈ 4 分钟

GATK简介

GATK全称叫做: Genome Analysis Toolkit. 是Broad Institute开发的用于二代重测序数据分析的一款软件.
目前主要用于人类的WGS以及WES基因测试流程, 具体流程介绍可以看官网的最佳实践

GATK3版本之前, 一直都是单机版本, 性能一直是瓶颈点, 做完一个WGS的流程大约需要3天时间. 因此在GATK4以后的版本之中, 引入Spark做分布式性能优化, GATK4.0版本可以讲整个WGS测序流程的时间压缩在半天之内, 性能提高将近10倍有余.

但是, 目前所有标注有Spark加速的工具都是BETA Tool, 虽然就我们测试来看敏感度和准确性都和单机版本没有太大区别, 但是由于整理功能开发阶段, 工具接口可能会调整, 因此如果想应用到生产系统上的话, 也请慎重选择.

WGS: Whole Genome Sequencing 全基因组测序

WES: Whole Exome Sequencing 全外显子测测序

WGS流程简介

在GATK的最佳实践里面, 有流程介绍, 也有样例程序供大家参考

但是如果大家之前没有接触过WGS的话, 看官网的介绍还是有点晕. 推荐看一下碱基矿工的GATK4.0和全基因组数据分析实践

好了, 言归正传, 我在这儿简单总结一下WGS的流程:

  1. 获取数据 — 脱机数据转化成FastQ格式
  2. 数据质控 — 使用Fastqc工具过滤掉低质量的数据
  3. 比对排序 — 使用Bwa + samtools工具对FastQ进行比对排序, 并将格式转化为Bam格式
  4. 碱基去重 — 使用GATK的MarkDuplicates工具完成该步骤
  5. 碱基矫正 — 使用GATK的BQSR工具完成该步骤
  6. 变异检测 — 使用GATK的HaplotypeCaller工具完成该步骤
  7. 变异控制 — 使用GATK的VQSR工具完成该步骤

最后我们实现的一个功能是将原始的FastQ个数的数据, 转化为VCF格式的数据, 完成整个WGS的流程.

VCF(Variant Call Format): 是种文本格式变异数据的格式, 可以直接用文本编辑器查看里面的字段及数据, 字段含义如文所示

使用Spark进行WGS分析流程

流程介绍

使用Spark进行WGS的流程如下图所示:

  1. 首先FastQ格式的原始数据,通过FastqToSam工具转化为UBam格式
  2. 接着使用BwaSpark方法进行比对, 输出经过比对的Bam格式数据
  3. 最后通过ReadsPipelineSpark进行变异检测,并将变异点输出为VCF格式

下面讲详细介绍每个工具的用法及命令

FastqToSam

官网简介

1
2
3
4
5
Converts a FASTQ file to an unaligned BAM or SAM file.

Output read records will contain the original base calls and quality scores will be translated depending on the base quality score encoding: FastqSanger, FastqSolexa and FastqIllumina.

There are also arguments to provide values for SAM header and read attributes that are not present in FASTQ (e.g see RG or SM below).

这个工具的作用就是做好格式转化, 并对BAM格式进行排序. 这个工具是单机的, 无法使用Spark加速. 官方工具转化一个NA1278 30X的样本大致需要3-4个小时. 我们团队这个工具进行了优化, 使得转化时间提高到1.5小时,大大降低了样本转化的时间

执行命令

1
/gatk/gatk FastqToSam -F1 NA12878_1.fastq.gz -F2 NA12878_2.fastq.gz -O NA12878_unaligned.bam -SM SM1 -PL illumina -R hg19.fa

FastqToSam文档

BwaSpark

简介

1
Align reads to a given reference using BWA on Spark

这个工具本质上是使用HDFS分片的能力, 让Spark对BWA软件分布化. Spark的每个Task都比对一个块大小的uBam. 每个块大小由参数--bam-partition-size指定, 默认值是使用Hadoop默认的块大小.
通过Spark的分布式话, 原有的比对时间从4小时,可以降低到1个小时左右.

执行命令

1
/gatk/gatk BwaSpark -I hdfs://hacluster/gatk/NA12878_unaligned.bam -O hdfs://hacluster/gatk/NA12878_aligned.bam -R hdfs://hacluster/gatk/hg19.2bit --spark-runner SPARK --spark-master yarn-cluster

BwaSpark文档

ReadsPipelineSpark

简介

1
Takes unaligned or aligned reads and runs BWA (if specified), MarkDuplicates, BQSR, and HaplotypeCaller to generate a VCF file of variants

在之前的GATK版本之中, 每个命令都是单独的, 只能通过自己编写脚本的方式讲这些工具集串行起来.
而在GATK4.0之中, 想用ReadsPipelineSpark这一个工具, 将整个变异流程全部统一起来, 未来变异检测只要执行这个工具即可.

实际上ReadsPipelineSpark的流程能够包括BwaSpark的步骤, 但为什么还要特别分开呢? 原因在于: 目前ReadsPipelineSpark的实现之中没有缓存必要的RDD, 导致重计算, 整个性能没有分开计算好. 因此目前分成了两个步骤, 等社区解决了这个问题之后, 可以讲这两个步骤合并. 这也符合社区对于这个工具的定位.

执行命令

1
/gatk/gatk ReadsPipelineSpark -I hdfs://hacluster/gatk/NA12878_aligned.bam -O hdfs://hacluster/gatk/NA12878.vcf -R hdfs://hacluster/gatk/hg19.2bit --known-sites hdfs://hacluster/gatk_ref/dbsnp.vcf --spark-runner SPARK --spark-master yarn-cluster

ReadsPipelineSpark文档

总结

GATK4.0 在今年初引入了Spark进行分布式化性能优化, 整个WGS流程的性能由原有的天级别降低到小时级别, 性能得到极大优化.

后续将从这个流程输出,分析一下目前使用Spark之中的性能瓶颈点以及源码级的分析

附录

数据集获取

Input

输入的数据可以有两种, 一种是FASTQ格式, 另外一种是BAM格式

NA12878是人类基因测试里面最常用的实验数据

reference

1
2
3
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.64.alt

knownsites

1
2
3
4
5
6
7
8
## 获取SNP的Knownsites: 1000G_phase1.snps.high_confidence
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi

## 获取Indel的Knownsites: Mills_and_1000G_gold_standard.indels
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi