构建生信分析流程是生物信息学从业人员必备的技能之一,对该项能力的评估常常是各大公司招录人员的参考项目之一。
在进行 ngsjs 项目时,我做了一张示意图来表示一些高通量测序数据分析项目重现性的要点(图一)。
ngsjs: A set of command line tools, NGS data analysis workflows, and R shiny plugins/R markdown document for exploring next-generation sequencing data.
一个好的生物信息分析流程可以让你事倍功半,有效减负,同时也有利于他人重复你的数据分析结果。
图一 高通量测序数据分析项目重现性的要点
其中,使用统一的管道(pipeline)、工作流程(workflow)就是其中最重要的一环。
根据生信信息学数据分析流程(管道、工作流程序)构建的风格和方式,大致有以下几大流派(注1):
注1: 围棋比赛中常用“流派”来形容各具特色的棋手/派系,参见资料,深度学习派代表 AlphaGo 现在在围棋届已经一骑绝尘。
脚本语言流的主要是通过简单的脚本语言(如 shell,R,Python,Perl)运行各类命令行脚本/程序。常见的几种工作模式:
前两种(1 和 2)是大多数生物信息学初学者(不具备封装和打包能力)最早开始接触生信分析流程的方式。后两种(3 和 4)是专业人员开发新工具、新流程的必备技能。
这里使用 htslib.org 所给WGS/WES Mapping to Variant Calls (v1.0)
作为工作模式1,2
的示例(已略去注释):
# mapping
bwa index <ref.fa>
bwa mem -R '@RG\tID:foo\tSM:bar\tLB:library1' <ref.fa> <read1.fa> <read1.fa> > lane.sam
samtools fixmate -O bam <lane.sam> <lane_fixmate.bam>
samtools sort -O bam -o <lane_sorted.bam> -T </tmp/lane_temp> <lane_fixmate.sam>
# Improvement
java -Xmx2g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R <ref.fa> -I <lane.bam> -o <lane.intervals> --known <bundle/b38/Mills1000G.b38.vcf>
java -Xmx4g -jar GenomeAnalysisTK.jar -T IndelRealigner -R <ref.fa> -I <lane.bam> -targetIntervals <lane.intervals> --known <bundle/b38/Mills1000G.b38.vcf> -o <lane_realigned.bam>
java -Xmx4g -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R <ref.fa> -knownSites >bundle/b38/dbsnp_142.b38.vcf> -I <lane.bam> -o <lane_recal.table>
java -Xmx2g -jar GenomeAnalysisTK.jar -T PrintReads -R <ref.fa> -I <lane.bam> --BSQR <lane_recal.table> -o <lane_recal.bam>
java -Xmx2g -jar MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT INPUT=<lane_1.bam> INPUT=<lane_2.bam> INPUT=<lane_3.bam> OUTPUT=<library.bam>
samtools merge <sample.bam> <library1.bam> <library2.bam> <library3.bam>
samtools index <sample.bam>
java -Xmx2g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R <ref.fa> -I <sample.bam> -o <sample.intervals> --known >bundle/b38/Mills1000G.b38.vcf>
java -Xmx4g -jar GenomeAnalysisTK.jar -T IndelRealigner -R <ref.fa> -I <sample.bam> -targetIntervals <sample.intervals> --known >bundle/b38/Mills1000G.b38.vcf> -o <sample_realigned.bam>
samtools index <sample_realigned.bam>
# Variant Calling
bcftools mpileup -Ou -f <ref.fa> <sample1.bam> <sample2.bam> <sample3.bam> | bcftools call -vmO z -o <study.vcf.gz>
bcftools mpileup -Ob -o <study.bcf> -f <ref.fa> <sample1.bam> <sample2.bam> <sample3.bam>
bcftools call -vmO z -o <study.vcf.gz> <study.bcf>
tabix -p vcf <study.vcf.gz>
bcftools stats -F <ref.fa> -s - <study.vcf.gz> > <study.vcf.gz.stats>
mkdir plots
plot-vcfstats -p plots/ <study.vcf.gz.stats>
bcftools filter -O z -o <study_filtered..vcf.gz> -s LOWQUAL -i'%QUAL>10' <study.vcf.gz>
工作模式3
和4
是开发生物信息工具的标准方式。很多生物信息学分析工具本身就是涉及到各类复杂计算和命令的集成式流程,如:
使用和开发这类工具/流程的主要原因:
我目前主要是 R 语言、Python 写命令行程序、函数、R 包/模块,同时用 CRAN、PyPI 以及 GitHub 分发。
同时,因为 R 语言目前还没有提供一个原生机制直接部署命令行可执行程序(Python、Node包均提供),我现在做了两手准备:
rbin
函数、以及 ngsjs 增加rbin
命令行程序一键收集 R 包中inst/bin
下面的文件。目前的体验来看,JavaScript 提供的 npm 和 yarn 包管理工具速度非常快和方便,很适合 R 语言用户同时使用(只需要会写一个 package.json 文件即可)。
Common Workflow language 语言流是近几年兴起的,专门用于数据分析流程构建的一类语言/工具。
这类语言/工具最核心的部分:定义每一个计算过程(脚本)的输入和输出,然后通过连接这些输入和输出,构成数据分析流程(图二,图三)(如 Galaxy, wdl,cromwell,nextflow,snakemake,bpipe 等)。
图二 CWL流程示意图一
图三 CWL 流程示意图二
示例项目:
使用和开发这类工具的主要原因:
扩展阅读:
Makefile 流主要是基于软件开发工具 Makefile 的 rule、target 语法运行流程。在 snakemake 工具出现之后(使得数据分析流程支持 CWL),使用Makefile
式 Rule 文件构建生物信息学分析流程的用户迅速增加。
pyflow-ATACseq 项目提供的 ATAC-seq 数据分析流程:
图五 ATAC-seq Snakemake 示例流程图
snakemake 示例文件:
rule targets:
input:
"plots/dataset1.pdf",
"plots/dataset2.pdf"
rule plot:
input:
"raw/{dataset}.csv"
output:
"plots/{dataset}.pdf"
shell:
"somecommand {input} {output}"
配置文件流(和 CWL 不冲突)主要是基于 JSON、YAML、TOML 等类型的配置文件,然后开发相应的解析器解析和执行流程。如 Galaxy、华为公司最近开源的 Kubegene(基于谷歌开发并开源的容器调度技术 kubernetes)、bashful 的流程文件。
很多计算机软件自动测试流程和构建工具也主要基于配置文件来构建和执行:如 circleci、travis。
这里给出一个基于配置文件的工具示例(图六):
图六 bashful 执行输出
bashful 输入文件格式及部分字段:
config:
show-failure-report: false
show-summary-errors: true
max-parallel-commands: 6
show-task-times: true
x-reference-data:
all-apps: &app-names
- some-lib-4
- utilities-lib
- important-lib
- some-app1
- some-app3
tasks:
- name: Cloning Repos
parallel-tasks:
- name: "Cloning <replace>"
cmd: example/scripts/random-worker.sh 2 <replace>
ignore-failure: true
for-each: *app-names
- name: Building Repos
parallel-tasks:
- name: "Building <replace>"
cmd: example/scripts/random-worker.sh 1 <replace>
ignore-failure: true
for-each: *app-names
复杂度较高的生物信息学流程/工具一般至少会提供一个配置文件来管理参数。用户目前也大多接受使用配置文件统一管理变量。
命令行参数也常常结合配置文件同时使用,这么做的主要原因:
Jupyter notebook 和 R markdown 分别由 Python 语言和 R 语言社区贡献,均可以用于整合文档、代码、以及代码的输出,构建动态、交互式文档和报告系统。
这两个工具已经风靡全世界的数据科学社区,同时也占据了生物信息分析流程中的下游统计分析、建模、以及可视化。
这两个工具兴起的主要原因:
Jupyter notebook 示例:
图七 Jupyter notebook
R markdown 示例:
图八 Jupyter notebook
以 R 语言为例,在一个 R 包开发过程中,常常集成 R markdown 文件来动态更新文档、教程和项目主页。
比如其中我开发的两个项目 configr、BioInstaller:
图九 configr 说明文档
图十 BioInstaller 项目主页
相关的 R 包:
我在这里设想了一个 R markdown 的应用场景:
软件和科学社区一直会有新的工具、思想、范式出现,生物信息学数据分析流程也不例外,我在这篇文章中所列的几种方式只能大致涵盖目前比较主流的几种方式。
还有一些”非主流“流程构建方式:
output (multi_omics) ~ input(2 * DNA@wgs@fq + 2 * DNA@wes@fq + DNA@chip_seq + 2 * RNA@rnaseq@fq + n * RNA@scrna + ...)
output (paper) ~ input(patient@survive + patient@habits_and_customs + mutation@proteinpaint +
—END—