使用snakemake 分析流程初步学习（一）

Snakemake是一种高效的工作流管理工具，常用于生物数据分析，如 ChIP-seq 数据处理。它通过定义规则（rule）来自动化分析流程，支持任务的并行执行与依赖管理，确保输入输出的正确性。Snakemake能有效简化复杂的 ChIP-seq 数据处理步骤，如原始数据的质量控制、对齐、峰值调用等。通过自动化流程，减少人为错误并提高分析效率。

snakemake官方教程：

https://snakemake.readthedocs.io/en/stable/tutorial/tutorial.html

先照着官方教程无脑运行一遍：

一、准备工作

1.下载示例数据

mkdir snakemake-tutorialcd snakemake-tutorial

# 下载数据

$ curl -L https://api.github.com/repos/snakemake/snakemake-tutorial-data/tarball -o snakemake-tutorial-data.tar.gz

# 解压

tar --wildcards -xf snakemake-tutorial-data.tar.gz --strip 1 "*/data" "*/environment.yaml"

2. 搭建环境

使用mamba来创建一个名为snakemake-tutorial的新的 conda 环境，并且指定了一个文件environment.yaml作为环境创建的基础。这个 YAML 文件包含了环境所需的所有包和依赖项的列表。

# 安装mamba

conda install -n base -c conda-forge mamba

# 使用mamba安装各种需要的软件

mamba env create --name snakemake-tutorial --file environment.yaml

3.环境搭建完成

$ conda activate snakemake-tutorial

$ snakemake --help#usage: snakemake [-h] [--dry-run] [--profile PROFILE] [--workflow-profile WORKFLOW_PROFILE] [--cache [RULE ...]] [--snakefile FILE]              

#  至此,完成了环境的搭建二、snakemake demo

新建一个名为 Snakefile 的文件，把下面每部分的内容依次写入到这个文件中。

step1: mapping reads

# 定义一个名为 'bwa_map' 的规则rule bwa_map:   # 指定输入文件，这里是参考基因组序列文件和样本的fastq格式的reads文件   input:       "data/genome.fa",  # 参考基因组序列文件的路径       "data/samples/A.fastq"  # 样本A的reads文件的路径   # 指定输出文件，这里是映射后的BAM文件   output:       "mapped_reads/A.bam"  # 输出的BAM文件的路径   # 使用shell命令来执行映射过程   shell:       # 这是一个管道命令，将BWA和Samtools命令串联起来       "bwa mem {input} | samtools view -Sb - > {output}"  # 执行BWA mem命令并将输出通过管道传递给Samtools       # {input} 是一个变量，代表输入文件的路径       # -Sb 参数告诉Samtools以BAM格式输出序列和质量分数       # - 是标准输入的表示，意味着Samtools将从BWA的输出中读取数据       # > 是重定向操作符，将Samtools的输出写入到指定的BAM文件中       # {output} 是一个变量，代表输出文件的路径

# ____________________________________________________

# 执行前面定义的流程 snakemake --cores 5 mapped_reads/A.bam

snakemake -np  mapped_reads/A.bam# -n :dryrun,将模拟执行工作流，但不会真正运行任何规则的脚本或命令# -p

Step 2: 泛化 read mapping 的规则

rule bwa_map:   input:       "data/genome.fa",       "data/samples/{sample}.fastq"   output:       "mapped_reads/{sample}.bam"   shell:       "bwa mem {input} | samtools view -Sb - > {output}"

#  将{sample}作为通配符

snakemake --cores 5 mapped_reads/{A,B}.bam

Step 3: 对 read 排序

rule samtools_sort:   input:       "mapped_reads/{sample}.bam"   output:       "sorted_reads/{sample}.bam"   shell:       "samtools sort -T sorted_reads/{wildcards.sample} "       "-O bam {input} > {output}"

# 运行snakemake --cores 5 sorted_reads/{A,B}.bam

Step 4: 建立索引 & 可视化

rule samtools_index:   input:       "sorted_reads/{sample}.bam"   output:       "sorted_reads/{sample}.bam.bai"   shell:       "samtools index {input}"

# 运行snakemake --dag sorted_reads/{A,B}.bam.bai | dot -Tsvg > dag.svg

Step 5: Calling genomic variants

rule bcftools_call:   input:       fa="data/genome.fa",       bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),       bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)   output:       "calls/all.vcf"   shell:       "bcftools mpileup -f {input.fa} {input.bam} | "       "bcftools call -mv - > {output}"

# 运行snakemake --cores 5 calls/all.vcf

snakemake -np calls/all.vcf

Step 6:使用自定义的脚本

可以添加自定义的脚本

新建一个绘图的脚本scripts/plot-quals.py

import matplotlibmatplotlib.use("Agg")import matplotlib.pyplot as pltfrom pysam import VariantFile

quals = [record.qual for record in VariantFile(snakemake.input[0])]plt.hist(quals)plt.savefig(snakemake.output[0])

定义画图的rule

Step 7:使用 "all" 定义 target rule

在工作流的顶部设置一个规则，该规则all将所有通常所需的目标文件作为input

至此，整个流程就算完成了。总体看一下：

到Snakefile的目录下，执行流程即可:

snakemake -n  ##  dryrun 测试下先

snakemake --cores 5