You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

如何从Python对象生成Snakemake通配符?基因组流程开发

用Snakemake通配符替代动态规则实现参考序列下载

我正在使用Snakemake开发基因组分析流程,随着输入输出类型日益多样,希望借助Python对象让代码保持清晰且易于扩展。目前我通过Python循环生成动态规则,想将其转换为Snakemake通配符的实现方式,但不清楚具体方法。


一、Python参考序列类定义

class Reference:
    def __init__(self, name, species, source, genome_seq, genome_seq_url, transcript_seq, transcript_seq_url, annotation_gtf, annotation_gtf_url, annotation_gff, annotation_gff_url) -> None:
        
        self.name = name
        self.species = species
        self.source = source
        self.genome_seq = genome_seq
        self.genome_seq_url = genome_seq_url
        self.transcript_seq = transcript_seq
        self.transcript_seq_url = transcript_seq_url
        self.annotation_gtf = annotation_gtf
        self.annotation_gtf_url = annotation_gtf_url
        self.annotation_gff = annotation_gff
        self.annotation_gff_url = annotation_gff_url

二、参考序列配置CSV文件

name,species,source,genome_seq,genome_seq_url,transcript_seq,transcript_seq_url,annotation_gtf,annotation_gtf_url,annotation_gff,annotation_gff_url
BDGP6_46,FruitFly,Ensembl,Drosophila_melanogaster.BDGP6.46.dna.toplevel.fa.gz,https://ftp.ensembl.org/pub/release-111/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.46.dna.toplevel.fa.gz,Drosophila_melanogaster.BDGP6.46.cdna.all.fa.gz,https://ftp.ensembl.org/pub/release-111/fasta/drosophila_melanogaster/cdna/Drosophila_melanogaster.BDGP6.46.cdna.all.fa.gz,Drosophila_melanogaster.BDGP6.46.111.gtf.gz,https://ftp.ensembl.org/pub/release-111/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.46.111.gtf.gz,Drosophila_melanogaster.BDGP6.46.111.gff3.gz,https://ftp.ensembl.org/pub/release-111/gff3/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.46.111.gff3.gz

三、原始Snakefile代码(含动态规则生成)

import csv
import pathlib

def get_references(references_path:str) -> dict:
    refs_table = dict()
    
    with open(references_path, 'r') as file:
        reader = csv.DictReader(file)
        for row in reader:
            ref_data = Reference(
                row['name'],
                row['species'],
                row['source'],
                row['genome_seq'],
                row['genome_seq_url'],
                row['transcript_seq'],
                row['transcript_seq_url'],
                row['annotation_gtf'],
                row['annotation_gtf_url'],
                row['annotation_gff'],
                row['annotation_gff_url']
            )

            refs_table[row['name']] = ref_data

        return refs_table
    
references_table = get_references('references.csv')

rule all:
    input:
        genome_seq     = expand("../resources/references/{ref_name}/{genome_seq}",         zip, 
                                genome_seq=[references_table[ref].genome_seq for ref in references_table.keys()], 
                                ref_name=[references_table[ref].name for ref in references_table.keys()]),

        transcript_seq = expand("../resources/references/{ref_name}/{transcript_seq}", zip, 
                                transcript_seq=[references_table[ref].transcript_seq for ref in references_table],
                                ref_name=[references_table[ref].name for ref in references_table]),

        annotation_gtf = expand("../resources/references/{ref_name}/{annotation_gtf}", zip, 
                                annotation_gtf=[references_table[ref].annotation_gtf for ref in references_table], 
                                ref_name=[references_table[ref].name for ref in references_table]),

        annotation_gff = expand("../resources/references/{ref_name}/{annotation_gff}", zip, 
                                annotation_gff=[references_table[ref].annotation_gff for ref in references_table.keys()], 
                                ref_name=[references_table[ref].name for ref in references_table.keys()]),

# 动态生成规则的代码
for ref_name, ref in references_table.items(): 

    pathlib.Path(f"../resources/references/{ref_name}/").mkdir(parents=True, exist_ok=True) 
    pathlib.Path(f"../logs/download/refs/").mkdir(parents=True, exist_ok=True) 
    pathlib.Path(f"../times/download/refs/").mkdir(parents=True, exist_ok=True) 

    genome_seq     = f"../resources/references/{ref_name}/{ref.genome_seq}"
    transcript_seq = f"../resources/references/{ref_name}/{ref.transcript_seq}"
    annotation_gtf = f"../resources/references/{ref_name}/{ref.annotation_gtf}"
    annotation_gff = f"../resources/references/{ref_name}/{ref.annotation_gff}"

    log_file  = f"../logs/download/refs/{ref_name}.txt"
    time_file = f"../times/download/refs/{ref_name}.txt"

    genome_seq_url     = ref.genome_seq_url
    transcript_seq_url = ref.transcript_seq_url
    annotation_gtf_url = ref.annotation_gtf_url
    annotation_gff_url = ref.annotation_gff_url

    rule_name = f"download_reference_{ref_name}"

    rule:
        name : rule_name
        output: 
            genome_seq = genome_seq,
            transcript_seq = transcript_seq,
            annotation_gtf = annotation_gtf, 
            annotation_gff = annotation_gff
        params:
            genome_seq_url     = genome_seq_url,
            transcript_seq_url = transcript_seq_url,
            annotation_gtf_url = annotation_gtf_url,
            annotation_gff_url = annotation_gff_url,
        log:
            log_file = log_file
        benchmark:
            time_file
        container:
            "dockers/general_image"
        threads: 
            1
        message:
            "Downloading {params.genome_seq_url} and {params.transcript_seq_url} and {params.annotation_gtf_url} and {params.annotation_gff_url}"
        shell:
            """ 
            wget {params.genome_seq_url} -O {output.genome_seq}           &> {log.log_file}
            wget {params.transcript_seq_url} -O {output.transcript_seq}   &> {log.log_file}
            wget {params.annotation_gtf_url} -O {output.annotation_gtf}   &> {log.log_file}
            wget {params.annotation_gff_url} -O {output.annotation_gff}   &> {log.log_file}
            """

四、改造后的通配符实现方案

核心逻辑

用{ref_name}作为通配符,通过全局的references_table动态获取对应参考序列的属性,替代循环生成规则的冗余逻辑,贴合Snakemake的设计范式。

完整改造代码

import csv
import pathlib

class Reference:
    def __init__(self, name, species, source, genome_seq, genome_seq_url, transcript_seq, transcript_seq_url, annotation_gtf, annotation_gtf_url, annotation_gff, annotation_gff_url) -> None:
        
        self.name = name
        self.species = species
        self.source = source
        self.genome_seq = genome_seq
        self.genome_seq_url = genome_seq_url
        self.transcript_seq = transcript_seq
        self.transcript_seq_url = transcript_seq_url
        self.annotation_gtf = annotation_gtf
        self.annotation_gtf_url = annotation_gtf_url
        self.annotation_gff = annotation_gff
        self.annotation_gff_url = annotation_gff_url

def get_references(references_path:str) -> dict:
    refs_table = dict()
    
    with open(references_path, 'r') as file:
        reader = csv.DictReader(file)
        for row in reader:
            ref_data = Reference(
                row['name'],
                row['species'],
                row['source'],
                row['genome_seq'],
                row['genome_seq_url'],
                row['transcript_seq'],
                row['transcript_seq_url'],
                row['annotation_gtf'],
                row['annotation_gtf_url'],
                row['annotation_gff'],
                row['annotation_gff_url']
            )

            refs_table[row['name']] = ref_data

        return refs_table
    
references_table = get_references('references.csv')
REF_NAMES = list(references_table.keys())

# 预创建所有需要的目录
for ref_name in REF_NAMES:
    pathlib.Path(f"../resources/references/{ref_name}/").mkdir(parents=True, exist_ok=True)
pathlib.Path(f"../logs/download/refs/").mkdir(parents=True, exist_ok=True)
pathlib.Path(f"../times/download/refs/").mkdir(parents=True, exist_ok=True)

rule all:
    input:
        expand("../resources/references/{ref_name}/{genome_seq}", 
               ref_name=REF_NAMES,
               genome_seq=lambda wildcards: references_table[wildcards.ref_name].genome_seq),
        expand("../resources/references/{ref_name}/{transcript_seq}", 
               ref_name=REF_NAMES,
               transcript_seq=lambda wildcards: references_table[wildcards.ref_name].transcript_seq),
        expand("../resources/references/{ref_name}/{annotation_gtf}", 
               ref_name=REF_NAMES,
               annotation_gtf=lambda wildcards: references_table[wildcards.ref_name].annotation_gtf),
        expand("../resources/references/{ref_name}/{annotation_gff}", 
               ref_name=REF_NAMES,
               annotation_gff=lambda wildcards: references_table[wildcards.ref_name].annotation_gff),

rule download_reference:
    output:
        genome_seq = "../resources/references/{ref_name}/{genome_seq}",
        transcript_seq = "../resources/references/{ref_name}/{transcript_seq}",
        annotation_gtf = "../resources/references/{ref_name}/{annotation_gtf}",
        annotation_gff = "../resources/references/{ref_name}/{annotation_gff}",
    params:
        genome_seq_url=lambda wildcards: references_table[wildcards.ref_name].genome_seq_url,
        transcript_seq_url=lambda wildcards: references_table[wildcards.ref_name].transcript_seq_url,
        annotation_gtf_url=lambda wildcards: references_table[wildcards.ref_name].annotation_gtf_url,
        annotation_gff_url=lambda wildcards: references_table[wildcards.ref_name].annotation_gff_url,
        genome_seq=lambda wildcards: references_table[wildcards.ref_name].genome_seq,
        transcript_seq=lambda wildcards: references_table[wildcards.ref_name].transcript_seq,
        annotation_gtf=lambda wildcards: references_table[wildcards.ref_name].annotation_gtf,
        annotation_gff=lambda wildcards: references_table[wildcards.ref_name].annotation_gff,
    log:
        "../logs/download/refs/{ref_name}.txt"
    benchmark:
        "../times/download/refs/{ref_name}.txt"
    container:
        "dockers/general_image"
    threads: 1
    message:
        "Downloading reference files for {wildcards.ref_name}"
    shell:
        """
        # 可选:若预创建目录可省略此步骤
        mkdir -p ../resources/references/{wildcards.ref_name}/ ../logs/download/refs/ ../times/download/refs/
        
        wget {params.genome_seq_url} -O {output.genome_seq} &> {log}
        wget {params.transcript_seq_url} -O {output.transcript_seq} &> {log}
        wget {params.annotation_gtf_url} -O {output.annotation_gtf} &> {log}
        wget {params.annotation_gff_url} -O {output.annotation_gff} &> {log}
        """

关键改动说明

  1. 简化rule all的expand逻辑:不再使用zip,通过lambda wildcards从references_table中动态匹配ref_name对应的文件名,代码更简洁。
  2. 单规则复用替代动态生成:创建一个带{ref_name}通配符的download_reference规则,所有参考序列的下载逻辑复用同一规则,消除冗余代码。
  3. 动态参数与路径绑定:通过lambda wildcards从全局配置中提取对应参考序列的URL和文件名,实现配置与业务逻辑分离。
  4. 统一预创建目录:在规则执行前批量创建所需目录,也可选择在shell命令中添加mkdir -p确保目录存在。

内容的提问来源于stack exchange,提问作者Zingo

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.06.27 09:35:55