prasopes.py 9.53 KB
Newer Older
Yan's avatar
Yan committed
1
2
3
4
5
6
7
#!/usr/bin/env python3

from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
from matplotlib.widgets import SpanSelector
from PyQt5 import QtCore
from PyQt5 import QtWidgets
8
9
10
#from rawFin2 import load_raw
from rawParse import load_raw
#from scipy.signal import find_peaks
Yan's avatar
Yan committed
11
12
13
14
15
16
17
import sys
import matplotlib
import numpy as np
import logging
matplotlib.use("Qt5Agg")


18
def zoom_factory(ax, base_scale, plot=None, an=None):
Yan's avatar
Yan committed
19
    def zoom_fun(event, chart, ann):
Yan's avatar
Yan committed
20
21
        if event.inaxes == ax:
            if event.button == 'up':
22
                #zoom in
Yan's avatar
Yan committed
23
24
                scale_factor = 1/base_scale
            elif event.button == 'down':
25
                #zoom out
Yan's avatar
Yan committed
26
27
                scale_factor = base_scale
            else:
28
                #should not happen
Yan's avatar
Yan committed
29
30
                scale_factor = 1
                print(event.button)
31
32
33

            if event.key == 'shift':
                data = event.ydata
Yan's avatar
Yan committed
34
35
                new_top = data + (ax.get_ylim()[1] - data)*scale_factor
                ax.set_ylim([new_top*-0.01,new_top])
36
37
            else:
                data = event.xdata
Yan's avatar
Yan committed
38
39
40
41
42
43
44
                x_left = data - ax.get_xlim()[0]
                x_right = ax.get_xlim()[1] - data
                ax.set_xlim([data - x_left*scale_factor,
                            data + x_right*scale_factor])
            if ann is not None:
                if type(chart['x']) != type(None):
                    ann_spec(event.inaxes, chart['x'], chart['y'], ann)
45
            ax.figure.canvas.draw()
Yan's avatar
Yan committed
46
47
48

    fig = ax.get_figure()  # get the figure of interest
    # attach the call back
49
    fig.canvas.mpl_connect('scroll_event', lambda event: zoom_fun(event, plot, an))
Yan's avatar
Yan committed
50
51
52
53
54

    # return the function
    return zoom_fun


Yan's avatar
Yan committed
55
def pan_factory(axis, plot=None, an=None):
56
    def pan_fun(event, ax):
Yan's avatar
Yan committed
57
58
59
60
        if event.dblclick and event.inaxes == ax:
            ax.get_figure()
            ax.autoscale(True)
            ax.set_ylim(ax.get_ylim()[1]*-0.01,ax.get_ylim()[1]*1.1)
Yan's avatar
Yan committed
61
            ann_spec(spectrum, plot['x'], plot['y'], an)
Yan's avatar
Yan committed
62
63
            ax.figure.canvas.draw()
        elif event.button == 1 and event.inaxes == ax:
Yan's avatar
Yan committed
64
            ax.start_pan(event.x, event.y, event.button)
65
            id_drag = fig.canvas.mpl_connect('motion_notify_event', lambda action: drag_fun(action, ax))
Yan's avatar
Yan committed
66
            id_release = fig.canvas.mpl_connect('button_release_event',
67
                                                lambda action: drag_end(action, id_drag, id_release, plot, an, ax))
Yan's avatar
Yan committed
68

69
    def drag_fun(event, ax):
Yan's avatar
Yan committed
70
71
72
        ax.drag_pan(1, 'x', event.x, event.y)
        ax.figure.canvas.draw()

Yan's avatar
Yan committed
73
    def drag_end(event, id_drag, id_release, chart, ann, ax):
Yan's avatar
Yan committed
74
75
76
        if event.button == 1:
            fig.canvas.mpl_disconnect(id_drag)
            fig.canvas.mpl_disconnect(id_release)
Yan's avatar
Yan committed
77
78
79
            if ann is not None:
                if type(chart['x']) != type(None):
                    ann_spec(ax, chart['x'], chart['y'], ann)
80
            ax.figure.canvas.draw()
Yan's avatar
Yan committed
81
82
    fig = axis.get_figure()
    fig.canvas.mpl_connect('button_press_event', lambda action: pan_fun(action, axis))
Yan's avatar
Yan committed
83

Yan's avatar
Yan committed
84
85
86
def pick_masses(x_min, x_max, ax, an, plot):
    ax.set_xlim(x_min, x_max)
    ann_spec(spectrum, plot['x'], plot['y'], an)
Yan's avatar
Yan committed
87
    ax.figure.canvas.draw()
Yan's avatar
Yan committed
88

Yan's avatar
Yan committed
89
def pick_times(x_min, x_max, mpl_spectrum, data_set, mpl_chromatogram, an, ms_s):
Yan's avatar
Yan committed
90
91
    start_scan = None
    end_scan = None 
Yan's avatar
Yan committed
92
93
    for i, j in enumerate(data_set['chrom_dat'][0, :]):
        if j > x_min and start_scan is None:
Yan's avatar
Yan committed
94
            start_scan = i
Yan's avatar
Yan committed
95
        if j > x_max and end_scan is None:
Yan's avatar
Yan committed
96
            end_scan = i
Yan's avatar
Yan committed
97
    mpl_spectrum.clear()
98
    an.clear()
Yan's avatar
Yan committed
99
100
101
102
103
    ms_s['x'] = data_set['masses']
    ms_s['y'] = np.mean(data_set['matrix'][start_scan:end_scan], axis=0)
    spectrum_plot(spectrum, ms['x'], ms['y'], an)
    mpl_chromatogram.clear()
    chrom_plot(mpl_chromatogram, data_set['chrom_dat'][0, :], data_set['chrom_dat'][1, :])
104
105
106
    mpl_chromatogram.plot(data_set['chrom_dat'][0, start_scan:end_scan],
                          data_set['chrom_dat'][1, start_scan:end_scan],
                          'b.')
Yan's avatar
Yan committed
107
108

def populate(mpl_chromatogram, mpl_spectrum, data_set, time_sel, an, ms_s):
109
    an.clear()
Yan's avatar
Yan committed
110
    def update_span_selector(chrom, spec, d_set, ms_spec):
111
112
113
114
115
116
        time_selector = SpanSelector(chrom, 
                                     lambda x_min, x_max: pick_times(x_min, x_max, spec, d_set, chrom, an, ms_spec),
                                     'horizontal',
                                     useblit=True,
                                     rectprops=dict(alpha=0.15, facecolor='purple'),
                                     button=3)
Yan's avatar
Yan committed
117
118
119
120
121
122
123
124
125
126
127
        return time_selector
    mpl_spectrum.clear()
    mpl_chromatogram.clear()
    ms_s['x'] = data_set['masses']
    ms_s['y'] = np.mean(data_set['matrix'], axis=0)
    spectrum_plot(mpl_spectrum, ms_s['x'], ms_s['y'], an)
    chrom_plot(mpl_chromatogram, data_set['chrom_dat'][0, :], data_set['chrom_dat'][1, :])
    time_sel[0] = update_span_selector(mpl_chromatogram, mpl_spectrum, data_set, ms_s)
    mpl_chromatogram.figure.canvas.draw()
    mpl_spectrum.figure.canvas.draw()

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
def ann_spec(ms_spec, mass, intensity, an, ann_limit=0.01):
    def sub_peaks(peakz, x_value, y_value, x, y, coef_x=10, coef_y=10):
        sp = sorted(peakz,
                key=lambda h: y_value[h],
                reverse=True)
        sub=[]
        rx = x/coef_x
        ry = y/coef_y
        for g in sp:
            add = True
            yv = y_value[g]
            xv = x_value[g]
            for u in sub:
                if abs(yv - y_value[u]) < ry and abs(xv - x_value[u]) < rx:
                    add = False
            if add:
                sub.append(g)
        return sub
146

Yan's avatar
Yan committed
147
    borders=ms_spec.get_xlim()
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    wstart = np.argmax(mass>borders[0])
    wend = np.argmax(mass>borders[1])
    if wend == 0:
        wend = np.argmax(mass)
    #peaks = find_peaks(intensity[wstart:wend], height=ms_spec.get_ylim()[1] / 100)
    peaks = []
    l = ms_spec.get_ylim()[1] * ann_limit
    q = intensity
    for i,j in enumerate(q[wstart:wend], start=wstart):
        #wend selects the first element which not to include.
        #Thus q[i+1] always exists as it points to that last element.
        if j > l and j > q[i-1] and j > q[i+1]:
            peaks.append(i)
    s_peaks = sub_peaks(peaks, mass, intensity,
                        ms_spec.get_xlim()[1] - ms_spec.get_xlim()[0],
                        ms_spec.get_ylim()[1] - ms_spec.get_ylim()[0])
164
165
166
167
    
    for j in an:
        j.remove()
    an.clear()
Yan's avatar
Yan committed
168

Yan's avatar
Yan committed
169
    for i in s_peaks:
170
171
172
        an.append(ms_spec.annotate('{0:.2f}'.format(mass[i]),
                  xy=(mass[i], intensity[i]),
                  textcoords='data'))
Yan's avatar
Yan committed
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191

def spectrum_plot(spc, mass, intensity, an):
    spc.plot(mass, intensity)
    spc.set_title("Spectrum:", loc ="right")
    spc.set_xlabel("m/z")
    spc.set_ylabel("ion count")
    spc.set_ylim(spc.get_ylim()[1] * -0.01, )
    spc.ticklabel_format(scilimits=(0, 0), axis='y')
    ann_spec(spc, mass, intensity, an)

def chrom_plot(chrom, times, tic):
    chrom.plot(times, tic)
    chrom.set_ylim(chrom.get_ylim()[1] * -0.011, chrom.get_ylim()[1] * 1.1)
    chrom.set_title("Chromatogram:", loc ="right")
    chrom.set_xlabel("time (min)")
    chrom.set_ylabel("total ion count")
    chrom.ticklabel_format(scilimits=(0, 0), axis='y')

def open_file(chrom, spc, d_set, time_sel, an, ms_s):
Yan's avatar
Yan committed
192
193
194
    filename=QtWidgets.QFileDialog.getOpenFileName(caption = "Open spectrum", filter="Finnigan RAW files (*.raw, *.RAW)")[0]
    if filename is '':
        return
Yan's avatar
Yan committed
195
196
    d_set['chrom_dat'], d_set['masses'], d_set['matrix'] = load_raw(filename)
    populate(chrom, spc, d_set, time_sel, an, ms_s)
Yan's avatar
Yan committed
197
198
199
200


if __name__=="__main__":
    #ds for data_set
Yan's avatar
Yan committed
201
    ds = dict(chrom_dat=None,masses=None,matrix=None)
202
203
    #mass spectrometry set
    ms = dict(x=None,y=None)
Yan's avatar
Yan committed
204

Yan's avatar
Yan committed
205
    p_logger = logging.getLogger('parseLogger')
Yan's avatar
Yan committed
206
    logging.basicConfig()
Yan's avatar
Yan committed
207
    p_logger.setLevel("WARN")
Yan's avatar
Yan committed
208
209
210
211
212
213
214
215
216
217
218

    graph = Figure(figsize=(5,4),dpi=100)
    graph.patch.set_facecolor("None")

    chromatogram = graph.add_subplot(211,facecolor=(1,1,1,0.8))
    spectrum = graph.add_subplot(212,facecolor=(1,1,1,0.8))
    graph.tight_layout()

    mpl_canvas = FigureCanvas(graph)
    mpl_canvas.setStyleSheet("background-color:transparent;")
    mpl_canvas.setAutoFillBackground(False)
219
    mpl_canvas.setFocusPolicy( QtCore.Qt.ClickFocus )
Yan's avatar
Yan committed
220
221
222
223
    
    timeSelector=[None]
    annotation=[]

224
    pan_factory(chromatogram)
Yan's avatar
Yan committed
225
    zoom_factory(chromatogram, 1.15)
226
227
    pan_factory(spectrum, ms, annotation)
    zoom_factory(spectrum, 1.15, ms, annotation)
Yan's avatar
Yan committed
228
    mass_selector = SpanSelector(spectrum, lambda x_min, x_max: pick_masses(x_min, x_max, spectrum, annotation, ms), 'horizontal', useblit=True, rectprops=dict(alpha=0.15, facecolor='purple'), button=3)
Yan's avatar
Yan committed
229
230
231
232
233
234

    app = QtWidgets.QApplication(sys.argv)
    main_window = QtWidgets.QMainWindow()

    file_menu = QtWidgets.QMenu('&File',main_window)
    main_window.menuBar().addMenu(file_menu)
Yan's avatar
Yan committed
235
    file_menu.addAction('&Open..', lambda: open_file(chromatogram, spectrum, ds, timeSelector, annotation, ms), QtCore.Qt.CTRL + QtCore.Qt.Key_O)
Yan's avatar
Yan committed
236
237
238
239
240
241
242
    file_menu.addAction('&Quit', main_window.close, QtCore.Qt.CTRL + QtCore.Qt.Key_Q)

    main_widget = QtWidgets.QWidget(main_window)
    main_window.setCentralWidget(main_widget)

    layout = QtWidgets.QVBoxLayout(main_widget)
    layout.addWidget(mpl_canvas)
243
    mpl_canvas.setFocus()
Yan's avatar
Yan committed
244

Yan's avatar
Yan committed
245
246
247
248

    if len(sys.argv) == 2:
        raw_file=sys.argv[1]
        ds['chrom_dat'],ds['masses'],ds['matrix'] = load_raw(raw_file)
249
        populate(chromatogram, spectrum, ds, timeSelector,annotation, ms)
Yan's avatar
Yan committed
250
251
252
    else:
        spectrum_plot(spectrum, [0], [0], annotation)
        chrom_plot(chromatogram, [0], [0])
Yan's avatar
Yan committed
253
254
255
256

    main_window.show()
    sys.exit(app.exec_())