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

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

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

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

软件获取：到官网sliverworkspace.com免费下载个人版，最新版本 2.0.277363

1.手动命令行运行

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

QC 完成后，然后运行下一个步骤：

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

如下面代码:

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 肿瘤突变分析流程，加上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等等，这样后续的步骤可以作为一个列表来表示：

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}等等，包括

• 样本信息
• 患者信息
• 分析结果
• 用药信息
• 引用文章链接
• 审核签名 等等

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

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

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

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

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

• 首先这些脚本都是生信开发人员编写的，绝大多数没有测试，从单元测试、集成测
  试、功能测试、压力测试、稳定性测试都没有，一旦，项目复杂度上升，这些脚本/
  工具的代码质量堪忧，很多公司都是一边运行一边调试。
• 其次基于命令行脚本的运行环境，没有友好的交互界面，对于使用者要求过高，难以普及大范围推广。对于使用者的要求基本上就是一个生信开发人员的要求：熟悉Linux 操作系统，熟悉各种常用分析软件和工具，能够从脚本错误输出中判断出原因并调试解决。

以上缺陷也是笔者开发“自动图形化生物信息分析系统”的初衷，后续内容生信流程的图形化实现请继续关注本系列文章。

您可以下载PPT或加QQ群:853718264讨论