使用 `snakemake` 的 MD 任务

MD task with `snakemake`

我想为分子动力学模拟创建一个非常简单的管道。该程序 (Amber) 只需要 3 个文件作为输入,并生成很多文件,其中一些我将来会用到。所以我的管道非常简单:

  1. 检查 *.in*.prmtop*.rst 是否在文件夹中(我保证它只有一个文件用于任何这些扩展),如果这些文件不存在则发出警告
  2. 运行 shell 命令(基于所有输入文件的名称)
  3. 检查是否生成了 *.outmdenmdinfo*.nc

就是这样。这是我处理的程序的标准方法。一个文件夹,一项任务,基于文件用途而非内容的简短文件名。

我写了一个简单的管道:

rule all:
    input: '{inp}.out'

rule amber:
    input:
        '{inp}.in',
        '{top}.prmtop',
        '{coord}.rst'
    output:
        '{inp}.out',
        'mden',
        'mdinfo',
        '{inp}.nc'
    shell:
        'pmemd.cuda'
        ' -O'
        ' -i {inp}.in'
        ' -o {inp}.out'
        ' -p {top}.prmtop'
        ' -c {coord}.rst'
        ' -r {inp}.rst'
        ' -x {inp}.nc'
        ' -ref {coord}.rst'

而且它不起作用。

  1. all 规则中的所有输入必须是明确的。 (为什么?为什么它不能是 regex 或通配符表达式?如果我在文件夹中看到 *.out 并且 shell 脚本的状态代码是 0,就这样,工作完成了)
  2. 我必须在 output 中使用 input 中的所有通配符,但我只想在 shell 或其他规则中使用一些通配符
  3. 我一定不希望得到像 mden 这样具有潜在“非唯一”名称的文件,因为它可能会随着另一项任务而改变,但我知道,这将只是一项任务,而且它是我的 MD 程序如何工作的直接方式(是的,我知道 Ambers 的 -e-inf 键,但它过于复杂的简单任务)。

所以,我想决定是否值得为此使用 snakemake。这是一个非常简单的任务,但我已经花了几个小时,我看到了很多文档,很多例子,但我不能应用到我的案例中。 snakemake 看起来正是我需要的,但是我不能用这个框架笼统地表达简单的任务,我不想指定明确的文件名,因为我会失去灵活性,我想 运行 数百个自动完成的简单任务,只是输入的文件会有所不同。我确定我还没有想出如何处理这个框架。也许你可以告诉我我应该怎么做?谢谢!

希望这会让您朝着正确的方向前进。

如果我没理解错的话,snakemake 的输入是一个包含 amber 输入文件的文件夹。您知道此文件夹包含一个 .in 个文件、一个 .prmtop 个文件和一个 .rst 个文件,但您不知道这些文件的全名。

如果你想让 snakemake 在单个输入文件夹上 运行,那么你根本不需要通配符,下面的脚本就可以了。

import glob
import os

input_folder = config['amber_folder']

# We don't know the name of input file. We only know it ends in '.in'
inp = glob.glob(os.path.join(input_folder, '*.in'))
assert len(inp) == 1
inp = inp[0]

name = os.path.splitext(os.path.basename(inp))[0]
output_folder = name + '_results'

out = os.path.join(output_folder, name + '.out')

rule all:
    input:
        out

rule amber:
    input:
        inp= inp,
        top= glob.glob(os.path.join(input_folder, '*.prmtop')),
        rst= glob.glob(os.path.join(input_folder, '*.rst')),
    output:
        out= out,
        nc= os.path.join(output_folder, name + '.nc'),
        mden= os.path.join(output_folder, 'mden'),
        mdinfo= os.path.join(output_folder, 'mdinfo'),
    shell:
        r"""
        pmemd.cuda \
         -O \
         -i {input.inp} \
         -o {output.out} \
         -p {input.top} \
         -c {input.rst} \
         -r {input.rst} \
         -x {output.nc} \
         -ref {input.rst}
        """

执行方式:

snakemake -j 1 -C amber_folder='your-input-folder'

如果你有很多输入文件夹,你可以编写一个 for 循环来执行上面的命令,但可能更好的方法是将输入列表传递给 snakemake 并让它处理它们。