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"
#include "Uniform.h"
 
class TTNeutralGenesFH;
class TTNeutralGenesSH;
class TTNtrlPhenotyperFH;
class TTNOhtaStats;
class TProtoNeutralGenes;
class TTNeutralGenes_byte;
// ----------------------------------------------------------------------------------------
 
//                                 N E U T R A L   T R A I T   (abstract base)
 
// ----------------------------------------------------------------------------------------
class TTNeutralGenes : public TTrait
{
protected:
  TProtoNeutralGenes* _myProto;
  const trait_t _type;
 
public:
 
  TTNeutralGenes ()
    : _myProto(0), _type(NTRL) { }
 
  TTNeutralGenes(const TTNeutralGenes& T)
    : _myProto(T._myProto), _type(NTRL) { }
 
  virtual ~TTNeutralGenes () {}
 
  void            set_proto          (TProtoNeutralGenes* proto) {_myProto = proto;}
 
  virtual TTNeutralGenes& operator= (const TTrait& T) = 0;
  virtual bool    operator== (const TTrait& T);
  virtual bool    operator!= (const TTrait& T);
  virtual void    init                 ( ) = 0;
  virtual void    init_sequence        ( ) = 0;
  virtual void    reset                ( ) = 0;
  virtual void*   set_trait            (void* value)           {return NULL;}
  virtual void**  get_sequence         ( )  const = 0;
  virtual unsigned int    get_allele   (int loc, int all) const = 0;
  virtual double  get_allele_value     (int loc, int all)  const = 0;
  virtual void    set_allele_value     (unsigned int locus, unsigned int allele, double value) = 0;
  virtual void    set_sequence         (void** seq) = 0;
  virtual trait_t get_type             ( )  const              {return _type;}
  virtual void    set_value            ( )                     { }
  virtual void*   getValue             ( )      const              {return 0;}
  virtual void    inherit              (const TTrait* mother, const TTrait* father) = 0;
  virtual void    mutate               ( ) = 0;
  virtual void    show_up              ( ) = 0;
  virtual TTNeutralGenes*  clone       ( ) = 0;
 
  //implements StorableComponent
  virtual void    store_data       (BinaryStorageBuffer* saver) = 0;
  virtual bool    retrieve_data    (BinaryStorageBuffer* reader) = 0;
};
// ----------------------------------------------------------------------------------------
 
//                                 N E U T R A L   T R A I T   (byte encoding)
 
// ----------------------------------------------------------------------------------------
class TTNeutralGenes_byte : public TTNeutralGenes
{
private:
  void (TTNeutralGenes_byte::* _mutate_func_ptr) (void);
  void (TProtoNeutralGenes::* _inherit_func_ptr) (sex_t, unsigned char*, const unsigned char**);
 
  unsigned char** _sequence; // ordered as [ploidy][locus]
 
public:
 
  TTNeutralGenes_byte ()
    : TTNeutralGenes(), _mutate_func_ptr(0), _inherit_func_ptr(0), _sequence(0) { }
 
  TTNeutralGenes_byte(const TTNeutralGenes_byte& T)
    : TTNeutralGenes(T), _mutate_func_ptr(T._mutate_func_ptr), _inherit_func_ptr(T._inherit_func_ptr),
    _sequence(0) { }
 
  virtual ~TTNeutralGenes_byte ();
 
  void            set_mut_func_ptr     (void(TTNeutralGenes_byte::* theFunc)(void))
  {_mutate_func_ptr = theFunc;}
  void            set_inherit_func_ptr (void(TProtoNeutralGenes::* theFunc)(sex_t, unsigned char*, const unsigned char**))
  {_inherit_func_ptr = theFunc;}
  void            set_allele           (unsigned int loc, unsigned int al, unsigned char val) {_sequence[al][loc] = val;}
 
  void            mutate_SSM           ( );
  void            mutate_KAM           ( );
  void            mutate_2all          ( );
  void            mutate_NULL          ( ) { }
  virtual TTNeutralGenes& operator= (const TTrait& T);
  virtual void    init                 ( );
  virtual void    init_sequence        ( );
  virtual void    reset                ( );
  virtual void**  get_sequence         ( )  const              {return (void**)_sequence;}
  virtual unsigned int    get_allele   (int loc, int all) const;
  virtual double  get_allele_value     (int loc, int all)  const ;
  virtual void    set_allele_value     (unsigned int locus, unsigned int allele, double value) ;
  virtual void    set_sequence         (void** seq);
  virtual void    inherit              (const TTrait* mother, const TTrait* father);
  virtual void    mutate               ( )                      {(this->*_mutate_func_ptr) ();}
  virtual void    show_up              ( );
  virtual TTNeutralGenes*  clone       ( )                      {return new TTNeutralGenes_byte(*this);}
 
  //implements StorableComponent
  virtual void    store_data       (BinaryStorageBuffer* saver);
  virtual bool    retrieve_data    (BinaryStorageBuffer* reader);
};  
// ----------------------------------------------------------------------------------------
 
//                     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 int  _2L; //genome size (ploidy*num_locus)
  unsigned int  _init_model;
  double _recomb_rate;
  void (TTNeutralGenes_byte::* _mutate_func_ptr) (void);
  void (TProtoNeutralGenes::* _inherit_func_ptr) (sex_t, unsigned char*, const unsigned char**);
 
  vector< TTNeutralGenesFH* > _writers;
  TTNOhtaStats* _ohtaStats;
  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;}
  double           get_mut_rate         ( )            {return _mut_rate;}
  unsigned int     get_genome_size      ( )            {return _2L;}
  unsigned int     get_num_mutations    ( )            {return RAND::Binomial( _mut_rate, _2L );}
  unsigned int     get_init_model       ( )            {return _init_model;}
  
  TTNeutralGenesSH* get_stater () {return _stats;}
 
  void    inherit_low                  (sex_t SEX, unsigned char* seq, const unsigned char** parent);
  void    inherit_free                 (sex_t SEX, unsigned char* seq, const unsigned char** parent);
 
  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;}
 
  virtual int                      get_phenotype_dimension () {return 0;}
 
  virtual int                      get_allele_number       () {return _allele_num;}
  virtual   int                    get_locus_number        () {return _locus_num;}
 
  virtual bool setParameters();
 
  virtual void loadFileServices ( FileServices* loader );
 
  virtual void loadStatServices ( StatServices* loader );
 
  virtual bool resetParameterFromSource (std::string param, SimComponent* cmpt); // {return false;}
 
  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_varcompWC();
  
  void setOutputOption(string opt);
  void set_write_fct( void (TTNeutralGenesFH::* fct_ptr) () ) {write_fct = fct_ptr;}
  
};
 
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTNOhtaStats : public TraitFileHandler<TProtoNeutralGenes> {
 
  TMatrix _pairwiseCombs;
 
public:
  TTNOhtaStats(TProtoNeutralGenes* T) : TraitFileHandler<TProtoNeutralGenes>(T,".nohta") {}
 
  virtual ~TTNOhtaStats () {}
  virtual void FHwrite ();
  virtual void FHread (string& filename) {}
};
 
 
// ----------------------------------------------------------------------------------------
 
//                     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;
  bool _is_diallelic_bitstring;
  
  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), _is_diallelic_bitstring(false), _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   setFstat_bitstring       (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_bitstring (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               (const TTrait* ind1, const TTrait* 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