搜索
查看: 130|回复: 0

图形化开放式生信分析系统开发 - 3 生信分析pipeline的进化

[复制链接]

4

主题

4

帖子

103

积分

注册会员

Rank: 2

积分
103
发表于 2019-9-9 12:35:55 | 显示全部楼层 |阅读模式
本帖最后由 豆浆包子 于 2019-9-19 12:59 编辑

自动图形化开放式生信分析系统开发
-3 生信分析pipeline的进化
接上两篇内容,本文主要讲述工作中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其实就是基于文件的工作流啊。
比如其中一个步骤:

输入fastq文件    ->       运行软件fastQC       ->       输出QC结果html文件

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

输入fastq文件    ->       运行软件bwa       ->       输出sam文件
file:////tmp/wps-sliver/ksohtml/wps5xkU0K.pngfile:////tmp/wps-sliver/ksohtml/wpsE4YZ0Y.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、 脚本连续运行

随着熟练程度提高,生信分析上用到的软件/工具也熟悉起来了,但是问题也暴露出来了,简单的一套GATKBestPractice肿瘤突变分析流程,加上CNV,SV分析从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等等,这样后续的步骤可以作为一个列表来表示:

bwa mem -t 8 -M -R  \
"@RG\tID:0bdd6f55\tLB:5fba\tPL:Illumina\tPU:3102\tSMSN"  \  
$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

继续完善
l 如何判断这一步是否真正完成了,运行过程有没有错误。如果有错误,停止后续步骤运行:这里首先想到的是,运行结束后,判断预期的输出文件是否存在,文件大小是否大于0,有些软件即使运行错误也会创建一个大小为0的文件。
l 比如计算这一步骤运行需要多少时间。在命令行shell 前面加上time

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

3、 一个脚本shell文件运行整个分析流程
上面的内容解决了shell脚本连续运行的问题,但是还有一些遗留问题可以改进:
l 输入文件如果指定一个目录是否更好一些? 如: $data
l 输出文件如果指定一个目录是否更好一些? 如: $result
l 运行的软件/工具/脚本路径使用变量替代,这样便于升级维护,升级时候只需要修改该变量的值就可以了。如:$bwa  $samtools  $gatk
l 运行过程中引用的reference文件,数据库文件的路径也用变量替代,升级版本的时候只需要修改变量的路径就可以了,这样便于升级维护  如 $hg19 (hg19.fa)
l 运行中的重要参数,一些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\tSMSN"  \  
$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的文件完善一下,在分析之前录入样本信息,具体样本信息的记录和操作。

运行分析之前,用SampleReport字段表示分析状态,扫描脚本根据SampleReport字段是否为空判断,该样本编号SampleNumber对应的文件是否已经分析过。分析开始后,更新SampleReport字段为当前日期,分析完成后,再更新为分析完成时的日期。
分析报告,首先我们准备一个分析报告模板,将需要填充的字段,用变量的形式表示,如${sn},${sampleReport}等等,包括

A样本信息
B患者信息
C分析结果
等等。

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

6、 然而这样就够了么?

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

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

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

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


产品完整PPT下载


欢迎加入QQ群讨论:853718264

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:基因组组装 求助python代码或解题思路
下一篇:生成genefuse 的fusion.csv
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 09:35 , Processed in 0.033993 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.