I currently have an MGF file containing MS2 spectral data (QE_2706_229_sequest_high_conf.mgf). The file template is here, as well as a snippet of example:BEGIN IONSTITLE=File3249 Spectrum10594 scans: 11084PEPMASS=499.59366 927079.3CHARGE=3+RTINSECONDS=1710SCANS=11084104.053180 3866.360000110.071530 178805.000000111.068610 1869.210000111.074780 10738.600000112.087240 13117.900000113.071150 7148.790000114.102690 4146.490000115.086840 11835.600000116.070850 6230.980000... ...END IONSThis unannotated spectral file contains thousands of these entries, the total file size is ~150 MB.I then have a series of text files which I need to parse. Each file is similar to the format above, with the first column being read into a NumPy array. Then the unannotated spectra file is parsed for each entry until a matching array is found from the annotated text files input.(Filename 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.0Once a match is found, an MGF annotated file is written that then contains the full entry information in the unannotated file, but with a line that specifies the filename of the annotated text file that matched that particular entry. The output is below:BEGIN IONSSEQ=GRPGPVAGHHQMPRTITLE=File3249 Spectrum10594 scans: 11084PEPMASS=499.59366 927079.3... ...END IONSThere may be a much more computationally inexpensive way to parse. Given 2,000 annotated files to search through, with the above large unannotated file, parsing currently takes ~ 12 hrs on a 2.6 GHz quad-core Intel Haswell CPU.import numpy as npimport subprocess as spimport sysfrom pyteomics import mgf, auxiliarydef main():    pep_files = []    if (len(sys.argv) > 0):        spec_in = sys.argv[1]    else:        print 'Not enough Command Line Arguments!'    path = '/DeNovo/QE_2706_229_sequest_high_conf.mgf'    print spec_in    pep_files.append(spec_in)    for ann_spectra in pep_files:        seq = ann_spectra[:ann_spectra.find('-') - 1]        print seq        a = np.genfromtxt(ann_spectra, dtype=float, invalid_raise=False, usemask=False, filling_values=0.0, usecols=(0))                      b = np.delete(a, 0)        entries = []        with mgf.read(path) as reader:            for spectrum in reader:                   if np.array_equal(b, spectrum['m/z array']):                    entries.append(spectrum)                    file_name = 'DeNovo/good_training_seq/{}.mgf'.format(ann_spectra[:-4])                    with open(file_name, 'wb') as mgf_out:                        for entry in entries:                            mgf_out.write('BEGIN IONS')                            mgf_out.write('\nSEQ={}'.format(seq))                            mgf_out.write('\nTITLE={}'.format(entry['params']['title']))                            mgf_out.write('\nPEPMASS={} {}'.format(entry['params']['pepmass'][0], entry['params']['pepmass'][1]))                            mgf_out.write('\nCHARGE={}'.format(entry['params']['charge']))                            mgf_out.write('\nRTINSECONDS={}'.format(str(entry['params']['rtinseconds'])))                            mgf_out.write('\nSCANS={}'.format(entry['params']['scans']))                            mgf_out.write('\n')                            p = np.vstack([entry['m/z array'], entry['intensity array']])                            output = p.T                            np.savetxt(mgf_out, output, delimiter=' ', fmt='%f')                            mgf_out.write('END IONS')if __name__ == '__main__':    main()The Python script is used with command line arguments with the following bash script:for f in *.txt ; do python2.7 mgf_parser.py $f; doneThis was used to be able to only parse a given number of files at a time. Suggestions on any more efficient parsing methods?
Faster way to parse file to array, compare to array in second file, write final file
python;performance;bioinformatics
Disclaimer : I have very limited skill/knowledge for anything related to np. Also, I haven't run your code to identify bottleneck which should be the very first thing to do before trying to perform optimisations.Common ways to make things faster are :not to do something if you don't need to.not to do something multiple times if once is enough.In the first category, your import subprocess as sp and from pyteomics import auxiliary are not required : you can get rid of this. This shouldn't change anything from a performance point of view but it's always good to make things easier.In the second category, it seems like you could open filename once for each ann_spectra.Also, for each entry, you could retrieve entry['params'] only once. Similarly, you could call mgf_out once.This is the code I have at this stage :def main():    pep_files = []    if (len(sys.argv) > 0):        spec_in = sys.argv[1]    else:        print 'Not enough Command Line Arguments!'    path = '/DeNovo/QE_2706_229_sequest_high_conf.mgf'    print spec_in    pep_files.append(spec_in)    for ann_spectra in pep_files:        seq = ann_spectra[:ann_spectra.find('-') - 1]        print seq        a = np.genfromtxt(ann_spectra, dtype=float, invalid_raise=False, usemask=False, filling_values=0.0, usecols=(0))                      b = np.delete(a, 0)        entries = []        file_name = 'DeNovo/good_training_seq/{}.mgf'.format(ann_spectra[:-4])        with open(file_name, 'wb') as mgf_out, mgf.read(path) as reader:            for spectrum in reader:                   if np.array_equal(b, spectrum['m/z array']):                    entries.append(spectrum)                    for entry in entries:                        param = entry['params']                        mgf_out.write('BEGIN IONS\nSEQ={}\nTITLE={}\nPEPMASS={} {}\nCHARGE={}\nRTINSECONDS={}\nSCANS={}\n').format(                            seq,                            param['title'],                            param['pepmass'][0],                            param['pepmass'][1],                            param['charge'],                            str(param['rtinseconds']),                            param['scans'])                        p = np.vstack([entry['m/z array'], entry['intensity array']])                        output = p.T                        np.savetxt(mgf_out, output, delimiter=' ', fmt='%f')                        mgf_out.write('END IONS')Now, I have realised something that I find quite confusing : you keep adding elements to entries and then you loop over it : it will loop once the first time, twice the second time, n times the n_th times. Is this something we really need/want to do ?Finally, you are doing things in an un-usual way : you have pep_files = [], then pep_files.append(spec_in). This could easily be written pep_files = [spec_in]. Then, you do for ann_spectra in pep_files:, you might as well get rid of pep_files, ann_spectra and the loop and work directly with spec_in (and maybe you should exit when the argument is not provided).