Commit d9682638 authored by Jens Behrmann's avatar Jens Behrmann

Add model of isotope patterns

parent 5a97b1c4
function [skeleton, gPattern, estNumCAtoms] = isotopicPattern(m,mzVec,gaussianSigma,showFlag)
% Producing an isotopic pattern for a desired mass and a given mzVector
% (from the measurement). The pattern is determined by assuming human
% tissue and average numbers of amino acids and C13 compared to C12 atoms,
% with a small correction for N- and other unlikely stable isotopes.
%
% Minitest: isotopicPattern(1198.7,600:0.16:4000,0.2,true)
%
% Input:
% m - Mass value of interest for which the isotopic pattern
% will be determined (doesn't have to occur in mzVec)
% mzVec - Vector containing the mass grid from a specific
% measurement
% gaussianSigma - Sigma of gaussian function that transforms skeleton into
% gPattern (default set to
% showFlag - If set to true skeleton and gPattern are shown in a
% figure (default set to false)
%
% Output:
% skeleton - Vector containing the estimated isotopic pattern
% (first entry is relative height of monoisotopic peak,
% second entry is relative height of mon. peak + 1 Da and
% so on)
% gPattern - Using skeleton and mzVec to produce an estimate what one
% sees for a given isotopic skeleton pattern, mass grid and
% gaussianSigma when assuming a gaussian function of
% "mass inaccuracy", i.e. the result is a (normalized)
% vector of the same length as mzVec
% estNumCAtoms - Estimated number of C-atoms of for the given mass (as
% floating number)
narginchk(2,4);
if nargin < 4
showFlag = false;
if nargin < 3
gaussianSigma = 0.2;
end
end
% For the following compare the sources (tables of amino acids):
% https://www.promega.com/~/media/files/resources/technical%20references/amino%20acid%20abbreviations%20and%20molecular%20weights.pdf
% http://webext.pasteur.fr/tekaia/aafreq.html
% [and there hs - homo sapiens]
% http://www.imgt.org/IMGTeducation/Aide-memoire/_UK/aminoacids/formuleAA/
aminoAcids = [... % Proportion in homo sapiens/ #C-atoms/ molecular weight
... % for the 20 amnio acids in their common sorting
6.97, 3, 89;
5.68, 6, 174;
3.58, 4, 132;
4.68, 4, 133;
2.31, 3, 121;
4.71, 5, 146;
6.95, 5, 147;
6.70, 2, 75;
2.64, 6, 155;
4.35, 6, 131;
9.87, 6, 131;
5.72, 6, 146;
2.20, 5, 149;
3.64, 9, 165;
6.37, 5, 115;
8.25, 3, 105;
5.34, 4, 119;
1.27,11, 204;
2.63, 9, 181;
6.00, 5, 117];
weightedMean = @(x,w)(1/sum(w)*(w(:))'*x(:));
weights = aminoAcids(:,1);
mAvgAminoAcid = weightedMean(aminoAcids(:,3),weights);
avgNumCperAA = weightedMean(aminoAcids(:,2),weights);
mH2O = 18; % One peptide bond produces one H2O as a side product
% (subtracting for mass estimate)
estNumCAtoms = (m - mH2O)/(mAvgAminoAcid - mH2O)*avgNumCperAA;
thr = 0.005; % Threshold to cut off very unlikely, insignificant isotopes
pC13 = 0.0107; % Proportion of carbon 13 from natural, stable
% C12 and C13 isotope occurence
isotopeMaxNum = 20; % Determines how many Daltons to the right from the
% monoisotopic mass is being looked at for significant
% isotopes
skeleton = binopdf(0:isotopeMaxNum,round(estNumCAtoms),pC13);
skeleton(skeleton<thr) = 0;
% L1-normalizing
l1normed = @(vec)(vec/sum(abs(vec)));
skeleton = l1normed(skeleton(:));
%% Add gaussian shape
gaussFun = @(vecSkeleton, sigma, x)(...
sum(repmat((vecSkeleton(:))',length(x),1).*exp(-((repmat(x(:), ...
1,length(vecSkeleton) )-repmat(0:isotopeMaxNum, length(x), 1))./...
(sqrt(2)*sigma)).^2),2));
% Scalar constant doesn't matter because of later rescaling to |*|_1 = 1
gPattern = zeros(1,length(mzVec));
% Estimating maxPeakWidth
if max(diff(mzVec)) > 1 % Assuming linear measurement mode
maxPeakWidthInDa = 15;
else % assuming reflector mode
maxPeakWidthInDa = 1.5;
end
[~,mHelp1] = min(abs(mzVec - ( m-maxPeakWidthInDa) ));
[~,mHelp2] = min(abs(mzVec - ( m+length(skeleton)-1+maxPeakWidthInDa) ));
gPattern(mHelp1:mHelp2) = gaussFun(skeleton, gaussianSigma, ...
mzVec(mHelp1:mHelp2)-m);
% L1-normalizing
gPattern(mHelp1:mHelp2) = l1normed(gPattern(mHelp1:mHelp2));
%% Show
if showFlag
subplot(1,2,1);
stem(skeleton);
subplot(1,2,2);
plot(mzVec(mHelp1:mHelp2),gPattern(mHelp1:mHelp2));
xlim([mzVec(mHelp1) mzVec(mHelp2)]);
end
end
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment