Snakemake Workflows
Basic Rule Structure
# Snakefile
rule all:
input:
expand("results/{sample}_counts.txt", sample=SAMPLES)
rule align:
input:
r1 = "data/{sample}_R1.fq.gz",
r2 = "data/{sample}_R2.fq.gz",
index = "ref/genome.fa"
output:
bam = "aligned/{sample}.bam"
threads: 8
shell:
"bwa mem -t {threads} {input.index} {input.r1} {input.r2} | "
"samtools sort -@ {threads} -o {output.bam}"
Config File
# config.yaml
samples:
- sample1
- sample2
- sample3
reference: "ref/genome.fa"
threads: 8
# Snakefile
configfile: "config.yaml"
SAMPLES = config["samples"]
REF = config["reference"]
rule all:
input:
expand("results/{sample}.bam", sample=SAMPLES)
Wildcards and Expand
# Define samples
SAMPLES = ["A", "B", "C"]
CHROMOSOMES = [str(i) for i in range(1, 23)] + ["X", "Y"]
# Expand for all combinations
rule all:
input:
expand("results/{sample}_{chrom}.vcf", sample=SAMPLES, chrom=CHROMOSOMES)
# Access wildcard in rule
rule process:
input:
"data/{sample}.bam"
output:
"results/{sample}.vcf"
params:
sample = lambda wildcards: wildcards.sample
shell:
"bcftools call -s {params.sample} {input} > {output}"
Python Integration
rule analyze:
input:
counts = "data/{sample}_counts.txt"
output:
results = "results/{sample}_de.csv"
run:
import pandas as pd
counts = pd.read_csv(input.counts, sep='\t')
# Process data
results = counts.groupby('gene').sum()
results.to_csv(output.results)
Conda Environments
rule fastqc:
input:
"data/{sample}.fq.gz"
output:
"qc/{sample}_fastqc.html"
conda:
"envs/qc.yaml"
shell:
"fastqc {input} -o qc/"
# envs/qc.yaml
channels:
- bioconda
- conda-forge
dependencies:
- fastqc=0.12.1
- multiqc=1.14
Container Support
rule align:
input:
r1 = "data/{sample}_R1.fq.gz"
output:
bam = "aligned/{sample}.bam"
container:
"docker://biocontainers/bwa:v0.7.17"
shell:
"bwa mem {input} | samtools sort -o {output}"
# Run with Singularity
snakemake --use-singularity --singularity-args "-B /data"
Resource Management
rule memory_intensive:
input:
"data/{sample}.bam"
output:
"results/{sample}.vcf"
threads: 4
resources:
mem_mb = 16000,
time = "4:00:00",
gpu = 1
shell:
"variant_caller --threads {threads} {input} > {output}"
Cluster Execution
# cluster.yaml
__default__:
partition: normal
time: "1:00:00"
mem: 4G
align:
partition: normal
time: "8:00:00"
mem: 32G
cpus-per-task: 8
# SLURM execution
snakemake --cluster "sbatch --partition={cluster.partition} \
--time={cluster.time} --mem={cluster.mem}" \
--cluster-config cluster.yaml --jobs 100
# Or use profile
snakemake --profile slurm
Checkpoints
# For rules that determine downstream files dynamically
checkpoint split_samples:
input:
"data/all_samples.txt"
output:
directory("split/")
shell:
"split_script.py {input} {output}"
def get_split_files(wildcards):
checkpoint_output = checkpoints.split_samples.get(**wildcards).output[0]
return expand("split/{sample}.txt", sample=glob_wildcards("split/{sample}.txt").sample)
rule aggregate:
input:
get_split_files
output:
"results/aggregated.txt"
shell:
"cat {input} > {output}"
Modular Workflows
# Snakefile
include: "rules/qc.smk"
include: "rules/align.smk"
include: "rules/call.smk"
rule all:
input:
rules.qc_all.input,
rules.call_all.input
# rules/qc.smk
rule fastqc:
input: "data/{sample}.fq.gz"
output: "qc/{sample}_fastqc.html"
shell: "fastqc {input} -o qc/"
rule qc_all:
input: expand("qc/{sample}_fastqc.html", sample=SAMPLES)
Temporary and Protected Files
rule align:
input:
"data/{sample}.fq.gz"
output:
bam = temp("aligned/{sample}.unsorted.bam"), # Auto-deleted
sorted = protected("aligned/{sample}.bam") # Write-protected
shell:
"bwa mem {input} > {output.bam} && "
"samtools sort {output.bam} -o {output.sorted}"
Benchmarking
rule heavy_computation:
input:
"data/{sample}.bam"
output:
"results/{sample}.txt"
benchmark:
"benchmarks/{sample}.tsv"
shell:
"heavy_tool {input} > {output}"
Logging
rule process:
input:
"data/{sample}.bam"
output:
"results/{sample}.txt"
log:
"logs/{sample}.log"
shell:
"(command {input} > {output}) 2> {log}"
Run Commands
# Dry run (show what would be executed)
snakemake -n
# Run with 8 cores
snakemake --cores 8
# Force rerun of specific rule
snakemake --forcerun align
# Run until specific target
snakemake results/sample1.vcf
# Generate DAG visualization
snakemake --dag | dot -Tpdf > dag.pdf
# Generate report
snakemake --report report.html
Complete RNA-seq Workflow
configfile: "config.yaml"
SAMPLES = config["samples"]
rule all:
input:
"results/multiqc_report.html",
"results/deseq2_results.csv"
rule fastp:
input:
r1 = "data/{sample}_R1.fq.gz",
r2 = "data/{sample}_R2.fq.gz"
output:
r1 = "trimmed/{sample}_R1.fq.gz",
r2 = "trimmed/{sample}_R2.fq.gz",
json = "qc/{sample}_fastp.json"
threads: 4
shell:
"fastp -i {input.r1} -I {input.r2} -o {output.r1} -O {output.r2} "
"--json {output.json} --thread {threads}"
rule salmon_quant:
input:
r1 = "trimmed/{sample}_R1.fq.gz",
r2 = "trimmed/{sample}_R2.fq.gz",
index = config["salmon_index"]
output:
directory("salmon/{sample}")
threads: 8
shell:
"salmon quant -i {input.index} -l A -1 {input.r1} -2 {input.r2} "
"-o {output} --threads {threads}"
rule multiqc:
input:
expand("qc/{sample}_fastp.json", sample=SAMPLES),
expand("salmon/{sample}", sample=SAMPLES)
output:
"results/multiqc_report.html"
shell:
"multiqc qc/ salmon/ -o results/"
rule deseq2:
input:
quants = expand("salmon/{sample}", sample=SAMPLES),
metadata = config["metadata"]
output:
"results/deseq2_results.csv"
script:
"scripts/deseq2.R"
Related Skills
- workflow-management/nextflow-pipelines - Nextflow alternative
- workflows/rnaseq-to-de - RNA-seq workflow
- read-qc/fastp-workflow - QC in pipelines
majiayu000