Как использовать список в табличной конфигурации Snakemake, для описания блоков секвенирования для биоинформационного конвейера

Question

Как использовать список в табличной конфигурации 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}
                    """

Answer 1

Пожалуйста, посмотрите на один из лучших рабочих процессов в https :: //github.com/snakemake-workflows.

Например, dna-seq-gatk-многоадресного вызова рабочий процесс определяет два файла tsv, один для единиц измерения и один для образцов.Это позволяет (а) добавлять больше атрибутов к выборкам и (б) иметь несколько единиц на выборку.

Answer 2

Что об этом ниже?

Я исключил файл Samples.yaml, так как считаю его необязательным с учетом вашего листа с образцами.

В правиле samtools_merge_bam вы собираете все файлы unit-bam, использующие один и тот же SampleSM.Эти файлы unit-bam создаются в map_paired_end, где лямбда-выражение собирает файлы fastq для каждого модуля.

Обратите внимание, что я удалил файлы unit-bam из правила all, поскольку (я думаю) это всего лишь промежуточные файлы, и их можно пометить как временные, используя флаг 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}
        """