一个真实的生物信息学流程示例
如何把shell脚本改造成bpipe流程
在本文中，我们将构建一个用于call突变的Bpipe流程。在这个流程中，我们会依次执行以下四个步骤：
- 用bwa进行比对
- 用samtools进行建索引和排序
- 用Picard移除PCR重复
- 用samtools找突变
要使流程运行成功，首先得安装上述软件，并加入环境变量方便调用。软件安装有点麻烦，也可以先阅读完整个章节再去实践操作。假设参考基因组名为reference.fa，并已经用bwa构建好索引。需要比对的reads位于本地目录中，并命名为s_1.txt。
shell脚本
为了展示如何将shell脚本转换为Bpipe脚本，我们先来写个简单的 shell 脚本吧：
#!/bin/bash
bwa aln -I -t 8 reference.fa s_1.fq > out.sai
bwa samse reference.fa out.sai s_1.fq > out.sam
samtools view -bSu out.sam | samtools sort - out.sorted
java -Xmx1g -jar /usr/local/picard-tools/MarkDuplicates.jar \
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000\
METRICS_FILE=out.metrics \
REMOVE_DUPLICATES=true \
ASSUME_SORTED=true \
VALIDATION_STRINGENCY=LENIENT \
INPUT=out.sorted.bam \
OUTPUT=out.dedupe.bam
samtools index out.dedupe.bam
samtools mpileup -uf reference.fa out.dedupe.bam | bcftools view -bvcg - > out.bcf
Step 1 - 将命令转换成Bpipe阶段
首先，我们需要做的就是确定哪些命令在逻辑上属于一起，并声明它们为 Bpipe 阶段。在每个阶段中，我们将原始的shell 命令转换为Bpipe exec语句：
align = {
exec "bwa aln -I -t 8 reference.fa s_1.fq > out.sai"
exec "bwa samse reference.fa out.sai s_1.fq > out.sam"
}
sort = {
exec "samtools view -bSu out.sam | samtools sort - out.sorted"
}
dedupe = {
exec """
java -Xmx1g -jar /usr/local/picard-tools/MarkDuplicates.jar
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
METRICS_FILE=out.metrics
REMOVE_DUPLICATES=true
ASSUME_SORTED=true
VALIDATION_STRINGENCY=LENIENT
INPUT=out.sorted.bam
OUTPUT=out.dedupe.bam
"""
}
index = {
exec "samtools index out.dedupe.bam"
}
call_variants = {
exec "samtools mpileup -uf reference.fa out.dedupe.bam | bcftools view -bvcg - > out.bcf"
}
Note：在原始的shell脚本中，运行MarkDuplicates的多行命令每行的末尾都有反斜杠\，但Bpipe本身就能理解多行命令，所以我们不需要在Bpipe的版本中加入反斜杠。
Step 2 - 连接不同阶段并运行
在第 1 步中，我们虽然定义了不同的阶段，但并没有定义如何将他们串联。所以我们需要在最后一行加入以下的命令将不同的阶段连接在一起：
Bpipe.run {
align + sort + dedupe + index + call_variants
}
到目前为止，虽然我们只增加或改了一些字符，但的确已有了一个全功能的 Bpipe 脚本，它已经比原来的 shell 脚本强大了不少。将上述文本保存为pipeline.txt，就可以调用 Bpipe 运行了：
bpipe run pipeline.txt
输出如下所示：
### ==============================================================================================
| Starting Pipeline at 2011-10-06 |
### ==============================================================================================
### ======================================== Stage align =========================================
[17bp reads: max_diff = 2
[bwa_aln](bwa_aln]) 38bp reads: max_diff = 3
...
虽然这个比对流程看起来是在前台运行，但实际上 Bpipe 已经在后台启动了它，你可以按 Ctrl + C，Bpipe 会提示运行选项：
Pipeline job running as process 33202. Terminate? (y/n):
如果你回答no，即使退出，工作也会继续运行。我们也可以使用bpipe jobs命令来查看正在运行的任务，或使用bpipe log命令继续查看输出的日志文件。
Step 3 - 定义输入输出变量
到目前为止，虽然我们的Bpipe流程已经能够正常运作了，但还有些简陋，比如，我们没有根据Bpipe的输入和输出变量来定义我们的输入和输出。为了添加输入和输出变量，我们需要用$input和$output来替换输入和输出。我们还可以借此机会设置一些全局变量比如参考基因组和Picard的路径，使这个脚本更加通用：
REFERENCE="reference.fa"
PICARD_HOME="/usr/local/picard-tools/"
align = {
exec "bwa aln -I -t 8 $REFERENCE $input > ${input}.sai"
exec "bwa samse $REFERENCE ${input}.sai $input > $output"
}
sort = {
exec "samtools view -bSu $input | samtools sort -o - - > $output"
}
dedupe = {
exec """
java -Xmx1g -jar $PICARD_HOME/MarkDuplicates.jar
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
METRICS_FILE=out.metrics
REMOVE_DUPLICATES=true
ASSUME_SORTED=true
VALIDATION_STRINGENCY=LENIENT
INPUT=$input
OUTPUT=$output
"""
}
index = {
exec "samtools index $input"
forward input
}
call_variants = {
exec "samtools mpileup -uf $REFERENCE $input | bcftools view -bvcg - > $output"
}
Bpipe.run {
align + sort + dedupe + index + call_variants
}
现在，这个Bpipe流程在工作方式上有了更大的灵活性！比如你要删除中间去重（dedupe）这个阶段，那就直接在最后面的Bpipe.run语句中删除dedupe即可。Bpipe会自动安排好输入文件名和输出文件名！
由于第一个阶段现在需要一个输入参数，因此在运行 bpipe 时必须提供一个输入，如下所示：
bpipe run pipeline.txt s_1.fq
Note：
- 1）第一个align的命令中
$input的变量是用${input}.sai的格式而不是$input.sai。这是因为在Groovy中".“是一个特殊操作符，用来表示变量的属性。而在bash中如果遇到这种情况，解决方法就是用大括号括起来，和周围的字符区别开来。 - 2）你可能已经注意到了在index阶段最后有一个forward input语句。这是因为index阶段的输出（BAM文件的索引）并不是后续阶段的正确输入（正确输入是BAM文件）。“forward”命令允许阶段建议哪些输出可能是下一个管道阶段要使用的最佳默认输入。在这种情况下，程序建议输入（BAM文件）应传递到下一个管道阶段作为默认输入。在下一阶段，确实想使用BAM文件索引，它可以使用from语句或使用Bpipe对输入文件扩展名的支持（例如，将其称为“$input.bai”）。
Step 4 - 重命名输出
我们的脚本还有一个缺陷：输出文件的后缀不是很令人满意。因为在 Bpipe 流程中，每一个阶段的后缀默认是该阶段名，比如默认情况下，第一个align阶段的输出文件名是s_1.fq.align。所以我们还需要在$input和$output变量中添加文件扩展名来自定义扩展名：
PICARD_HOME="/usr/local/picard-tools/"
REFERENCE="reference.fa"
align = {
exec """
bwa aln $REFERENCE $input > $output.sai;
bwa samse $REFERENCE $output.sai $input > $output.sam
"""
}
sort = {
exec "samtools view -bSu $input.sam | samtools sort -o - - > $output.bam"
}
dedupe = {
exec """
java -Xmx1g -jar $PICARD_HOME/MarkDuplicates.jar
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
METRICS_FILE=out.metrics
REMOVE_DUPLICATES=true
ASSUME_SORTED=true
VALIDATION_STRINGENCY=LENIENT
INPUT=$input.bam
OUTPUT=$output.bam
"""
}
index = {
exec "samtools index $input.bam"
}
call_variants = {
exec "samtools mpileup -uf $REFERENCE $input.bam | bcftools view -bvcg - > $output.vcf"
}
run {
align + sort + dedupe + index + call_variants
}
在这里，我们在输入和输出变量中添加了文件扩展名，这些扩展名告诉Bpipe每个阶段正在进行哪种操作。Bpipe将操作分为两类：
- 转换操作（Transform）是根据输入生成一种新型文件，即它们将文件扩展名修改为新文件。
- 过滤器操作（Filter）会在不更改文件类型的情况下修改文件，即它们使文件扩展名保持不变，但在名称正文中添加了一个组件。
Bpipe通过查看文件扩展名来推断正在发生的事情。Bpipe允许您根据需要对阶段操作的类型进行申请-请参见“转换”或“过滤器”注释，但是在大多数情况下，Bpipe会通过文件扩展名来推断正在发生的情况。真正重要的是我们没有规定了文件的全称，我们只是帮助它了解文件名是如何改变的。这意味着我们的流程是灵活的，不用使用硬编码的文件名，而使用可识别的明智名称。