首頁(yè) 猿問(wèn) Snakemake...

Snakemake 從輸入文件名派生多個(gè)變量

Python

慕斯王 2023-10-31 14:19:41

我在從輸入文件名派生變量時(shí)遇到問(wèn)題 - 特別是如果您想根據(jù)分隔符進(jìn)行拆分。我嘗試了不同的方法（我無(wú)法開(kāi)始工作），到目前為止唯一有效的方法最終失敗了，因?yàn)樗鼘ふ易兞康乃锌赡茏兓ㄒ约耙虼瞬淮嬖诘妮斎胛募Ｎ业膯?wèn)題 - 輸入文件按以下模式命名： 18AR1376_S57_R2_001.fastq.gz我一開(kāi)始對(duì)變量的初步定義：SAMPLES, = glob_wildcards("../run_links/{sample}_R1_001.fastq.gz")但這最終導(dǎo)致我的文件全部被命名18AR1376_S57，我想刪除 _S57 （指的是示例表 ID）。我在搜索時(shí)發(fā)現(xiàn)的一種可行的方法是：SAMPLES,SHEETID, = glob_wildcards("../run_links/{sample}_{SHEETID}_R1.001.fastq.gz"}但它會(huì)查找樣本和sheetid的所有可能組合，因此會(huì)查找不存在的輸入文件。然后我嘗試了一種更基本的 python 方法：SAMPLES, = glob_wildcards("../run_links/{sample}_R1_001.fastq.gz")ID,=expand([{sample}.split("_")[0]], sample=SAMPLES)``但這根本不起作用然后我嘗試保留原來(lái)的通配符 glob SAMPLES, = glob_wildcards("../run_links/{sample}_R1_001.fastq.gz") 但定義新的輸出文件名稱（基于我在另一個(gè)論壇中找到的說(shuō)明） - 但這給了我一個(gè)我無(wú)法弄清楚的語(yǔ)法錯(cuò)誤。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 """所以，有兩個(gè)選擇，工作流程開(kāi)始：將樣本拆分為樣本和樣本表 ID定義新的輸出名稱并在_{sample} 的分隔符上使用 split有人對(duì)如何解決這個(gè)問(wèn)題有建議嗎？謝謝

查看完整描述

1 回答

開(kāi)滿天機(jī)

TA貢獻(xiàn)1786條經(jīng)驗(yàn) 獲得超13個(gè)贊

glob_wildcards我將使用一個(gè)簡(jiǎn)單的 python 字典來(lái)定義附加到 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-]+)_(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您的規(guī)則中未使用的部分trim。

反對(duì) 回復(fù) 2023-10-31