在snakemake中通配符有问题

编程问答 2022-06-07

问题描述

我遇到通配符无法转换为假定值的麻烦。这是Snakefile：

import pandas as pd

configfile: "config.json"
experiments = pd.read_csv(config["experiments"],sep = '\t')
experiments['Name'] = [filename.split('/')[-1].split('_R' if ',' in filename else '.fa')[0] for filename in experiments['Files']]
name2sample = {experiments.iloc[i]['Name'] : experiments.iloc[i]['Sample'] for i in range(len(experiments))}
mg_experiments = experiments[experiments["Data type"] == 'dna']

def preprocess_input(wildcards):
    # get files with matching names
    df = experiments.loc[experiments['Name'] == wildcards.name,'Files']
    # get first value (in case multiple) and split on commas
    return df.iloc[0].split(',')

def join_reads_input(wildcards):
    df = mg_experiments.loc[mg_experiments['Sample'] == wildcards.sample,'Files']
    names = [filename.split('/')[-1].split('_R' if ',' in filename else '.fa')[0] for filename in df]
    return ['{}/Preprocess/Trimmomatic/quality_trimmed_{}{}.fq'.format(config["output"],name,fr) for name in names
        for files in df for fr in (['_forward_paired','_reverse_paired'] if ',' in files else [''])]

rule all:
    input:
        expand("{output}/Annotation/uniprotinfo.tsv",output = config["output"],sample = experiments["Sample"]),expand("{output}/Annotation/{sample}/protein2cog.tsv",expand("{output}/Preprocess/Trimmomatic/quality_trimmed_{name}{fr}.fq",fr = (['_forward_paired','_reverse_paired'] if experiments["Files"].str.contains(',').tolist() else ''),name = experiments['Name'])

rule preprocess:
    input:
        preprocess_input
    output:
        expand("{{output}}/Preprocess/Trimmomatic/quality_trimmed_{{name}}{fr}.fq",').tolist() else ''))
    threads:
        config["threads"]
    run:
        shell("python preprocess.py -i {reads} -t {threads} -o {output}/Preprocess -adaptdir MOSCA/Databases/illumina_adapters -rrnadbs MOSCA/Databases/rRNA_databases -d {data_type}",data_type = experiments.loc[experiments['Name'] == wildcards.name]["Data type"].iloc[0],reads = ",".join(input))

rule join_reads:
    input:
        join_reads_input
    output:
        expand("{output}/Assembly/{{sample}}/{{sample}}{fr}.fastq",fr = (['_forward','_reverse'] if experiments["Files"].str.contains(',').tolist() else ''))
    run:
        for file in input:
            print(file)
            if 'forward' in file:
                shell("touch {output}/Assembly/{wildcards.sample}/{wildcards.sample}_forward.fastq; cat {file} >> {output}/Assembly/{wildcards.sample}/{wildcards.sample}_forward.fastq",output = config["output"])
            elif 'reverse' in file:
                shell("touch {output}/Assembly/{wildcards.sample}/{wildcards.sample}_reverse.fastq; cat {file} >> {output}/Assembly/{wildcards.sample}/{wildcards.sample}_reverse.fastq",output = config["output"])
            else:
                shell("touch {output}/Assembly/{wildcards.sample}/{wildcards.sample}.fastq; cat {file} >> {output}/Assembly/{wildcards.sample}/{wildcards.sample}.fastq",output = config["output"])

rule assembly:
    input:
        expand("{output}/Assembly/{{sample}}/{{sample}}{fr}.fastq",').tolist() else ''))
    output:
        expand("{output}/Assembly/{{sample}}/contigs.fasta",output = config["output"])
    threads:
        config["threads"]
    run:
        reads = ",".join(input)
        shell("python assembly.py -r {reads} -t {threads} -o {output}/Assembly/{{sample}} -a {assembler}",assembler = config["assembler"])

由于我的笨拙，这可能会非常令人困惑。 rule preprocess运行预处理脚本，rule join_reads将通过示例（在下面的Preprocess/Trimmomatic/quality_trimmed文件中定义）获得的读取（experiments部分）汇总在一起，因此可以将它们一起提交组装。这是配置文件：

{
  "output": "output","threads": 14,"experiments": "experiments.tsv","assembler": "Metaspades"
}

这是experiments.tsv文件：

Files   Sample  Data type   Condition
path/to/mg_R1.fastq,path/to/mg_R2.fastq Sample  dna
path/to/a/0.01/mt_0.01a_R1.fastq,path/to/a/0.01/mt_0.01a_R2.fastq   Sample  mrna    c1
path/to/b/0.01/mt_0.01b_R1.fastq,path/to/b/0.01/mt_0.01b_R2.fastq   Sample  mrna    c1
path/to/c/0.01/mt_0.01c_R1.fastq,path/to/c/0.01/mt_0.01c_R2.fastq   Sample  mrna    c1
path/to/a/1/mt_1a_R1.fastq,path/to/a/1/mt_1a_R2.fastq   Sample  mrna    c2
path/to/b/1/mt_1b_R1.fastq,path/to/b/1/mt_1b_R2.fastq   Sample  mrna    c2
path/to/c/1/mt_1c_R1.fastq,path/to/c/1/mt_1c_R2.fastq   Sample  mrna    c2
path/to/a/100/mt_100a_R1.fastq,path/to/a/100/mt_100a_R2.fastq   Sample  mrna    c3
path/to/b/100/mt_100b_R1.fastq,path/to/b/100/mt_100b_R2.fastq   Sample  mrna    c3
path/to/c/100/mt_100c_R1.fastq,path/to/c/100/mt_100c_R2.fastq   Sample  mrna    c3

这里的问题是：猫报告了MissingOutputException，因为它找不到文件output/Assembly/{wildcards.sample}_forward.fastq（反之亦然）。这意味着wildcards.sample没有转换为“ Sample”，我不明白为什么。但是，尽管cat规则停止了必须再次执行的工作流程，但它仍然设法正确生成文件。从那里可以顺利进行，因为汇编规则已经有其输入文件。

为什么wildcards.sample不能转换为“ Sample”？

解决方法

这里有很多东西。我认为对于您的特定问题，当您使用关键字参数进行外壳处理时，它会阻止snakemake格式化其余的通配符。将// iterate till the minimum of two lists while(index1 < a.size() && index2 < b.size()) { if (a.get( index1 ) > b.get( index2 )) { union.add(i++,b.get( index2++ )); } else { union.add(i++,a.get( index1++ )); } } //add the elements of the bigger list to the union while(index1 < a.size()) { union.add(i++,a.get( index1++ )); } while(index2 < b.size()){ union.add(i++,b.get( index2++ )); } 更改为{output}/Assembly/{wildcards.sample}/{wildcards.sample}_reverse.fastq，并将sample作为参数传递给shell。

其他建议：

  • 将输入函数置于要应用的规则之上。
  • 在输入/输出中使用多个规则而不是复杂的扩展。您可以有两个规则使用相同的输出文件，一个规则使用成对的输入，另一个规则使用不成对的输入。
  • 如果您有只调用shell的run指令，请将其替换为shell。您可以将{output}/Assembly/{sample}/{sample}_reverse.fastq逻辑捕获到params指令中。您可以直接将config [assembler]放入shell格式的令牌中，例如reads=','.join(input)。
  • 在扩展中使用shell: python assembly.py ... -a {config[assembler]}，而不是使用allow_missing来转义通配符格式。
  • {{}}会将文件追加到输出中，即使输出不存在（不需要触摸）。
  • 尝试使您的行数少于100个字符，以便它们能正确显示在stackoverflow或github上。

我认为您可以对逻辑进行很多简化，但是我对您的工具了解不足，无法推荐更多细节。