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

Question

我在从输入文件名派生变量时遇到了问题，尤其是如果您想根据分隔符进行拆分。我尝试了不同的方法(我无法使用(，但到目前为止唯一有效的方法最终失败了，因为它在寻找变量的所有可能变化(因此也在寻找不存在的输入文件(。

我的问题-输入文件以以下模式命名：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"}
但它会查找示例和sheetid的所有可能组合，因此会查找不存在的输入文件。

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

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

但根本不起作用

然后我试着保留我原来的通配符globSAMPLES, = 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
2. 定义新的输出名称并对｛sample｝的分隔符_使用拆分

有人知道如何处理这个问题吗？感谢

谢谢！

Answer 1

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

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-]+)_(Sd+)_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>
"""

N。B： 我删除了您的规则trim中未使用的params部分。