对于 Snakemake 工作流程,我需要操作许多 BAM 文件中的标签,并且希望通过脚本( using the Snakemake script: directive )对它们进行管道处理。我这样做的具体方法是使用 pysam stream processing .
infile = pysam.AlignmentFile("-", "rb")
outfile = pysam.AlignmentFile("-", "wb", template=infile)
for s in infile:
(flowcell, lane) = s.query_name.split(':')[0:2]
rg_id = ".".join([flowcell, lane])
s.set_tag('RG',rg_id,'Z')
outfile.write(s)
这个脚本可以独立运行,但我无法弄清楚如何通过snakemake script
指令集成它。
我更喜欢这种方式来最大限度地减少 IO 和 RAM 使用。
编辑:采用直接加载来修复 RG 标签。
# parameters passed from snakemake
bam_file = snakemake.input[0]
fixed_bam_file = snakemake.output[0]
bamfile = pysam.AlignmentFile(bam_file, "rb")
fixed_bamfile = pysam.AlignmentFile(fixed_bam_file, "wb", template = bamfile)
for i, read in enumerate(bamfile.fetch()):
(flowcell, lane) = read.query_name.split(':')[0:2]
rg_id = ".".join([flowcell, lane])
read.set_tag('RG',rg_id,'Z')
fixed_bamfile.write(read)
if not (i % 100000):
print("Updated the read group for {} reads to {}".format(i, rg_id))
bamfile.close()
fixed_bamfile.close()
编辑:Snakemakes run:
和 shell:
指令设置 workdir:
目录,而 script:
指令相对于执行 Snakefile 的目录进行操作(保持一切整洁)。因此,将流处理器放在 script:
下会出现问题。
最佳答案
使用 shell
而不是 script
指令:
rule all:
input:
expand('{sample}_edited.bam'), sample=['a', 'b', 'c']
rule somename:
input:
'{sample}.bam'
output:
'{sample}_edited.bam'
shell:
'''
cat {input} > python edit_bam.py > {output}
'''
关于Snakemake 脚本访问 stdin/stdout 以进行流处理,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/54206452/