xmotif.h

/**************************************************************************
 * Copyright (c) 2005-2011 T. M. Murali                                   *
 * Copyright (c) 2008-2011 Naveed Massjouni                               *
 * Copyright (c) 2004 Greg Grothaus                                       *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm 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.                                    *
 *                                                                        *
 * Biorithm 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 Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#include<string>
#include<iostream>
#include<fstream>
#include<vector>
#include<map>
#include<sstream>
#include<cstdlib>
#include<ctime>
#include<set>
#include<cmath>
#include<cfloat>
#include<algorithm>


// just for defining MyNT.
#include "enrichment.h"
// for ltqnorm().
#include "ltqnorm.h"


#include "motifDataset.h"
#include "itemset.h"

using namespace std;

#ifndef XMOTIFDATASETGREG
#define XMOTIFDATASETGREG

enum XmotifDiscriminantMethod { DISCRIMINANT_ALL_SAMPLES = 1, DISCRIMINANT_ONE_CLASS_SAMPLES,
                                //ALL_SAMPLES_WEIGHTED, 
                                DISCRIMINANT_ONE_CLASS_SAMPLES_WEIGHTED };

enum XmotifSampleExpansionMethod { ADD_SAMPLE_MATCH_GENE_RANGE = 1, ADD_SAMPLE_MATCH_GENE_RANGE_ONE_CLASS,
                             ADD_SAMPLE_SEARCH, ADD_SAMPLE_SEARCH_ONE_CLASS,
                                   ADD_SAMPLE_SEARCH_ONLY_ADD, ADD_SAMPLE_SEARCH_ONLY_ADD_ONE_CLASS,
ADD_SAMPLE_NONE};

enum XmotifGeneSelectionMethod { GENE_SELECTION_MAX = 1, GENE_SELECTION_STOUFFER,
                                 // used when i already know the genes
                                 // (e.g., when i read pre-computed
                                 // xMotifs from a file) and all i
                                 // want to do is compute stats for
                                 // the genes. THIS OPTION SHOULD
                                 // NEVER BE SET FROM THE COMMAND
                                 // LINE.
                                 GENE_SELECTION_TRIVIAL};
                             
enum XmotifPValueMethod {PVALUE_RANK, PVALUE_TTEST};

struct XMotifMethods
{
private:
  XmotifDiscriminantMethod _discSelectionMethod;
  XmotifGeneSelectionMethod _geneSelectionMethod;
  XmotifSampleExpansionMethod _sampleExpansionMethod;
  XmotifPValueMethod _pValueMethod;
public:
  XMotifMethods(XmotifDiscriminantMethod _dsm, XmotifGeneSelectionMethod _gsm,
                XmotifSampleExpansionMethod _asm, XmotifPValueMethod _pvm)
      : _discSelectionMethod(_dsm), _geneSelectionMethod(_gsm),
            _sampleExpansionMethod(_asm), _pValueMethod(_pvm)
    {}
  
  virtual ~XMotifMethods()
    {}

  XmotifDiscriminantMethod getDiscriminantSelectionMethod() const
    {
      return(_discSelectionMethod);
    }
  XmotifSampleExpansionMethod getSampleExpansionMethod() const
    {
      return(_sampleExpansionMethod);
    }
  XmotifGeneSelectionMethod getGeneSelectionMethod() const
    {
      return(_geneSelectionMethod);
    }
  XmotifPValueMethod getPValueMethod() const
    {
      return(_pValueMethod);
    }
};

  

class xmotif
{
private:
  XMotifMethods methods;

  struct _ProbeInterval
  {
  private:
    int probeIndex;
    double low;
    double high;
    double pval;
    double zscore;
    
  public:
    _ProbeInterval(int i, double l, double h, double p)
        : probeIndex(i), low(l), high(h), pval(p)
      {
        zscore = ltqnorm(1 - p);
      }

    double getHigh() const
      {
        return(high);
      }
    double getLow() const
      {
        return(low);
      }
    int getProbeIndex() const
      {
        return(probeIndex);
      }
    
    double getPValue() const
      {
        return(pval);
      }

    double getZScore() const
      {
        return(zscore);
      }
    
    
  };
  
  struct compareProbeIntervals
    : public binary_function< const xmotif::_ProbeInterval&, const xmotif::_ProbeInterval&, bool> 
  {
    bool operator()(const _ProbeInterval& __x, const _ProbeInterval& __y) const
      {
        return(__x.getZScore() > __y.getZScore());
      }
  };

        double pCutoff;

        vector<vector<int> > prevVal, nextVal, prevCountVal, nextCountVal;
        vector<int> numData;    //number of data points for each gene
        
        vector<vector<vector<double> > > countTables;
        double countPVal(int,int,int);
        double countPVal2(int,int,int);

        // Calculates p value using t test method.
        double pValueTTest(const set<int> &discriminatingSet,
                const set<int> &otherSet, int probeIndex);

        void calculateCountTables();
        void calculateFactTables();
        void calculateWeightsTables();
        vector<double> logFactTable;
        double logProbeSize;

  // Use gsm to decide which method to use for computing matching probes.
        ItemsetRange computeMatchingProbes(const set<int> &discriminatingSet,
                const set<int> &otherSet);


  ItemsetRange computeMatchingProbesKS(set<int>& D);

  // Compute the probes that match a discriminant by computing for
  // each probe the p-value of the interval spanned in that probe by
  // the samples in the discriminant and keeping only probes where the
  // p-value is less than a threshold. 
  ItemsetRange computeMatchingProbesCount(const set<int> &discriminatingSet,
    const set<int> &otherSet);

  // Compute the probes that match a discriminant by computing for
  // each probe by computing per-probe p-values as in
  // computeMatchingProbesCount and then using the subset of probes
  // that maximise Stouffer's Z.
  ItemsetRange computeMatchingProbesStouffersZ(
        const set<int>& discriminatingSet, const set<int>& otherSet);

  // Given a discriminant and the probes that match that discriminant,
  // simply compute per-probe p-values. This method is useful for
  // filling up pre-computed xMofits that you read from a file, when
  // these xMotifs only have row and column information.
  ItemsetRange computeMatchingProbesTrivial(set<int>& D, set< int > &probes);

  bool relaxRange;


  /* compute the range of probe j in the discriminant D. return false
     if the discriminant is empty or has other bad properties. return
     the indices and actual endpoints of the range in minIndex,
     maxIndex, minExpression, and maxExpression. return the number of
     points in the discriminant that are in the range in numPoints and
     the total number of points in the dataset that are within the
     range in totalCount. */
  bool computeRangeInDiscriminant(int j, const set<int> &D,
                                  int &minIndex, int &maxIndex,
                                  double &minExpression, double &maxExpression,
                                  int &numPoints, int &totalCount);
  

        set<int> randomSubset(set<int>&,int);
        set<int> randomSubset(vector<int>, int);
        set<int> randomSubset(vector<int>);
        set<int> randomSubset(vector<int>, vector<double>, double);
        set<int> randomSubset(int,int);
        set<int> randomSubset(int);
        int randInt(int,int);
        double randDouble(double,double);

        // Returns the discriminating set in discriminatingSet.
        // otherSet will contain indexes of all samples in all classes
        // except chosenClass.
        void getDiscriminatingSet(int discriminantMethod, int &chosenClass,
                set<int> &discriminatingSet, set<int> &otherSet);

        map<set<int>,bool> alreadyFoundDiscriminant;
        vector<vector<double> > classWeights,nextWeights;
        vector<double> normalize(vector<double>);
        vector<double> updateWeights(vector<double>,set<int>);
        
        bool expandSamples(set<int> &discriminatingSet, set<int> &otherSet,
                ItemsetRange &, int);
        bool expandSamplesMatchRange(set<int> &, ItemsetRange &, int, bool);
        bool expandSamplesGradientDescent(set<int> &discriminatingSet,
                set<int> &otherSet, ItemsetRange &, int, bool);
        bool expandSamplesGradientDescentKeepOriginal(set<int> &, set<int> &,
                ItemsetRange &, int, bool);
        itemset buildItemset(set<int>&, ItemsetRange&, int);

        motifDataset* myData;
        set<itemset> motifs;
  // true if i read pre-computed motifs from a file.
  bool _motifsRead;
  
public:
        xmotif(motifDataset* in, const XMotifMethods& methods);

  void changePval(double);

  void computeMotifs(int ns, void (*processFunc)(const itemset&));


  bool getMotifsRead() const
    {
      return(_motifsRead);
    }
  
  void printMotifs(ostream* out);

  void readMotifs(string infile);
  
  set<itemset> getMotifs();
};

#endif