有3个不同程序、2个输入目录和一个预定义的.txt文件的Snakemake



这将是一个有点苛刻的问题,但到目前为止,由sankemake doku主演的《死亡》并没有产生预期的结果,我希望有更有经验的人能帮助我度过难关。

让我描述一下我所拥有的:

1:外壳上的三个程序,执行方式如下:

./program_1 $(cat file/from/dir_1/foo.fa) $(cat file/from/dir_2/bar.fa) > output.txt
./program_2 $(cat file/from/dir_1/foo.fa) $(cat file/from/dir_2/bar.fa) > output.txt
./program_3 $(cat file/from/dir_1/foo.fa) $(cat file/from/dir_2/bar.fa) > output.txt

是的,cat是必要的,它们使用文件中的字符串,而不是文件本身。


2:两个目录看起来像这样:

hsa-let-7f-5p.fa
hsa-let-7g-5p.fa
hsa-miR-100-3p.fa
hsa-miR-100-5p.fa
hsa-miR-101-2-5p.fa

(更多(

NM_000044.fa
NM_000059.fa
NM_000075.fa
NM_000088.fa
NM_000103.fa

(更多(


3:具有所有所需输出组合的.txt文件:

hsa-miR-29b-3p__NM_138473_programm_1.txt
hsa-miR-29b-3p__NM_138473_programm_2.txt
hsa-miR-29b-3p__NM_138473_programm_3.txt
hsa-miR-545-3p__NM_002332_programm_1.txt
hsa-miR-545-3p__NM_002332_programm_2.txt
hsa-miR-545-3p__NM_002332_programm_3.txt

(更多(


并不是两个目录中的所有文件都被使用,有些文件被多次使用,我不希望所有的组合,只希望指定一次。输出文件应该是独立的,并根据上面指定的.txt命名。最终的sankefile应该在集群上很好地并行化。


工作流程非常简单:

1.读取输出文件名
2.从两个目录中检索组合所需的输入
3.在所有3个程序上运行它们,并输出预定义的txt。通过将程序stdout管道传输到一个文件中。

完成


但是。。。如何告诉Snakemake?有人告诉我这很方便,但到目前为止运气不好。如果有问题,请提问。提前感谢(:



我根据我所知调整了下面的代码。我现在正在从文件中读出;程序";在顶部;dirs";在";规则一";拥有完整的路径。我还切换了shell命令中的两个输入,因为我很困惑,它们是这样的。是的,我知道文件名是次优的。

out = []
f = open("snakemake_output_small.txt", "r")
for line in f:
out.append(line.replace('n', ''))   
# Get distinct filenames
hsa = set([x.split('__')[0] for x in out])
nm = [x.split('__')[1] for x in out]
nm = set([re.sub('_program_.*', '', x) for x in nm])
program = ['full/path/to/program_1', 'full/path/to/program_2', 'full/path/to/program_3']
# Force snakemake to use exactly these wildcard values
# i.e. do not match by regex
wildcard_constraints:
hsa= '|'.join([re.escape(x) for x in hsa]),
nm= '|'.join([re.escape(x) for x in nm]),
program= '|'.join([re.escape(x) for x in program]),

rule all:
input:
out,

rule one:
input:
hsa= 'full/path/to/dir1/{hsa}.fa',
nm= 'full/path/to/dir2/{nm}.fa',
output:
'{hsa}__{nm}_{program}.txt',
shell:
r"""
{wildcards.program} {input.nm} {input.hsa} > {output}
"""

现在它告诉我:

Building DAG of jobs...
MissingInputException in line 21 of Snakefile:
Missing input files for rule all:
hsa-miR-545-3p__NM_002332_program_3.txt
hsa-miR-29b-3p__NM_138473_program_1.txt
hsa-miR-29b-3p__NM_138473_program_1.txt
hsa-miR-545-3p__NM_002332_program_2.txt
hsa-miR-29b-3p__NM_138473_program_3.txt
hsa-miR-545-3p__NM_002332_program_2.txt

这是100%的用户错误,但我缺少什么。如果我不是完全失明,输入目录是正确的。


我也忘了提到,其中两个程序有自己的参数,必须包括在内。程序2的A-d。以及用于程序3的-P此处/a/parameter/文件。因此shell命令可能需要分开。Sry这一切都很混乱。

像以前一样,如果有问题,请提问。

我更喜欢将输出文件名以及输入和程序的组合放在csv文件中,用panda读取它,并使用数据帧来指导工作流程。在这里,我解析输出文件名以提取相关的输入,这在我看来是模糊的。

为了容纳传递给每个程序的选项,请使用字典,或者更好的是,使用yaml配置文件,并使用program通配符查询该字典:

# Read this list from file:
out= ['hsa-miR-29b-3p__NM_138473_program_1.txt',
'hsa-miR-29b-3p__NM_138473_program_2.txt',
'hsa-miR-29b-3p__NM_138473_program_3.txt',
'hsa-miR-545-3p__NM_002332_program_1.txt',
'hsa-miR-545-3p__NM_002332_program_2.txt',
'hsa-miR-545-3p__NM_002332_program_3.txt']
hsa = set([x.split('__')[0] for x in out])
nm = [x.split('__')[1] for x in out]
nm = set([re.sub('_program_.*', '', x) for x in nm])
# Better to put this in a yaml file and read it with `--configfile progs.yaml`
programs = {'program_1': 
{'path': '/path/to/program_1', 'opts': '-d foo -P bar'}, 
'program_2': 
{'path': '/path/to/program_2', 'opts': '-P spam eggs'}, 
'program_3': 
{'path': '/path/to/program_3', 'opts': ''}
}
wildcard_constraints:
hsa= '|'.join([re.escape(x) for x in hsa]),
nm= '|'.join([re.escape(x) for x in nm]),
program= '|'.join([re.escape(x) for x in programs.keys()]),

rule all:
input:
out,

rule one:
input:
hsa= 'dir1/{hsa}.fa',
nm= 'dir2/{nm}.fa',
output:
'{hsa}__{nm}_{program}.txt',
params:
path= lambda wc: programs[wc.program]['path'],
opts= lambda wc: programs[wc.program]['opts'],
shell:
r"""
{params.path} {params.opts} {input.hsa} {input.nm} > {output}
"""

相关内容

最新更新