snakemake：规则的可选输入

Question

我想知道是否有办法在规则中加入可选输入。一个示例案例是排除未配对的读取以进行对齐（或仅具有未配对的读取）。 pseudo 规则示例：

rule hisat2_align:
    input:
        rU: lambda wildcards: ('-U '+ read_files[wildcards.reads]['unpaired']) if wildcards.read_type=='trimmed' else '',
        r1: lambda wildcards: '-1 '+ read_files[wildcards.reads]['R1'],
        r2: lambda wildcards: '-2 '+ read_files[wildcards.reads]['R2']
    output:
        'aligned.sam'
    params:
        idx: 'index_prefix',
        extra: ''
    shell:
        'hisat2 {params.extra} -x {params.idx} {input.rU} {input.r1} {input.r2}'

在这里，没有修剪读取（rU=''）会导致丢失输入文件错误。我可以通过调整后的 input/shell 语句的重复规则或通过 params 处理输入来解决这个问题（我相信还有其他方法）。我正在尝试巧妙地处理这个问题，以便可以通过 snakemake 包装器（目前是自定义包装器）运行完成此步骤。

我见过的最接近的例子是 https://groups.google.com/d/msg/snakemake/qX7RfXDTDe4/XTMOoJpMAAAJ 和约翰内斯的回答。但是我们有一个条件分配（例如input: 'a' if condition else 'b'）而不是一个可选的。

任何 help/guidance 将不胜感激。

ps。可选输入也可以帮助处理不同数量的 hisat2 索引（如此处所述：https://snakemake-wrappers.readthedocs.io/en/stable/wrappers/hisat2.html）。

编辑

澄清潜在的输入：

1）单独使用单端读取，在rU中声明。示例的读取文件可能是

sample1_single_1.fastq.gz
sample1_single_2.fastq.gz

在这种情况下 r1 和 r2 可能是空列表或根本没有在规则中声明。

2) 使用双端读取并在r1和r2中声明它们。示例的读取文件可能是

sample1_paired_1_R1.fastq.gz
sample1_paired_1_R2.fastq.gz
sample1_paired_2_R1.fastq.gz
sample1_paired_2_R2.fastq.gz

在这种情况下，`rU`` 可能是空列表或根本没有在规则中声明。

3) 一起使用成对和单端读取（例如，trimmomatic 的输出，其中一些对被破坏）。示例的读取文件可能是

sample1_paired_1_R1.fastq.gz
sample1_paired_1_R2.fastq.gz
sample1_paired_2_R1.fastq.gz
sample1_paired_2_R2.fastq.gz
sample1_unpaired_1_R1.fastq.gz
sample1_unpaired_1_R2.fastq.gz
sample1_unpaired_2_R1.fastq.gz
sample1_unpaired_2_R2.fastq.gz

作为解决方案。我最终使用了@timofeyprodanov 方法。我没有意识到一个空列表可以用于此。感谢所有的回答和评论！

Answer 1

一种解决方案是通过输出文件名或路径传递结束信息。像下面这样的东西应该适用于 the existing wrapper:

def get_fastq_reads(wcs):
    if wcs.endedness == 'PE':  # Paired-end
        return ["reads/{wcs.sample}.1.fastq.gz", "reads/{wcs.sample}.2.fastq.gz"]

    if wcs.endedness == 'SE':  # Single-end
        return ["reads/{wcs.sample}.fastq.gz"]

    raise(ValueError("Unrecognized wildcard value for 'endedness': %s" % wcs.endedness))

rule hisat2:
    input:
      reads=get_fastq_reads
    output:
      "mapped/{sample}.{endedness}.bam"
    log:                                # optional
      "logs/hisat2/{sample}.{endedness}.log"
    params:                             # idx is required, extra is optional
      idx="genome",
      extra="--min-intronlen 1000"
    wildcard_constraints:
        endedness="(SE|PE)"
    threads: 8                          # optional, defaults to 1
    wrapper:
      "0.27.1/bio/hisat2"

有了这条规则，就可以将 reads/tardigrade.fastq.gz 映射到

> snakemake mapped/tardigrade.SE.bam

或reads/tardigrade.{1,2}.fastq.gz与

> snakemake mapped/tardigrade.PE.bam

索引注释

我被the note on the index files搞糊涂了，觉得可能是错的。 HISAT2 不接受该参数的文件，而是接受所有索引文件共有的单个前缀，因此应该只有一个参数。文档中的示例 idx="genome.fa" 具有误导性。构建 HISAT2 附带的玩具参考 (22_20-21M.fa) 得到的索引是 22_20-21M_snp.{1..8}.ht2，在这种情况下，人们将使用 idx="22_20-21M_snp".

Answer 2

我认为mere的做法，即让输出的文件名包含结束信息，是Snakemake中最自然的做法。另一种要求您有重复规则但没有条件语句的方法是使用 ruleorder 指令，例如ruleorder: align_pe > align_se。如果可选输入存在，则将使用更高优先级的规则。

Answer 3

我通常使用 expand 和空列表或非空列表来做到这一点：

rule a:
    input:
        expand('filename', proxy=[] if no_input else [None])

snakemake：规则的可选输入

snakemake: optional input for rules

python

bioinformatics

snakemake

编辑

索引注释