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

Question

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

我使用Snakemake表格(使用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}
                    """

Answer 1

下面这个怎么样？

我已经排除了Samples.yaml，因为我认为根据您的样本表，它是不必要的。

在规则samtools_merge_bam中，您收集共享同一SampleSM的所有单元bam文件。这些单元bam文件是在map_paired_end中创建的，其中的lambda表达式收集每个单元的fastq文件。

还要注意，我已经从所有规则中删除了unit-bam文件，因为(我认为)这些文件只是中间文件，可以使用temp()标志将它们标记为临时文件。

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}
        """

Answer 2

请查看https://github.com/snakemake-workflows的其中一个最佳实践工作流程

例如，dna-seq-gatk-variant-calling工作流程定义了两个tsv文件，一个用于单位，另一个用于样本。这允许您(a)向样本添加更多属性，以及(b)每个样本具有多个单位。