# -------------------------------------------------------------------------
# Copyright (C) 2005-2013 Martin Strohalm
# 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
# ----