Snakemake 使用输入文件名的一部分作为输出答案

【问题标题】：Snakemake use part of input file name for the outputSnakemake 使用输入文件名的一部分作为输出
【发布时间】：2021-10-23 04:47:39
【问题描述】：

我有多个 .fastq 格式的数据文件，名称类似于

start-Number_NAME-Sample_info-Machine_info-Lab_info-end_R1.fastq.gz（和 R2 类似）所以两个例子可以是

start-1000_John-brain_1-hiseq_2500-LAB_KA-end_R1.fastq.gz
start-1000_John-brain_1-hiseq_2500-LAB_KA-end_R2.fastq.gz

start-1200_Smith-Liver_5-Novaseq_6000-LAB_RH-end_R1.fastq.gz
start-1200_Smith-Liver_5-Novaseq_6000-LAB_RH-end_R2.fastq.gz

它们都有相同的结构，但是它们的数量、名称、样品、机器和实验室信息各不相同，这让我很难制定一个涵盖所有这些的蛇形规则。

我的目标是在这些输入文件上使用工具，例如与 bwa mem 对齐。但是我想创建一个简短的输出名称。这样输入文件将是全名，而输出名称将仅包含 Number_NAME 信息，因此我可以创建 1000_John.bam 和 1200_Smith.bam 之类的文件

我正在使用 Snakemake 创建此管道，我尝试了多种方法但无法使其正常工作。

第一个选项：通过此链接 (Snakemake, how to change output filename when using wildcards)

import pathlib
import glob
import os
indir = pathlib.Path("FASTQ/chr1/")
paths = indir.glob("start-*R?.fastq.gz")
SAMPLES = set([x.stem.split("-")[1] for x in paths]) # ID's
print(SAMPLES)

rule all:
    input:
        expand("output/{sample}_mapped.bam",sample=SAMPLES)

def find_fastq(wildcards):
    fastqs = [str(x) for x in indir.glob(f"{wildcards.sample}*.fastq.gz")]
    return sorted(fastqs)

rule bwa:
    input:
        fastqs = find_fastq
    output:
        mapped = "output/{sample}_mapped.bam"
    params:
        ref = "ref.fa"
    shell:
        "bwa mem {params.ref} {input.fastqs} {input.fastqs} | samtools sort > {output.mapped}"

我得到了正确的 Number_NAME id 集，但是 bwa 规则给了我

规则 bwa 中的错误：工作编号：54 输出：输出/1000_John_mapped.bam 贝壳： bwa 内存参考.fa | samtools 排序 > 输出/1000_John_mapped.bam （其中一个以非零退出代码退出的命令；请注意，snakemake 使用 bash 严格模式！）

所以 fastq 不作为 bwa mem 命令的输入。我已经尝试了几种不同的迭代路径（使用 * 或？）和 SAMPLES 的列表理解。

我还尝试为 R1 和 R2 创建两个 find_fastq 函数，但它不起作用。即使我只将一个文件解析为 bwa mem 它仍应作为单端对齐运行，因此使用 bwa mem {params.ref} {input.fastqs}

第二个选项：

我还尝试简单地创建两个单独的列表，其中包含 ID =["1000_John","1200_Smith"] 和 INFO = ["brain_1-hiseq_2500-LAB_KA-","Liver_5-Novaseq_6000-LAB_RH"] 但是那么如果我在我的 bwa 规则中使用扩展

f1 = expand("start-{sample}_{info}-end_R1.fastq.gz",sample=ID,info=INFO),
f2 = expand("start-{sample}_{info}-end_R2.fastq.gz",sample=ID,info=INFO)

它失败了，因为扩展是两组的乘积。而不仅仅是跨两个列表的相同索引的比较。 SetA = ["A","B"] SetB = ["C","D"] -> 展开 AC、AD、BC、BD，我想要的是 AC 和 BD。

第三个选项：创建一个包含所有不同信息的列表，然后将切片信息解析为扩展函数。

fastqdir = glob.glob('FASTQ/chr1/'+'*_R1.fastq.gz')
print(fastqdir)
SAMPLES = [i.split("-",2) for i in fastqdir] 
SAMPLES = [item for sublist in SAMPLES for item in sublist] #creating a list with all elements
#['start','1000_John','brain_1-hiseq_2500-LAB_KA-end_R1.fastq.gz','start','1200_Smith','Liver_5-Novaseq_6000-LAB_RH-end_R1.fastq.gz']
print(SAMPLES[1::3]) #all ID's

f1 = expand("start-{sample}_{info}",sample=SAMPLES[1::3],info=SAMPLES[2::3]),
f2 = expand("start-{sample}_{info}",sample=SAMPLES[1::3],info=SAMPLES[2::3])

但同样没有用。

质量检查：

有没有办法只使用扩展函数而不是乘积对两个列表的同一索引中的元素进行成对比较？
或者任何人都可以帮助使用路径的第一个选项以及如何在提供全名作为输入的同时正确使用重命名这一想法。其他建议也很感激

【问题讨论】：

标签： python alignment snakemake fastq

【解决方案1】：

有多种方法可以解决此类问题。您的第二种方法并不遥远，除了 expand() 使用 itertools 函数 product 并因此生成两个通配符列表的所有组合。但是，您可以改为使用zip，这将为您提供所需的组合（请参见此处：https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#the-expand-function）。那将是：

expand("start-{sample}_{info}-end_R1.fastq.gz", zip, sample=ID,info=INFO)

我建议在这里考虑一个样本表。您可以设计一个简单的 tsv 文件，其中一列带有您所需的样本名称，另一（两）列包含 fastq 文件的确切路径和名称。可以在 Snakefile 中读取此工作表，例如 pandas 数据帧，并允许您使用输入函数获取每个样本的 fastq 文件。

此选项还允许您轻松地将新样本添加到实验中，而无需更新 Snakefile 中的任何列表。

【讨论】：