生信分析云平台产品开发 - 3 生信分析pipeline的进化

接上两篇内容,本文主要讲述工作中NGS从科研进入医学临床领域,工作中接触到生信流程,以及最终在自动图形化开放式生信分析系统开发中生信workflow设计实现的过程。 1.手动命令行运行 2. 脚本连续运行 3.一个脚本 shell 文件运行整个分析流程 4. 自动扫描文件并运行脚本 5. 带报告的自动扫描并触发运行脚本


接上两篇内容,本文主要讲述工作中NGS从科研进入医学临床领域,工作中接触到生信流程,以及最终在自动图形化开放式生信分析系统开发中生信workflow设计实现的过程。

接触二代测序,生信分析,那真是打开了一个新世界的大门,各种名次术语满天飞,搞的头晕脑胀。什么“什么是高通量测序/NGS”、Sanger法测序(一代测序)、外显子测序(whole exon sequencing)、mRNA测序 (RNA-seq)、SNP/SNV(单核苷酸位点变异)、INDEL (基因组小片段插入)、copy number variation (CNV)基因组拷贝数变异、structure variation (SV)基因组结构变异等等。

百度了各种相关的分析软件和文件格式,什么fastq,fastq,bam,vcf等等。下面分阶段描述生信分析流程升级/进化的过程:

1.手动命令行运行

经过几个月接触,自学、爬坑,慢慢搞清楚了部分内容,在似懂非懂之间开始了生信流程分析,终于有一天明白过来,这所谓的pipeline其实就是基于文件的工作流啊。比如其中一个步骤:

attachments-2019-09-Uxmum0ND5d8314ec7daf6.png

QC 完成后,然后运行下一个步骤:

attachments-2019-09-lPkzzeCS5d8314fdcf0ea.png


运行模式,一个输入或者多个输入文件,通过软件分析/计算得到一个或者多个输出文件。然后输出文件部分或者全部作为下一个步骤的输入文件。这时候手动分析的话,只能手动的一个一个输入命令,完成每一个步骤,直到得到最后结果。


如下面代码:

bwa mem -t 8 -M -R \ 
"@RG\tID:0bdd6f55\tLB:5fba\tPL:Illumina\tPU:3102\tSM: B1701" \
B1701_R1.fq.gz B1701_R2.fq.gz | samtools view -bS - > B1701.bam

gatk ReorderSam \
-R /opt/ref/hg19.fa \
-I B1701.bam \
-O B1701_reordered.bam

gatk SortSam \
-I B1701_reordered.bam \
-O B1701_sorted.bam \
-SO coordinate

2. 脚本连续运行

随着熟练程度提高,生信分析上用到的软件/工具也熟悉起来了,但是问题也暴露出来了,简单的一套 GATK Best Practice 肿瘤突变分析流程,加上CNVSV 分析从 fastq 文件开始到最后得到过滤的 vcf 结果,一共有 30 多个步骤。自己一条一条输入次数多了就开始烦躁了。

这时候自然会考虑,如何减少手动输入,将这些脚本自动化。

脚本自动运行:当然这需要一点编程基础了。其实总的来看,每一个步骤的输入和输出可以根据最开始的输入文件来判断。例如 B1701_R1.fastq.gz,bwa map 之后得到B1701_R1.bam,所以只需要获得最初的文件前缀,作为 SampleNumber 字段,后续的中间输出,最终的输出文件都以这个 SampleNumber 为前缀,以扩展名作为区分。这时候脚本就可以连续运行了。以 shell 为例:总的脚本运行:

workrun.sh B1701_R1.fastq.gz B1701_R2.fastq.gz

脚本的第一步,就是获取输入文件:B1701_R1.fastq.gz B1701_R2.fastq.gz经过匹配计算,可以得到 B1701 作为 SampleNumber,并保存在变量$SN 中。后续的输出都以$SN.bam $SN_sortted.bam $SN_marked.bam 等等,这样后续的步骤可以作为一个列表来表示:

export SN=1701

bwa mem -t 8 -M -R \
"@RG\tID:0bdd6f55\tLB:5fba\tPL:Illumina\tPU:3102\tSM:$SN" \
$SN_R1.fq.gz $SN_R2.fq.gz | samtools view -bS - >$SN.bam

gatk ReorderSam \
-R /opt/ref/hg19.fa \
-I /opt/result/$SN.bam \
-O $SN_reordered.bam

gatk SortSam \
-I $SN_reordered.bam \
-O $SN_sorted.bam \
-SO coordinate

运行脚本之前使用 B1701 替换变量$SN 得到要运行的真实的 shell 命令

bwa mem -t 8 -M -R \ 
"@RG\tID:0bdd6f55\tLB:5fba\tPL:Illumina\tPU:3102\tSM: B1701" \
B1701_R1.fq.gz B1701_R2.fq.gz | samtools view -bS - > B1701.bam

gatk ReorderSam \
-R /opt/ref/hg19.fa \
-I B1701.bam \
-O B1701_reordered.bam

gatk SortSam \
-I B1701_reordered.bam \
-O B1701_sorted.bam \
-SO coordinate

继续完善:

  • 如何判断这一步是否真正完成了,运行过程有没有错误。如果有错误,停止后续步骤运行:这里首先想到的是,运行结束后,判断预期的输出文件是否存在,文件大小是否大于 0,有些软件即使运行错误也会创建一个大小为 0 的文件。

  • 比如计算这一步骤运行需要多少时间。在命令行 shell 前面加上 time

time gatk SortSam \ 
-I B1701_reordered.bam \
-O B1701_sorted.bam \
-SO coordinate

3.一个脚本 shell 文件运行整个分析流程


上面的内容解决了 shell 脚本连续运行的问题,但是还有一些遗留问题可以改进:


  • 输入文件如果指定一个目录是否更好一些? 如: $data

  • 输出文件如果指定一个目录是否更好一些? 如: $result

  • 运行的软件/工具/脚本路径使用变量替代,这样便于升级维护,升级时候只需要修改该变量的值就可以了。如:$bwa $samtools $gatk

  • 运行过程中引用的 reference 文件,数据库文件的路径也用变量替代,升级版本的时候只需要修改变量的路径就可以了,这样便于升级维护 如 $hg19 (hg19.fa)

  • 运行中的重要参数,一些 cutoff 值,配置的运行资源 如: $threads


这样经过以上替换,前面的 shell 脚本就替换为:

export SN=B1701
export data=/opt/data
export result=/opt/result
export bwa=/opt/tools/bwa
export samtools=/opt/tools/samtools
export bwa=/opt/tools/gatk
export hg19=/opt/ref/hg19.fa
export threads=8

time $bwa mem -t $threads -M -R \
"@RG\tID:0bdd6f55\tLB:5fba\tPL:Illumina\tPU:3102\tSM:$SN" \
$data/$SN_R1.fq.gz $data/$SN_R2.fq.gz \
| $samtools view -bS - >$result/$SN.bam

time $gatk ReorderSam \
-R $hg19 \
-I $result/$SN.bam \
-O $result/$SN_reordered.bam

time $gatk SortSam \
-I $result/$SN_reordered.bam \
-O $result/$SN_sorted.bam \
-SO coordinate

这时候已经将整套流程简单精简为一个 shell 脚本,如命名为 workrun.sh,每次运行整套流程之前,将变量$SN 的值修改为需要的值就可以了。如果要升级软件、升级 reference 文件版本,修改 shell 脚本相应变量值即可。

到这里就结束了么?还能继续改进么?请继续往下看。


4. 自动扫描文件并运行脚本


前面我们通过变量定义两个目录$data$result 分别来表示,分析流程的输入文件目录$data 和分析输出文件目录$result,这时候如果我们写一个脚本,按照一定周期判断$data 目录下是否有符合要求的文件,如果有文件符合要求,就运行前面的 workrun.sh 启动分析流程。待整个分析流程结束后,将$SN对应的 SampleNumber 值写入一个文件,下次扫描判断文件对应的 SampleNumber 是否已经分析过。

5. 带报告的自动扫描并触发运行脚本

前面已经实现了自动扫描并分析文件,这时候我们需要将保存$SN 的文件完善一下,在分析之前录入样本信息,具体样本信息的记录和操作。

参照文章:自动化图形生物信息分析系统开发-2 样本信息处理

运行分析之前,用 SampleReport字段表示分析状态,扫描脚本根据 SampleReport字段是否为空判断,该样本编号 SampleNumber对应的文件是否已经分析过。分析开始后,更新SampleReport字段为当前日期,分析完成后,再更新为分析完成时的日期。

分析报告,首先我们准备一个分析报告模板,将需要填充的字段,用变量的形式表示,如${sn}${sampleReport}等等,包括

  1. 样本信息
  2. 患者信息
  3. 分析结果
  4. 用药信息
  5. 引用文章链接
  6. 审核签名
  7. 等等

等分析结束后,从样本保存文件,和分析流程最终输出文件中获取数据并填充,得到整个分析报告。像这些数据处理过程,使用 shell 就有些吃力了,我这里使用 python 改写了上面的脚本,并实现了对数据处理,报告填充功能。


到这里,基本上就达到绝大多数公司的生信自动化分析水平了


6. 然而到这里就足够了么?

这里讲的生信的应用领域是医学临床领域,然而上述水平到这里最多也就是“工具”、“脚本”的水平,真要应用于临床,作为一个 “医疗产品”来要求,还有相当远的的距离。毕竟医学是严肃的事情,直接影响到人的健康和生命,希望各位生信大佬理解。

从“软件工程”的角度,上述内容也远远达不到一个软件产品的标准:

  • 首先这些脚本都是生信开发人员编写的,绝大多数没有测试,从单元测试、集成测试、功能测试、压力测试、稳定性测试都没有,一旦,项目复杂度上升,这些脚本/工具的代码质量堪忧,很多公司都是一边运行一边调试。

  • 其次基于命令行脚本的运行环境,没有友好的交互界面,对于使用者要求过高,难以普及大范围推广。对于使用者的要求基本上就是一个生信开发人员的要求:熟悉Linux 操作系统,熟悉各种常用分析软件和工具,能够从脚本错误输出中判断出原因并调试解决。


后续内容生信流程的图形化实现请继续关注本系列文章。


产品完整PPT下载


欢迎回复讨论或者加入QQ群:853718264

  • 发表于 2019-09-09 18:19
  • 阅读 ( 113 )
  • 分类:软件工具

0 条评论

请先 登录 后评论
不写代码的码农
豆浆包子

6 篇文章

作家榜 »

  1. 合肥国肽生物 113 文章
  2. 祝让飞 104 文章
  3. 刘永鑫 64 文章
  4. 生信分析流 47 文章
  5. SXR 44 文章
  6. 调研图 38 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章