我有一个用 Python 编写的巨大管道,它使用非常大的 .gz 文件(压缩约 14GB),但需要更好的方法将某些行发送到外部软件(来自 blast-legacy/2.2.26 的 formatdb)。我有一个很久以前有人为我写的 Perl 脚本,它的速度非常快,但我需要在 Python 中做同样的事情,因为管道的其余部分是用 Python 编写的,我必须保持这种方式。Perl 脚本使用两个文件句柄,一个用于保存 .gz 文件的 zcat,另一个用于存储软件需要的行(每 4 行中的 2 行)并将其用作输入。它涉及生物信息学,但不需要经验。该文件为 fastq 格式,软件需要 fasta 格式。每 4 行是一个 fastq 记录,取第 1 行和第 3 行并在第 1 行的开头添加“>”,这是 formatdb 软件将用于每条记录的 fasta 等价物。
perl 脚本如下:
#!/usr/bin/perl
my $SRG = $ARGV[0]; # reads.fastq.gz
open($fh, sprintf("zcat %s |", $SRG)) or die "Broken gunzip $!\n";
# -i: input -n: db name -p: program
open ($fh2, "| formatdb -i stdin -n $SRG -p F") or die "no piping formatdb!, $!\n";
#Fastq => Fasta sub
my $localcounter = 0;
while (my $line = <$fh>){
if ($. % 4==1){
print $fh2 "\>" . substr($line, 1);
$localcounter++;
}
elsif ($localcounter == 1){
print $fh2 "$line";
$localcounter = 0;
}
else{
}
}
close $fh;
close $fh2;
exit;
它真的很好用。我怎么能在 Python 中做同样的事情?我喜欢 Perl 如何使用这些文件句柄,但我不确定如何在不创建实际文件的情况下在 Python 中做到这一点。我能想到的就是 gzip.open 文件并将我需要的每条记录的两行写入一个新文件并将其与“formatdb”一起使用,但这太慢了。有任何想法吗?我需要将它放入 python 管道中,所以我不能只依赖 perl 脚本,而且我还想知道一般如何做到这一点。我假设我需要使用某种形式的子流程模块。
这是我的 Python 代码,但它再次变慢并且速度是这里的问题(巨大的文件):
#!/usr/bin/env python
import gzip
from Bio import SeqIO # can recognize fasta/fastq records
import subprocess as sp
import os,sys
filename = sys.argv[1] # reads.fastq.gz
tempFile = filename + ".temp.fasta"
outFile = open(tempFile, "w")
handle = gzip.open(filename, "r")
# parses each fastq record
# r.id and r.seq are the 1st and 3rd lines of each record
for r in SeqIO.parse(handle, "fastq"):
outFile.write(">" + str(r.id) + "\n")
outFile.write(str(r.seq) + "\n")
handle.close()
outFile.close()
cmd = 'formatdb -i ' + str(tempFile) + ' -n ' + filename + ' -p F '
sp.call(cmd, shell=True)
cmd = 'rm ' + tempFile
sp.call(cmd, shell=True)