Использование модуля общей памяти STAR в Snakemake для последовательных задач выравнивания

#python #shared-memory #snakemake

Вопрос:

Я пишу конвейер Snakemake для обработки последовательности scRNAseq, в котором в качестве инструмента выравнивания используется ЗВЕЗДА.

Загрузка индекса генома в память для каждого задания выравнивания занимает очень много времени. Поскольку в моем распоряжении много памяти, я считаю, что можно использовать модуль «общая память» от STAR. С его помощью я могу загрузить геном в память и хранить его там, пока все работы не будут выполнены.

С помощью оболочки этого очень легко достичь, например, здесь.

Но чинить его в конвейере змеиных змеев кажется нетривиальной задачей, особенно если модуль «Общая память» STAR не требует никаких входных или выходных данных. Например, я пытался:

 # Not full pipeline
rule umi_tools_extract:
    input:
        read1="{sample}_R1.fq.gz",
        read2="{sample}_R2.fq.gz",
        whitelist="whitelist.txt"
    output:
        "{sample}_R1_extracted.fq.gz"
    shell:
        """
        umi_tools extract --bc-pattern=CCCCNNNN 
                          --stdin {input.read1} 
                          --read2-in {input.read2} 
                          --stdout {output} 
                          --read2-stdout --filter-cell-barcode --error-correct-cell 
                          --whitelist={input.whitelist}
        """

rule STAR:
    input:
        fq="{sample}_R1_extracted.fq.gz",
        genomeDir=directory("/path/to/genome")
    output:
        "{sample}_Aligned.sortedByCoord.out.bam"
    threads:32
    shell:
        """
        STAR --runThreadN {threads} 
             --genomeLoad LoadAndKeep 
             --genomeDir {input.genomeDir} 
             --readFilesIn {input.fq} 
             --readFilesCommand zcat 
             --limitBAMsortRAM 20000000000 
             --outFilterMultimapNmax 1 
             --outFilterType BySJout 
             --outSAMstrandField intronMotif 
             --outFilterIntronMotifs RemoveNoncanonical 
             --outFilterMismatchNmax 6 
             --outSAMtype BAM SortedByCoordinate 
             --outFileNamePrefix {wildcards.sample}_
        """
... ...

rule STAR_unload:
    input:
        genomeDir=dirctory("/path/to/genome")
#    output:  ## Put something here??
#        ""
    shell:
        """
        STAR --genomeLoad Remove 
             --genomeDir {input.genomeDir} 
        """
 

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

Любая помощь будет признательна!

Комментарии:

1. Этот вопрос, вероятно, следует перенести на bioinformatics.stackexchange.com !

2. Привет @KonradRudolph Перед отправкой вопроса я провел поиск по ключевым словам с помощью «змеиного змея с общей памятью STAR» как в SO, так и в Bioinformatics StackExchange и обнаружил, что у SO больше связанных действий, поэтому я разместил его здесь 🙂 Спасибо за ваше предложение!

3. Такой поиск, вероятно, вводит в заблуждение: переполнение стека содержит значительно больше контента, чем сайт биоинформатики, но (а) большая часть этого контента является исторической, и (б) ваш вопрос утонет в море других тем, в то время как Биоинформатика уделяет больше внимания таким вопросам. Таким образом, вероятность получить хороший ответ намного выше у последнего, даже несмотря на то, что он имеет значительно меньшую общую активность. В любом случае, рад, что ты нашел решение.

Ответ №1:

модуль STAR «общая память» не требует никакого ввода или вывода чего-либо

Я не знаком со ЗВЕЗДОЙ и опцией общей памяти, но проблему с вводом/выводом можно легко решить, используя фиктивные файлы флагов, которые сигнализируют о выполнении шага. Напр.:

 rule index_genome:
    input:
        fa= 'genome.fa',
    output:
        done= touch('index.done'),
    shell:
        r"""
        STAR index ... {input.fa} 
        """

rule load_genome:
    input:
        'index.done',
    output:
        touch('loading.done'),
    shell:
        r"""
        STAR --genomeLoad LoadAndExit --genomeDir ...
        """

rule align:
     input:
         fq= 'reads.fastq',
         idx= 'loading.done',
     output:
         ...
     shell:
         r"""
         STAR ...
         """

rule unload_genome:
    # Delete the loading.done flag file otherwise subsequent runs of the pipeline 
    # will fail to load the genome again if STAR alignment is needed.
    input:
        # Generic function that aggregates all alignment outputs
        bams= get_files(),
        idx= 'loading.done',
    output:
        "logs/STARunload_Log.out"
    shell:
        r"""
        STAR --genomeLoad Remove 
             --genomeDir {input.genomeDir} 
             --outFileNamePrefix logs/STARunload_
        rm {input.idx}
        """"
 

Ответ №2:

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

Используя --outFileNamePrefix флаг STAR , я могу создать фиктивный файл журнала во время выгрузки генома. Но это не так чисто, как ответ Log.final.out дариобера , так как ЗВЕЗДА автоматически создает три файла Log.out Log.progress.out .

 # Unload STAR genome index
rule STAR_unload:
input:
    # Generic function that aggregates all alignment outputs
    bams=get_files(),
    genomeDir=directory("/path/to/index")
output:
    "logs/STARunload_Log.out"
shell:
    """
    STAR --genomeLoad Remove 
         --genomeDir {input.genomeDir} 
         --outFileNamePrefix logs/STARunload_
    """

rule some_downstream_step:
input:
    bam="Aligned.sortedByCoord.out.bam"
    dummy="logs/STARunload_Log.out"
output:
    "genes_assigned.bam"
shell:
    """
    ... ...
    """
 

Комментарии:

1. Обратите внимание, что при выгрузке генома вы также должны удалить файл loading.done флага. Я отредактировал свой ответ, чтобы добавить правило unload_genome .

2. Да, я полагаю, что мог бы также отметить загрузку. сделано, файл помечен как временный? т. е. temp(touch('loading.done')) .