前幾天柒巫，幫師妹分析轉(zhuǎn)錄組數(shù)據(jù)刷袍，由于她只需要FPKM，并且我剛好在學(xué)snakemake妓布，就試著用 snakemake 寫了一個簡單的從 raw data 到 FPKM 的流程姻蚓。

snakemake簡介：

• snakemake是一款基于python3的軟件
• snakemake的主要功能是快速搭建分析流程
• 官方網(wǎng)址

snakemake優(yōu)勢：

• 支持并行運算
• 支持斷點運算
• 支持流程控制
• 支持內(nèi)存控制
• 支持CPU核心控制
• 支持運行時間控制
• 支持。匣沼。狰挡。。释涛。加叁。

anaconda 安裝 snakmake

• 查看環(huán)境
  conda info --envs
• 創(chuàng)建 python3 環(huán)境
  conda create -n py3 python=3
• 激活 py3 環(huán)境
  conda activate py3
• 安裝 snakemake
  conda install snakemake -y
• 退出環(huán)境
  conda deactivate

raw2fpkm pipeline (trim_galore -> hisat2 -> samtools -> stringtie)

• raw data filter by trim-galore
• clean data mapping with hisat2
• mapping result sort as BAM
• quantify with stringtie

raw2fpkm.py 腳本，如下：

id={"874","898","908","970","971","972","974","976","978"}
hisat2_index="~/reference/G.barbadense/zju_v1.1/index/hisat2/zju"
annotation_gff="~/reference/G.barbadense/zju_v1.1/Hai7124_V1.1.Gene.gff"

rule all:
    input:
        expand("03.align/SRR8089{id}_hisat_sorted.bam"),
        expand("04.quantify/SRR8089{id}.exp.gtf"),
        expand("04.quantify/SRR8089{id}_fpkm.txt")

rule trim_galore:
    input:
        "01.raw/SRR8089{id}_1.fastq.gz",
        "01.raw/SRR8089{id}_2.fastq.gz"
    output:
        "02.clean/SRR8089{id}_1_trim.fastq.gz",
        "02.clean/SRR8089{id}_2_trim.fastq.gz"
    log:
        "02.clean/SRR8089{id}_trim.log"
    shell:
        "trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 5 \
        -o {output[0]} {output[1]} {input[0]} {input[0]} > {log} 2>&1"

rule hisat2:
    input:
        "02.clean/SRR8089{id}_1_trim.fastq.gz",
        "02.clean/SRR8089{id}_2_trim.fastq.gz"
    output:
        "03.align/SRR8089{id}_hisat.sam"
    log:
        "03.align/SRR8089{id}_hisat.log"
    shell:
        "hisat2 -p 30 -x {hisat2_index} -1 {input[0]} -2 {input[1]} -S {output[0]} > {log} 2>&1"

rule bam2sam:
    input:
        "03.align/SRR8089{id}_hisat.sam"
    output:
        "03.align/SRR8089{id}_hisat.bam"
    shell:
        "samtools view -b -S -h -@ 20 {input[0]} > {output[0]}"

rule bam_sort:
    input:
        "03.align/SRR8089{id}_hisat.bam"
    output:
        "03.align/SRR8089{id}_hisat_sorted.bam"
    shell:
        "samtools sort -@ 20 -o {output[1]} {input[0]} "

rule bam_index:
    input:
        "03.align/SRR8089{id}_hisat_sorted.bam"
    output:
        "03.align/SRR8089{id}_hisat_sorted.bam.bai"
    shell:
        "samtools index {input[0]} > {output[1]}"

rule stringtie:
    input:
        "03.align/SRR8089{id}_hisat_sorted.bam"
    output:
        "04.quantify/SRR8089{id}.gtf",
        "04.quantify/SRR8089{id}.exp.txt"
    log:
        "04.quantify/SRR8089{id}_quantify.log"
    shell:
        "stringtie -e -B -p 27 -G {annotation_gff} -o {output[0]} -A {output[1]} {input[0]} > {log} 2>&1"

rule fpkm:
    input:
        "04.quantify/SRR8089{id}.exp.txt"
    output:
        "04.quantify/SRR8089{id}_fpkm.txt"
    shell:
        "sed 's//_/g' {input[0]} | awk '{print $1'\t'$9}' > {output[0]}"

snakemake 的運行：
$ snakemake -s raw2fpkm.py -n -p # 查看是否有邏輯錯誤唇撬，并打印每一步的命令行
$ snakemake -s raw2fpkm.py --dag | dot -Tpdf > raw2fpkm_dag.pdf # 生成流程的邏輯拓撲圖

1554707436(1).jpg

生信初學(xué)，拋磚引玉窖认，snakmake更具體的用法請查看官網(wǎng)豫柬。

參考來源：
https://www.bilibili.com/video/av45832590
https://www.bilibili.com/video/av15908415
https://www.bilibili.com/video/av47588837/