Snakemake 从每个样本目录问题输入 fastq 文件用于宏基因组学分析

问题描述

我正在开发一个新的 snakemake 宏基因组学管道来修剪 fastq 文件，并通过 kraken 运行它们。每个样本都有一个包含正向和反向读取的目录。

Sample_1/r1_paired.fq.gz
Sample_1/r2_paired.fq.gz
Sample_2/r1_paired.fq.gz
Sample_2/r2_paired.fq.gz

我提供了一个用户可以上传的样本表，其中包含样本名称和读取名称。我使用pandas 来解析样本表并提供snakefile 所需的名称。这是我的蛇文件。

 #Extract sample names from CSV
import pandas as pd
import os
df = pd.read_csv("sample_table_test.csv")
print(df)
samples = df.library.to_list()
print("Samples being processed:",samples)
R1 = df.r1_file.to_list()
R2 = df.r2_file.to_list()
print(R1,R2)


rule all:
    input:
        expand("{sample}.bracken",sample=samples),#Trimmomatic to trim paired end reads
rule trim_reads:
    input:
        "{sample}/{R1}","{sample}/{R2}",output:
        "{sample}/{R1}_1_trim_paired.fq.gz","{sample}/{R2}_2_trim_paired.fq.gz",conda:
        "env.yaml",shell:
        "trimmomatic PE -threads 8 {input} {input} {output} {output} SLIDINGWINDOW:4:30 LEADING:2 TRAILING:2 MINLEN:50"

#Kraken2 to bin reads and assign taxonomy
rule kraken2:
    input:
        "{sample}/{R1}_1_trim_paired.fq.gz",output:
        "{sample}_report.txt","{sample}_kraken_cseqs#.fq",shell:
        "kraken2 --gzip-compressed --paired --classified-out {output} {input} {input} --db database/minikraken2_v1_8GB/ --report {sample}_report.txt --threads 1"

#Bracken estimates abundance of a species within a sample
rule bracken:
    input:
        "{sample}_report.txt",output:
        "{sample}.bracken",shell:
        "bracken -d database/minikraken2_v1_8GB/ -i {input} -o {output} -r 150"

我收到以下错误，并一直在努力寻找更好的方法来编写我的蛇文件以避免此问题。非常感谢这里的任何帮助。

WildcardError in line 19 of /Metagenomics/Metagenomics/snakemake/Snakefile:
Wildcards in input files cannot be determined from output files:
'R1'

谢谢！

解决方法

问题出在您的 rule kraken2：

rule kraken2:
    input:
        "{sample}/{R1}_1_trim_paired.fq.gz","{sample}/{R2}_2_trim_paired.fq.gz",output:
        "{sample}_report.txt","{sample}_kraken_cseqs#.fq",

规则中的所有通配符都应从 output 部分确定。每个规则的逻辑是它提供某些文件作为可能的输出。在您的情况下，规则提供文件 "{sample}_report.txt" 和 "{sample}_kraken_cseqs#.fq"，其中 {sample} 成为一个自由级别，并替换为将模式解析为文件名的特定值。现在 Snakemake 可以确定此规则的输入，但前提是它拥有所有信息。好的，{sample} 的值是根据输出定义的，但是 {R1} 和 {R2} 的值是什么？

您有多种选择。第一个是在 output: 中的某处定义这些值。看起来这不是你的情况。第二个选项是全局定义这些值（您可能正在尝试这样做）：

R1 = df.r1_file.to_list()
R2 = df.r2_file.to_list()

在这种情况下，{R1} 和 {R2} 不应是通配符，而是 expand 函数的参数：

rule kraken2:
    input:
        expand("{{sample}}/{R1}_1_trim_paired.fq.gz",R1=R1),expand("{{sample}}/{R1}_1_trim_paired.fq.gz",R2=R2)
    output:
        "{sample}_report.txt",

甚至更好：

expand("{{sample}}/{R}_1_trim_paired.fq.gz",R=R1+R2)

请注意，通配符 {sample} 现在必须放在双大括号中才能与 expand 函数的参数区分开来。

还有其他选项，例如从其他 vildcards 的值中解析 {R1} 的值，例如 lambda wildcards: ...，但我想这不是您所需要的。