snakemake PICARD 合并 bam 文件

Question

我是使用 snakemake 的新手，在使用 PICARD MergeSamFiles 将 bam 文件合并为一个 bam 文件时遇到问题。我想将 1_sorted.bam 2_sorted.bam ...10_sorted.bam 合并到一个带有目录名的 bam 文件中。

import snakemake.io 
import os.path

PICARD="/data/src/picard.jar"
(SAMPLES,)=glob_wildcards("bam/{sample}_sorted.bam")
NAME=os.path.dirname

def bam_inputs(wildcards):
    files = expand("bam/{sample}_sorted.bam", sample=SAMPLES)
    INPUT = "I="+files 
    return INPUT

rule all:
    input: "bam/{NAME}.bam"

rule merge_bams:
    input: bam_inputs
    output: "bam/{NAME}.bam"
    params: mrkdup_jar="/data/src/picard.jar"
    shell: "java -Xmx16G -jar {params.mrkdup_jar} MergeSamFiles \
    {input} \
    O={output} \
    SORT_ORDER=coordinate \
    ASSUME_SORTED=false \
    USE_THREADING=true"

错误：

Building DAG of jobs...
WildcardError in line 12 of /data/data/Samples/snakemake-example/WGS-test/step3.smk:
Wildcards in input files cannot be determined from output files:
'NAME'

我不知道如何将所有 bam 文件合并为一个，也不知道如何将目录名称设置为最终 bam 文件的变量。请指教。

更新：

import snakemake.io

PICARD="/data/src/picard.jar"
(SAMPLES,)=glob_wildcards("bam/{sample}_sorted.bam")
#NAME=os.path.dirname
NAME="test"

rule all:
    input: "bam/{name}.bam".format(name=NAME)

rule merge_bams:
    input: expand("bam/{sample}_sorted.bam",sample=SAMPLES)
    output: "bam/{name}.bam".format(name=NAME)
    params: mrkdup_jar="/data/src/picard.jar"
    shell: """java -Xmx16G -jar {params.mrkdup_jar} MergeSamFiles \
    {"I=" + input} \
    O={output} \
    SORT_ORDER=coordinate \
    ASSUME_SORTED=false \
    USE_THREADING=true """

ERROR:

RuleException in line 11 of /data/data/Samples/snakemake-example/WGS-test/step3.smk:
NameError: The name '"I=" + input' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}}

MergeSamFiles \
I= sub1_sorted.bam I=sub2_sorted.bam I=sub3_sorted.bam \
O= sub.bam \
SORT_ORDER=coordinate \
        ASSUME_SORTED=false \
        USE_THREADING=true

Answer 1

让我们看看rule all。您需要向snakemake 展示您实际希望构建为目标的文件。没有通配符：只是一些明确的东西。你说应该是bam文件，有目录名？

rule all:
    input: f"bam/{NAME}.bam"

请注意，我使用 f 字符串将 {NAME} 从通配符转换为来自变量 NAME 的精确字符串值。您可以选择任何其他方式来执行此操作，即"bam/{name}.bam".format(name=NAME)

接下来，请记住，现在“all”规则中的{NAME} 和“merge_bams”规则中的{NAME} 是不同的实体，因此它们没有任何共同之处。此外，通配符不一定等于您在第 6 行定义的 NAME 变量。我会以不同的方式调用通配符以避免误解。

还有一件事：我不确定你在 bam_inputs 函数中做了什么：

INPUT = "I="+files 

expand 函数的结果应该足以指定 merge_bams 规则的输入。如果您需要为列表中的每个文件添加“I=”，请尝试在 shell: 部分中正确执行：

rule merge_bams:
    input: bam_inputs
    output: "bam/{NAME}.bam"
    params: mrkdup_jar="/data/src/picard.jar"
    shell: f"""java -Xmx16G -jar {{params.mrkdup_jar}} MergeSamFiles 
        {" ".join(["I=" + s for s in input])} 
        O={{output}} 
        SORT_ORDER=coordinate 
        ASSUME_SORTED=false 
        USE_THREADING=true"""