Snakemake - файл параметров рассматривается как подстановочный знак - PullRequest
Новая
       39

Snakemake - файл параметров рассматривается как подстановочный знак

1 голос
/ 06 февраля 2020

Я написал конвейер в Снейкмейке. Это конвейер ATA C -seq (конвейер биоинформатики для анализа данных геномики из эксперимента c). В основном, до этапа объединения с выравниванием я использую подстановочный знак {sample_id}, чтобы затем переключиться на подстановочный знак {sample} (объединение двух или более sample_ids в один образец).

рабочий DAG здесь (для простоты показан только один образец; оранжевый и синий {sample_id} s объединены в один зеленый {sample}

Все правило выглядит следующим образом : 

configfile: "config.yaml"
SAMPLES_DICT = dict()

with open(config['SAMPLE_SHEET'], "r+") as fil:
    next(fil)
    for lin in fil.readlines():
        row = lin.strip("\n").split("\t")
        sample_id = row[0]
        sample_name = row[1]
        if sample_name in SAMPLES_DICT.keys():
            SAMPLES_DICT[sample_name].append(sample_id)
        else:
            SAMPLES_DICT[sample_name] = [sample_id]

SAMPLES = list(SAMPLES_DICT.keys())
SAMPLE_IDS = [sample_id for sample in SAMPLES_DICT.values() for sample_id in sample]
rule all:
    input:
        # FASTQC output for RAW reads
        expand(os.path.join(config['FASTQC'], '{sample_id}_R{read}_fastqc.zip'),
               sample_id = SAMPLE_IDS,
               read = ['1', '2']),

        # Trimming
        expand(os.path.join(config['TRIMMED'],
                            '{sample_id}_R{read}_val_{read}.fq.gz'),
               sample_id = SAMPLE_IDS,
               read = ['1', '2']),

        # Alignment
        expand(os.path.join(config['ALIGNMENT'], '{sample_id}_sorted.bam'),
               sample_id = SAMPLE_IDS),

        # Merging
        expand(os.path.join(config['ALIGNMENT'], '{sample}_sorted_merged.bam'),
               sample = SAMPLES),

        # Marking Duplicates
        expand(os.path.join(config['ALIGNMENT'], '{sample}_sorted_md.bam'),
               sample = SAMPLES),

        # Filtering
        expand(os.path.join(config['FILTERED'],
                            '{sample}.bam'),
               sample = SAMPLES),
        expand(os.path.join(config['FILTERED'],
                            '{sample}.bam.bai'),
               sample = SAMPLES),

        # multiqc report
        "multiqc_report.html"

    message:
        '\n#################### ATAC-seq pipeline #####################\n'
        'Running all necessary rules to produce complete output.\n'
        '############################################################'

Я знаю, что это слишком грязно, я должен оставить только необходимые биты, но здесь мое понимание змейского дела не получается, потому что я не знаю, что я должен оставить и что Я должен удалить.

Это работает, насколько мне известно, точно так, как я хочу.

Однако я добавил правило: 

rule hmmratac:
    input:
        bam = os.path.join(config['FILTERED'], '{sample}.bam'),
        index = os.path.join(config['FILTERED'], '{sample}.bam.bai')
    output:
        model = os.path.join(config['HMMRATAC'], '{sample}.model'),
        gappedPeak = os.path.join(config['HMMRATAC'], '{sample}_peaks.gappedPeak'),
        summits = os.path.join(config['HMMRATAC'], '{sample}_summits.bed'),
        states = os.path.join(config['HMMRATAC'], '{sample}.bedgraph'),
        logs = os.path.join(config['HMMRATAC'], '{sample}.log'),
        sample_name = '{sample}'
    log:
        os.path.join(config['LOGS'], 'hmmratac', '{sample}.log')
    params:
        genomes = config['GENOMES'],
        blacklisted = config['BLACKLIST']
    resources:
        mem_mb = 32000
    message:
        '\n######################### Peak calling ########################\n'
        'Peak calling for {output.sample_name}\n.'
        '############################################################'
    shell:
        'HMMRATAC -Xms2g -Xmx{resources.mem_mb}m '
        '--bam {input.bam} --index {input.index} '
        '--genome {params.genome} --blacklist {params.blacklisted} '
        '--output {output.sample_name} --bedgraph true &> {log}'

И в правило all, после фильтрации, перед multiq c я добавил: 

    # Peak calling
    expand(os.path.join(config['HMMRATAC'], '{sample}.model'),
           sample = SAMPLES),

Соответствующие фрагменты config.yaml: 

# Path to blacklisted regions
BLACKLIST: "/mnt/data/.../hg38.blacklist.bed"

# Path to chromosome sizes
GENOMES: "/mnt/data/.../hg38_sizes.genome"

# Path to filtered alignment
FILTERED: "alignment/filtered"

# Path to peaks
HMMRATAC: "peaks/hmmratac"

Это ошибка * Я получаю (Это продолжается для каждый ввод и вывод правила). * Технически это предупреждение, но оно останавливает выполнение snakemake, поэтому я называю это ошибкой. 

File path alignment/filtered//mnt/data/.../hg38.blacklist.bed.bam contains double '/'. This is likely unintended. It can also lead to inconsistent results of the file-matching approach used by Snakemake.
WARNING:snakemake.logging:File path alignment/filtered//mnt/data/.../hg38.blacklist.bed.bam contains double '/'. This is likely unintended. It can also lead to inconsistent results of the file-matching approach used by Snakemake.

На самом деле это не ... - я просто не сделал Я чувствую себя в безопасности, предоставляя здесь абсолютный путь.

В течение пары дней я боролся с этой ошибкой. Черновая документация, слушал вступление. Я понимаю, что приведенное выше описание далеко от совершенства (оно огромно b c Я даже не знаю, как это сделать, чтобы привести минимальный воспроизводимый пример ...), но я в отчаянии и надеюсь, что вы можете быть терпеливы со мной .

Буду очень признателен за любые предложения относительно того, как его погуглить, где искать ошибку.

1 Ответ

0 голосов
/ 07 февраля 2020

Технически это предупреждение, но оно останавливает выполнение snakemake, поэтому я называю это ошибкой.

Было бы полезно опубликовать логи от snakemake, чтобы увидеть, завершается ли snakemake завершением ошибка и если да, то какая ошибка.

Однако, в дополнение к предложению Eri c C. использовать wildcards.sample вместо {sample} в качестве имени файла, я думаю, что это довольно подозрительно :

alignment/filtered//mnt/data/.../hg38.blacklist.bed.bam

/mnt/ обычно находится в root файловой системы, и вы добавляете к нему относительный путь (alignment/filtered). Вы уверены, что это правильно?

...