Mar*_*que 6 python subprocess samtools
假设我从shell运行以下命令
{
samtools view -HS header.sam; # command1
samtools view input.bam 1:1-50000000; # command2
} | samtools view -bS - > output.bam # command3
Run Code Online (Sandbox Code Playgroud)
对于那些不熟悉samtools视图的人(因为这是stackoverflow).这实际上是在创建一个具有新标头的新bam文件.bam文件通常是大型压缩文件,因此即使在某些情况下通过文件也可能非常耗时.一种替代方法是进行command2,然后使用samtools reheader来切换标头.这会两次通过大文件.上面的命令一次性通过bam,这对于较大的bam文件是有用的(即使在压缩时它们也会大于20GB - WGS).
我的问题是如何使用subprocess在python中实现这种类型的命令.
我有以下内容:
fh_bam = open('output.bam', 'w')
params_0 = [ "samtools", "view", "-HS", "header.sam" ]
params_1 = [ "samtools", "view", "input.bam", "1:1-50000000"]
params_2 = [ "samtools", "view", "-bS", "-" ]
sub_0 = subprocess.Popen(params_0, stderr=subprocess.PIPE, stdout=subprocess.PIPE)
sub_1 = subprocess.Popen(params_1, stderr=subprocess.PIPE, stdout=subprocess.PIPE)
### SOMEHOW APPEND sub_1.stdout to sub_0.stdout
sub_2 = subprocess.Popen(params_2, stdin=appended.stdout, stdout=fh_bam)
Run Code Online (Sandbox Code Playgroud)
任何帮助是极大的赞赏.谢谢.
我认为由于涉及的文件大小,连接内存中前两个子进程的输出是不可行的。我建议将前两个子进程的输出包装在一个类似的文件中。看起来您只需要 read 方法,因为 popen 只会从其标准输入文件中读取,而不是查找或写入。下面的代码假设从 read 返回一个空字符串足以表明流位于 EOF
class concat(object):
def __init__(self, f1, f2):
self.f1 = f1
self.f2 = f2
def read(self, *args):
ret = self.f1.read(*args)
if ret == '':
ret = self.f2.read(*args)
return ret
fh_bam = open('output.bam', 'w')
params_0 = [ "samtools", "view", "-HS", "header.sam" ]
params_1 = [ "samtools", "view", "input.bam", "1:1-50000000"]
params_2 = [ "samtools", "view", "-bS", "-" ]
sub_0 = subprocess.Popen(params_0, stderr=subprocess.PIPE, stdout=subprocess.PIPE)
sub_1 = subprocess.Popen(params_1, stderr=subprocess.PIPE, stdout=subprocess.PIPE)
### Somehow Append sub_1.stdout to sub_0.stdout
sub_2 = subprocess.Popen(params_2, stdin=concat(sub_0.stdout, sub_1.stdout), stdout=fh_bam)
Run Code Online (Sandbox Code Playgroud)
澄清一下,f1.read将阻塞,并且仅''在管道关闭/EOF'd 时返回。concat.read仅在发生这种情况后尝试读取 from f2,因此 fromf1和 的输出f2不会交织在一起。当然,f1重复读取末尾会产生轻微的开销,这可以通过设置标志变量来指示从哪个文件读取来避免。但我怀疑它会显着提高性能。