如何用matlab编写分段函数_如何用 wdl 编写生信流程

生信分析离不开各种各样的流程。一个好的流程可以帮我们快速的完成各种各样的分析。然而，如何让我们的流程有更高的复用性，可移植性，稳定性，对于流程的使用而言，也非常的重要。 现状

对于如何进行生信分析，相信是八仙过海，各显神通

• 如果实现一个简单文件的文件，可能执行某条语句或运行某个软件即可。

• 如果满足一套复杂的要求，可能需要多条语句，用makefile，bash，或其他脚本语言串连成一个模块即可

• 如果减少一套分析方案的手动执行过程，则需要将多个模块串成pipeline或workflow…

workflow搭建就像将各个零件按设计图纸组装成机器一样，图纸不同，机器不同。但是如果各部分零件的规格和接口都不一样，图纸再简单，也是个复杂的工程，甚至只能一次次调整，耗时又耗力。目前生信分析的模块就是这样，没有统一的构建标准。

除了对流程进行定义缺乏标准外，如何调度和监控的任务，实现任务的并发、资源调度也是一件让人抓狂的事情。WDL是board推出的一套便于human readable and writable的流程编写规则，其主要的特点在于更适合于 生信人员(而非纯IT人员)来书写流程。

例子

我们先看一个hello_world的案例

task hello {
  input {
    String pattern
    File in
  }

  command {
    egrep '${pattern}' '${in}'
  }

  runtime {
    docker: "broadinstitute/my_image"
  }

  output {
    Array[String] matches = read_lines(stdout())
  }
}

workflow wf {
  call hello
}

具体介绍前，从上面的例子我们可以总结WDL 的一些特点，

(1)分层清晰，即使不了解WDL也能明白每部分的作用；

(2)明确定义输入和输出，方便模块间引用。

WDL组件

WDL核心是两大部分，task 和 workflow，二者组成结构相似，具体如下：

PART1

task

task是对模块的定义，把通常用的一些命令组合到一起，形成一个复用程度高的集合。以下部分wdl中task包括的元素：

input

为该task的输入模块。如果该模块有外部输入，则是必须提供该部分。

在input内部，主要是进行变量的声明。变量的声明格式为 Type(？+) name=value ,value可以不存在，由外来传入。

其中Type的类型有String、Int、Float、Boolean、File、Array、Pair、Map，以及自定义的Struct类型。

？表示变量可不用传入信息，+只能跟在Aarry后面，表示数组至少有一个元素。

当定义为File类型时，会自动检查文件是否存在，不存在会报错。有以下例子

Int? i                     # An integer value, i is optional
Float f = 27.3             # A floating point number
Boolean b = true           # A boolean true/false
String s = "hello, world"  # A string value
File f = "path/to/file"    # A file
Array[Int] a1

声明内部也支持各种表达式的，例如“+-*/”等，这一块跟所有的编程语言很类似，不再赘述。

input之外的声明

input的声明是为了参数外传，对于task里的其他中间变量可以写在input外面。

这一部分是可选，写到input里面和外面好像效果都是一样的。

command

具体执行的命令，最基本的是shell，也支持其他语言的嵌套，软件的运行命令等。shell的所有语法，基本在这个代码块中都适用。例如hello_world中的

egrep '${pattern}' '${in}'

其中{pattern} 是变量的引用，也可以用~｛pattern｝来引用。官方建议用~｛pattern｝来引用。

runtime

配置该task的运行的环境，一般有docker，memory，cpu等。这一块用来设置运行的环境，资源等，需要根据不同的运行平台进行设置。也可以是变量，由外部来导入。

output

输出模块，如果有输出，则一定要申明，其他模块在调用该模块的结果时可以直接引用，比如调用例子的结果方式是：hello.matches

声明的方式同input中一样，也支持表达式，例如“+-*/”等，如${name}+"txt"；也支持变量的内插，例如"${name}.txt"，当输出类型是File时，建议value值用双引号引起来，这样在下一个模块调用时才是绝对路径。"${name}.txt"，当输出类型是File时，建议value值用双引号引起来，这样在下一个模块调用时才是绝对路径。

output中有个glob，可以glob某个pattern的文件，返回一个array。在某些情况下会很有用。

在output中还支持很多函数，例如read_int()，详情见3中函数部分。

meta

可选，例如存放作者和版本信息等。

parameter_meta

可选，相当于-h 里面的信息。可以使用wdl-aid来生成wdl的文档，方便分享。

PART2

workflow

将定义的task有机的组合到一起，就形成了我们常说的流程。流程通常是DAG(有向无环图)的瀑布模式。

一个wdl里面只能有一个workflow定义，否则会报错。workflow的组成和task很相似，input，output，meta等定义方式一致，不再赘述。需要说明的是，workflow的output可以直接引用task的结果，比如：

output {
    File outfile1=task1.result
    File outfile2=task2.result
}

workflow与task最大的不同是没有command，而是有call，即调用task任务，同时为实现多样的DAG关系也存在特殊定义，比如scatter，if等，具体说明如下。

call

call是对task的直接调用，可以直接使用input的参数，也可以调用其他task的结果。正是因为call的顺序和输入输出的先后依赖关系，才构成了workflow的主要DAG关系。示例如下。

call my_task as my_task_alias2 {
    input: 
        threshold=threshold,
        file1=task1.output1,
        file2=task2.output1
  }

my_task是定义好的task名称，as后面是调用的时候的名称，如果同一个模块进行多次调用，需要给不同的别名，以便来区分每次调用后的输出。大括号里面的input，是将workflow的变量传参到task里, 在workflow里面变量直接调用即可，不需要$和双引号。

scatter

一个流程除了串联，还有并发，这个可以用scatter结合Array来实现。这个特性适用于我们有多个样品或多个文件，并发运行相同的分析，从而有效的节省时间，让DAG更简洁。示例如下：

scatter(i in integers) {
  call task1{input: num=i}
  call task2{input: num=task1.output}
}

scatter 和for循环不同，它是一种并发关系。在此例中，integers是一个array，调度器会对array里面的元素进行并发处理，而不是一个个循环处理。

if

对变量进行判断，选择性的运行某个task。注意这个语法没有else。

if(number!=1) { call task1 {input: num=i} }

其他

import

使用import 可以导入task.wdl 或者sub-workflow.wdl，被workflow直接调用。

import "~/wdl/hello_word.wdl" as hello

sub-workflow

workflow本身可以被另一个workflow import进去，以sub-workflow的形式被调用，可以很方便的作为一个整体去复用。

函数

wdl提供了许多标准的函数来帮助实现简单的读写功能。常用有：

• stdout()：标准输出流

• stderr()：标准错误流

• read_lines()：按行读取，返回array

• 其他的见说明

定义struct

定义一个数据结构，有点像面向对象的类,一般放在workflow或者task之外。例如：

struct Person {
    String name
    Int age
    File input1
}

声明：Person a

调用：a.name

使用

WDL可以用cromwell(java) 或者widdler(python)运行，在下载cromwell(https://github.com/broadinstitute/cromwell/releases)时，可以同时下载womtool。womtool可以根据WDL内需要的的input变量生成一个json版的输入文件模板。

输入文件

使用womtools 生成

java -jar womtool-48.jar inputs test.wdl 

即可生成json文件，根据提示的File，Int，Array输入不同的类型的数据：

{
  "task1.inputfile": "new.txt",
  "task1.age": "18",
  "task1.word": "Say hello"
  "task1.samples": ["jia","yi"],
  "task1.info": [{"name":"jia","math":"89"},{"name":"yi","math":"93"}]
}

cromwell简单运行

java -Dconfig.file=test.config -jar cromwell-48.jar run hello_world.wdl -i input.json -o workflow.json

各种配置具体可以查看cromwell官网(https://cromwell.readthedocs.io/en/develop/)

-Dconfig.file 流程配置文件(可选)

可以设置一系列配置信息，比如打开断点续跑功能Call Cache。

call-caching {
  enabled = true
  invalidate-bad-cache-results = true
}

run *.wdl 一种运行模式

run模式用于本地测试一个独立的WDL, 不用配置任何信息，简单粗暴，直接运行WDL。另一种模式是server，服务器引擎模式，依赖于环境配置。

-i 输入文件

上面提到的WDL需要的各种参数形成的json文件。

-o workflow 配置文件

可以配置workflow的一下信息，比如最后的目录整理。

{
  "final_workflow_outputs_dir": ""~/Result",
  "use_relative_output_paths": true,
  "final_workflow_log_dir": "~/Log",
  "write_to_cache": true,
  "read_from_cache": true
}

输出结果

输出分为层级，里面包含输入输出，运行脚本，日志，等各种信息记录。

ref

• https://github.com/openwdl/wdl/blob/master/versions/1.0/SPEC.md#arrayint-rangeint

作者：童蒙审稿：同尘

往期回顾

SpatialDB引发的10x技术原理的探究 图形掰弯利器--circlize 单细胞RNA系列专题之一：单细胞RNA测序中质控之重要细节(下篇) 单细胞RNA系列专题之一：单细胞RNA测序中质控之重要细节(上篇) 大火的单细胞分析软件Seurat，不能正常安装，怎么办？ 最新版SRA数据上传操作说明 ATAC-seq分析干货-2 ATAC-seq分析干货-1