使用目录或目录中的所有文件作为蛇形的输入

如何解决使用目录或目录中的所有文件作为蛇形的输入

我是蛇形的新手。我想在snakemake中使用目录或目录中的所有文件作为输入。例如，两个不同编号的目录。 bam 文件，

--M1
    M1-1.bam
    M1-2.bam
--M2
    M2-3.bam
    M2-5.bam

我只想合并M1-1.bam，M1-2.bam到M1.bam； M2-3.bam、M2-5.bam 至 M2.bam；我尝试使用wildcards和expand后跟this和this，代码如下，

config.yaml

SAMPLES:
  M1:
    - 1
    - 2
  M2:
    - 3
    - 5
rawdata: path/to/rawdata
outpath: path/to/output
reference: path/to/reference

蛇形文件

configfile:"config.yaml"
SAMPLES=config["SAMPLES"]
REFERENCE=config["reference"]
RAWDATA=config["rawdata"]
OUTPATH=config["outpath"]

ALL_INPUT = []
for key,values in SAMPLES.items():
    ALL_INPUT.append(f"Map/bwa/merge/{key}.bam")
    ALL_INPUT.append(f"Map/bwa/sort/{key}.sort.bam")
    ALL_INPUT.append(f"Map/bwa/dup/{key}.sort.rmdup.bam")
    ALL_INPUT.append(f"Map/bwa/dup/{key}.sort.rmdup.matrix")
    ALL_INPUT.append(f"SNV/Mutect2/result/{key}.vcf.gz")
    ALL_INPUT.append(f"Map/bwa/result/{key}")
    for value in values:
        ALL_INPUT.append(f"Map/bwa/result/{key}/{key}-{value}.bam")
        for num in {1,2}:
            ALL_INPUT.append(f"QC/fastp/{key}/{key}-{value}.R{num}.fastq.gz")

rule all:
    input:
        expand("{outpath}/{all_input}",all_input=ALL_INPUT,outpath=OUTPATH)
        
rule fastp:
    input:
        r1= RAWDATA + "/{key}-{value}.R1.fastq.gz",r2= RAWDATA + "/{key}-{value}.R2.fastq.gz"
    output:
        a1="{outpath}/QC/fastp/{key}/{key}-{value}.R1.fastq.gz",a2="{outpath}/QC/fastp/{key}/{key}-{value}.R2.fastq.gz"
    params:
        prefix="{outpath}/QC/fastp/{key}/{key}-{value}"
    shell:
        """
        fastp -i {input.r1} -I {input.r2} -o {output.a1} -O {output.a2} -j {params.prefix}.json -h {params.prefix}.html
        """

rule bwa:
    input:
        a1="{outpath}/QC/fastp/{key}/{key}-{value}.R1.fastq.gz",a2="{outpath}/QC/fastp/{key}/{key}-{value}.R2.fastq.gz"
    output:
        o1="{outpath}/Map/bwa/result/{key}/{key}-{value}.bam"
    params:
        mem="4000",rg="@RG\\tID:{key}\\tPL:ILLUMINA\\tSM:{key}"
    shell:
        """
        bwa mem -t {threads} -M -R '{params.rg}' {REFERENCE} {input.a1} {input.a2} | samtools view -b -o {output.o1}
        """

## get sample index from raw fastq
key_ids,value_ids = glob_wildcards(RAWDATA + "/{key}-{value}.R1.fastq.gz")
# remove duplicate sample name,and this is useful when there is only one sample input
key_ids = list(set(key_ids))

rule merge:
    input:
        expand("{outpath}/Map/bwa/result/{key}/{key}-{value}.bam",outpath=OUTPATH,key=key_ids,value=value_ids)
    output:
        "{outpath}/Map/bwa/merge/{key}.bam"
    shell:
        """
        samtools merge {output} {input}
        """

合并命令中的{input}将是

M1-1.bam M1-2.bam M1-3.bam M1-5.bam M2-1.bam M2-2.bam M2-3.bam M2-5.bam

实际上，对于M1样本，{input}应该是M1-1.bam M1-2.bam；对于 M2，M2-3.bam M2-5.bam。我也读过 this，但我不知道是否有很多目录，每个目录都有不同的文件。

然后我尝试使用目录作为输入，对于 merge rule

rule mergebam:
    input:
        "{outpath}/Map/bwa/result/{key}"
    output:
        "{outpath}/Map/bwa/merge/{key}.bam"
    log:
        "{outpath}/Map/bwa/log/{key}.merge.bam.log"
    shell:
        """
        samtools merge {output} `ls {input}/*bam` > {log} 2>&1
        """

但这给了我MissingInputException error

Missing input files for rule merge:
/{outpath}/Map/bwa/result/M1

任何想法将不胜感激。

解决方法

我还没有完全解析你的问题，但无论如何我都会试一试......在规则 merge 中你有：

expand("{outpath}/Map/bwa/result/{key}/{key}-{value}.bam",outpath=OUTPATH,key=key_ids,value=value_ids)

这意味着您收集了outpath、key 和 value 的所有组合。

大概您想要value 内的所有组合，每个outpath和key。所以使用：

expand("{{outpath}}/Map/bwa/result/{{key}}/{{key}}-{value}.bam",value=value_ids)

,

如果您将 config.yaml 更改为以下内容，是否可以使用 expand 使实现更容易？

SAMPLES:
  M1:
    - M1-1
    - M2-2
  M2:
    - M2-3
    - M2-5