专栏首页图形化开放式生信分析系统开发图形化开放式生信分析系统开发 - 3 生信分析流程的进化
原创

图形化开放式生信分析系统开发 - 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 肿瘤突变分析流程,加上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 等等,这样后续的步骤可以作为一个列表来表示:

  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. 然而到这里就足够了么?

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

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

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 满分室间质评之GATK Somatic SNV+Indel+CNV+SV(下)性能优化

    #此处是原先Manta分析SV的步骤一,生成runWorkflow.py,因为这一不步速度很快,所以串行执行 rm -f ${result}/${sn}/r...

    SliverWorkspace
  • 使用Gatk Germline spns-indels Pipeline分析遗传病(耳聋)

    这次没有拿到遗传病的室间质评的数据,直接从NCBI上找一些数据来分析。NCBI上搜索deaf,点击第一条搜索结果,最后几经跳转找到数据下载页面:https://...

    SliverWorkspace
  • 图形化开放式生信分析系统开发 - 1 需求分析及技术实现

    从三年前开始,工作的原因接触到了NGS(二代测序)技术和相关的生信分析,在公司技术到临床应用转化过程中遇到一系列问题,在问题中挣扎、解决问题的过程中逐渐有了开发...

    SliverWorkspace
  • Rabbitmq分布式事务

    今天小编带来的是分享课题是分布式事务。小编是在一家O2O公司做程序员,今天就以公司的一则订单业务来作为分享课题的场景。

    黑洞代码
  • MySQL常见的图形化工具

    MySQL作为一款非常流行的、开源的关系型数据库,应用非常广泛。因为MySQL开源的缘故,图形化管理维护工众多,除了系统自带的命令行管理工具之外,还有许多其他的...

    java乐园
  • [Linux]shell基础教程1-变量、字符串、数组、注释

    这是bash的一个特殊参数,但是也可以用在其他的shell中,比如sh、zsh、 tcsh 或者dash。使用echo命令可以查看正在使用的shell名称。

    祥知道
  • JS示例01-鼠标移入移出事件

    1、匿名函数 2、鼠标事件 3、document.getElementById() 4、window.onload 5、行间事件提取

    专注APP开发
  • python 魔法方法连载三 __setattr()__

    上一篇我们讲到了当实例对象获取一个不存在的属性时,只要你重载了__getattr()__,就能定制错误提示。有了__getattr__(), 那必须有 __se...

    周辰晨
  • Maven系列第8篇:大型Maven项目,快速按需任意构建,必备神技能!相知恨晚!

    整个maven系列的内容前后是有依赖的,如果之前没有接触过maven,建议从第一篇看起,本文尾部有maven完整系列的连接。

    路人甲Java
  • Xcode真机测试could not find developer disk image解决方法

            在使用Xcode进行真机调试的时候,有时根据真机的系统不同,会出现could not find developer disk image 错误,...

    珲少

扫码关注云+社区

领取腾讯云代金券