在 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.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:
"""
... ...
"""
我正在编写一个用于 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.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:
"""
... ...
"""