如何在 Snakemake 表格配置中使用列表，用于描述生物信息学管道的测序单元

如何在 Snakemake 表格配置中使用列表。

我使用 Snakemake Tabular（与 BWA mem 映射）配置来描述我的测序单元（在单独的行上测序的文库）。在分析的下一阶段，我必须合并测序单元（映射的 .bed 文件）并获取合并的 .bam 文件（每个样本一个）。现在我使用 YAML 配置来描述哪些单元属于哪些样本。但我希望为此目的使用表格配置，

我不清楚如何从制表符分隔文件的单元格中编写和调用列表（包含示例信息）。

这就是我的单位表格配置的样子：

Unit    SampleSM    LineID  PlatformPL  LibraryLB   RawFileR1   RawFileR2
sample_001.lane_L1  sample_001  lane_L1 ILLUMINA    sample_001  /user/data/sample_001.lane_L1.R1.fastq.gz   /user/data/sample_001.lane_L1.R2.fastq.gz
sample_001.lane_L2  sample_001  lane_L2 ILLUMINA    sample_001  /user/data/sample_001.lane_L2.R1.fastq.gz   /user/data/sample_001.lane_L2.R2.fastq.gz
sample_001.lane_L8  sample_001  lane_L8 ILLUMINA    sample_001  /user/data/sample_001.lane_L8.R1.fastq.gz   /user/data/sample_001.lane_L8.R2.fastq.gz
sample_002.lane_L1  sample_002  lane_L1 ILLUMINA    sample_002  /user/data/sample_002.lane_L1.R1.fastq.gz   /user/data/sample_002.lane_L1.R2.fastq.gz
sample_002.lane_L2  sample_002  lane_L2 ILLUMINA    sample_002  /user/data/sample_002.lane_L2.R1.fastq.gz   /user/data/sample_002.lane_L2.R2.fastq.gz

这是我的示例 YAML 配置的样子：

samples:
 "sample_001": ["sample_001.lane_L1", "sample_001.lane_L2", "sample_001.lane_L8"]
 "sample_002": ["sample_002.lane_L1", "sample_002.lane_L2"]

我的 Snakemake 代码：

import pandas as pd
import os

workdir: "/user/data/snakemake/"

configfile: "Samples.yaml"

units_table = pd.read_table("Units.tsv").set_index("Unit", drop=False)

rule all:
    input:
        expand('map_folder/{unit}.bam', unit=units_table.Unit),
        expand('merge_bam_folder/{sample}.bam', sample=config["samples"]),

rule map_paired_end:
    input:
        r1 = lambda wildcards: expand(units_table.RawFileR1[wildcards.unit]),
        r2 = lambda wildcards: expand(units_table.RawFileR2[wildcards.unit])
    output:
        bam = 'map_folder/{unit}.bam'
    params: 
        bai = 'map_folder/{unit}.bam.bai',
        ref='/user/data/human_g1k_v37.fasta.gz',
        SampleSM = lambda wildcards: units_table.SampleSM[wildcards.unit],
        LineID = lambda wildcards: units_table.LineID[wildcards.unit],
        PlatformPL = lambda wildcards: units_table.PlatformPL[wildcards.unit],
        LibraryLB = lambda wildcards: units_table.LibraryLB[wildcards.unit]
    threads:
        16  
    shell:
            r"""
                    seqtk mergepe {input.r1} {input.r2}\
                    | bwa mem -M -t {threads} -v 3 \
                    {params.ref} - \
                    -R "@RG\tID:{params.LineID}\tSM:{params.SampleSM}\tPL:{params.PlatformPL}\tLB:{params.LibraryLB}"\
                    | samtools view -u -Sb - \
                    | samtools sort - -m 4G -o {output.bam} 

                    samtools index {output.bam}
                    """

rule samtools_merge_bam:
    input:  
        lambda wildcards: expand('map_folder/{file}.bam', file=config['samples'][wildcards.sample])
    output:
        bam = 'merge_bam_folder/{sample}.bam'
    threads:
        1
    shell:  
                    r"""
                    samtools merge {output.bam} {input}

                    samtools index {output.bam}
                    """

下面这个呢？

我已排除 Samples.yaml，因为我认为鉴于您的示例表，没有必要这样做。

In rule samtools_merge_bam您收集共享同一 SampleSM 的所有 unit-bam 文件。这些unit-bam文件创建于map_paired_end其中 lambda 表达式收集每个单元的 fastq 文件。

另请注意，我已从所有规则中删除了 unit-bam 文件，因为（我认为）这些只是中间文件，可以使用以下命令将它们标记为临时文件：temp() https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#protected-and-temporary-files flag.

import pandas as pd
import os

workdir: "/output/dir" 

units_table = pd.read_table("Units.tsv")
samples= list(units_table.SampleSM.unique())

rule all:
    input:
        expand('merge_bam_folder/{sample}.bam', sample= samples),

rule map_paired_end:
    input:
        r1 = lambda wildcards: units_table.RawFileR1[units_table.Unit == wildcards.unit],
        r2 = lambda wildcards: units_table.RawFileR2[units_table.Unit == wildcards.unit],
    output:
        bam = 'map_folder/{unit}.bam'
    params: 
        bai = 'map_folder/{unit}.bam.bai',
        ref='/user/data/human_g1k_v37.fasta.gz',
        SampleSM = lambda wildcards: list(units_table.SampleSM[units_table.Unit == wildcards.unit]),
        LineID = lambda wildcards: list(units_table.LineID[units_table.Unit == wildcards.unit]),
        PlatformPL = lambda wildcards: list(units_table.PlatformPL[units_table.Unit == wildcards.unit]),
        LibraryLB = lambda wildcards: list(units_table.LibraryLB[units_table.Unit == wildcards.unit]),
    threads:
        16  
    shell:
        r"""
        seqtk mergepe {input.r1} {input.r2}\
        | bwa mem -M -t {threads} -v 3 \
        {params.ref} - \
        -R "@RG\tID:{params.LineID}\tSM:{params.SampleSM}\tPL:{params.PlatformPL}\tLB:{params.LibraryLB}"\
        | samtools view -u -Sb - \
        | samtools sort - -m 4G -o {output.bam} 

        samtools index {output.bam}
        """

rule samtools_merge_bam:
    input:  
        lambda wildcards: expand('map_folder/{unit}.bam',
            unit= units_table.Unit[units_table.SampleSM == wildcards.sample])
    output:
        bam = 'merge_bam_folder/{sample}.bam'
    threads:
        1
    shell:  
        r"""
        samtools merge {output.bam} {input}

        samtools index {output.bam}
        """