search

eugenemel / maven

Maven GUI: Metabolomics Analysis and Visualization Engine
https://github.com/eugenemel/maven/releases
GNU General Public License v3.0
19 stars 9 forks source link

Review correct for natural 13C abundance in isotopes computation #382

Closed PMSeitzer closed 3 years ago

PMSeitzer commented 3 years ago

This work also has applicability to #380

Ensure that the definition makes sense here.

May want to generalize this to reflect natural abundance of all isotopes

PMSeitzer commented 3 years ago

Note that as of completion of #371, this is now in PeakGroup::pullIsotopes().

PMSeitzer commented 3 years ago

See Also #387

PMSeitzer commented 3 years ago

Implementation is found in MassCalculator::computeIsotopes(), in mzMassCalculator.cpp.

The critical computations are all in this snippet:

    const double abC12 = 0.9893;
    const double abC13 = 0.0107;
    const double abN14 = 0.9963620;
    const double abN15=  0.0036420;
    const double abS32=  0.949926;
    const double abS34=  0.042524;
    const double abH  =  0.999885;
    const double abH2  =  0.00011570;
...
    const double D_Delta = 2.01410177811-1.00782503224;
    const double C_Delta = 13.00335483521-12.0;
    const double N_Delta = 15.0001088989-14.00307400446;
    const double S_Delta = 33.96786701-31.9720711744;
...
        isotopes[i].abundance=
                 mzUtils::nchoosek(CatomCount,c)*pow(abC12,CatomCount-c)*pow(abC13,c)
               * mzUtils::nchoosek(NatomCount,n)*pow(abN14,NatomCount-n)*pow(abN15,n)
               * mzUtils::nchoosek(SatomCount,s)*pow(abS32,SatomCount-s)*pow(abS34,s)
               * mzUtils::nchoosek(HatomCount,d)*pow(abH,HatomCount-d)  *pow(abH2,d);

where c is the number of 13C, n is the number of 15N, etc

PMSeitzer commented 3 years ago

It might be worth refining this computation to be more accurate (e.g., more significant digits). Similarly, will want to document this somehow.

PMSeitzer commented 3 years ago

This formula appears to be roughly accurate, as verified with a number of compounds checked against: https://www.sisweb.com/mstools/isotope.htm

The only thing that is slightly confusing is that if isotopes are not checked in the settings, they will not be included in the computation. For instance, if only 13C is checked, the above calculator gives very similar results for C10 and C100. However, C10S10 varies significantly, b/c there is a nontrivial amount of 34S that is not being accounted for.

It's worth mentioning that at this time, only C, H, S, and N are handled by the isotopes calculator.