如何从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} """
关键改动说明
- 简化
rule all的expand逻辑:不再使用zip,通过lambda wildcards从references_table中动态匹配ref_name对应的文件名,代码更简洁。 - 单规则复用替代动态生成:创建一个带
{ref_name}通配符的download_reference规则,所有参考序列的下载逻辑复用同一规则,消除冗余代码。 - 动态参数与路径绑定:通过
lambda wildcards从全局配置中提取对应参考序列的URL和文件名,实现配置与业务逻辑分离。 - 统一预创建目录:在规则执行前批量创建所需目录,也可选择在shell命令中添加
mkdir -p确保目录存在。
内容的提问来源于stack exchange,提问作者Zingo
相关产品推荐
相关产品推荐

