如何让Snakemake RNA-seq管道按样本类型自动切换子流程

Got it, let's break down how to solve this cleanly—no need to duplicate rules or run separate pipelines. The key is to leverage Snakemake's dynamic input/output functions and conditional logic to adapt to each sample's type, either via a config file or automatic file detection.

Option 1: Explicit Sample Type in Config File #

This is straightforward if you know your sample types upfront. Start by defining sample types in your config.yaml:

samples:
  sample_A: SE
  sample_B: PE
  sample_C: SE

Step 1: Write Helper Functions for Input/Output #

Create reusable functions to fetch the correct files for SE/PE samples. Add these to your Snakefile:

def get_raw_input(wildcards):
    sample_type = config["samples"][wildcards.sample]
    if sample_type == "SE":
        return f"data/raw/{wildcards.sample}.fastq.gz"
    else:
        return [
            f"data/raw/{wildcards.sample}_R1.fastq.gz",
            f"data/raw/{wildcards.sample}_R2.fastq.gz"
        ]

def get_qc_output(wildcards):
    sample_type = config["samples"][wildcards.sample]
    if sample_type == "SE":
        return f"data/qc/{wildcards.sample}_qc.fastq.gz"
    else:
        return [
            f"data/qc/{wildcards.sample}_R1_qc.fastq.gz",
            f"data/qc/{wildcards.sample}_R2_qc.fastq.gz"
        ]

Step 2: Build Universal Rules (No Duplication!) #

Instead of writing separate rules for SE/PE, use a single rule that adapts its command based on input count. For example, a quality control rule:

rule quality_control:
    input:
        get_raw_input
    output:
        get_qc_output
    params:
        cmd_flags = lambda wildcards, input, output: (
            f"-i {input} -o {output}" if len(input) == 1
            else f"-i1 {input[0]} -i2 {input[1]} -o1 {output[0]} -o2 {output[1]}"
        )
    shell:
        """
        somecommand {params.cmd_flags}
        """

Step 3: Merge to Single Output #

For the merging step, write another helper function to fetch QC outputs, then build a rule that handles both cases:

def get_merge_input(wildcards):
    sample_type = config["samples"][wildcards.sample]
    if sample_type == "SE":
        return f"data/qc/{wildcards.sample}_qc.fastq.gz"
    else:
        return [
            f"data/qc/{wildcards.sample}_R1_qc.fastq.gz",
            f"data/qc/{wildcards.sample}_R2_qc.fastq.gz"
        ]

rule merge_to_single_output:
    input:
        get_merge_input
    output:
        f"data/merged/{wildcards.sample}_final.bam"
    shell:
        """
        if [ {len(input)} -eq 1 ]; then
            # SE case: process single file to merged output
            star --readFilesIn {input} --outFileNamePrefix data/merged/{wildcards.sample}_
        else
            # PE case: process paired files to single merged BAM
            star --readFilesIn {input[0]} {input[1]} --outFileNamePrefix data/merged/{wildcards.sample}_
        fi
        mv data/merged/{wildcards.sample}_Aligned.sortedByCoord.out.bam {output}
        """

Step 4: Define Final Targets #

Use expand to generate targets for all samples:

rule all:
    input:
        expand("data/merged/{sample}_final.bam", sample=config["samples"].keys())

Option 2: Automatic Detection of Sample Type (No Config Needed) #

If you don't know sample types upfront, use a checkpoint to run the download first, then auto-detect file counts to determine SE/PE.

Step 1: Checkpoint for Data Download #

First, define a checkpoint to handle downloading (Snakemake will run this before resolving downstream rules):

checkpoint download_raw_data:
    output:
        directory("data/raw")
    shell:
        """
        # Your download command here (e.g., wget, sra-toolkit)
        # Ensure PE files follow a consistent naming pattern like {sample}_R1.fastq.gz
        """

Step 2: Auto-Detect Samples and Types #

Add functions to parse downloaded files and determine sample types:

import glob
import os

def get_all_samples():
    # Fetch output from the download checkpoint
    checkpoint_dir = checkpoints.download_raw_data.get_output()[0]
    # Extract unique sample names from file names
    samples = set()
    for file in glob.glob(os.path.join(checkpoint_dir, "*.fastq.gz")):
        # Adjust split logic to match your file naming convention
        sample_name = os.path.basename(file).split("_")[0]
        samples.add(sample_name)
    return sorted(samples)

def get_sample_files(wildcards):
    checkpoint_dir = checkpoints.download_raw_data.get_output()[0]
    return glob.glob(os.path.join(checkpoint_dir, f"{wildcards.sample}*.fastq.gz"))

def is_pe_sample(wildcards):
    return len(get_sample_files(wildcards)) == 2

Step 3: Adapt Rules to Auto-Detected Types #

Rewrite your rules to use these auto-detection functions. For example:

rule quality_control:
    input:
        get_sample_files
    output:
        lambda wildcards: (
            f"data/qc/{wildcards.sample}_qc.fastq.gz" if not is_pe_sample(wildcards)
            else [
                f"data/qc/{wildcards.sample}_R1_qc.fastq.gz",
                f"data/qc/{wildcards.sample}_R2_qc.fastq.gz"
            ]
        )
    params:
        cmd_flags = lambda wildcards, input, output: (
            f"-i {input} -o {output}" if len(input) == 1
            else f"-i1 {input[0]} -i2 {input[1]} -o1 {output[0]} -o2 {output[1]}"
        )
    shell:
        """
        somecommand {params.cmd_flags}
        """

Step 4: Final Targets #

Use the auto-detected sample list for your final rule:

rule all:
    input:
        expand("data/merged/{sample}_final.bam", sample=get_all_samples())

Key Tips to Avoid Headaches #

• Consistent Naming: Ensure PE files follow a predictable pattern (e.g., sample_R1.fastq.gz, sample_R2.fastq.gz) so detection functions work reliably.
• Reuse Logic: Keep conditional logic in helper functions/params instead of shell scripts where possible—this makes the pipeline easier to maintain.
• Test Incrementally: Test with one SE and one PE sample first to verify the logic works before scaling to full datasets.

内容的提问来源于stack exchange，提问作者athiebaut