通过脚本和snakemake实现自动化
到目前为止,我们已经完成了所有工作,并复制并粘贴了许多命令来完成所需的操作。这可行!但是也可能很耗时,并且更容易出错。接下来,我们将向你展示如何将所有这些命令放入Shell脚本中。
一个shell脚本是一个文本文件的完整的shell命令,运行时就如同你在命令行交互方式运行它们。
在这里,我们将创建一个从中获取并一次运行它们全部的命令。
编写shell脚本
让我们将质量控制过程中的所有命令放入一个脚本中。
我们称之为run_qc.sh
。该sh
在的结尾告诉你,这是一个bash脚本。
使用nano(编辑器)编辑文件run-qc.sh
,并将以下内容放在其中:
cd ${PROJECT}
mkdir -p quality
cd quality
ln -s ../data/*.fq.gz ./
printf "I see $(ls -1 *.fq.gz | wc -l) files here.\n"
fastqc *.fq.gz -t 4
multiqc .
现在,这是一个shell脚本,您可以使用一个命令来执行所有这些命令的去-试试跑下吧!:
cd ~/
bash run-qc.sh
重新运行shell脚本
假设您想重新运行脚本。你要怎么做?
好吧,请注意,quality
目录是在脚本开始创建的,所有内容都在该目录中执行。因此,要像这样删除quality
目录,
rm -fr quality
然后你可以运行
bash run-qc.sh
编写Shell脚本的一些技巧
1.使它可执行
您可以bash
通过一些magic来摆脱上面的命令部分:
把
#! /bin/bash
放在文件顶部,然后运行
chmod +x ~/run-qc.sh
您现在可以运行
./run-qc.sh
代替bash run-rnaseq.sh
。
您可能会想,好吧,为什么这很重要?好的,您可以对R脚本和Python脚本执行相同的操作(但是放在/usr/bin/env Rscript
或/usr/bin/env python
放在顶部,而不是/bin/bash
)。这基本上用脚本的编写语言来注释脚本,因此您不必自己了解或记住。
所以:这不是必须的,但这是一个很好的技巧。
您也可以始终通过指定或来强制脚本以特定语言运行。bash ``Rscript
2.另一个很好的补充:使它很好地报错
Shell脚本的一个怪异的方面是(默认情况下)即使有错误,它们也可以继续运行。这是不好的行为,我们应该将其关闭。
您可以通过重新运行上面的脚本而不删除目录来观察此行为rnaseq/
-该mkdir
命令将打印错误,因为目录仍然存在,但是每个shell脚本的一个很好的补充就是使它在第一个错误时失败。通过放
set -e
在顶部-告诉bash在第一个错误时退出,而不是勇敢地继续前进。
3.最后一个不错的补充:使shell脚本打印出它们正在运行的命令!
你可能会注意到,shell脚本为您提供了它的运行命令的输出,但不告诉你它的运行命令。
如果添加
set -x
在Shell脚本的顶部,然后重新运行它,
cd ~/
rm -fr quality
./run-qc.sh
然后您将看到正在运行的全部命令!
关于shell脚本的最后说明:
set -e
并且set -x
仅在shell脚本中起作用-它们是bash命令。您需要在Python和R中使用其他方法。
Snakemake自动化!
通过shell脚本实现自动化非常棒,但是这里存在一些问题。
首先,您必须每次都运行整个工作流程,并且每次都要重新计算所有内容。如果您运行的工作流需要4天,并且在最后更改了命令,则必须手动进入,然后运行依赖于已更改命令的内容。
其次,它是非常明确的,并且不是很通用。如果要在其他RNAseq数据集上运行,则必须更改许多命令。
snakemake是帮助解决这些问题的几种工作流程系统之一。(您可以在此处阅读文档。)让我们看一下!
首先,让我们激活我们的snakemake环境
source deactivate
source activate snake
我们将自动化相同的脚本进行修剪,但是使用snakemake。
rule all:
input:
"trim/TARA_135_SRF_5-20_rep1_1m_1.qc.fq.gz",
"trim/TARA_135_SRF_5-20_rep1_1m_2.qc.fq.gz"
rule trim_reads:
input:
r1="data/TARA_135_SRF_5-20_rep1_1m_1.fq.gz",
r2="data/TARA_135_SRF_5-20_rep1_1m_2.fq.gz",
adapters="trim/combined.fa"
ouput:
p1="trim/TARA_135_SRF_5-20_rep1_1m_1.qc.fq.gz",
p2="trim/TARA_135_SRF_5-20_rep1_1m_2.qc.fq.gz",
s1="trim/TARA_135_SRF_5-20_rep1_1m_1_s1.qc.fq.gz",
s2="trim/TARA_135_SRF_5-20_rep1_1m_2_s2.qc.fq.gz"
shell:'''
trimmomatic PE {input.r1} \
{input.r2} \
{output.p1} {output.s1} \
{output.p2} {output.s2} \
ILLUMINACLIP:{input.adapters}:2:40:15 \
LEADING:2 TRAILING:2 \
SLIDINGWINDOW:4:2 \
MINLEN:25
'''
现在我们可以运行
cd $PROJECT
snakemake
您应该看到“什么都没做”。
那是因为修剪的文件已经存在!
让我们修复一下:
rm trim/TARA_135_SRF_5-20_rep1*
现在,当您运行时snakemake
,您应该看到正在运行Trimmomatic。是的!然后,如果snakemake
再次运行,您将发现它不需要执行任何操作-所有文件都是“最新的”。
添加环境
在整个研讨会中,我们一直在使用conda环境。我们展示了您必须使用来在Bioconda课程中导出塔拉环境 conda env export -n tara -f $PROJECT/tara_conda_environment.yaml
。我们也可以在snakemake规则中使用此环境!
rule all:
input:
"trim/TARA_135_SRF_5-20_rep1_1m_1.qc.fq.gz",
"trim/TARA_135_SRF_5-20_rep1_1m_2.qc.fq.gz"
rule trim_reads:
input:
r1="data/TARA_135_SRF_5-20_rep1_1m_1.fq.gz",
r2="data/TARA_135_SRF_5-20_rep1_1m_2.fq.gz",
adapters="trim/combined.fa"
ouput:
p1="trim/TARA_135_SRF_5-20_rep1_1m_1.qc.fq.gz",
p2="trim/TARA_135_SRF_5-20_rep1_1m_2.qc.fq.gz",
s1="trim/TARA_135_SRF_5-20_rep1_1m_1_s1.qc.fq.gz",
s2="trim/TARA_135_SRF_5-20_rep1_1m_2_s2.qc.fq.gz"
conda: "tara_conda_environment.yaml"
shell:'''
trimmomatic PE {input.r1} \
{input.r2} \
{output.p1} {output.s1} \
{output.p2} {output.s2} \
ILLUMINACLIP:{input.adapters}:2:40:15 \
LEADING:2 TRAILING:2 \
SLIDINGWINDOW:4:2 \
MINLEN:25
'''
我们现在不打算在集群上运行它,因为它要求您能够下载内容,而我们不能执行此操作。但是,这是将来执行此操作的语法。
其他资源
今天,我们已经介绍了snakemake
的一些基础知识,但是,如果您需要其他教程,可以在这里添加一个。