python - 将文件解析为数组的更快方法，与第二个文件中的数组进行比较

Question

我目前有一个包含 MS2 光谱数据的 MGF 文件 (QE_2706_229_sequest_high_conf.mgf)。文件模板在下面的链接中，以及示例片段:

http://www.matrixscience.com/help/data_file_help.html

BEGIN IONS
TITLE=File3249 Spectrum10594 scans: 11084
PEPMASS=499.59366 927079.3
CHARGE=3+
RTINSECONDS=1710
SCANS=11084
104.053180 3866.360000
110.071530 178805.000000
111.068610 1869.210000
111.074780 10738.600000
112.087240 13117.900000
113.071150 7148.790000
114.102690 4146.490000
115.086840 11835.600000
116.070850 6230.980000
... ...
END IONS

这个未注释的光谱文件包含数千个这样的条目，总文件大小约为 150 MB。然后我有一系列需要解析的文本文件。每个文件都类似于上面的格式，第一列被读入一个 numpy 数组。然后为每个条目解析未注释的光谱文件，直到从注释的文本文件输入中找到匹配的数组。

(文件名 GRPGPVAGHHQMPR)

       m/z            i matches
 104.05318      3866.4
 110.07153    178805.4
 111.06861      1869.2
 111.07478     10738.6
 112.08724     13117.9
 113.07115      7148.8
 114.10269      4146.5
 115.08684     11835.6
 116.07085      6231.0

一旦找到匹配项，就会写入一个 MGF 注释文件，然后包含未注释文件中的完整条目信息，但其中一行指定与该特定条目匹配的注释文本文件的文件名。输出如下:

BEGIN IONS
SEQ=GRPGPVAGHHQMPR
TITLE=File3249 Spectrum10594 scans: 11084
PEPMASS=499.59366 927079.3
... ...
END IONS

可能有一种计算成本更低的解析方法。给定 2,000 个带注释的文件进行搜索，使用上述大型未注释文件，在 2.6 GHz 四核 Intel Haswell cpu 上解析当前需要大约 12 小时。

下面是工作代码:

import numpy as np
import sys
from pyteomics import mgf
from glob import glob

def main():
    """
    Usage: python mgf_parser
    """

    pep_files = glob('*.txt')
    mgf_file = 'QE_2706_229_sequest_high_conf.mgf'
    process(mgf_file, pep_files)

def process(mgf_file, pep_files):
    """Parses spectra from annotated text file. Converts m/z values to numpy array.

        If spectra array matches entry in MGF file, writes annotated MGF file.
    """

    ann_arrays = {}
    for ann_spectra in pep_files:
        a = np.genfromtxt(ann_spectra, dtype=float, invalid_raise=False, 
                          usemask=False, filling_values=0.0, usecols=(0))              
        b = np.delete(a, 0)
        ann_arrays[ann_spectra] = b

    with mgf.read(mgf_file) as reader:
        for spectrum in reader:
            for ann_spectra, array in ann_arrays.iteritems():       
                if np.array_equal(array, spectrum['m/z array']):
                    print '> Spectral match found for file {}.txt'.format(ann_spectra[:-4])
                    file_name = '{}.mgf'.format(ann_spectra[:-4])
                    spectrum['params']['seq'] = file_name[52:file_name.find('-') - 1]
                    mgf.write((spectrum,), file_name)        


if __name__ == '__main__':
    main()

这过去只能一次解析给定数量的文件。关于更有效的解析方法的建议？

Answer 1

我看到了改进的空间，因为您正在为每个小文件重复解析整个 MGF 文件。如果您重构代码使其仅被解析一次，您可能会获得不错的加速。

以下是我将如何调整您的代码，同时摆脱 bash 循环并使用 mgf.write函数，它可能比 np.savetxt 慢一点，但更易于使用:

from pyteomics import mgf
import sys
import numpy as np

def process(mgf_file, pep_files):
    ann_arrays = {}
    for ann_spectra in pep_files:
        a = np.genfromtxt(ann_spectra, invalid_raise=False,
                          filling_values=0.0, usecols=(0,))              
        b = np.delete(a, 0)
        ann_arrays[ann_spectra] = b

    with mgf.read(mgf_file) as reader:
        for spectrum in reader:
            for ann_spectra, array in ann_arrays.iteritems():
                if np.allclose(array, spectrum['m/z array']):
                # allclose may be better for floats than array_equal
                    file_name = 'DeNovo/good_training_seq/{}.mgf'.format(
                                                      ann_spectra[:-4])
                    spectrum['params']['seq'] = ann_spectra[
                                                      :ann_spectra.find('-') - 1]
                    mgf.write((spectrum,), file_name)    

if __name__ == '__main__':
     pep_files = sys.argv[1:]
     mgf_file = '/DeNovo/QE_2706_229_sequest_high_conf.mgf'
     process(mgf_file, pep_files)

然后为了实现与 bash 循环相同的效果，您可以将其称为

python2.7 mgf_parser.py *.txt

如果展开的参数列表太长，可以使用glob而不是依赖 bash 来扩展它:

from glob import iglob
pep_files = iglob(sys.argv[1])

并且这样调用它以防止被 bash 扩展:

python2.7 mgf_parser.py '*.txt'