/usr/share/mmass/mspy/mod_pattern.py is in mmass 5.5.0-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 | # -------------------------------------------------------------------------
# Copyright (C) 2005-2013 Martin Strohalm <www.mmass.org>
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# Complete text of GNU GPL can be found in the file LICENSE.TXT in the
# main directory of the program.
# -------------------------------------------------------------------------
# load libs
import math
import numpy
# load stopper
from mod_stopper import CHECK_FORCE_QUIT
# load blocks
import blocks
# load objects
import obj_compound
import obj_peaklist
# load modules
import calculations
import mod_basics
import mod_signal
import mod_peakpicking
# ISOTOPIC PATTERN FUNCTIONS
# --------------------------
def pattern(compound, fwhm=0.1, threshold=0.01, charge=0, agentFormula='H', agentCharge=1, real=True, model='gaussian'):
"""Calculate isotopic pattern for given compound.
compound (str or mspy.compound) - compound
fwhm (float) - gaussian peak width
threshold (float) - relative intensity threshold for isotopes (in %/100)
charge (int) - charge to be calculated
agentFormula (str or mspy.compound) - charging agent formula
agentCharge (int) - charging agent unit charge
real (bool) - get real peaks from calculated profile
model (gaussian, lorentzian, gausslorentzian) - peak shape function
"""
# check compound
if not isinstance(compound, obj_compound.compound):
compound = obj_compound.compound(compound)
# check agent formula
if agentFormula != 'e' and not isinstance(agentFormula, obj_compound.compound):
agentFormula = obj_compound.compound(agentFormula)
# add charging agent to compound
if charge and agentFormula != 'e':
formula = compound.formula()
for atom, count in agentFormula.composition().items():
formula += '%s%d' % (atom, count*(charge/agentCharge))
compound = obj_compound.compound(formula)
# get composition and check for negative atom counts
composition = compound.composition()
for atom in composition:
if composition[atom] < 0:
raise ValueError, 'Pattern cannot be calculated for this formula! --> ' + compound.formula()
# set internal thresholds
internalThreshold = threshold/100.
groupingWindow = fwhm/4.
# calculate pattern
finalPattern = []
for atom in composition:
# get isotopic profile for current atom or specified isotope only
atomCount = composition[atom]
atomPattern = []
match = mod_basics.ELEMENT_PATTERN.match(atom)
symbol, massNumber, tmp = match.groups()
if massNumber:
isotope = blocks.elements[symbol].isotopes[int(massNumber)]
atomPattern.append([isotope[0], 1.]) # [mass, abundance]
else:
for massNumber, isotope in blocks.elements[atom].isotopes.items():
if isotope[1] > 0.:
atomPattern.append(list(isotope)) # [mass, abundance]
# add atoms
for i in range(atomCount):
CHECK_FORCE_QUIT()
# if pattern is empty (first atom) add current atom pattern
if len(finalPattern) == 0:
finalPattern = _normalize(atomPattern)
continue
# add atom to each peak of final pattern
currentPattern = []
for patternIsotope in finalPattern:
# skip peak under relevant abundance threshold
if patternIsotope[1] < internalThreshold:
continue
# add each isotope of current atom to peak
for atomIsotope in atomPattern:
mass = patternIsotope[0] + atomIsotope[0]
abundance = patternIsotope[1] * atomIsotope[1]
currentPattern.append([mass, abundance])
# group isotopes and normalize pattern
finalPattern = _consolidate(currentPattern, groupingWindow)
finalPattern = _normalize(finalPattern)
# correct charge
if charge:
for i in range(len(finalPattern)):
finalPattern[i][0] = (finalPattern[i][0] - mod_basics.ELECTRON_MASS*charge) / abs(charge)
# group isotopes
finalPattern = _consolidate(finalPattern, groupingWindow)
# get real peaks from profile
if real:
prof = profile(finalPattern, fwhm=fwhm, points=100, model=model)
finalPattern = []
for isotope in mod_signal.maxima(prof):
finalPattern.append(isotope)
centroid = mod_signal.centroid(prof, isotope[0], isotope[1]*0.99)
if abs(centroid-isotope[0]) < fwhm/100.:
finalPattern[-1][0] = centroid
# normalize pattern
finalPattern = _normalize(finalPattern)
# discard peaks below threshold
filteredPeaks = []
for peak in finalPattern:
if peak[1] >= threshold:
filteredPeaks.append(list(peak))
finalPattern = filteredPeaks
return finalPattern
# ----
def gaussian(x, minY, maxY, fwhm=0.1, points=500):
"""Make Gaussian peak.
mz (float) - peak m/z value
minY (float) - min y-value
maxY (float) - max y-value
fwhm (float) - peak fwhm value
points (int) - number of points
"""
# make gaussian
return calculations.signal_gaussian(float(x), float(minY), float(maxY), float(fwhm), int(points))
# ----
def lorentzian(x, minY, maxY, fwhm=0.1, points=500):
"""Make Lorentzian peak.
mz (float) - peak m/z value
minY (float) - min y-value
maxY (float) - max y-value
fwhm (float) - peak fwhm value
points (int) - number of points
"""
# make gaussian
return calculations.signal_lorentzian(float(x), float(minY), float(maxY), float(fwhm), int(points))
# ----
def gausslorentzian(x, minY, maxY, fwhm=0.1, points=500):
"""Make half-Gaussian half-Lorentzian peak.
mz (float) - peak m/z value
minY (float) - min y-value
maxY (float) - max y-value
fwhm (float) - peak fwhm value
points (int) - number of points
"""
# make gaussian
return calculations.signal_gausslorentzian(float(x), float(minY), float(maxY), float(fwhm), int(points))
# ----
def profile(peaklist, fwhm=0.1, points=10, noise=0, raster=None, forceFwhm=False, model='gaussian'):
"""Make profile spectrum for given peaklist.
peaklist (mspy.peaklist) - peaklist
fwhm (float) - default peak fwhm
points (int) - default number of points per peak width (not used if raster is given)
noise (float) - random noise width
raster (1D numpy.array) - m/z raster
forceFwhm (bool) - use default fwhm for all peaks
model (gaussian, lorentzian, gausslorentzian) - peak shape function
"""
# check peaklist type
if not isinstance(peaklist, obj_peaklist.peaklist):
peaklist = obj_peaklist.peaklist(peaklist)
# check raster type
if raster != None and not isinstance(raster, numpy.ndarray):
raster = numpy.array(raster)
# get peaks
peaks = []
for peak in peaklist:
peaks.append([peak.mz, peak.intensity, peak.fwhm])
if forceFwhm or not peak.fwhm:
peaks[-1][2] = fwhm
# get model
shape = 0
if model == 'gaussian':
shape = 0
elif model == 'lorentzian':
shape = 1
elif model == 'gausslorentzian':
shape = 2
# make profile
if raster != None:
data = calculations.signal_profile_to_raster(numpy.array(peaks), raster, float(noise), shape)
else:
data = calculations.signal_profile(numpy.array(peaks), int(points), float(noise), shape)
# make baseline
baseline = []
for peak in peaklist:
if not baseline or baseline[-1][0] != peak.mz:
baseline.append([peak.mz, -peak.base])
# apply baseline
data = mod_signal.subbase(data, numpy.array(baseline))
return data
# ----
def matchpattern(signal, pattern, pickingHeight=0.75, baseline=None):
"""Compare signal with given isotopic pattern.
signal (numpy array) - signal data points
pattern (list of [mz,intens]) - theoretical pattern to compare
pickingHeight (float) - centroiding height
baseline (numpy array) - signal baseline
"""
# check signal type
if not isinstance(signal, numpy.ndarray):
raise TypeError, "Signal must be NumPy array!"
# check baseline type
if baseline != None and not isinstance(baseline, numpy.ndarray):
raise TypeError, "Baseline must be NumPy array!"
# check signal data
if len(signal) == 0:
return None
# get signal intensites for isotopes
peaklist = []
for isotope in pattern:
peak = mod_peakpicking.labelpeak(
signal = signal,
mz = isotope[0],
pickingHeight = pickingHeight,
baseline = baseline
)
if peak:
peaklist.append(peak.intensity)
else:
peaklist.append(0.0)
# normalize peaklist
basepeak = max(peaklist)
if basepeak:
peaklist = [p/basepeak for p in peaklist]
else:
return None
# get rms
rms = 0
for x, isotope in enumerate(pattern):
rms += (isotope[1] - peaklist[x])**2
if len(pattern) > 1:
rms = math.sqrt(rms/(len(pattern)-1))
return rms
# ----
def _consolidate(isotopes, window):
"""Group peaks within specified window.
isotopes: (list of [mass, abundance]) isotopes list
window: (float) grouping window
"""
if isinstance(isotopes, numpy.ndarray):
isotopes = isotopes.tolist()
isotopes.sort()
f = (window/1.66)*(window/1.66)
buff = []
buff.append(isotopes[0])
for current in isotopes[1:]:
previous = buff[-1]
if (previous[0] + window) >= current[0]:
mass = (previous[0]*previous[1] + current[0]*current[1]) / (previous[1] + current[1])
#ab1 = previous[1] * math.exp( - ((previous[0]-mass)*(previous[0]-mass)) / f )
#ab2 = current[1] * math.exp( - ((current[0]-mass)*(current[0]-mass)) / f )
#buff[-1] = [mass, ab1+ab2]
buff[-1] = [mass, previous[1] + current[1]]
else:
buff.append(current)
return buff
# ----
def _normalize(data):
"""Normalize data."""
# get maximum Y
maximum = data[0][1]
for item in data:
if item[1] > maximum:
maximum = item[1]
# normalize data data
for x in range(len(data)):
data[x][1] /= maximum
return data
# ----
|