snakemake 如何正确使用 glob_wilcards

Question

我有很多配对的 fastq 文件，我在运行 trim_galore 包之后遇到问题，因为它用 _1_val_1 和 _2_val_2 命名了 fastq 文件，因为例子： AD50_CTGATCGTA_1_val_1.fq.gz 和 AD50_CTGATCGTA_2_val_2.fq.gz.

我想继续 snakemake 并使用

import os
import snakemake.io
import glob

DIR="AD50"
(SAMPLES,READS,) = glob_wildcards(DIR+"{sample}_{read}.fq.gz")
READS=["1","2"]

rule all:
    input:
        expand(DIR+"{sample}_dedup_{read}.fq.gz",sample=SAMPLES,read=READS)

rule clumpify:
    input:
        r1=DIR+"{sample}_1_val_1.fq.gz",
        r2=DIR+"{sample}_2_val_2.fq.gz"
    output:
        r1out=DIR+"{sample}_dedup_1.fq.gz",
        r2out=DIR+"{sample}_dedup_2.fq.gz"
    shell:
        "clumpify.sh in={input.r1} in2={input.r2} out={output.r1out} out2={output.r2out} dedupe subs=0"

错误是：

Building DAG of jobs...
MissingInputException in line 13 of /home/peterchung/Desktop/Rerun-Result/clumpify.smk:
Missing input files for rule clumpify:
AD50/AD50_CTGATCGTA_2_val_2_val_2.fq.gz
AD50/AD50_CTGATCGTA_2_val_1_val_1.fq.gz

我厌倦了另一种方式，不知何故最接近的是它检测到丢失的输入，比如 AD50_CTGATCGTA_1_val_2.fq.gz 和 AD50_CTGATCGTA_2_val_1.fq.gz 不存在。

我不确定我使用的glob_wildcards函数是否正确，因为它有很多下划线。我累了：

 glob_wildcards(DIR+"{sample}_{read}_val_{read}.fq.gz")

但效果不佳。

Answer 1

Glob 通配符实际上是一个将正则表达式应用于目录列表的包装器。默认情况下，通配符会贪婪地匹配 .*。您需要更具体地指定您的通配符，并确保您的规则输入匹配相同的模式匹配。

查看您的示例：

AD50_CTGATCGTA_2_val_2.fq.gz
^ {sample}         ^ ^{read}

{sample} 通配符会消耗所有内容，直到正则表达式不再匹配为止，直到 val。 {read} 然后只剩下 2.

在 rule all 中，然后请求 {sample}_dedup_{read}.fq.gz，即 {AD50_CTGATCGTA_2_val}_dedup_{2}.fq.gz（留在花括号中以显示通配符所在的位置）。当它与 clumpify 匹配时，您请求作为输入： {AD50_CTGATCGTA_2_val}_2_val_2.fq.gz，这就是您丢失该文件的原因。

要解决问题，您有几种选择：

如果示例应包含 1_val 部分，则您需要更新 clumpify 的输入以匹配您现有的文件名（删除额外的 _2_val 等）
如果示例只应包含 AD50_CTGATCGTA，请构建一个更具体的文件名。考虑添加 wildcard constraints 来限制您的匹配，例如[\d]+ 供阅读。这似乎是你在最后一个例子中得到的结果。

最后，expand function by default applies the product of the supplied iterables. That's why you are getting AD50_CTGATCGTA_1_val_2.fq.gz, for example. You need to add zip as the second argument to override that default and match the order of wildcards returned from glob_wildcards. See here也是。

snakemake 如何正确使用 glob_wilcards

snakemake how to use glob_wilcards properly

python

wildcard

glob

snakemake