如何用Snakemake批量处理多物种ID文件生成对应.id文件？

需求说明 #

• 输入目录input_data下有多个命名为{species}_multicopy_Ids.tsv的文件，每文件首列是uniq_id（如OG开头），第二列是逗号分隔的protein_id
• 需将每个文件的protein_id拆分后，写入对应uniq_id命名的.id文件，每个物种的.id文件单独存放在out_dir/{species}目录下
• 已实现单文件处理脚本，但批量配置失败

现有脚本问题分析 #

1. Shell脚本问题：原脚本中用awk -F "\t" '{print $5}'取protein_id列，但实际该列是第二列，应取$2；且脚本直接将.id文件输出到当前目录，无法指定物种专属输出目录
2. Snakefile问题：
   • 输入文件名拼写不一致：glob匹配的是_multicopy_Ids.tsv，但rule里输入是_multicopy_GeneIDs.tsv
   • dynamic使用不当，Snakemake中动态输出更适合用directory配合明确路径
   • 未确保输出目录out_dir/{species}提前创建

修正后的解决方案 #

调整后的Shell脚本（00_scripts/myshellscript.sh） #

#!/bin/bash

# 接收两个参数：输入文件路径、输出目录路径
input_file="$1"
output_dir="$2"

# 创建输出目录（如果不存在）
mkdir -p "$output_dir"

# 跳过表头（输入文件无表头可删除此句）
tail -n +2 "$input_file" | while read -r line
do
    # 提取uniq_id（第一列）
    uniq_id=$(echo "$line" | awk -F "\t" '{print $1}')
    # 提取protein_id列（第二列）
    protein_ids=$(echo "$line" | awk -F "\t" '{print $2}')
    # 拆分并写入对应.id文件，输出到指定目录
    echo "$protein_ids" | tr -d ' ' | tr ',' '\n' > "${output_dir}/${uniq_id}.id"
done

最终Snakefile #

import os
from snakemake.io import glob_wildcards

# 定义目录
input_dir = "input_data"
output_root = "out_dir"

# 获取所有物种名
species_list = glob_wildcards(os.path.join(input_dir, "{species}_multicopy_Ids.tsv")).species

# 目标：所有物种的输出目录下的.id文件
rule all:
    input:
        expand(os.path.join(output_root, "{species}", "{uniq_id}.id"),
               species=species_list,
               uniq_id=get_uniq_ids(input_dir, species_list))

# 规则：处理单个物种的输入文件
rule process_species:
    input:
        os.path.join(input_dir, "{species}_multicopy_Ids.tsv")
    output:
        directory(os.path.join(output_root, "{species}"))
    params:
        out_dir=os.path.join(output_root, "{species}")
    shell:
        "bash 00_scripts/myshellscript.sh {input} {params.out_dir}"

# 辅助函数：获取所有物种的uniq_id列表
def get_uniq_ids(input_dir, species_list):
    uniq_ids = set()
    for species in species_list:
        input_file = os.path.join(input_dir, f"{species}_multicopy_Ids.tsv")
        with open(input_file, 'r') as f:
            next(f)  # 跳过表头
            for line in f:
                uniq_id = line.strip().split('\t')[0]
                uniq_ids.add(uniq_id)
    return sorted(uniq_ids)

关键说明 #

1. 辅助函数get_uniq_ids：遍历所有输入文件提取uniq_id，确保rule all能明确所有目标输出文件
2. directory输出：用directory标记物种输出目录，Snakemake会自动识别该目录下所有文件为规则输出
3. 参数传递：将输出目录作为参数传递给shell脚本，确保每个物种的文件写入对应专属目录

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