用Snakemake跑pipeline到底有多么优雅 Day 21 & 22 # 100天生信/数据科学自我挑战#

如何优雅的跑pipeline回顾

Day 17 - 19 我们提到写（跑）Pipline是一件可以很简单（bash/Python/Perl乃至是R里面一个for loop），也可以是一个很复杂的事儿。

举例来说假如我们开始分析一批包含5个个体的RNA测序数据

步骤包含，

1.  从NCBI SRA数据库下载5个RNA seq 源数据（fastq格式）

2.  QC这两个fastq格式的文件

3.  将fastq比对到参考基因组上生成BAM文件

Bash的代码如下

indexFile="INDEX_ref" # path to indexed refrence file


rsids=( SRR8992621 ERR4393948 )
for rsid in "${rsids[@]}"
do
    # 第一步根据rsid下载fastq文件
    fastq-dump --split-files $rsid -0 ./
    #第二步 QC fastq
    fastp --in1 $rsid"_1.fastq" --in2 $rsid"_2.fastq" --out1 $rsid"_trimed_1.fastq" --out2 $rsid"_trimed_2.fastq"
    # 第三步用hisat比对到参考基因组，生成sam，并同时用samtools将 sam转成bam
    hisat2 -x $indexFile --met-file $rsid".met" -1 $rsid"_trimed_1.fastq" -2 $rsid"_trimed_2.fastq" | samtools view -Sb -F 4 -o $rsid".bam"
done

前面我们说到，Bash代码存在很多问题

比如所有文件都放在在一个文件夹，

比如如果某种因为原因pipeline跑的时候断了，不能从断点处续跑，

无法在超算上并行加速，

无法删掉中间不需要的文件夹，

以来软件安装可能会让可移植性变差等等（详见下文分解）

Snakemake跑同样的流程有多么优雅

Snakemake是python写的一个pipeline manager，与Python 和conda（不熟悉Conda的小伙伴请从文末找conda详解链接）都无缝衔接，所以功能异常强大，让写和跑pipeline都变得异常简单。

如果用Snakemake跑上面的这个流程，我们只需要四个文件

第一个文件是配置文件1 config.yaml，config 主要是告诉Snakemake我们的一些参数和文件地址

sra_list: #需要下载的SRA文件列表
  /Users/yanjunzan/Documents/barley_en/data/test_RNA.txt
working_dir: # 数据路径
  /Users/yanjunzan/Documents/barley_en/data/
result_dir:# 结果路径
  /Users/yanjunzan/Documents/barley_en/results/
ref:# 参考基因组
  Barley_MorexV3_pseudomolecules.fasta
gff3:# 参考基因组注释文件
  Hv_Morex.pgsb.Jul2020.gff3
fastp:# 软件运行参数
  qualified_quality_phred: "15"

第二个文配置文件2 mapping.yaml，软件的下载地址和版本信息，这里我们用的是conda

channels: # 所需要软件的channel位置
  - conda-forge
  - bioconda
  - defaults
dependencies:# 所需要软件的名称和版本号
  - hisat2=2.1.0
  - fastp=0.19.5
  - sra-tools=2.11.0
  - samtools=1.13

第三个文件test_RNA.txt 包含我们的SRA ID的输入文件，每行一个个体，第一列是SRA ID

[12:02:57][yanjun@176-30-ultuna: /Volumes/4TZS/ytb/genotyping]$ cat /Users/yanjunzan/Documents/barley_en/data/test_RNA.txt
Run Bytes LibraryLayout
SRR8992621 500061793 PAIRED
ERR4393948 500153103 PAIRED

第四个文件 是我们的Snakemake代码文件，snakefile_SRR.sk

import os
import pandas as pd
configfile: "./config.yaml"

###########################
# define sample
###########################

df = pd.read_csv(config["sra_list"],sep="\t")
SAMPLES = df["Run"].tolist()


WORKING_DIR = config["result_dir"]
RESULT_DIR = config["result_dir"]

###########################
# Input functions for rules
###########################

def sample_is_single_end(sample):
    """This function checks if raeds are single or paired end"""
    if os.stat(sample).st_size < 100:
        return True
    else:
        return False



rule all:
    input:
        qcfile = expand(RESULT_DIR + "fastp/{sample}.html",sample=SAMPLES),
        fq1 = expand(WORKING_DIR + "trimmed/{sample}_R1_trimmed.fq.gz",sample = SAMPLES),
        fq2 = expand(WORKING_DIR + "trimmed/{sample}_R2_trimmed.fq.gz",sample = SAMPLES),
        bam = expand(WORKING_DIR + "mapped/{sample}.bam", sample = SAMPLES)
    message:
        "Job done! Removing temporary directory"

rule get_SRR_files:
    output:
        fw = temp(WORKING_DIR + "fastq/{sample}_1.fastq"),
        rev= temp(WORKING_DIR + "fastq/{sample}_2.fastq")
    conda:
        "./mapping.yaml"
    params:
       SRA = "{sample}",
       DIR = config["result_dir"]+"fastq/"
    message:
        "using fastq-dump to download SRA data files to {output.fw}."
    conda:
        "./mapping.yaml"
    shell:
        "touch {output.rev}; fastq-dump --split-files {params.SRA} -O {params.DIR}"


rule fastp:
    input:
        fw = WORKING_DIR + "fastq/{sample}_1.fastq",
        rev= WORKING_DIR + "fastq/{sample}_2.fastq"
    output:
        fq1  = WORKING_DIR + "trimmed/{sample}_R1_trimmed.fq.gz",
        fq2  = WORKING_DIR + "trimmed/{sample}_R2_trimmed.fq.gz",
        html = RESULT_DIR + "fastp/{sample}.html"
    conda:
        "./mapping.yaml"
    message:"trimming {wildcards.sample} reads to {output.fq1}"
    threads: 1
    log:
        RESULT_DIR + "fastp/{sample}.log.txt"
    params:
        qualified_quality_phred = config["fastp"]["qualified_quality_phred"]
    shell:
        "fastp --thread {threads}  --html {output.html} \
        --qualified_quality_phred {params.qualified_quality_phred} \
        --in1 {input.fw} --in2 {input.rev} --out1 {output.fq1} --out2 {output.fq2}; \
        2> {log}"

rule hisat_mapping:
    conda:
        "mapping.yaml"
    input:
        fq1  = WORKING_DIR + "trimmed/{sample}_R1_trimmed.fq.gz",
        fq2  = WORKING_DIR + "trimmed/{sample}_R2_trimmed.fq.gz",
        indexFiles = [config["working_dir"]+ "MorexV3/" + config["ref"].replace("fasta","") + str(i) + ".ht2" for i in range(1,9)]
    output:
        bams  = protected(WORKING_DIR + "mapped/{sample}.bam"),
        met   = RESULT_DIR + "logs/{sample}_met.txt"
    params:
        indexName = config["working_dir"]+ "MorexV3/" + config["ref"].replace(".fasta",""),
        sampleName = WORKING_DIR + "trimmed/{sample}_R2_trimmed.fq.gz"
    message:
        "mapping reads to genome to bam files {params.sampleName}."
    threads: 2
    shell:
        "hisat2 -p {threads} --met-file {output.met} -x {params.indexName} \
        -1 {input.fq1} -2 {input.fq2} | samtools view -Sb -F 4 -o {output.bams}"

我们只需要直接运行

snakemake -s snakefile_SRR.sk

就会形成如下的文件结构，和最终的BAM。

今天我们先不介绍Snakemake具体怎么写，先介绍下它的优点，明天我们开始介绍snakemake的语法规则以及如何从头写一个pipeline

Snakemake如何应对 For loop中的问题

Day 17-19我们提到用for loop跑pipeline有六七个难题。我们前面也说了在for loop里怎么解决，现在我们看看以上的Snakemake代码是如何解决这些问题的。

第一，我们所有文件都在一个文件夹，不合理。

在For loop里

我们应该把源fastq 文件、QC后的fastq文件，BAM等按照文件类型分门别类，这样容易管理。

解决方案也很简单，先建立几个文件夹，然后在每一步输出的时候把响应的文件写入即可。

在Snakemake里

我们不需要创建文件夹，只需要写出路径而已，如果这个路径已经存在，那么文件直接写入这个路径，如果这个路径不存在，Snakemake会直接创建这个路径对应的文件夹，然后把文件直接写入这个路径。

第二， 如果出于某种原因,跑了两个个体后，代码在第三个个体的某一步中断了，如何在中断处继续跑？

在for loop 

这个就不那么容易了，需要几行代码，但是也没那么难。

我们要在跑pipeline之前，先查查输出到了哪一个个体。

然后在和输入的SRA  id比对，只把没跑的SRA id作为input输入到for loop里。

在Snakemake里

Snakemake会自动check 每个个体跑到了那一步，直接从断掉的那一步开始跑。

第三，如果我们不是有2个个体，而是500个,如何记录pipline现在跑到第几个个体的第几步了？

在for loop里

我们需要在每一步都输出一个log，也就是每个个体跑loop中的每一步的时候标准输出，重定向到一个文件中。

将来如果某一个个体没有成功，直接去log里可以看到到底是什么出了问题。容易debug

在Snakemake里，

rule fastp我们加了输出重定向把结果放到log里（见代码 rule fastp），同样的方法可以把其余rule的输出放到log里

第四，如果我们不是有2个个体，而是500个，如何在超算上并行？

For loop是不行了。

我们需要用parallel或者其余编程语言里的并行方式把每一个SRA ID 发送到一个核心中。

Snakemake与超算兼容，

在跑的时候只需要加上snakemake -s Snakefile --cluster 'sbatch -t 60 --mem=2g -c 1 -j 10'

第五，如果我们有500个个体，所用的硬盘又有限，怎么办？

For loop里

上回我们说过，我们需要在第三步跑成功的时候，将1，2，的输出删掉，但是当第三部不成功的时候保留1，2的输出。

所以需要我们去查，如果bam成功输出，则删掉1，2的fastq。

如果bam不成功，则查QC是不是成功，如果成功删掉原始fastq，如果QC不成功则保留原始fastq

Snakemake里

细心的小伙伴可能发现了，如果我们需要删除某一步骤的output只需要加上temp，在下一步成功运行后自动会删掉加上temp的文件。

如果我们要保证某一步输出永远不被删掉可以加上protected（会把文件权限更改成只读）

第六， 其中涉及到的软件有三个 fastq-dump，fastp，hisat，试想如果我们写了一个pipeline 打算发表，那么别人每次用都得装所有依赖的软件，涉及到不同平台（win10，Mac，Linux），不同的软件版本，有时候可能移植性比较差，别人会很难重复你的过程。

 For loop里

我们在pipeline运行前，要写代码检查所有软件是不是已经安装。如果没有安装则写相应的代码去安装。

Snakemake里

我们配置文件2的目的就是告诉snakemake我们所有软件的版本号，以及下载地址。

每次运行这个sankemake conda都会到对应的channel去下载（不熟悉conda的请去我们的Day1-10，链接见文末），所以不论谁跑这个pipline只需有conda，有snakemake就够了。

现在大家应该知道Snakemake这类pipeline manager的好处了。

从明天开始我们介绍，到底怎么从头写一个snakemake pipeline，它运行的逻辑什么，除了以上提到的优点它还有哪些优点。

关于博主

大家好，我是山石，山西农业大学（本）、天津大学（硕）、瑞典乌普萨拉大学（博），瑞典农业科学大学（博士后）。

在生命科学领域学习工作了15年，目前为瑞典农业大学生物大数据分析Research Fellow。

从普通二本院校到世界百强大学，一路踉跄走到了求学路的尽头。虽未获得成功，但是心中的小火苗还没有熄灭。此刻怀着仅剩的执着，开启了职场终身学习的征程。

100天生信/数据科学自我挑战回顾

0. 成为更好的自己之 << 100天生信/数据科学自我挑战>>

1. Day 1-100天生信/数据科学自我挑战（如何制定一个强有力的入门进阶计划）

2. Day 2 # 100天生信/数据科学自我挑战# -- Conda 入门到精通之初识Conda

3.Day 3 Conda 分类、安装和配置# 100天生信/数据科学自我挑战#

4. Day 4 Conda 环境使用的一些坑 # 100天生信/数据科学自我挑战#

5. Day 5 (忙的一地鸡毛没有更新)

6. Day 6 Conda Channel配置 #100天生信/数据科学自我挑战#

7. Day 7 & 8 Conda进阶技能 # 100天生信/数据科学自我挑战#

8. Day 9 Conda 不同电脑同步应用环境

9. Day 16 Conda 无网络状态下在不同电脑同步应用环境

10.如何优雅的跑pipline（1）Day 17 & 18 # 100天生信/数据科学自我挑战#

11.用loop的方式跑pipeline 遇到问题和解决方案-Day 19 & 20 # 100天生信/数据科学自我挑战#

如果你还不熟悉 <<100天生信/数据科学自我挑战>>

不管你是想入门生信/数据科学的初学者，还是想进阶生信/数据科学技能的从业者，可能都有不知如何下手或者坚持了几天就因为各种原因放弃了的经历。

如果是这样，那么这个100天生信/数据科学自我挑战就是为你而设计的（超简单、纯免费，自虐升级，无广告和付费内容，详情见下文）。

大家好，我是山石，山西农业大学（本）、天津大学（硕）、瑞典乌普萨拉大学（博），瑞典农业科学大学（博士后）。

2020年年底以来我在自己的社交账号上开始分享从生命科学湿实验硕士，转型生物信息学的博士以后的学习科研和经历。

其实这个问题同样困惑了我很久。回顾学习的过程，就像上台阶一样，往往登上一个台阶之后需要很久才有可能进阶，也或许很久都止步不前（比如我现在）。

究其原因就是找不到合适的方法，有时候即便找到了合适的方法，很难在学习新东西和做Projects之间找到一个平衡，无法长久坚持下去，养成终身学习的习惯。

而这个行业的进步又很快，每周都有新的方法和软件出来，为了避免被后浪早早拍死，凉透了，我发起这个和大家一起进步的自我挑战。

 << 100天生信/数据科学自我挑战>>是什么

今天我发起一个100天生信/数据科学自我挑战。

这个挑战其实很简单，只包括两个内容。

第一，参与者每天至少花5分钟去学习生信或者数据科学的知识; 

第二把学习的过程分享到自己的社交媒体中，任何社交媒体都可以，比如微博，知乎，B站或者微信（记得加入话题# 100天生信/数据科学自我挑战# 哦）。

请大家转发这条动态，让更多的小伙伴加入进来共同进步。

 << 100天生信/数据科学自我挑战>>为什么

学习是一种习惯，终身学习的习惯一旦养成，将很难改变，日积月累简单的一个好习惯，会让大家终身受益。

英国的一份统计调查显示，养成一个习惯的平均时间是66天。我们每天学习至少5分钟，100天自我挑战。

看似不长，只要坚持下去我相信一半以上的人都能够养成终身学习的习惯。

在自媒体分享的原因有三个，

其一，给自己立一个flag，创造舆论压力来鞭策自己，

其二，把自己学过的内容归档，每周，每月，每个季度都有可以回头总结和升华，做到查漏补缺，慢慢形成知识体系。

其三，通过100天持续不断的分享记录自己的学习过程，能够遇到很多同伴，一起讨论相互帮助。

最后，这个学习过程和做Project的过程可能会让大家敲开心仪的公司的实习offer和心仪导师的升学offer。

我们抱怨行业内卷，想躺平，但是试想像一下如果一位申请者拿了过去365天每天不间断的学习笔记来求职，面试官会不会拒而不见。

一个人的处境只会因为这个人做了什么而改变，而不是这个人身上的标签。所以请大家转发这条动态，让更多的小伙伴加入进来共同进步。

<< 100天生信/数据科学自我挑战>>怎么做

学习的内容不需要很复杂。

如果你是一位初学者，可以是一些学科的基本概念和技术。

例如什么是二代测序，Python 编程基础，什么是回归分析，什么是聚类等等。我的自媒体有很多相关内容（生信基础、linux，一二三代测序、RNA seq、 Variant calling等等）。

如果你已经有一定的基础，期望进阶学习，可以用公共数据重复文章中的内容，做一个RNA seq，Single cell sequencing 的课题，也可以去kaggle，Github上去学习一个解决具体实践问题的小project。

挑战赛的参与方式非常简单，大家只需要根据自己的情况，列一个提纲，每天循序渐进的学习一点点，记录在自己的社交媒体上，如果遇到合适的伙伴，可以一起积极交流互动。

请大家转发这条动态，让更多的小伙伴加入进来共同进步。

对于我自己来说，过去的半年我一直在梳理自己学习过程的一些基本概念和基础知识。

我会继续完成生信，R，Python编程，统计，GWAS，多组学整合，基因组选育的内容，并在我的自媒体上更新。

之后我计划在学科前沿，Web Developing和 App Developing努努力。

如果能把自己科研中的一些实用结果做成App和网站，能让科研成果触及更多的同行。

Document Your Journey To A Better Version Of Yourself

我创建了一个微信群，不好意思在自己自媒体天天更新的小伙伴，欢迎入群哈（私信我15934076136，我拉大家入群）。

不愿意入群的小伙伴，我在知乎提了一个问题” 100天生信/数据科学自我挑战，记录成为更好的自己？”。

大家也可以匿名更新，希望我们有足够多的朋友在一起更新，相互见证彼此的成长。

最后再次请大家转发这条动态，让更多的小伙伴加入进来共同进步。

期待100天后大家都成为更好的自己, 有一群优秀的小伙伴在100天后等你。

数量遗传学、GWAS、GS多组学分析系列课程链接

      1. 用时七年搜集整理的数量遗传学全网最全课程都包含些什么内容？

      2. 数量遗传学百年峥嵘开启信息化、智慧化现代农业和精准医疗新篇章

      3. 生信、多组学分析的生物学、生物信息学、群体数量遗传学基础(1)

      4.遗传学、多组学分析课程之（四）--为什么说突变、重组、单倍型及其在群体中的扩散是分析中最最重要的概念

      5.一、二、三代测序与分子标记的检出(数量、多组学分析课程)

      6.一、二、三代测序区别与联系(数量、多组学分析课程)
      7.第一代测序（Sanger 测序）原理 (数量、多组学分析课程)

      8.二代测序（Illumina测序）原理 (数量、多组学分析课程)

      9.第三代测序（PacBio测序）原理 (数量、多组学分析课程)

10第三代测序（Nanopore测序）原理 (数量、多组学分析课程)
      11. 第11讲 - 10X测序原理 (数量、多组学分析课程)

      12.第12讲 - 基于PCR的分子标记检测原理 (数量、多组学分析课程)

      13.第13讲 - 基因芯片的分子标记检测原理 (数量、多组学分析课程)
      14.第14讲 涉及生物信息学涵盖知识点小节(数量、多组学分析课程)

      15 - 重测序、外显子测序、简化基因组测序(数量、多组学分析课程)

      16 - 细数表观组学(数量、多组学分析课程)
      17 - 极简甲基化测序 (数量、多组学分析课程)

      18-极简Chip seq (数量、多组学分析课程)

      19 - 极简ATAC seq (数量、多组学分析课程)

      20 - 生信必备技能之极简Linux (数量、多组学分析课程)

      21 - 生信必备技能之Linux文件操作和系统资源查看

     22 - 生信必备Linux技能之软件运行与查阅帮助(数量、多组学分析课程)

     23 - 生信必备Linux技能之 Vim文档代码编辑 (数量、多组学分析课程)

     24 Linux常用命令总结- rysnc、grep 和find等-生信必备Linux技能

     24-2   零基础如何快速入门生信分析并且告别调包侠成长攻略

26 生信必备理论技能之--认识测序文库和测序过程

     27- 生信实战技能之—Fastq与FastQC

     28- 生信实战技能之—序列比对与BAM SAM 文件

     29- 生信实战技能之— Variant Calling 常见试验设计和注意细节

     30- Variant Calling （2）BAM文件质控

     31- Variant Calling （3）变异检出和质控

     32- RNA测序（1）RNA-seq类别、方法与应用场景

    33- RNA-seq建库、测序下机数据分析思路方法#数量多组学分析系列课程

数量遗传学前言讲座

1. 60年人工选择与三代遗传学家的探索-弗吉尼亚鸡体重双向选择系

2. 四篇文章NC+ NG +PG+Genes阐释数量遗传经典问题

高校求职启发

      1.2018-2020国内找高校教职受到的几点启发

      2.高校招聘启事晦涩词背后对应的考核待遇与支持

      3-4.高高校招聘启发系列之三四-材料准备和简历投递之后

      5.高校招聘启发分享系列五-拿到offer之后

科研Tips

1. 1000块自制表型组半自动检测设备

2. 没经费也能干大项目-50块钱测序一个人的基因组