overlaps.cpp

Go to the documentation of this file.

/* overlaps.cpp :  Mass spectrum Overlap finding algorithms
 * Copyright (C) 2020  Daniel Haley
 * 
 * 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.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include "atomprobe/isotopes/overlaps.h"
#include "atomprobe/helper/misc.h"
#include "../helper/helpFuncs.h"

#include <vector>
#include <map>
#include <utility>
#include <string>
#include <algorithm>
#include <set>

#define EQ_TOLV(f,g,h) (fabs( (f) - (g)) < (h))


namespace AtomProbe
{ 

using std::vector;
using std::pair;
using std::map;
using std::set;
using std::string;

//DEBUG ONLY
using std::cerr;
using std::endl;


//Find the masses that overlap with one another 
void findOverlapsFromSpectra(const vector<vector<pair<float,float> > > &massDistributions,
            float massDelta,const set<unsigned int> &skipIDs,
            vector<pair<size_t ,size_t> > &overlapIdx,
            vector<pair<float,float> > &overlapMasses)

{
    
    //Loop through each "fingerprint" spectrum 
    for(unsigned int ui=0;ui<massDistributions.size();ui++)
    {
        //Skip any pair where we dont understand one ion
        if(skipIDs.find(ui) !=skipIDs.end())
            continue;
        //Loop through each peak
        for(unsigned int uj=0;uj<massDistributions[ui].size();uj++)
        {
            if(skipIDs.find(uj) !=skipIDs.end())
                continue;
            float m1;
            m1=massDistributions[ui][uj].first;

            //Loop over other fingerprints 
            for(unsigned int uk=ui+1;uk<massDistributions.size();uk++)
            {
                //Loop over other fingerprint's peaks
                for(unsigned int ul=0;ul<massDistributions[uk].size();ul++)
                {
                    //If within tolerance, then record
                    float m2;
                    m2=massDistributions[uk][ul].first;
                    if(EQ_TOLV(m1,m2,massDelta))
                    {
                        //Mass distributions ui and uk overlap
                        overlapIdx.push_back(std::make_pair(ui,uk));
                        overlapMasses.push_back(std::make_pair(m1,m2));
                    }
                
                }
            }

        }

    }
}

//Find all overlaps from the ions listed in the given rangefile
// using the mass delta, and assuming they occupy charge states from
// 1->maxCharge
void findOverlaps(const AbundanceData &natData,
            const RangeFile &rng, float massDelta, unsigned int maxCharge,
            vector<pair< size_t,size_t> > &overlapIdx, vector<pair<float,float> > &overlapMasses )
{
    ASSERT(maxCharge > 0);
    ASSERT(massDelta >0.0f);
    //Minimum acceptable probability before we consider a given mass peak to not be observable
    const float PEAK_MIN_PROBABILITY=10.0/1e6; //10ppm
    const float PEAK_GROUP_DIST = 0.01; // Small m/z value in Da to group isotopes, e.g. those in Mo2
    
    vector<pair<unsigned int ,vector<size_t> > > overlaps;

    //Create a list of the ions in the rangefile,
    // decomposed so that the natural abundance table
    // can look up their "fingerprint" theoretical intensities
    vector<vector<pair<float,float> > > massDistributions;
    set<unsigned int> skipIDs;

    for(unsigned int ui=0;ui<rng.getNumIons(); ui++)
    {

        //FIXME: Better error handling?
        //Obtain the fragments for this ion. 
        // abort routine if we can't do this
        vector<pair<string,size_t> > fragments;
        vector<pair<float,float> > nullDist;
        if(!rng.decomposeIonById(ui,fragments))
        {
            skipIDs.insert(ui);
            //FIXME: Hack, insert one distribution for each charge
            for(unsigned int uj=0;uj<maxCharge;uj++)
                massDistributions.push_back(nullDist);
            continue;
        }


        //Convert to element idx (in table)/frequency pairing
        vector<pair<size_t,unsigned int> > symbolIdxFreq;
        natData.getSymbolIndices(fragments,symbolIdxFreq);

        //check that each fragment is a real element
        bool skip;
        skip=false;
        for(unsigned int uj=0;uj<symbolIdxFreq.size();uj++)
        {
            if(symbolIdxFreq[uj].first == (size_t)-1)
            {
                skip=true;
                break;
            }

        }

        //If any fragment is not a real element, skip this
        if(skip)
        {
            skipIDs.insert(ui);
            //FIXME: Hack, insert one empty distribution for each charge
            for(auto uj=0u;uj<maxCharge;uj++)
                massDistributions.push_back(nullDist);
            continue;
        }

        vector<size_t> elementIdx,freq;
        elementIdx.resize(symbolIdxFreq.size());
        freq.resize(symbolIdxFreq.size());
        for(unsigned int uj=0;uj<symbolIdxFreq.size();uj++)
        {
            elementIdx[uj]=symbolIdxFreq[uj].first;
            freq[uj] = symbolIdxFreq[uj].second;
        }

        //Generate the isotope 
        vector<pair<float,float> > massDist;
        natData.generateGroupedIsotopeDist(elementIdx,freq,massDist,PEAK_GROUP_DIST);


        //Kill any pair with a probability below our threshold
        vector<bool> killVec(massDist.size(),false);
        for(unsigned int uj=0;uj<massDist.size();uj++)
            killVec[uj]=(massDist[uj].second < PEAK_MIN_PROBABILITY);


        vectorMultiErase(massDist,killVec); 

        //Sum number of electrons
        unsigned int nElectrons=0;
        for(auto uj=0;uj<elementIdx.size(); uj++)
            nElectrons+=natData.getAtomicNumber(elementIdx[uj]);
        
        //Generate one mass distribution per charge state,
        // up to the maximum atomic number (number of electrons)
        unsigned int maxMultiplicity=std::min(maxCharge,nElectrons);
            

        for(unsigned int uj=1;uj<=maxMultiplicity;uj++)
        {
            vector<pair<float,float> > tmpDist;
            tmpDist.clear();

            if(uj <=nElectrons)
            {
                tmpDist=massDist;
                for(unsigned int uk=0;uk<tmpDist.size();uk++)
                    tmpDist[uk].first/=uj;
            }
            
            massDistributions.push_back(tmpDist);
        }

    }


    //Now, we have the mass distributions, check to see if any of them
    // are within tolerance of one another
    ASSERT(massDistributions.size() == rng.getNumIons()*maxCharge);

    overlapMasses.clear();
    overlapIdx.clear();

    findOverlapsFromSpectra(massDistributions,massDelta,skipIDs, overlapIdx,overlapMasses);
}


void findOverlaps(const AbundanceData &natData,
            const vector<string> &ionNames, float massDelta, unsigned int maxCharge,
            vector<pair< size_t,size_t> > &overlapIdx, vector<pair<float,float> > &overlapMasses )
{
    RangeFile rng;
    RGBf rgb;

    //Fake a rangefile
    for(unsigned int ui=0;ui<ionNames.size();ui++)
    {
        rng.addIon(ionNames[ui],ionNames[ui],rgb);
    }

    //call find overlaps
    findOverlaps(natData, rng, massDelta,maxCharge,overlapIdx,overlapMasses);

}
#ifdef DEBUG
bool testOverlapDetect()
{
    RangeFile rng;
    RGBf c;
    c.red=c.green=c.blue=1.0f;

    unsigned int ionIdTi,ionIdO;
    ionIdTi =rng.addIon("Ti","Ti",c);
    ionIdO =rng.addIon("O","O",c);
    rng.addRange(23.5,24.5,ionIdTi);
    rng.addRange(13.5,14.5,ionIdO);
    AbundanceData natData;
    if(natData.open("naturalAbundance.xml"))
        return false;

    vector<pair< size_t,size_t> > overlapIdx;
    vector<pair<float,float> > overlapMasses;

    //Ti overlaps with O
    const unsigned int MAX_CHARGESTATE=3;
    findOverlaps(natData,rng, 0.1,MAX_CHARGESTATE,overlapIdx,overlapMasses);

    TEST(overlapIdx.size() >=1,"Overlaps exist");

    return true;
    
}
#endif
}