如何在Snakemake中自动识别并合并技术/生物学重复样本

问题描述 #

我正在编写Snakemake工作流，想要通过sample.tsv文件自动判断样本是技术重复还是生物学重复，并在对应阶段执行重复样本合并规则。

样本表格式 #

sample.tsv的结构如下：

|sample    | unit_bio | unit_tech | fq1 | fq2 |                                                                                                                      
|----------|----------|-----------|-----|-----|
| bCalAnn1 |        1 |         1 | /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_1_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_1_R2.fastq.gz |
| bCalAnn1 |        1 |         2 | /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_2_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn1_1_2_R2.fastq.gz |
| bCalAnn2 |        1 |         1 | /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_1_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_1_R2.fastq.gz |
| bCalAnn2 |        1 |         2 | /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_2_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn2_1_2_R2.fastq.gz |
| bCalAnn2 |        2 |         1 | /home/assembly_downstream/data/arima_HiC/bCalAnn2_2_1_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn2_2_1_R2.fastq.gz |
| bCalAnn2 |        3 |         1 | /home/assembly_downstream/data/arima_HiC/bCalAnn2_3_1_R1.fastq.gz | /home/assembly_downstream/data/arima_HiC/bCalAnn2_3_1_R2.fastq.gz |

当前工作流代码 #

import pandas as pd
import os
import yaml

configfile: "config.yaml"

samples = pd.read_table(config["samples"], dtype=str)

rule all:
    input:
        expand(config["arima_mapping"] + "final/{sample}_{unit_bio}_{unit_tech}.bam", zip, 
            sample=samples["sample"], unit_bio=samples["unit_bio"], unit_tech=samples["unit_tech"])


..
some rules
..

rule add_read_groups:
    input:
        config["arima_mapping"] + "paired/{sample}_{unit_bio}_{unit_tech}.bam"
    output:
        config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam"
    params:
        platform = "ILLUMINA",
        sampleName = "{sample}",
        library = "{sample}",
        platform_unit ="None"
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/paired_read_groups/{sample}_{unit_bio}_{unit_tech}.log"
    shell:
        "picard AddOrReplaceReadGroups I={input} O={output} SM={params.sampleName} LB={params.library} PU={params.platform_unit} PL={params.platform} 2> {log}"

rule merge_tech_repl:
    input:
        config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam" 
    output:
        config["arima_mapping"] + "merge_tech_repl/{sample}_{unit_bio}_{unit_tech}.bam"
    params:
        val_string = "SILENT" 
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/merged_tech_repl/{sample}_{unit_bio}_{unit_tech}.log"
    threads:
        2 #verwendet nur maximal 2
    shell:
        "picard MergeSamFiles -I {input} -O {output} --ASSUME_SORTED true --USE_THREADING true --VALIDATION_STRINGENCY {params.val_string} 2> {log}"


rule mark_duplicates:
    input:
        config["arima_mapping"] + "merge_tech_repl/{sample}_{unit_bio}_{unit_tech}.bam" if config["tech_repl"] else config["arima_mapping"] + "paired_read_groups/{sample}_{unit_bio}_{unit_tech}.bam"
    output:
        bam = config["arima_mapping"] + "final/{sample}_{unit_bio}_{unit_tech}.bam",
        metric = config["arima_mapping"] + "final/metric_{sample}_{unit_bio}_{unit_tech}.txt"
    #params:
    conda:
        "../envs/arima_mapping.yaml"
    log:
        config["logs"] + "arima_mapping/mark_duplicates/{sample}_{unit_bio}_{unit_tech}.log"
    shell:
        "picard MarkDuplicates I={input} O={output.bam} M={output.metric} 2> {log}"

现有方案的问题 #

目前我在配置文件中用一个布尔值tech_repl控制mark_duplicates规则的输入来源（是从merge_tech_repl还是add_read_groups获取），但这种方式不够灵活：实际场景中可能部分样本有任意数量的重复，其他样本没有重复，无法针对性处理。

我想要实现的逻辑是：

• 检查sample.tsv，当样本名+unit_bio编号相同，但unit_tech编号不同时，合并这些技术重复样本；
• 无重复的样本直接跳过合并规则；
• 后续还要实现生物学重复的类似逻辑（unit_bio不同但样本名相同的情况）。

初步尝试的问题 #

我尝试用输入函数来动态获取重复样本，但目前的写法无法返回所有匹配的重复样本列表，而是逐个返回，不符合合并工具的输入要求；同时也不知道如何让无重复的样本跳过合并规则：

input_function(wildcards):
        return expand("{sample}_{unit_bio}_{i}.bam", sample = wildcards.sample,
                                                    unit_bio = wildcards.unit_bio,
                                                    i = samples["sample"].str.count(wildcards.sample))

rule tech_duplicate_check:
    input:
        input_function #(that returns a list of 2-n duplicates, where n could be different for each sample)
    output:
        {sample}_{unit_bio}.bam
    shell:
        MergeTechDupl_tool {input} # input is a list

解决方案 #

1. 预处理样本表，生成分组信息 #

在工作流开头对样本表进行分组，标记哪些样本组需要合并技术重复：

import pandas as pd
import os
import yaml

configfile: "config.yaml"

samples = pd.read_table(config["samples"], dtype=str)

# 生成技术重复分组：按sample和unit_bio分组，收集对应的unit_tech
tech_repl_groups = samples.groupby(["sample", "unit_bio"])["unit_tech"].apply(list).reset_index()
# 标记需要合并的组（unit_tech数量>1）
tech_repl_groups["needs_merge"] = tech_repl_groups["unit_tech"].apply(lambda x: len(x) > 1)

# 生成最终需要的输出文件列表：区分是否需要合并
final_outputs = []
for _, row in tech_repl_groups.iterrows():
    sample = row["sample"]
    unit_bio = row["unit_bio"]
    if row["needs_merge"]:
        # 合并后输出文件不带unit_tech
        final_outputs.append(f"{config['arima_mapping']}final/{sample}_{unit_bio}.bam")
    else:
        # 无重复，沿用原文件名格式
        unit_tech = row["unit_tech"][0]
        final_outputs.append(f"{config['arima_mapping']}final/{sample}_{unit_bio}_{unit_tech}.bam")

2. 更新rule all #

rule all:
    input:
        final_outputs

3. 重写merge_tech_repl规则 #

使用输入函数动态获取该组下的所有技术重复文件，并调整输出文件名（去掉unit_tech）：

rule merge_tech_repl:
    input:
        lambda wildcards: expand(
            f"{config['arima_mapping']}paired_read_groups/{{sample}}_{wildcards.unit_bio}_{unit_tech}.bam",
            sample=wildcards.sample,
            unit_tech=tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]["unit_tech"].iloc[0]
        )
    output:
        f"{config['arima_mapping']}merge_tech_repl/{{sample}}_{{unit_bio}}.bam"
    params:
        val_string = "SILENT" 
    conda:
        "../envs/arima_mapping.yaml"
    log:
        f"{config['logs']}arima_mapping/merged_tech_repl/{{sample}}_{{unit_bio}}.log"
    threads: 2
    shell:
        "picard MergeSamFiles {' '.join([f'-I {f}' for f in input])} -O {output} --ASSUME_SORTED true --USE_THREADING true --VALIDATION_STRINGENCY {params.val_string} 2> {log}"

4. 调整mark_duplicates规则的输入逻辑 #

用输入函数判断当前样本组是否需要合并，自动选择输入来源：

def get_mark_duplicates_input(wildcards):
    # 检查当前sample+unit_bio是否需要合并技术重复
    group = tech_repl_groups[
        (tech_repl_groups["sample"] == wildcards.sample) & 
        (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
    ]
    if group["needs_merge"].iloc[0]:
        return f"{config['arima_mapping']}merge_tech_repl/{wildcards.sample}_{wildcards.unit_bio}.bam"
    else:
        unit_tech = group["unit_tech"].iloc[0][0]
        return f"{config['arima_mapping']}paired_read_groups/{wildcards.sample}_{wildcards.unit_bio}_{unit_tech}.bam"

rule mark_duplicates:
    input:
        get_mark_duplicates_input
    output:
        bam = lambda wildcards: f"{config['arima_mapping']}final/{wildcards.sample}_{wildcards.unit_bio}.bam" if tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]["needs_merge"].iloc[0] else f"{config['arima_mapping']}final/{wildcards.sample}_{wildcards.unit_bio}_{tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]['unit_tech'].iloc[0][0]}.bam",
        metric = lambda wildcards: f"{config['arima_mapping']}final/metric_{wildcards.sample}_{wildcards.unit_bio}.txt" if tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]["needs_merge"].iloc[0] else f"{config['arima_mapping']}final/metric_{wildcards.sample}_{wildcards.unit_bio}_{tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]['unit_tech'].iloc[0][0]}.txt"
    conda:
        "../envs/arima_mapping.yaml"
    log:
        lambda wildcards: f"{config['logs']}arima_mapping/mark_duplicates/{wildcards.sample}_{wildcards.unit_bio}.log" if tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]["needs_merge"].iloc[0] else f"{config['logs']}arima_mapping/mark_duplicates/{wildcards.sample}_{wildcards.unit_bio}_{tech_repl_groups[
                (tech_repl_groups["sample"] == wildcards.sample) & 
                (tech_repl_groups["unit_bio"] == wildcards.unit_bio)
            ]['unit_tech'].iloc[0][0]}.log"
    shell:
        "picard MarkDuplicates I={input} O={output.bam} M={output.metric} 2> {log}"

5. 生物学重复合并的扩展思路 #

生物学重复的逻辑类似，只需按sample分组，收集对应的unit_bio，标记需要合并的组，然后新增merge_bio_repl规则，调整后续规则的输入输出即可。

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