snakemake:规则的可选输入

2024-05-29 06:46:08 发布

您现在位置:Python中文网/ 问答频道 /正文

我想知道是否有一种方法可以在规则中有可选的输入。 一个例子是排除未成对读取以进行对齐(或仅具有未成对读取)。伪规则示例:

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来处理输入(我相信还有其他方法)。我正努力处理这个问题,这样这个步骤就可以通过一个蛇形包装器(目前是一个定制的)来运行。在

我见过的最接近的例子是https://groups.google.com/d/msg/snakemake/qX7RfXDTDe4/XTMOoJpMAAAJ 以及约翰尼斯的回答。但是我们有一个条件赋值(例如input: 'a' if condition else 'b'),而不是可选的。在

如有任何帮助/指导,我们将不胜感激。在

可选输入也可以帮助处理不同数量的hisat2索引(如这里所述:https://snakemake-wrappers.readthedocs.io/en/stable/wrappers/hisat2.html)。在

编辑

为了澄清潜在的输入:

1)单独使用单端读取并在rU中声明它们。读取示例的文件可能是

^{pr2}$

在本例中,r1和{}可能是空列表,或者根本没有在规则中声明。在

2)使用成对的结束读取并在r1r2中声明它们。读取示例的文件可能是

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方法。我不知道空名单可以用来做这个。谢谢你的回答和评论!在


Tags: 示例readinput规则rufastqr2gz
3条回答

一种解决方案是通过输出文件名或路径传递endedness信息。下面的内容应该与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

^{pr2}$

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"。在

我认为mere的方法是Snakemake中最自然的方法,即让输出文件名包含关于受欢迎程度的信息。另一种需要重复规则而不是条件语句的方法是使用ruleorder指令,例如ruleorder: align_pe > align_se。如果存在可选输入,则将使用更高优先级规则。在

我通常使用extend与空或非空列表一起使用:

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

相关问题 更多 >

    热门问题