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

接上两篇内容，本文主要讲述工作中NGS从科研进入医学临床领域，工作中接触到生信流程，以及最终在实现的过程。

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

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

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 等等，这样后续的步骤可以作为一个列表来表示：

  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 脚本就替换为：

  SN=B1701
  data=/opt/data
  result=/opt/result
  bwa=/opt/tools/bwa
  samtools=/opt/tools/samtools
  bwa=/opt/tools/gatk
  hg19=/opt/ref/hg19.fa
  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. 然而到这里就足够了么?

医学&临床应用是严肃的事情

**这里讲的生信的应用领域是医学临床领域，然而上述水平到这里最多也就是“工具”、“脚本”的水平，真要应用于临床，作为一个 “医疗产品”来要求，还有相当远的的距离。