如何解决在snakemake中通配符有问题
我遇到通配符无法转换为假定值的麻烦。这是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上。
我认为您可以对逻辑进行很多简化,但是我对您的工具了解不足,无法推荐更多细节。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。