在 Snakemake 中使用 STAR 共享内存模块进行顺序对齐任务

Using STAR shared memory module in Snakemake for sequential alignment tasks

我正在编写一个用于 scRNAseq 序列处理的 Snakemake 管道,它使用 STAR 作为比对工具。

为每个比对作业将基因组索引加载到内存中非常耗时。由于我有很多内存可供使用,我认为使用 STAR 的“共享内存”模块是可行的。有了它,我可以将基因组加载到内存中并保存在那里,直到完成所有工作。

使用shell这很容易实现,例如here

但是在 Snakemake 管道中连接它似乎是一项不平凡的任务,尤其是 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

如有任何帮助,我们将不胜感激!

the STAR "shared memory" module doesn't require any input or output anything

我不熟悉 STAR 和共享内存选项,但是关于 input/output 的问题可以使用表示步骤已完成的虚拟标志文件轻松解决。例如:

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

我找到了一个类似于@dariober 解决方案的解决方法,使用虚拟文件。

使用 STAR 中的 --outFileNamePrefix 标志,我可以在基因组卸载期间创建一个虚拟日志文件。但这不像 dariober 的回答那么干净,因为 STAR 会自动创建三个文件 Log.final.outLog.outLog.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:
    """
    ... ...
    """