ttneutralgenes.h

Go to the documentation of this file.

 
#ifndef TTNEUTRALGENES_H
#define TTNEUTRALGENES_H
 
#include <cmath>
#include <vector>
#include "ttrait_with_map.h"
#include "types.h"
#include "lifecycleevent.h"
#include "filehandler.h"
#include "stathandler.h"
#include "datatable.h"
#include "metapop.h"
#include "binarystoragebuffer.h"
 
class TTNeutralGenesFH;
class TTNeutralGenesSH;
class TTNtrlPhenotyperFH;
class TProtoNeutralGenes;
// ----------------------------------------------------------------------------------------
 
//                                 N E U T R A L   T R A I T
 
// ----------------------------------------------------------------------------------------
class TTNeutralGenes : public TTrait 
{
private:
  TProtoNeutralGenes* _myProto;
  //parameters:
  unsigned int _allele_num;
  unsigned int _locus_num;
  unsigned int _ploidy;
  double _mut_rate;
  double _recomb_rate;
  unsigned int _mut_model;
  unsigned int _2L; //genome size (ploidy*num_locus)
  unsigned short _init_model;
  void (TTNeutralGenes::* _mutate_func_ptr) (void);
  void (TProtoNeutralGenes::* _inherit_func_ptr) (sex_t, unsigned char*, unsigned char**);
 
  //_sequence:
  unsigned char** _sequence; // ordered as [ploidy][locus]
 
  const trait_t _type;
   
public:
  
  TTNeutralGenes () 
    : _myProto(0),_allele_num(0),_locus_num(0),_ploidy(2),_mut_rate(0),_recomb_rate(0),_mut_model(0),
    _2L(0), _init_model(0), _mutate_func_ptr(0), _inherit_func_ptr(0), _sequence(0), _type(NTRL) { }
  
  TTNeutralGenes(const TTNeutralGenes& T) 
    : _myProto(T._myProto),
    _allele_num(T._allele_num), _locus_num(T._locus_num), _ploidy(2), _mut_rate(T._mut_rate),
    _recomb_rate(T._recomb_rate),_mut_model(T._mut_model), _2L(T._2L),
    _init_model(T._init_model),_mutate_func_ptr(T._mutate_func_ptr), _inherit_func_ptr(T._inherit_func_ptr), 
    _sequence(0), _type(NTRL) { }
  
  virtual ~TTNeutralGenes ();
    
  unsigned int    get_ploidy           ( )                                                 {return _ploidy;}
  unsigned int    get_locus_num        ( )                     {return _locus_num;}
  unsigned int    get_allele_num       ( )                     {return _allele_num;}
  
  void            set_proto          (TProtoNeutralGenes* proto) {_myProto = proto;}
  void            set_locus_num        (int value)             {_locus_num = value;}
  void            set_allele_num       (int value)             {_allele_num = value;}
  void            set_mut_rate         (double value)          {_mut_rate = value;}
  void            set_2L               (unsigned int val)      {_2L = val;}
  void            set_recomb_rate      (double value)          {_recomb_rate = value;}
  void            set_mut_model        (int value)             {_mut_model = value;}
  void            set_init_model       (unsigned short val)    {_init_model = val;}
  void            set_mut_func_ptr     (void(TTNeutralGenes::* theFunc)(void))
  {_mutate_func_ptr = theFunc;}
  void            set_inherit_func_ptr (void(TProtoNeutralGenes::* theFunc)(sex_t, unsigned char*, unsigned char**))
  {_inherit_func_ptr = theFunc;}
  void            set_allele           (unsigned int loc, unsigned int al, unsigned char val) {_sequence[al][loc] = val;}
  
 // static void     set_recomb_template           (unsigned int size);
  
  void            mutate_SSM           ( );
  void            mutate_KAM           ( );
  void            mutate_2all          ( );
  void            mutate_NULL          ( ) { }
  virtual TTNeutralGenes& operator= (const TTrait& T);
  virtual bool    operator== (const TTrait& T);
  virtual bool    operator!= (const TTrait& T);
  virtual void    init                 ( );
  virtual void    init_sequence        ( );
  virtual void    reset                ( );
  virtual void*   set_trait            (void* value)           {return NULL;}
  inline virtual void**  get_sequence  ( )  const              {return (void**)_sequence;}
  virtual double  get_allele_value     (int loc, int all)  
    {return ( !(loc<(int)_locus_num) || !(all<(int)_ploidy) ? 0 : (double)_sequence[all][loc]);}
  virtual void    set_allele_value (unsigned int locus, unsigned int allele, double value)
    {assert(locus < _locus_num && allele < 2); _sequence[allele][locus] = (unsigned char)value;}
  virtual void    set_sequence         (void** seq);
  inline virtual trait_t get_type      ( )  const              {return _type;}
  virtual void    set_value            ( )                     { }
  virtual void*   getValue             ( )      const              {return 0;}
  virtual void    inherit              (TTrait* mother, TTrait* father);
  virtual void    mutate               ( )                      {(this->*_mutate_func_ptr) ();}
  virtual void    show_up              ( );
  virtual TTNeutralGenes*  clone       ( )                      {return new TTNeutralGenes(*this);}
  
  //implements StorableComponent
  virtual void    store_data       (BinaryStorageBuffer* saver)  
  {
    for(unsigned int i = 0; i < _locus_num; ++i)
      for(unsigned int j = 0; j < _ploidy; ++j)
        saver->store(&_sequence[j][i], 1);
  }
  virtual bool    retrieve_data    (BinaryStorageBuffer* reader)
  { 
    for(unsigned int i = 0; i < _locus_num; ++i)
      for(unsigned int j = 0; j < _ploidy; ++j)
        reader->read(&_sequence[j][i], 1);
    return true;
  }
};  
// ----------------------------------------------------------------------------------------
 
//                     N E U T R A L   T R A I T   P R O T O T Y P E
 
// ----------------------------------------------------------------------------------------
class TProtoNeutralGenes : public TTProtoWithMap {
 
private:
  unsigned int _allele_num;
  unsigned int _locus_num;
  unsigned int _ploidy;
  double _mut_rate;
  int _mut_model;
  unsigned short _init_model;
  double _recomb_rate;
  void (TTNeutralGenes::* _mutate_func_ptr) (void);
  void (TProtoNeutralGenes::* _inherit_func_ptr) (sex_t, unsigned char*, unsigned char**);
  
  vector< TTNeutralGenesFH* > _writers;
  TTNeutralGenesSH* _stats;
  const trait_t _type;
  
public:
  
  TProtoNeutralGenes ( );
  
  TProtoNeutralGenes(const TProtoNeutralGenes& T);
  
  virtual ~TProtoNeutralGenes ( );
  
  unsigned int     get_ploidy           ( )                                        {return _ploidy;}
  unsigned int     get_locus_num        ( )            {return _locus_num;}
  unsigned int     get_allele_num       ( )            {return _allele_num;}
  
  TTNeutralGenesSH* get_stater () {return _stats;}
 
  void    inherit_low                  (sex_t SEX, unsigned char* seq, unsigned char** parent);
  void    inherit_free                 (sex_t SEX, unsigned char* seq, unsigned char** parent);
 
  //implementation of TraitPrototype:
  virtual void                     init  (){}
  
  virtual void                     reset (){TTProtoWithMap::reset();}
  
  virtual TTNeutralGenes*          hatch ();
  
  virtual TProtoNeutralGenes*      clone () {return new TProtoNeutralGenes(*this);}
  
  virtual   trait_t                 get_type () const {return _type;}
  //implementation of SimComponent:
  virtual bool setParameters();
  
  virtual void loadFileServices ( FileServices* loader );
  
  virtual void loadStatServices ( StatServices* loader );
  
  virtual bool resetParameterFromSource (std::string param, SimComponent* cmpt); // {return false;}
 
  //implementation of StorableComponent:
  virtual void store_data    ( BinaryStorageBuffer* saver ) {saver->store(&_locus_num,sizeof(int));}
  
  virtual bool retrieve_data ( BinaryStorageBuffer* reader );
};
 
// ----------------------------------------------------------------------------------------
 
//                     N E U T R A L   T R A I T   F I L E   H A N D L E R
 
// ----------------------------------------------------------------------------------------
 
class TTNeutralGenesFH: public TraitFileHandler<TProtoNeutralGenes> {
  
  string _output_option;
  void (TTNeutralGenesFH::* write_fct) ();
  
public:
        
  TTNeutralGenesFH ( TProtoNeutralGenes* TP )
    : TraitFileHandler<TProtoNeutralGenes> (TP, ".txt"), write_fct(0)
  { }
  
  virtual ~TTNeutralGenesFH ( ) { }
 
  virtual void FHwrite  ();
  
  virtual void FHread (string& filename);
 
  void write_TAB();
  void write_patch_TAB(Patch* patch, sex_t SEX, age_idx AGE, ofstream &FH);
  void write_PLINK ();
  void print_PLINK_PED(ofstream& FH, age_idx Ax, Patch *patch);
  void write_PLINK_BED (ofstream &BED);
  void write_GENEPOP();
  void write_patch_GENEPOP(Patch* patch, sex_t SEX, age_idx AGE, ofstream &FH, unsigned int digits);
  void write_FSTAT();
  void write_patch_FSTAT(Patch* patch, sex_t SEX, age_idx AGE, ofstream &FH, unsigned int digits);
  void write_Fst_i();
//  void write_freq();
  void write_varcompWC();
  
  void setOutputOption(string opt);
  void set_write_fct( void (TTNeutralGenesFH::* fct_ptr) () ) {write_fct = fct_ptr;}
  
};
 
// ----------------------------------------------------------------------------------------
 
//                     N E U T R A L   T R A I T   S T A T   H A N D L E R
 
// ----------------------------------------------------------------------------------------
 
class TTNeutralGenesSH: public TraitStatHandler<TProtoNeutralGenes, TTNeutralGenesSH> {
  
  DataTable< unsigned int > _alleleCountTable;
  DataTable< double > _alleleFreqTable;
  DataTable< double > _heteroTable;
  TMatrix _globalAlleleFreq;
  unsigned int _table_set_gen, _table_set_age, _table_set_repl;
  
  double Theta_FF, Theta_MM, Theta_FM;
  double _mean_theta, _mean_alpha;
  TMatrix *_coa_matrix;
  
  double _sib_prop[4];
  double _sib_coa[4];
  
  double _ho, _hs, _ht, _hsnei, _htnei, _nb_all_local, _nb_all_global,
    _fst, _fis, _fit, _fix_loc_local, _fix_loc_global;
  double _fst_WH;
  double _fst_WC, _fis_WC, _fit_WC;
  double *_fst_WC_loc, *_fis_WC_loc, *_fit_WC_loc;
 
  double _fst_W1, _fst_W2;
  TMatrix *_fst_matrix;
  
  //Nei's genetic distance:
  TMatrix *_D;
  double _meanD;
  
public:
  
        TTNeutralGenesSH (TProtoNeutralGenes* TP) 
    : TraitStatHandler<TProtoNeutralGenes, TTNeutralGenesSH>(TP), _table_set_gen(0), _table_set_age(0),
    _table_set_repl(0), _coa_matrix(0), _fst_WC_loc(0), _fis_WC_loc(0), _fit_WC_loc(0), _fst_matrix(0),
        _D(0)
        { }
  
  virtual ~TTNeutralGenesSH ( ) 
  {
    if(_coa_matrix != NULL) delete _coa_matrix;
    if(_fst_matrix != NULL) delete _fst_matrix;
    if(_fst_WC_loc) delete[]_fst_WC_loc;
    if(_fis_WC_loc) delete[]_fis_WC_loc;
    if(_fit_WC_loc) delete[]_fit_WC_loc;
    if(_D != NULL) delete _D;
  }
  
  virtual void init ( ) ;
  
  virtual bool setStatRecorders (std::string& token);
  
  void    setFreqRecorders (age_t AGE);
  void    setFreqRecordersPerPatch (age_t AGE);
  void    setFstatRecorders (age_t AGE);
  void    setFstat2Recorders (age_t AGE);
  void    setFstatWCRecorders (age_t AGE);
  void    setCoaMatrixRecorders (age_t AGE, unsigned char dim);
  void    setFstMatrixRecorders (age_t AGE, unsigned char dim);
  void    setNeiGeneticDistanceRecorders (age_t AGE, bool pairwise);
  void    setDxyRecorders (age_t AGE, bool patchwise);
  
  void setAdultAlleleFreq     () {setAlleleTables(ADULTS);}
  void setOffspringAlleleFreq () {setAlleleTables(OFFSPRG);}\
  void setHeterozygosity      (age_t AGE);
  void setAdultHeterozygosity     () {setHeterozygosity(ADULTS);}
  void setOffspringHeterozygosity () {setHeterozygosity(OFFSPRG);}
  
  double getGlobalAlleleFreq (unsigned int loc, unsigned int all) {
    return _globalAlleleFreq.get(loc, all);
  }
  
  double getHeterozygosity (unsigned int loc) {
    double het = 0;
    for(unsigned int i = 0; i < _heteroTable.getNumGroups(); ++i )
      het += _heteroTable.get(i, loc, 0);
    return het/_heteroTable.getNumGroups(); //mean per patch heterozygosity
  }
  
  void   setAlleleTables          (age_t AGE);
  void   setHeteroTable           (age_t AGE);
  void   allocateTables           (unsigned int loci, unsigned int all);
  
  DataTable<double>* getAlleleFreqTable () {return &_alleleFreqTable;}
  
  DataTable< unsigned int >* getAlleleCountTable () {return &_alleleCountTable;}
  
  DataTable<double>*   getHeteroTable      () {return &_heteroTable;}
  
  TMatrix* getGlobalFreqs () {return &_globalAlleleFreq;}
  
  void   setFstMatrix             (age_t AGE, unsigned char dim);
  void   setAdultsFstMatrix       ()               {setFstMatrix(ADULTS, 3);}
  void   setAdultsFstWithin       ()               {setFstMatrix(ADULTS, 1);}
  void   setAdultsFstBetween      ()               {setFstMatrix(ADULTS, 2);}
  void   setOffsprgFstMatrix      ()               {setFstMatrix(OFFSPRG, 3);}
  void   setOffsprgFstWithin      ()               {setFstMatrix(OFFSPRG, 1);}
  void   setOffsprgFstBetween     ()               {setFstMatrix(OFFSPRG, 2);}
  double getWeightedFst           ()               {return _fst_WH;}
  double getFst_ij                (unsigned int i)
  {
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_fst_matrix->getNbCols()) + 1 );
    return _fst_matrix->get(i/scale, i%scale);
  }
  void   setFst_li(unsigned int N, unsigned int L, double **array);
//  /**Computes raw Fst following the original definition (=var(p)/p_bar(1 - p_bar)).*/
//  double getFstWright(unsigned int i) {if(i == 1) return _fst_W1; else return _fst_W2;}
  
  void   setFstat                 (age_t AGE);
  void   setOffsprgFstat          ()               {setFstat(OFFSPRG);}
  void   setAdultsFstat           ()               {setFstat(ADULTS);}
  double setHo                    (age_idx age_pos);
  double setHs                    (age_idx age_pos);
  double setHt                    (age_idx age_pos);
  double getHsnei                 ()               {return _hsnei;}
  double getHtnei                 ()               {return _htnei;}
  double getHo                    ()               {return _ho;}
  double getHs                    ()               {return _hs;}
  double getHt                    ()               {return _ht;}
  double getFst                   ()               {return _fst;}
  double getFis                   ()               {return _fis;}
  double getFit                   ()               {return _fit;}
  void   setFstat2                 (age_t AGE);
  void   setOffsprgFstat2          ()               {setFstat2(OFFSPRG);}
  void   setAdultsFstat2           ()               {setFstat2(ADULTS);}
  deque<double> setHo2                    (age_idx age_pos);
  deque<double> setHs2                    (age_idx age_pos);
  deque<double> setHt2                    (age_idx age_pos);
  
  void   setFstatWeirCockerham    (age_t AGE);
  void   setFstatWeirCockerham_MS   (age_t AGE);
  void   setOffspringFstatWeirCockerham()          {setFstatWeirCockerham(OFFSPRG);}
  void   setAdultsFstatWeirCockerham()             {setFstatWeirCockerham(ADULTS);}
  double getFstWC                 ()               {return _fst_WC;}
  double getFisWC                 ()               {return _fis_WC;}
  double getFitWC                 ()               {return _fit_WC;}
  void   setLociDivCounter        (age_t AGE);
  double getNbAllLocal            ()               {return _nb_all_local;}
  double getNbAllGlobal           ()               {return _nb_all_global;}
  double getFixLocLocal           ()               {return _fix_loc_local;}
  double getFixLocGlobal          ()               {return _fix_loc_global;}
 
  double Coancestry               (void** ind1, void** ind2, unsigned int nb_locus);
  void   setCoaMatrix             (age_idx age_pos, unsigned char dim);
  void   setAdultsCoaMatrix       ()               {setCoaMatrix(ADLTx, 3);}
  void   setOffsprgCoaMatrix      ()               {setCoaMatrix(OFFSx, 3);}
  void   setAdultsCoaWithin       ()               {setCoaMatrix(ADLTx, 1);}
  void   setOffsprgCoaWithin      ()               {setCoaMatrix(OFFSx, 1);}
  void   setAdultsCoaBetween      ()               {setCoaMatrix(ADLTx, 2);}
  void   setOffsprgCoaBetween     ()               {setCoaMatrix(OFFSx, 2);}
  void   setAdults_Theta          ();
    
  double getCoa                   (unsigned int i)
  {
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_coa_matrix->getNbCols()) + 1 );
    return _coa_matrix->get(i/scale, i%scale);
  }
  double getMeanTheta             ()           {return _mean_theta;}
  double getMeanAlpha             ()           {return _mean_alpha;}
  double getTheta_FF              ()                   {return Theta_FF;}
  double getTheta_MM              ()                   {return Theta_MM;}
  double getTheta_FM              ()                   {return Theta_FM;}
  void   setSibStats              ();
  void   setSibCoa                (Individual *I1, Individual *I2);
  double getSibProportions        (unsigned int i) {return _sib_prop[i];}
  double getSibCoaMeans           (unsigned int i) {return _sib_coa[i];}
  
  void    setAdltNeiGeneticDistance    ( ) {setNeiGeneticDistance(ADULTS);}
  void    setOffsprgNeiGeneticDistance ( ) {setNeiGeneticDistance(OFFSPRG);}
  void    setNeiGeneticDistance        (age_t AGE);
  double  getNeiGeneticDistance        (unsigned int i)  
  { 
    unsigned int scale = (unsigned int)pow( 10.0, (int)log10((float)_D->getNbCols()) + 1 );
    return _D->get(i/scale,i%scale);
  }
  double  getMeanNeiGeneticDistance    ( ) {return _meanD;}
  
  double getDxyOffspringPerPatch  (unsigned int patch1, unsigned patch2) {return getDxyPerPatch(OFFSx, patch1, patch2);}
  double getDxyAdultPerPatch      (unsigned int patch1, unsigned patch2) {return getDxyPerPatch(ADLTx, patch1, patch2);}
  double getDxyPerPatch           (age_idx age, unsigned int patch1, unsigned patch2);
  double getDxy                   (unsigned int age_class);
};
 
// ----------------------------------------------------------------------------------------
 
//                     N E U T R A L   T R A I T   F R E Q  I N I T I A L I Z E R
 
// ----------------------------------------------------------------------------------------
 
//LCE_NtrlInit
//
class LCE_NtrlInit : public virtual LifeCycleEvent
{
 
  TMatrix _init_freq;
  unsigned int _nLoci;
 
public:
 
  LCE_NtrlInit ( );
 
  virtual ~LCE_NtrlInit ( ) { }
 
  virtual void execute ();
 
  void init_allele_freq(sex_t SEX, age_idx AGE, unsigned int size, unsigned int deme, double *values);
 
  virtual LifeCycleEvent* clone ( ) {return new LCE_NtrlInit();}
 
  virtual bool setParameters ();
 
  //SimComponent implementation:
  virtual void loadFileServices ( FileServices* loader ) {}
  virtual void loadStatServices ( StatServices* loader ) {}
  virtual bool resetParameterFromSource (std::string param, SimComponent* cmpt) {return false;}
  virtual age_t removeAgeClass ( ) {return 0;}
  virtual age_t addAgeClass ( ) {return 0;}
  virtual age_t requiredAgeClass () {return 0;}
 
};
 
#endif //TTNEUTRALGENES_H