snakemake从输入文件名派生多个变量

如何解决snakemake从输入文件名派生多个变量

我在从输入文件名中导出变量时遇到问题-尤其是如果您要基于分隔符进行拆分时。我尝试了不同的方法（我无法使用），但到目前为止，唯一可行的方法最终还是失败了，因为它正在寻找变量的所有可能的变体（因此不存在输入文件）。

我的问题-输入文件以以下模式命名： 18AR1376_S57_R2_001.fastq.gz

一开始我对变量的初始定义：
SAMPLES,= glob_wildcards("../run_links/{sample}_R1_001.fastq.gz")

但是最后我的文件都被命名为18AR1376_S57，我想删除_S57（指的是样本工作表ID）。

我在搜索时发现的一种可行的方法是：
SAMPLES,SHEETID,= glob_wildcards("../run_links/{sample}_{SHEETID}_R1.001.fastq.gz"}
但它会查找sample和sheetid的所有可能组合，因此会查找不存在的输入文件。

然后我尝试了一种更基本的python方法：

SAMPLES,= glob_wildcards("../run_links/{sample}_R1_001.fastq.gz")
ID,=expand([{sample}.split("_")[0]],sample=SAMPLES)``

但这根本不起作用

然后我尝试保留原始的通配符 SAMPLES,= glob_wildcards("../run_links/{sample}_R1_001.fastq.gz") 但要定义新的输出文件名（基于我在另一个论坛上找到的说明）-但这给了我一个我不知道的语法错误。

rule trim:
    input:
        R1 = "../run_links/{sample}_R1_001.fastq.gz",R2 = "../run_links/{sample}_R2_001.fastq.gz"
    params:
        prefix=lambda wildcards,output: 
    output:
        R1_Pair = ["./output/trimmed/{}_R1_trim_PE.fastq".format({sample}.split("_")[0])],#or the below version
        R1_Sing = ["./output/trimmed/{}_R1_trim_SE.fastq".format(a) for a in {sample}.split("_")],R2_Pair = ["./output/trimmed/{}_R2_trim_PE.fastq".format(a) for a in {sample}.split("_")],R2_Sing = ["./output/trimmed/{}_R2_trim_SE.fastq".format(a) for a in {sample}.split("_")]
    resources:
        cpus=8
    log:
        "../logs/trim_{sample}.log"
    conda:
        "envs/trim.yaml"
    shell:
        """
        trimmomatic PE -trimlog {log} -threads {resources.cpus} {input.R1} {input.R2} {output.R1_Pair} {output.R1_Sing} {output.R2_Pair} {output.R2_Sing} HEADCROP:10 ILLUMINACLIP:scripts/Truseq-Adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50
        """

所以，有两个选择，

1. 工作流程的开始：将样品分为样品ID和样品表ID
2. 定义新的输出名称，并在{sample}的定界符_上使用split

有人对此有技巧吗？ 谢谢

谢谢！

解决方法

我将使用一个简单的python字典来定义附加到fastq文件的样本名称，而不是使用glob_wildcards：

import os
import re

d = dict()

fastqPath = "."
for fastqF in [f for f in os.listdir(fastqPath) if(re.match(r"^[\w-]+_R1_001\.fastq\.gz",f))]:
        s = re.search(r"(^[\w-]+)_(S\d+)_R1_(001.fastq.gz)",fastqF)
        samplename = s.group(1)
        fastqFfile = os.path.join(fastqPath,fastqF)
        fastqRfile = os.path.join(fastqPath,s.group(1) + "_" + s.group(2) + "_R2_" + s.group(3))
        if(os.path.isfile(fastqRfile)):
                d[samplename] = {"read1":os.path.abspath(fastqFfile),"read2":os.path.abspath(fastqRfile)}

然后，fastq输入文件非常易于使用：

rule all:
        input:
                expand("output/trimmed/{sample}_R1_trim_PE.fastq",sample=d)

rule trim:
        input:
                R1 = lambda wildcards: d[wildcards.sample]["read1"],R2 = lambda wildcards: d[wildcards.sample]["read2"]
        output:
                R1_Pair="output/trimmed/{sample}_R1_trim_PE.fastq",R1_Sing="output/trimmed/{sample}_R1_trim_SE.fastq",R2_Pair="output/trimmed/{sample}_R2_trim_PE.fastq",R2_Sing="output/trimmed/{sample}_R2_trim_SE.fastq"
        resources:
                cpus=8
        log:
                "../logs/trim_{sample}.log"
        conda:
                "envs/trim.yaml"
        shell:
                """
                trimmomatic PE -trimlog {log} -threads {resources.cpus} {input.R1} {input.R2} {output.R1_Pair} {output.R1_Sing} {output.R2_Pair} {output.R2_Sing} HEADCROP:10 I>
                """

注意：我删除了规则params中未使用的trim部分。