/****************************************************************************
* Copyright (C) 2008 Peter Mortensen and Matthias Mann *
* This file is part of MSQuant. *
* *
* MSQuant is distributed under the terms of *
* the GNU General Public License. See src/COPYING.TXT or *
* <http://www.gnu.org/licenses/gpl.txt> for details. *
* *
* MSQuant 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 2 of the *
* License, or (at your option) any later version. *
* *
* MSQuant 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. *
* *
* You should have received a copy of the GNU General *
* Public License along with MSQuant; if not, write to *
* the Free Software Foundation, Inc., 59 Temple *
* Place, Suite 330, Boston, MA 02111-1307 USA *
* *
* Purpose: Calculates isotope distributions for peptides and others, *
* according to the algorithm given in:. *
* *
* [42]H. Kubinyi *
* Calculation of isotope distributions in mass spectrometry. *
* A trivial solution for a non-trivial problem, Anal. Chim. Acta *
* *
****************************************************************************/
/****************************************************************************
* CEBI *
* Software Development Group *
* Peter Mortensen *
* E-mail: NUKESPAMMERSdrmortensen@get2netZZZZZZ.dk *
* WWW: http://www.cebi.sdu.dk/ *
* *
* Program for post-processing of result from search in mass *
* spectrometric data. *
* *
* FILENAME: isotopeDistribution.cs *
* TYPE: CSHARP *
* *
* CREATED: PM 2008-01-24 Vrs 1.0. *
* UPDATED: PM 2008-xx-xx *
* *
* *
****************************************************************************/
//Notes:
// 1. File theoisotope.vb from DTAsuperCharge contains the same
// kind of computation, but for averagine.
//
// 2. Other general from DTAsuperCharge (helper classes for finding
// the invalid masses for ion series, e.g. b/y/<ECD>):
//
// a. PILperiodicTable.vb
//
// b. PILmolecule.vb. It seems to be mostly for progressively
// adding amino acids (for ion series) and (efficiently)
// compute the mass. E.g. newMoleculeByAddition(), used
// from crawlCombinations()/theoPeptideMasses.vb.
//
// "... peptides are built up by using newMoleculeByAddition() for
// each residue ..."
//
// c. In peakReductor.vb amino acids are specified by composition, e.g.:
//
// addPeptideByComp(someTheoMasses, "M", 5, 9, 1, 1, 1) 'M
//
// addPeptideByComp() uses Class theoPeptideMasses.
//
// addPeptideByComp() is internal to class peakReductor/
// peakReductor.vb and uses addAmass() of theoPeptideMasses and
// creates an instance of PILmolecule to pass to addAmass().
//
//
// Conclusion:
// We can use PILmolecule and newMoleculeByAddition for peptides.
//
// To do:
// 1. Convert the DTASC PILmolecule to C#.
// 2. Write some unit tests.
// 3. Let PILpeptide.cs have knowlegde of composition for amino
// acids. Use PILmolecule.
// 4. Extend PILpeptide and PILmolecule to return
// composition (to use in QuantWindow.vb for passing
// to findIsotopeDistribution, class isotopeDistribution ).
//
//
// d.
using System;
using System.Collections.Generic;
using System.Text;
// Use this?
/****************************************************************************
* <placeholder for header> *
****************************************************************************/
namespace massSpectrometryBase
{
/****************************************************************************
* <placeholder for header> *
****************************************************************************/
public class isotopeDistribution
{
//No instance for now. But for efficiency we may want to keep an
//instance around.
//Note: instead of space for 4 fixed and nPeak we could
//simply have a list of isotope abundances...
struct atomIsotopeStruct
{
public List<double> abundance3;
//Convenience initialisation/"constructor".
public atomIsotopeStruct(
double aIso1,
double aIso2,
double aIso3,
double aIso4
)
{
abundance3 = new List<double>(4);
abundance3.Add(aIso1);
abundance3.Add(aIso2);
abundance3.Add(aIso3);
abundance3.Add(aIso4);
}
} //struct atomIsotopeStruct
/****************************************************************************
* <placeholder for header> *
****************************************************************************/
public static List<double> findIsotopeDistribution(
int aCarbons,
int aHydrogens,
int aNitrogens,
int aOxygens,
int aSulfurs)
{
List<double> toReturn = new List<double>();
int Natom = 5; //Number of atoms in the list.
double prec = 0.000001; //Pruning threshold factor.
List<int> composition = new List<int>(Natom);
//List<double> composition = new List<double>();
//List<double> composition = new List<double>( 5);
composition.Add( aCarbons);
composition.Add( aHydrogens);
composition.Add( aNitrogens);
composition.Add( aOxygens);
composition.Add( aSulfurs);
//nPeak: number of heaviest isotope.
List<int> nPeak = new List<int>(Natom);
//List<int> nPeak = new List<int>(5);
//List<int> nPeak = new List<int>();
nPeak.Add( 2);
nPeak.Add( 2);
nPeak.Add( 2);
nPeak.Add( 3);
nPeak.Add( 4);
atomIsotopeStruct carbon =
new atomIsotopeStruct(98.900, 1.100, 0.000, 0.000);
atomIsotopeStruct hydrogen =
new atomIsotopeStruct(99.985, 0.015, 0.000, 0.000);
atomIsotopeStruct nitrogen =
new atomIsotopeStruct(99.630, 0.370, 0.000, 0.000);
atomIsotopeStruct oxygen =
new atomIsotopeStruct(99.760, 0.040, 0.200, 0.000);
atomIsotopeStruct sulfur =
new atomIsotopeStruct(95.020, 0.750, 4.210, 0.020);
List<atomIsotopeStruct> abund = new List<atomIsotopeStruct>(Natom);
abund.Add( carbon);
abund.Add( hydrogen);
abund.Add( nitrogen);
abund.Add( oxygen);
abund.Add( sulfur);
//This is where the result will be...
List<double> cPattern = new List<double>(1);
cPattern.Add( 1.0000); //Calculated pattern.
int p = 0; //Initialise = low
int q = 0; //Initialise = high
for ( int atom = 0; atom < Natom; atom++)
{
int numAtoms = composition[atom]; // composition(atom): number of atoms in the molecule.
for (int i = 0; i < numAtoms; i++)
{
//Used after the loop below.
List<double> dummy = new List<double>(100);
//For now: straight conversion. However it is not very efficient.
for (int j = 0; j < 100; j++)
{
dummy.Add( 0.0);
}
for (int k = p; k <= q; k++)
{
int len2 = nPeak[atom]; //nPeak: number of heaviest isotope.
for (int atomIsotope = 0; atomIsotope < len2; atomIsotope++)
{
//Future: reduce to: k +atomIsotope
int index1 = (k + 1) + (atomIsotope + 1) - 1;
index1--; //Because we use zero based indexes.
double ab = abund[atom].abundance3[atomIsotope];
//abund: abundancies of isotopes of the atom type atom.
dummy[index1] += cPattern[k] * ab;
} //for atomIsotope
} //for k
q += nPeak[atom] - 1; //q: new number of heaviest isotope of the pattern.
double max = 0.0;
//Find the largest value of dummy(k).
for (int k = p; k <= q; k++)
{
if ( dummy[k] > max)
{
max = dummy[k];
}
} //for k
//Normalize dummy(k) to 1, so max is 1.00.
for (int k = p; k <= q; k++)
{
dummy[k] /= max;
} //for k
//Eliminate small peaks at left side.
for (int k = p; k <= q; k++)
{
if ( dummy[k] > prec)
{
p = k;
k = q; //Loop exit...
}
else
{
int peter2 = 2; //For breakpoints.
}
} //for k
//Eliminate small peaks at right side.
for (int k = q; k >= p; k--)
{
if ( dummy[k] > prec)
{
q = k;
k = p; //Loop exit...
}
else
{
int peter2 = 2; //For breakpoints.
}
} //for k
int newLen = q + 1;
//Not very efficient as we create a new one
//every time. We could reuse the old, adding new
//elements at the end if needed.
cPattern = new List<double>(newLen); //Erase old isotope peak pattern.
for (int j = 0; j < newLen; j++)
{
cPattern.Add( 0.0);
}
//Create new isotope pattern.
for (int k = p; k <= q; k++)
{
cPattern[k] = dummy[k];
} //for k
} //for i, next atom of the same atom type
} //for atom, next atom type
toReturn = cPattern;
return toReturn;
} //findIsotopeDistribution()
} //class isotopeDistribution
} //namespace massSpectrometryBase
Generated by script codePublish.pl at 2009-01-05T15:20:59.