Snakemake流程featureCounts规则WildcardError报错求助

Snakemake 通配符的解析逻辑是从待生成的输出文件名，反向推导输入文件里通配符的具体取值：

• 其余规则的输出路径本身包含{sample}占位符，比如samtools_sort的输出为aligned/{sample}.sorted.bam，Snakemake可以直接从待生成的文件名中提取到样本名，匹配对应输入文件，因此不会触发报错。
• featureCounts规则的输出是固定文件名raw_Counts，路径中没有{sample}占位符，Snakemake无法反推输入路径中{sample}的取值，因此抛出通配符错误。

另外该规则的输入逻辑本身存在设计错误：featureCounts需要一次性传入所有样本的排序后BAM文件，才能生成整合的原始计数矩阵，并非单样本逐次运行的规则，即使通配符可以正常解析，原有写法也无法得到正确结果。

1. 修正featureCounts规则 #

移除输入路径中的单样本通配符，用expand生成所有样本的BAM文件列表，同时将GTF路径改为解压后的文本格式GTF（featureCounts不支持直接读取gz压缩的注释文件）：

rule featureCounts:
    input:
        samples=expand("aligned/{sample}.sorted.bam", sample=SAMPLE),   
        gtf="genome/Homo_sapiens.GRCh38.106.gtf"
    output:
        "raw_Counts"
    threads:
        16
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output} {input.samples}"

2. 修复基因组下载规则的逻辑错误 #

原有get_genome_gtf、get_genome_fa规则的shell命令为分行独立字符串，Snakemake不会自动拼接执行，会导致cd genome命令不生效，文件下载到错误路径，同时需要将解压后的GTF、FA文件明确声明为输出，避免规则重复运行：

rule get_genome_gtf:
    "下载人GRCh38版本基因组注释文件"
    output:
        gtf_gz = "genome/Homo_sapiens.GRCh38.106.gtf.gz",
        gtf = "genome/Homo_sapiens.GRCh38.106.gtf"
    shell:
        "cd genome && wget ftp://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz && gunzip -k Homo_sapiens.GRCh38.106.gtf.gz"

rule get_genome_fa:
    "下载人GRCh38版本基因组序列文件"
    output:
        fa_gz = "genome/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz",
        fa = "genome/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa"
    shell:
        "cd genome && wget ftp://ftp.ensembl.org/pub/release-106/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz && gunzip -k Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz"

同步修改rule all中对应GTF的目标路径，将原来的"genome/Homo_sapiens.GRCh38.106.gtf.gz"替换为"genome/Homo_sapiens.GRCh38.106.gtf"。

3. 修复比对规则的管道错误 #

原有HISAT2_align规则的shell命令为分行独立字符串，管道符无法跨字符串生效，会导致比对结果无法正常传给samtools生成BAM，合并为单条完整命令即可：

rule HISAT2_align:
    input:
        read1=rules.FASTP.output.trimmed1,
        read2=rules.FASTP.output.trimmed2,
        index=rules.HISAT2_index.output
    output:
        bam="aligned/{sample}.bam",
        metrics="logs/{sample}_HISATmetrics.txt"
    threads: 16
    shell:
        "hisat2 --threads {threads} -x index -1 {input.read1} -2 {input.read2} 2> {output.metrics} | samtools view -Sbh -o {output.bam} -"

4. 可选优化 #

glob_wildcards提取的样本名列表可能存在重复值，建议在初始化后加一步去重，避免expand生成重复路径：

(SAMPLE,FRR) = glob_wildcards("rawReads/{sample}_{frr}.fastq.gz")
SAMPLE = list(set(SAMPLE))

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