执行snakemake进行迭代映射答案

【问题标题】：Execute snakemake for iterative mapping执行snakemake进行迭代映射
【发布时间】：2021-07-05 03:01:26
【问题描述】：

我想使用 Snakemake 对我的样本进行迭代映射。我事先不知道样本需要重新映射多少次。对于某些样品，它可能是 2-3 次，对于其他样品，可能是 10 次。如果我理解正确，Snakemake 不能使用 while 循环，但也许某种检查点是可能的？

基本上我想在这个循环中做的是每次调用我的 Illumina 读取的 fasta 序列，直到它不再改变。（这是使用 bowtie2 > samtools view > samtools mpileup > bcftools call > bcftools view > bcftools index > bcftools 共识完成的）

我之前编写了一个 bash 脚本来执行此操作，但 Snakemake 确实可以加快此过程。在 bash 脚本中，我使用了我编写的 R 脚本，用于计算新旧 fasta 文件之间的差异数量。如果这个数字= 0，那么它可以停止循环，如果它不一样，那么它必须重新运行上面的步骤。理想情况下，它应该至少有 3 个循环，最多 15 个循环。

如果有人可以提供帮助，那就太好了！提前致谢

【问题讨论】：

我不确定 snakemake 是否适合执行迭代逻辑，但如果是的话，它可能需要在给定迭代中发现的差异数量在文件名中进行编码。最后一轮将生成名称中带有“0”的文件。一个问题是规则应该能够根据其输出的命名方式提前知道其输入应如何命名，但“0”代不一定是“1”代的子代。
如果事先确定了迭代次数（比如15），那么你可以决定文件名中的“世代数”从14到0，而不考虑实际的差异数，但是如果在“0”代之前达到实际零，您将有无用的步骤。
确实，我确实设置了最大值（否则在某些情况下这永远不会结束）。目前我通过使用 for 循环（for in {0..15}）构建它，如果循环至少运行了 5 次并且没有变化，则中断。

标签： loops iteration mapping bioinformatics snakemake

【解决方案1】：

由于您的步骤序列bowtie2 > samtools view > samtools mpileup ... 必须连续完成，即对于给定的输入 fastq，它不能并行化，我想我会按照您在 bash 脚本中的描述进行操作。类似的东西，主要是伪代码：

rule all:
    input:
        ['a.consensus.fastq', 'b.consensus.fastq', ...]

rule iterative_mapping:
    input:
        fq= '{sample}.fastq',
    output:
        fq= '{sample}.consensus.fastq',
    run:
        i = 0
        while i < 15 and fastq has changed:
            shell("bowtie2 > samtools view ...")
            shell("check_fastq_has_changed.R ...")

这将为不同的{sample}s 并行执行迭代映射。有意义吗？

【讨论】：

确实，现在我这样解决了，但我希望有另一种方法。不幸的是，我觉得这会很困难。
@LauraVP 我很想知道来自现有工作流管理器或完全假设的其他策略。对我来说，上面的解决方案看起来很不错，因为你使用 python 代码来实现逻辑，所以你有很多控制和可读性（比 bash 中的等价物更好）。
@dariober 我想一个“理想”的解决方案会对每个 shell 步骤都有一个单独的规则，以便从 snakemake 作业管理中受益。