蛇形如何正确使用 glob

问题描述

我有很多成对的 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")

但效果不佳。

解决方法

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

通过你的例子：

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

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

在规则 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 默认应用提供的可迭代对象的乘积。例如，这就是您得到 AD50_CTGATCGTA_1_val_2.fq.gz 的原因。您需要添加 zip 作为第二个参数以覆盖该默认值并匹配从 glob_wildcards 返回的通配符的顺序。另请参阅here。

glob python snakemake wildcard

蛇形如何正确使用 glob_wilcards

问题描述

解决方法