TTQOhtaStats Class Reference

TTQOhtaStats. More...

#include <ttquanti.h>

Inheritance diagram for TTQOhtaStats:

Collaboration diagram for TTQOhtaStats:

Public Member Functions
	TTQOhtaStats (TProtoQuanti *T)
virtual	~TTQOhtaStats ()
virtual void	FHwrite ()
virtual void	FHread (string &filename)
Public Member Functions inherited from TraitFileHandler< TProtoQuanti >
	TraitFileHandler (TProtoQuanti *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, TProtoQuanti *trait_proto)
virtual void	set_multi (bool rpl_per, bool gen_per, int rpl_occ, TMatrix *Occ, string path, TProtoQuanti *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< TProtoQuanti >
TProtoQuanti *	_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

TTQOhtaStats.

Constructor & Destructor Documentation

TTQOhtaStats()

TTQOhtaStats::TTQOhtaStats ( TProtoQuanti *  T )

inline

: TraitFileHandler<TProtoQuanti>(T,".qohta") {}

~TTQOhtaStats()

virtual TTQOhtaStats::~TTQOhtaStats ( )

inlinevirtual

{}

Member Function Documentation

FHread()

virtual void TTQOhtaStats::FHread ( string &  filename )

inlinevirtual

Implements FileHandler.

{}

FHwrite()

void TTQOhtaStats::FHwrite ( )

virtual

Implements TraitFileHandler< TProtoQuanti >.

{
 
  //  Output Ohta stats :
  //    Only if pop alive
  //    Only for diallelic mutation models
  //    Only for a single trait
  //    Only for multi-population scenarios
  if ((!_pop->isAlive()) || (_FHLinkedTrait->get_allele_model() > 2)
        || (_FHLinkedTrait->get_num_traits() > 1) || (_pop->getPatchNbr() < 2)) {
    error("TTQOhtaStats::quanti trait and population are not compatible with Ohta stats (min. 2 patches, one di-allelic trait), file will not be written.\n");
    return ;
  }
 
  // sets the pop ptr to a sub sampled pop, if sub sampling happened
  if(!get_service()) fatal("TTQuantiFH: link to file services not set!!\n");
  Metapop* pop =  get_service()->getSampledPop();
  int patchNbr = pop->getPatchNbr();
  unsigned int nb_trait = _FHLinkedTrait->get_num_traits();
  unsigned int nb_locus = _FHLinkedTrait->get_num_locus();
  TTQuanti_diallelic* trait;
 
 
  _pairwiseCombs.copy(nChooseKVec(nb_locus, 2));
 
  unsigned int num_comb = _pairwiseCombs.nrows();
 
  ostringstream REC;
  Individual* ind;
  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("TTQOhtaStats::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 > meanAlleleFreq(4,0.0); // allele counter at the two loci,"A" = [0], "a" = [1], "B" = [2], "b" = [3]
  vector< double > meanHapFreq(4,0.0);  // genotype counter
 
  vector< vector< double > > alleleFreq = 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);
 
    meanAlleleFreq.assign(4, 0.0);
    meanHapFreq.assign(4,0.0);
 
    for(int patch = 0; patch < patchNbr; patch++) {
 
      current_patch = pop->getPatch(patch);
 
      alleleFreq[patch].assign(4,0.0);
      hapFreq[patch].assign(4,0.0);
 
      if(patchSizes[patch] == 0) continue;
 
      // @TODO adapt Ohta stat output for any age class not only adults
      //  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 = dynamic_cast<TTQuanti_diallelic*>( ind->getTrait(_FHLinkedTraitIndex));
 
          //  Homologous copy 1
          //  locus 1;  "A" = 0, "a" = 1
          if ( !trait->get_allele_bit(loc1, 0) ) { //reference allele is allele '0'
            ++alleleFreq[patch][0];
            a1 = 0;
          }
          else {
            ++alleleFreq[patch][1];
            a1 = 1;
          }
 
          //  locus 2
          if (!trait->get_allele_bit(loc2, 0)) {
            ++alleleFreq[patch][2];
            a2 = 0;
          }
          else {
            ++alleleFreq[patch][3];
            a2 = 1;
          }
 
          ++hapFreq[patch][ twoLocHapMap[a1][a2] ];
 
          //  Homologous copy 2
          //  locus 1
          if (!trait->get_allele_bit(loc1, 1)) {
            ++alleleFreq[patch][0];
            a1 = 0;
          }
          else {
            ++alleleFreq[patch][1];
            a1 = 1;
          }
 
          //  locus 2
          if (!trait->get_allele_bit(loc2, 1)) {
            ++alleleFreq[patch][2];
            a2 = 0;
          }
          else {
            ++alleleFreq[patch][3];
            a2 = 1;
          }
 
          ++hapFreq[patch][ twoLocHapMap[a1][a2] ];
        } // all individuals
      } // two sexes
 
 
      for (size_t geno = 0; geno < 4; geno++) {
        alleleFreq[patch][geno] /= patchSizes[patch];
        meanAlleleFreq[ geno ] += alleleFreq[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++)
      meanAlleleFreq[geno] /= extant_patches;
 
    for (size_t hap = 0; hap < 4; hap++)
      meanHapFreq[hap] /= extant_patches;
 
    //  eliminate pairs with fixed alleles
    if ( !(meanAlleleFreq[0]*meanAlleleFreq[1]) && !(meanAlleleFreq[2] * meanAlleleFreq[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] -
                            (alleleFreq[patch][ reverseHapMap[hap][0] ] * alleleFreq[patch][ reverseHapMap[hap][1]] ), 2);
 
          //  Ohta (1982), Eq. 11
          Dst[pcomb] += pow((alleleFreq[patch][reverseHapMap[hap][0]] * alleleFreq[patch][reverseHapMap[hap][1]]) -
                            (meanAlleleFreq[reverseHapMap[hap][0]] * meanAlleleFreq[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] -
                             (meanAlleleFreq[reverseHapMap[hap][0]] * meanAlleleFreq[reverseHapMap[hap][1]]), 2);
        }
 
        double denom = alleleFreq[patch][0] * alleleFreq[patch][1] * alleleFreq[patch][2] * alleleFreq[patch][3];
 
        if (denom == 0.0)
          rSquare[patch][pcomb] = 0;
        else
          rSquare[patch][pcomb] = pow(hapFreq[patch][0] - alleleFreq[patch][0]*alleleFreq[patch][2], 2)   // p(AB) - p(A)*p(B)
                            / 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 < 4; geno++)
//        cout << geno << " : " << alleleFreq[patch][geno] << endl;
//      for (size_t hap = 0; hap < 4; hap++)
//        cout << hap << " : " << hapFreq[patch][hap] << endl;
//      for (size_t hap = 0; hap < 4; hap++)
//        cout << "Exp_" << reverseHapMap[hap][0] << reverseHapMap[hap][1] << " : "
//        << (alleleFreq[patch][reverseHapMap[hap][0]] * alleleFreq[patch][reverseHapMap[hap][1]]) << endl;
//    }
//    cout << "*****************" << endl;
//    for (size_t geno = 0; geno < 4; geno++)
//      cout << geno << " : " << meanAlleleFreq[geno] << endl;
//    for (size_t hap = 0; hap < 4; hap++)
//      cout << hap << " : " << meanHapFreq[hap] << endl;
//    for (size_t hap = 0; hap < 4; hap++)
//      cout << "Exp_" << reverseHapMap[hap][0] << reverseHapMap[hap][1] << " : "
//      << (meanAlleleFreq[reverseHapMap[hap][0]] * meanAlleleFreq[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 quanti could not open output file: \"%s\"\n",filename.c_str());
 
  #ifdef _DEBUG_
  message("TTQOhtaStats::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< TProtoQuanti >::_FHLinkedTrait, TraitFileHandler< TProtoQuanti >::_FHLinkedTraitIndex, _pairwiseCombs, FileHandler::_pop, ADLTx, Patch::assign(), TMatrix::copy(), error(), fatal(), FEM, Patch::get(), TMatrix::get(), TTQuanti_diallelic::get_allele_bit(), TProtoQuanti::get_allele_model(), FileHandler::get_filename(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), FileHandler::get_service(), Metapop::getPatch(), Metapop::getPatchNbr(), FileServices::getSampledPop(), Individual::getTrait(), Metapop::isAlive(), MAL, message(), nChooseKVec(), TMatrix::nrows(), and Patch::size().

Member Data Documentation

_pairwiseCombs

TMatrix TTQOhtaStats::_pairwiseCombs

private

Referenced by FHwrite().

---

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

• ttquanti.h
• ttquanti.cc