python-interface/mspy/mod_pattern.py at master · hide5stm/python-interface · GitHub

Repository navigation

History

351 lines (277 loc) · 11.1 KB

Raw

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

# -------------------------------------------------------------------------

# 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

# ----

python-interface/mspy/mod_pattern.py at master · hide5stm/python-interface · GitHub

Repository navigation

Expand file tree

Breadcrumbs

mod_pattern.py

Latest commit

History

Breadcrumbs

mod_pattern.py

File metadata and controls