-
Notifications
You must be signed in to change notification settings - Fork 0
/
Spectrum Alignment
62 lines (61 loc) · 2.56 KB
/
Spectrum Alignment
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
from urllib.request import urlretrieve
from pyopenms import *
gh = "https://raw.githubusercontent.com/OpenMS/pyopenms-extra/master"
urlretrieve (gh + "/src/data/YIC(Carbamidomethyl)DNQDTISSK.mzML", "observed.mzML")
exp = MSExperiment()
MzMLFile().load("observed.mzML", exp)
spectra = exp.getSpectra()
observed_spectrum = spectra[0]
tsg = TheoreticalSpectrumGenerator()
theo_spectrum = MSSpectrum()
p = tsg.getParameters()
p.setValue("add_y_ions", "true")
p.setValue("add_b_ions", "true")
p.setValue("add_metainfo", "true")
tsg.setParameters(p)
peptide = AASequence.fromString("YIC(Carbamidomethyl)DNQDTISSK")
tsg.getSpectrum(theo_spectrum, peptide, 1, 2)
import numpy as np
from matplotlib import pyplot as plt
def mirror_plot(obs_mz, obs_int, theo_mz, theo_int, title):
obs_int = [element / max(obs_int) for element in obs_int] # relative intenstiy
theo_int = [element * -1 for element in theo_int] # invert the intensity for the mirror plot
plt.figure(figsize=(12,8))
plt.bar(obs_mz, obs_int, width = 3.0)
plt.bar(theo_mz, theo_int, width = 3.0)
plt.title(title)
plt.ylabel('intensity')
plt.xlabel('m/z')
obs_mz, obs_int = observed_spectrum.get_peaks()
print(min(obs_mz)) # 212.012451171875
print(max(obs_mz)) # 795.2837524414062
theo_mz, theo_int = [], []
for mz, intensity in zip(*theo_spectrum.get_peaks()):
if mz >= 200.0 and mz <= 800.0:
theo_mz.append(mz)
theo_int.append(intensity)
title = 'Observed vs theoretical spectrum'
mirror_plot(obs_mz, obs_int, theo_mz, theo_int, title)
alignment = []
spa = SpectrumAlignment()
p = spa.getParameters()
p.setValue("tolerance", 0.5)
p.setValue("is_relative_tolerance", "false")
spa.setParameters(p)
spa.getSpectrumAlignment(alignment, theo_spectrum, observed_spectrum)
print("Number of matched peaks: " + str(len(alignment)))
print("ion\ttheo. m/z\tobserved m/z")
for theo_idx, obs_idx in alignment:
ion_name = theo_spectrum.getStringDataArrays()[0][theo_idx].decode()
ion_charge = theo_spectrum.getIntegerDataArrays()[0][theo_idx]
print(ion_name + "\t" + str(ion_charge) + "\t"
+ str(theo_spectrum[theo_idx].getMZ())
+ "\t" + str(observed_spectrum[obs_idx].getMZ()))
theo_mz, theo_int, obs_mz, obs_int = [], [], [], []
for theo_idx, obs_idx in alignment:
theo_mz.append(theo_spectrum[theo_idx].getMZ())
theo_int.append(theo_spectrum[theo_idx].getIntensity())
obs_mz.append(observed_spectrum[obs_idx].getMZ())
obs_int.append(observed_spectrum[obs_idx].getIntensity())
title = 'Observed vs theoretical spectrum (aligned)'
mirror_plot(obs_mz, obs_int, theo_mz, theo_int, title)