TTNOhtaStats Class Reference

TTNOhtaStats. More...

#include <ttneutralgenes.h>

Inheritance diagram for TTNOhtaStats:

Collaboration diagram for TTNOhtaStats:

Public Member Functions
	TTNOhtaStats (TProtoNeutralGenes *T)
virtual	~TTNOhtaStats ()
virtual void	FHwrite ()
virtual void	FHread (string &filename)
Public Member Functions inherited from TraitFileHandler< TProtoNeutralGenes >
	TraitFileHandler (TProtoNeutralGenes *trait_proto, const char *ext)
virtual	~TraitFileHandler ()
virtual void	FHread (string &filename)=0
virtual void	set (bool rpl_per, bool gen_per, int rpl_occ, int gen_occ, int rank, string path, TProtoNeutralGenes *trait_proto)
virtual void	set_multi (bool rpl_per, bool gen_per, int rpl_occ, TMatrix *Occ, string path, TProtoNeutralGenes *trait_proto)
Public Member Functions inherited from FileHandler
	FileHandler (const char *ext)
virtual	~FileHandler ()
virtual void	init ()
	Called by notifier during simulation setup, performs file checking. More...
virtual vector< string >	ifExist ()
	Checks if any file associated with the current file name already exists on disk. More...
virtual void	set (bool rpl_per, bool gen_per, int rpl_occ, int gen_occ, int rank, string path)
	Sets the hanlder parameters. More...
virtual void	set_multi (bool rpl_per, bool gen_per, int rpl_occ, TMatrix *Occ, string path)
virtual void	update ()
	Updates the inner replicate and generation counters and calls FHwrite if needed by the the periodicity of the file. More...
Metapop *	get_pop_ptr ()
	Returns the pointer to the current metapop through the FileServices interface. More...
void	set_pop_ptr (Metapop *pop_ptr)
FileServices *	get_service ()
	Returns pointer to the FileServices. More...
void	set_service (FileServices *srv)
std::string &	get_path ()
void	set_path ()
std::string &	get_extension ()
void	set_extension (const char *ext)
std::string &	get_filename ()
	Builds and returns the current file name depending on the periodicity of the file. More...
bool	get_isInputHandler ()
void	set_isInputHandler (bool val)
bool	get_isReplicatePeriodic ()
void	set_isReplicatePeriodic (bool val)
unsigned int	get_ReplicateOccurrence ()
void	set_ReplicateOccurrence (unsigned int val)
bool	get_isGenerationPeriodic ()
void	set_isGenerationPeriodic (bool val)
unsigned int	get_GenerationOccurrence ()
void	set_GenerationOccurrence (unsigned int val)
unsigned int	get_ExecRank ()
	unused yet... More...
void	set_ExecRank (int val)
TMatrix *	get_OccMatrix ()
void	set_OccMatrix (TMatrix *occ)
bool	get_isMasterExec ()
void	set_isMasterExec (bool is)
Public Member Functions inherited from Handler
virtual	~Handler ()

Private Attributes
TMatrix	_pairwiseCombs

Additional Inherited Members
Protected Attributes inherited from TraitFileHandler< TProtoNeutralGenes >
TProtoNeutralGenes *	_FHLinkedTrait
int	_FHLinkedTraitIndex
Protected Attributes inherited from FileHandler
Metapop *	_pop
	Pointer to the current metapop, set during initialization within the init function. More...

Detailed Description

TTNOhtaStats.

Constructor & Destructor Documentation

TTNOhtaStats()

TTNOhtaStats::TTNOhtaStats ( TProtoNeutralGenes *  T )

inline

: TraitFileHandler<TProtoNeutralGenes>(T,".nohta") {}

~TTNOhtaStats()

virtual TTNOhtaStats::~TTNOhtaStats ( )

inlinevirtual

{}

Member Function Documentation

FHread()

virtual void TTNOhtaStats::FHread ( string &  filename )

inlinevirtual

Implements FileHandler.

{}

FHwrite()

void TTNOhtaStats::FHwrite ( )

virtual

Implements TraitFileHandler< TProtoNeutralGenes >.

{
  Metapop* pop = get_pop_ptr();
  int patchNbr = pop->getPatchNbr();
  unsigned int nb_allele = _FHLinkedTrait->get_allele_num();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
 
  //  Output Ohta stats :
  //    Only if pop alive
  //    Only for diallelic mutation models
  //    Only for multi-population scenarios
  if ((!pop->isAlive()) || (nb_allele > 2) || (patchNbr < 2)) {
    error("TTNOhtaStats::neutral trait and population are not compatible with Ohta stats (min. 2 patches, di-allelic loci), file will not be written.\n");
    return;
  }
 
  _pairwiseCombs.copy(nChooseKVec(nb_locus, 2));
 
  unsigned int num_comb = _pairwiseCombs.nrows();
 
  ostringstream REC;
  Individual* ind;
  TTrait* trait;
  Patch* current_patch;
 
  vector<double> Dis(num_comb, 0.0),
                 Dst(num_comb, 0.0),
                 Disp(num_comb, 0.0),
                 Dstp(num_comb, 0.0);
 
  vector< vector<double> > rSquare =  vector< vector<double> > (patchNbr, vector<double>(num_comb, 0.0));
 
  vector<bool> NA(num_comb, false);
 
  int total_size = 0;
  int extant_patches = 0;
 
  vector<int> patchSizes(patchNbr, 0);
 
  // pre-record patch sizes and number of extant patches
  for(int patch = 0; patch < patchNbr; patch++) {
 
    current_patch = pop->getPatch(patch);
 
    patchSizes[patch] = 2 * (current_patch->size(FEM, ADLTx) + current_patch->size(MAL, ADLTx));
 
    if(patchSizes[patch]) ++extant_patches;  // count patches with individuals
 
    total_size += patchSizes[patch];
  }
 
#ifdef _DEBUG_
  message("TTNOhtaStats::FHwrite:: computing association stats of %i combinations\n", num_comb);
  fflush(stdout);
#endif
 
 
  unsigned int a1, a2;
  unsigned int twoLocHapMap[2][2] = {{0,1},{2,3}}; // genot: AB, Ab, aB, ab; A = 0, a = 1, B = 2, b = 3
  unsigned int reverseHapMap[4][2] = {{0,2},{0,3},{1,2},{1,3}}; // AB -> {0,2} ...
 
  vector< double > meanGenoFreq(4,0.0); // allele counter at the two loci,0 = [0], 1 = [1], 0 = [2], 1 = [3]
  vector< double > meanHapFreq(4,0.0);  // genotype counter
 
  vector< vector< double > > genoFreq = vector< vector< double > > (patchNbr, vector< double >(4,0.0));
  vector< vector< double > > hapFreq  = vector< vector< double > > (patchNbr, vector< double >(4,0.0));
 
  //  cycling through all pairwise locus combinations
  for (size_t pcomb = 0; pcomb < num_comb; pcomb++) {
 
    size_t loc1 = _pairwiseCombs.get(pcomb,0),
           loc2 = _pairwiseCombs.get(pcomb,1);
 
    unsigned int refLoc1, refLoc2;
 
    current_patch = pop->getPatch(0);
 
    ind = current_patch->get(FEM, ADLTx, 0);
 
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    //  reference alleles from the first female adult of the first patch
    refLoc1 = trait->get_allele(loc1, 0);
    refLoc2 = trait->get_allele(loc2, 0);
 
    meanGenoFreq.assign(4, 0.0);
    meanHapFreq.assign(4,0.0);
 
    for(int patch = 0; patch < patchNbr; patch++) {
 
      current_patch = pop->getPatch(patch);
 
      genoFreq[patch].assign(4,0.0);
      hapFreq[patch].assign(4,0.0);
 
      if(patchSizes[patch] == 0) continue;
 
      //  All adults in a patch
      for(unsigned int s =0; s < 2; ++s) { // two sexes
 
        for(unsigned int j = 0, size = current_patch->size(sex_t(s), ADLTx); j < size; j++) {
 
          ind = current_patch->get(sex_t(s), ADLTx, j);
 
          trait = ind->getTrait(_FHLinkedTraitIndex);
 
          //  Homologous copy 1
          //  locus 1
          if (trait->get_allele(loc1, 0) == refLoc1) {
            ++genoFreq[patch][0];
            a1 = 0;
          }
          else {
            ++genoFreq[patch][1];
            a1 = 1;
          }
 
          //  locus 2
          if (trait->get_allele(loc2, 0) == refLoc2) {
            ++genoFreq[patch][2];
            a2 = 0;
          }
          else {
            ++genoFreq[patch][3];
            a2 = 1;
          }
 
          ++hapFreq[patch][twoLocHapMap[a1][a2]];
 
          //  Homologous copy 2
          //  locus 1
          if (trait->get_allele(loc1, 1) == refLoc1) {
            ++genoFreq[patch][0];
            a1 = 0;
          }
          else {
            ++genoFreq[patch][1];
            a1 = 1;
          }
 
          //  locus 2
          if (trait->get_allele(loc2, 1) == refLoc2) {
            ++genoFreq[patch][2];
            a2 = 0;
          }
          else {
            ++genoFreq[patch][3];
            a2 = 1;
          }
 
          ++hapFreq[patch][twoLocHapMap[a1][a2]];
        } //  All individuals
      } // two sexes
 
      for (size_t geno = 0; geno < 4; geno++) {
        genoFreq[patch][geno] /= patchSizes[patch];
        meanGenoFreq[geno] += genoFreq[patch][geno];
      }
 
      for (size_t hap = 0; hap < 4; hap++) {
        hapFreq[patch][hap] /= patchSizes[patch];
        meanHapFreq[hap] += hapFreq[patch][hap];
      }
    } //  All patches
 
    for (size_t geno = 0; geno < 4; geno++)
       meanGenoFreq[geno] /= extant_patches;
 
     for (size_t hap = 0; hap < 4; hap++)
       meanHapFreq[hap] /= extant_patches;
 
    //  eliminate pairs with fixed alleles
    if ( !(meanGenoFreq[0]*meanGenoFreq[1]) && !(meanGenoFreq[2] * meanGenoFreq[3]))
      NA[pcomb] = true;
 
    if (!NA[pcomb]) {
 
      for(int patch = 0; patch < patchNbr; patch++) {
 
        if(!patchSizes[patch]) continue;
 
        for (size_t hap = 0; hap < 4; hap++) {
          //  Ohta (1982), Eq. 10
          Dis[pcomb] += pow(hapFreq[patch][hap] -
                            (genoFreq[patch][reverseHapMap[hap][0]] * genoFreq[patch][reverseHapMap[hap][1]]), 2);
 
          //  Ohta (1982), Eq. 11
          Dst[pcomb] += pow((genoFreq[patch][reverseHapMap[hap][0]] * genoFreq[patch][reverseHapMap[hap][1]]) -
                            (meanGenoFreq[reverseHapMap[hap][0]] * meanGenoFreq[reverseHapMap[hap][1]]), 2);
 
          //  Ohta (1982), Eq. 12
          Disp[pcomb] += pow(hapFreq[patch][hap] - meanHapFreq[hap], 2);
 
          //  Ohta (1982), Eq. 13
          Dstp[pcomb] += pow(meanHapFreq[hap] -
                             (meanGenoFreq[reverseHapMap[hap][0]] * meanGenoFreq[reverseHapMap[hap][1]]), 2);
        }
 
        double denom = genoFreq[patch][0] * genoFreq[patch][1] * genoFreq[patch][2] * genoFreq[patch][3];
 
        if (denom == 0.0)
          rSquare[patch][pcomb] = 0;
        else
          rSquare[patch][pcomb] = pow(hapFreq[patch][0] - genoFreq[patch][0]*genoFreq[patch][2], 2)
                            / denom;
      }
 
      //  Ohta (1982), Eq. 10-13
      Dis[pcomb] /= extant_patches;
      Dst[pcomb] /= extant_patches;
      Disp[pcomb] /= extant_patches;
      Dstp[pcomb] /= extant_patches;
 
    }
 
//    cout << "\nLocus " << loc1+1 << " and locus " << loc2+1 << endl;
//    cout << "\"A allele\" : " << refLoc1 << "\n\"B allele\" : " << refLoc2 << endl;
//    for(int patch = 0; patch < patchNbr; patch++) {
//      cout << "Patch " << patch+1 << endl;
//      for (size_t geno = 0; geno < Gtypes.size(); geno++)
//        cout << geno << " : " << genoFreq[patch][geno] << endl;
//      for (size_t hap = 0; hap < Htypes.size(); hap++)
//        cout << hap << " : " << hapFreq[patch][hap] << endl;
//      for (size_t hap = 0; hap < Htypes.size(); hap++)
//        cout << "Exp_" << reverseHapMap[hap][0] << reverseHapMap[hap][1] << " : " << (genoFreq[patch][reverseHapMap[hap][0]] * genoFreq[patch][reverseHapMap[hap][1]]) << endl;
//    }
//    cout << "*****************" << endl;
//    for (size_t geno = 0; geno < Gtypes.size(); geno++)
//      cout << geno << " : " << meanGenoFreq[geno] << endl;
//    for (size_t hap = 0; hap < Htypes.size(); hap++)
//      cout << hap << " : " << meanHapFreq[hap] << endl;
//    for (size_t hap = 0; hap < Htypes.size(); hap++)
//      cout << "Exp_" << reverseHapMap[hap][0] << reverseHapMap[hap][1] << " : " << (meanGenoFreq[reverseHapMap[hap][0]] * meanGenoFreq[reverseHapMap[hap][1]]) << endl;
//    cout << "*****************" << endl;
 
  } //  All pairwise locus combinations
 
  string filename = get_filename();
  ofstream FILE;
 
  //open file for writing
  FILE.open(filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
  if(!FILE) fatal("Trait neutral could not open output file: \"%s\"\n",filename.c_str());
 
#ifdef _DEBUG_
  message("TTNOhtaStats::FHwrite (%s)\n",filename.c_str());
  fflush(stdout);
#endif
 
  FILE << "loc1\tloc2\tDst\tDis\tDstp\tDisp";
  for(int patch = 0; patch < patchNbr; patch++)
    FILE << "\tr_" << patch+1;
  FILE << endl;
 
  for (size_t pcomb = 0; pcomb < num_comb; pcomb++) {
    if (!NA[pcomb]) {
      FILE << _pairwiseCombs.get(pcomb,0)+1 << "\t" << _pairwiseCombs.get(pcomb,1)+1 << "\t"
           << Dst[pcomb] << "\t" << Dis[pcomb] << "\t" << Dstp[pcomb] << "\t" << Disp[pcomb];
      for(int patch = 0; patch < patchNbr; patch++)
        FILE << "\t" << rSquare[patch][pcomb];
      FILE << endl;
    }
  }
 
  FILE.close();
  std::ios_base::sync_with_stdio(true); // reset
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, _pairwiseCombs, ADLTx, Patch::assign(), TMatrix::copy(), error(), fatal(), FEM, Patch::get(), TMatrix::get(), TTrait::get_allele(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_filename(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_pop_ptr(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), Metapop::isAlive(), MAL, message(), nChooseKVec(), TMatrix::nrows(), and Patch::size().

Member Data Documentation

_pairwiseCombs

TMatrix TTNOhtaStats::_pairwiseCombs

private

Referenced by FHwrite().

---

The documentation for this class was generated from the following files:

• ttneutralgenes.h
• ttneutralgenes.cc