-
-
Notifications
You must be signed in to change notification settings - Fork 307
/
MRMMapper.py
105 lines (85 loc) · 4.45 KB
/
MRMMapper.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
import copy, sys
import pyopenms
"""
python pyTOPP/MRMMapper.py --in ../source/TEST/TOPP/MRMMapping_input.chrom.mzML \
--tr ../source/TEST/TOPP/MRMMapping_input.TraML --out MRMMapper.tmp.out \
&& diff MRMMapper.tmp.out ../source/TEST/TOPP/MRMMapping_output.chrom.mzML
"""
def algorithm(chromatogram_map, targeted, precursor_tolerance, product_tolerance, allow_unmapped=True, allow_double_mappings=False):
output = copy.copy(chromatogram_map)
output.clear(False);
empty_chromats = []
output.setChromatograms(empty_chromats);
notmapped = 0
for chrom in chromatogram_map.getChromatograms():
mapped_already = False
for transition in targeted.getTransitions():
if (abs(chrom.getPrecursor().getMZ() - transition.getPrecursorMZ()) < precursor_tolerance and
abs(chrom.getProduct().getMZ() - transition.getProductMZ()) < product_tolerance):
if mapped_already:
this_peptide = targeted.getPeptideByRef(transition.getPeptideRef() ).sequence
other_peptide = chrom.getPrecursor().getMetaValue("peptide_sequence")
print "Found mapping of", chrom.getPrecursor().getMZ(), "/", chrom.getProduct().getMZ(), "to", transition.getPrecursorMZ(), "/",transition.getProductMZ()
print "Of peptide", this_peptide
print "But the chromatogram is already mapped to", other_peptide
if not allow_double_mappings: raise Exception("Cannot map twice")
mapped_already = True
precursor = chrom.getPrecursor();
peptide = targeted.getPeptideByRef(transition.getPeptideRef() )
precursor.setMetaValue("peptide_sequence", peptide.sequence)
chrom.setPrecursor(precursor)
chrom.setNativeID(transition.getNativeID())
if not mapped_already:
notmapped += 1
print "Did not find a mapping for chromatogram", chrom.getNativeID()
if not allow_unmapped: raise Exception("No mapping")
else:
output.addChromatogram(chrom)
if notmapped > 0:
print "Could not find mapping for", notmapped, "chromatogram(s)"
dp = pyopenms.DataProcessing()
# dp.setProcessingActions(ProcessingAction:::FORMAT_CONVERSION)
pa = pyopenms.ProcessingAction().FORMAT_CONVERSION
dp.setProcessingActions(set([pa]))
chromatograms = output.getChromatograms();
for chrom in chromatograms:
this_dp = chrom.getDataProcessing()
this_dp.append(dp)
chrom.setDataProcessing(this_dp)
output.setChromatograms(chromatograms);
return output
def main(options):
precursor_tolerance = options.precursor_tolerance
product_tolerance = options.product_tolerance
out = options.outfile
chromat_in = options.infile
traml_in = options.traml_in
# precursor_tolerance = 0.05
# product_tolerance = 0.05
# out = "/tmp/out.mzML"
# chromat_in = "../source/TEST/TOPP/MRMMapping_input.chrom.mzML"
# traml_in = "../source/TEST/TOPP/MRMMapping_input.TraML"
ff = pyopenms.MRMFeatureFinderScoring()
chromatogram_map = pyopenms.MSExperiment()
fh = pyopenms.FileHandler()
fh.loadExperiment(chromat_in, chromatogram_map)
targeted = pyopenms.TargetedExperiment();
tramlfile = pyopenms.TraMLFile();
tramlfile.load(traml_in, targeted);
output = algorithm(chromatogram_map, targeted, precursor_tolerance, product_tolerance)
pyopenms.MzMLFile().store(out, output);
def handle_args():
import argparse
usage = ""
usage += "\nMRMMapper maps measured chromatograms (mzML) and the transitions used (TraML)"
parser = argparse.ArgumentParser(description = usage )
parser.add_argument('--in', dest="infile", help = 'An input file containing chromatograms')
parser.add_argument("--tr", dest="traml_in", help="TraML input file containt the transitions")
parser.add_argument("--out", dest="outfile", help="Output file with annotated chromatograms")
parser.add_argument("--precursor_tolerance", dest="precursor_tolerance", default=0.1, help="Precursor tolerance when mapping (in Th)", metavar='0.1', type=float)
parser.add_argument("--product_tolerance", dest="product_tolerance", default=0.1, help="Product tolerance when mapping (in Th)", metavar='0.1', type=float)
args = parser.parse_args(sys.argv[1:])
return args
if __name__ == '__main__':
options = handle_args()
main(options)