专栏首页科技记者​宏转录组学习笔记(三)--通过脚本和snakemake实现自动化

​宏转录组学习笔记(三)--通过脚本和snakemake实现自动化

还是接上次的教程翻译,宏转录组学习笔记(二)宏转录组学习笔记(一)

通过脚本和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是帮助解决这些问题的几种工作流程系统之一。(您可以在此处阅读文档。)[1]让我们看一下!

首先,让我们激活我们的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的一些基础知识,但是,如果您需要其他教程,可以在这里[2]添加一个。

Reference

[1]https://snakemake.readthedocs.io/en/stable/

[2]https://ngs-docs.github.io/2018-cicese-metatranscriptomics/automation/

(教程结束!)

本文分享自微信公众号 - 科技记者(kejijizhe)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-03-29

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 宏转录组学习笔记(一)

    前面提到,已经有两家公司通过宏转录组(Metatranscriptomics)测序检测肠道微生物,面向消费者提供检测服务。对宏转录组充满了好奇,有这样的比方说,...

    用户1075469
  • 让windows 10 内置ubuntu(WSL)成为扩增子分析生产力

    先提示下,由于现在大部分电脑的win10版本是1903或者更低,wsl的性能相比2004版本的wsl2有一定差距。据说前者不是真正的linux内核,后者才是。又...

    用户1075469
  • 使用melonnpan通过扩增子或宏基因组测序数据有效预测微生物群落的代谢图谱

    热心肠研究院的这个介绍让我对这个软件产生了好奇,我决定学习一下这个软件的使用,看看它和picrust的区别在哪,picrust2刚刚发布,看看是棋逢对手还是略胜...

    用户1075469
  • 辞旧迎新

    1. 周五的组会,demo time, 想到过去一段时间加班加点的项目快要完工,心中是有点小开心的。想来想去,决定简单展示一下,顺便听听大家的意见。 于是...

    包子面试培训
  • 前端这条路怎么走,作为一名后端er,说说我的见解

    近期都游荡在各大群里看大家的讨论,经常看到关于程序员生涯的一些讨论,颇有感触,最近的国庆的确过得有些堕落,都没怎么更新,仔细相信还是应该分享点经验给大家的!想必...

    风间影月
  • 【JAVA进阶】之类型转换

    如果涉及到日期的转换,则需要在 BeanUtils之前使用注册方法转换一下日期,代码如下

    用户5640963
  • BAT加持的B端赋能时代,如何才能成就下一个马云?

    在乱世之中总是容易有人会崭露头角,移动互联网时代的乱战背景之下,同样会如此。无论在哪个领域里,我们总是能够看到在互联网巨头戒备森严的夹击之下突出重围的英雄们。滴...

    孟永辉
  • AI开发最大升级:Pandas与Scikit-Learn合并,新工作流程更简单强大!

    对于许多数据科学家来说,一个典型的工作流程是在Scikit-Learn进行机器学习之前,用Pandas进行探索性的数据分析。新版本的Scikit-Learn将会...

    新智元
  • 最新:iOS 13 适配

    DarkMode 主要从两个方面来适配,一是颜色,二是图片,适配的代码不是很多,接下来让我们一起来看看具体是怎么操作的吧。

    猿_人类
  • MySQL8.0基础教程 - 事务隔离级别解决之道

    隔离性是事务的基本特性之一,它可以防止数据库在并发处理时出现数据不一致的情况。最严格的情况下,我们可以采用串行化的方式来执行每一个事务,这就意味着事务之间是相互...

    JavaEdge

扫码关注云+社区

领取腾讯云代金券