优雅地处理可变数量的输入和输出文件

如何解决优雅地处理可变数量的输入和输出文件

您有一个管道，其中一些样本有一个输入文件并生成一个输出文件，而其他样本有两个输入文件并生成两个输出文件。生物信息学领域的典型案例是 NGS 测序文库，其中一些样本是双端测序，其他样本是单端测序。如果您需要修剪读取并对齐它们，您必须考虑输入/输出文件的可变数量。

处理这个问题最合适的方法是什么？我觉得使用检查点有点矫枉过正，因为在我看来检查点有点神秘，但我可能错了...

对于文件的输入/输出数量只能为 1 或 2 的情况，我会这样做：首先根据输入的数量在 if-else 指令中有一个 run文件。如果输入数为 1，则触摸第二个文件。对于以下规则，这次再次使用 if-else 检查第二个文件是否有 0 个字节。

这是一个示例（未测试但应该是正确的）：

import os

samples = {'S1': ['s1.R1.fastq'],'S2': ['s2.R1.fastq','s2.R2.fastq']}

rule all:
    input:
        expand('bam/{sample_id}.bam',sample_id= samples),rule cutadapt:
    input:
        fastq= lambda wc: samples[wc.sample_id],output:
        r1= 'cutadapt/{sample_id}.R1.fastq',r2= touch('cutadapt/{sample_id}.R2.fastq'),run:
        if len(input.fastq) == 1:
            # {output.r1} created,{output.r2} touched
            shell('cutadapt ... -o {output.r1} {input.fastq}')
        else:
            shell('cutadapt ... -o {output.r1} -p {output.r2} {input.fastq}')

rule align:
    input:
        r1= 'cutadapt/{sample_id}.R1.fastq',r2= 'cutadapt/{sample_id}.R2.fastq',output:
        'bam/{sample_id}.bam',run:
        if os.path.getsize(input.r2) == 0:
        # or 
        # if len(samples[wildcards.sample_id]) == 1:
            shell('hisat2 ... {input.r1} > {output.bam}')
        else:
            shell('hisat2 ... -1 {input.r1} -2 {input.r2} > {output.bam}')

这行得通，但我觉得人为地创建第二个文件只是为了让工作流程保持愉快是很尴尬的。有没有更好的解决方案？

解决方法

对于单端和双端情况应有两个单独的规则：

rule single_end:
    input:
        fastq = 's{n}.R1.fastq'
    output:
        r1 = 'cutadapt/S{n}.R1.fastq'
    shell:
        'cutadapt ... -o {output.r1} {input.fastq}'

rule pair_ends:
    input:
        r1 = 's{n}.R1.fastq',r2 = 's{n}.R2.fastq'
    output:
        r1 = 'cutadapt/S{n}.R1.fastq',r2 = 'cutadapt/S{n}.R2.fastq'
    shell:
        'cutadapt ... -o {output.r1} -p {output.r2} {input}'

现在有一个歧义：只要可以应用 pair_ends，也可以应用 single_end。这个问题可以用 ruleorder:

解决

ruleorder pair_ends > single_end

,

我同意这里的检查站是多余的。

您可以使用输入函数和参数函数组合任意输入和命令行。诀窍是 Form.resx 规则的可选 r2 输出，它可能存在也可能不存在，这使工作流的 DAG 复杂化。在这种情况下，使用 Snakemake 的 directories as output 功能可能是合适的：

cutadapt