Nemo  2.3.46
TTNeutralGenesFH Class Reference

A file handler to save the neutral markers genotypes in the FSTAT format (extended). More...

#include <ttneutralgenes.h>

+ Inheritance diagram for TTNeutralGenesFH:
+ Collaboration diagram for TTNeutralGenesFH:

Public Member Functions

 TTNeutralGenesFH (TProtoNeutralGenes *TP)
 
virtual ~TTNeutralGenesFH ()
 
virtual void FHwrite ()
 
virtual void FHread (string &filename)
 
void write_TAB ()
 
void write_GENEPOP ()
 
void write_FSTAT ()
 
void write_Fst_i ()
 
void write_freq ()
 
void write_varcompWC ()
 
void setOutputOption (string opt)
 
void set_write_fct (void(TTNeutralGenesFH::*fct_ptr)())
 
- Public Member Functions inherited from TraitFileHandler< TProtoNeutralGenes >
 TraitFileHandler (TProtoNeutralGenes *trait_proto, const char *ext)
 
virtual ~TraitFileHandler ()
 
virtual void set (bool rpl_per, bool gen_per, int rpl_occ, int gen_occ, int rank, 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 bool 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...
 
Metapopget_pop_ptr ()
 Returns the pointer to the current metapop through the FileServices interface. More...
 
void set_pop_ptr (Metapop *pop_ptr)
 
FileServicesget_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)
 
TMatrixget_OccMatrix ()
 
void set_OccMatrix (TMatrix *occ)
 
bool get_isMasterExec ()
 
void set_isMasterExec (bool is)
 
- Public Member Functions inherited from Handler
virtual ~Handler ()
 

Private Attributes

string _output_option
 
void(TTNeutralGenesFH::* write_fct )()
 

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

A file handler to save the neutral markers genotypes in the FSTAT format (extended).

By default, the file extension is ".dat" for the genotype file. It is changed to ".fsti" for the per-locus/per-patch Weir& Hill Fst's and ".freq" for the per-locus/per-patch allele frequencies. Also implements the FileHandler::FHread method to load a population's genotypes from an FSTAT file (using the 'source_pop' population parameter).

Constructor & Destructor Documentation

TTNeutralGenesFH::TTNeutralGenesFH ( TProtoNeutralGenes TP)
inline
virtual TTNeutralGenesFH::~TTNeutralGenesFH ( )
inlinevirtual
243 { }

Member Function Documentation

void TTNeutralGenesFH::FHread ( string &  filename)
virtual

reads in only FSTAT format

Implements TraitFileHandler< TProtoNeutralGenes >.

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FileHandler::_pop, Patch::add(), ADLTx, ADULTS, error(), fatal(), FEM, TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), IndFactory::makeNewIndividual(), OFFSx, TTrait::set_sequence(), and Individual::setPedigreeClass().

962 {
963  unsigned int digit, nloci_infile, nall_infile, npatch_infile;
964  // unsigned int ploidy = _FHLinkedTrait->get_ploidy();
965  unsigned int nb_all = _FHLinkedTrait->get_allele_num();
966  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
967  unsigned int patchNbr = _pop->getPatchNbr();
968 
969  bool is_extended = false;
970 
971  ifstream FILE(filename.c_str(),ios::in);
972 
973  if(!FILE) fatal("could not open FSTAT input file \"%s\"\n", filename.c_str());
974 
975  FILE >> npatch_infile >> nloci_infile >> nall_infile >> digit;
976 
977  if(npatch_infile != patchNbr) fatal("# of patch in FSTAT file differs from simulation settings\n");
978 
979  if(nloci_infile != nb_locus && nloci_infile != nb_locus +4) {
980 
981  fatal("# of loci in FSTAT file (%i loci) differs from simulation settings (%i loci)\n",nloci_infile, nb_locus);
982 
983  } else if (nloci_infile == nb_locus + 4) {
984  is_extended = true;
985  } else
986  is_extended = false;
987 
988  if(nall_infile != nb_all) fatal("# of alleles in FSTAT file differs from simulation settings\n");
989 
990  digit = (int) pow(10.0,(double)digit);
991  string loc_name;
992 
993  for(unsigned int i = 0; i < nloci_infile; ++i) {
994  FILE >> loc_name; //we don't do anything with this here
995  }
996 
997 
998  unsigned int pop, genot, all0, all1, age, sex, ped, origin;
999  age_idx agex;
1000  Individual *ind;
1001  unsigned char* seq[2];
1002  int lnbr = nloci_infile +2; //line counter
1003  seq[0] = new unsigned char[nb_locus];
1004  seq[1] = new unsigned char[nb_locus];
1005 
1006  while(FILE>>pop) {
1007  //check pop value
1008  if(pop > patchNbr) fatal("found an illegal pop identifier in FSTAT file (%i) at line %i",pop, lnbr);
1009  for(unsigned int i = 0; i < nb_locus; ++i) {
1010  FILE>>genot;
1011 
1012  all0 = genot/digit;
1013  all1 = genot%digit;
1014 
1015  if(all0 <= nb_all)
1016  seq[0][i] = all0 - 1;
1017  else {
1018  error("in FSTAT input file at line %i, locus %i : \
1019 first allele value %d is greater than the max value specified (%i)!\n", lnbr, i+1, all0, nb_all);
1020  fatal("Please check the input file.\n");
1021  }
1022 
1023  if(all1 <= nb_all)
1024  seq[1][i] = all1 - 1;
1025  else {
1026  error("in FSTAT input file at line %i, locus %i : \
1027 second allele value %i is greater than the max value specified (%i)!\n", lnbr, i+1, all1, nb_all);
1028  fatal("Please check the input file.\n");
1029  }
1030  }
1031 
1032  if(is_extended) FILE >> age >> sex >> ped >> origin;
1033  else {
1034  age = ADULTS;
1035  sex = FEM;
1036  ped = 0;
1037  origin = 0;
1038  }
1039 
1040  agex = (age == ADULTS ? ADLTx : OFFSx);
1041 
1042  ind = _pop->makeNewIndividual(0, 0, static_cast<sex_t> (sex), origin - 1);
1043  ind->setPedigreeClass((unsigned char)ped);
1044  ind->getTrait(_FHLinkedTraitIndex)->set_sequence((void**)seq);
1045  _pop->getPatch(pop-1)->add(static_cast<sex_t> (sex), agex, ind);
1046 
1047  lnbr++;
1048 
1049  }
1050 
1051  FILE.close();
1052 
1053  delete [] seq[0];
1054  delete [] seq[1];
1055 }
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
#define ADULTS
Adults age class flag (breeders).
Definition: types.h:54
age_idx
Array index of the age classes in the patch sizes and containers arrays.
Definition: types.h:41
TTrait * getTrait(IDX T)
Trait accessor.
Definition: individual.h:277
Individual * makeNewIndividual(Individual *mother, Individual *father, sex_t sex, unsigned short homepatch)
Creates an individual with pointers to parents, sex and home ID set but no genetic data...
Definition: indfactory.cc:152
void setPedigreeClass(Individual *mother, Individual *father)
Definition: individual.h:115
virtual void set_sequence(void **seq)=0
Called to set the sequence pointer to an existing trait.
void fatal(const char *str,...)
Definition: output.cc:90
Patch * getPatch(unsigned int i)
Patch accessor, return the ith+1 patch in the metapop.
Definition: metapop.h:251
unsigned int getPatchNbr()
Definition: metapop.h:270
int error(const char *str,...)
Definition: output.cc:73
Definition: types.h:42
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
Definition: types.h:42
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
This class contains traits along with other individual information (sex, pedigree, etc. ).
Definition: individual.h:49
void add(sex_t SEX, age_idx AGE, Individual *ind)
Adds an individual to the appropriate container, increments its size, eventually resizing it...
Definition: metapop.h:540
Definition: types.h:37
int _FHLinkedTraitIndex
Definition: filehandler.h:220
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
void TTNeutralGenesFH::FHwrite ( )
virtual

Implements TraitFileHandler< TProtoNeutralGenes >.

References FileHandler::_pop, Metapop::isAlive(), and write_fct.

596 {
597  if(!_pop->isAlive()) return;
598 
599  (this->*write_fct)();
600 }
bool isAlive()
Checks if the population still contains at least one individual in any sex or age class...
Definition: metapop.h:299
void(TTNeutralGenesFH::* write_fct)()
Definition: ttneutralgenes.h:235
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
void TTNeutralGenesFH::set_write_fct ( void(TTNeutralGenesFH::*)()  fct_ptr)
inline

References write_fct.

Referenced by TProtoNeutralGenes::loadFileServices().

257 {write_fct = fct_ptr;}
void(TTNeutralGenesFH::* write_fct)()
Definition: ttneutralgenes.h:235
void TTNeutralGenesFH::setOutputOption ( string  opt)

References _output_option.

Referenced by TProtoNeutralGenes::loadFileServices().

951 {
952 
953  if(opt != "1")
954  _output_option = opt;
955  else
956  _output_option = "locus"; //default
957 }
string _output_option
Definition: ttneutralgenes.h:233
void TTNeutralGenesFH::write_freq ( )

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, ADULTS, ALL, TTNeutralGenesSH::allocateTables(), fatal(), DataTable< T >::get(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), FileHandler::get_service(), TProtoNeutralGenes::get_stater(), TTNeutralGenesSH::getAlleleFreqTable(), Metapop::getCurrentAge(), FileServices::getGenerationReplicateFileName(), TTNeutralGenesSH::getHeteroTable(), Metapop::getPatchNbr(), message(), OFFSPRG, TTNeutralGenesSH::setAlleleTables(), TTNeutralGenesSH::setHeteroTable(), Metapop::size(), and warning().

1123 {
1124  unsigned int nb_allele = _FHLinkedTrait->get_allele_num();
1125  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1126  unsigned int patchNbr = _pop->getPatchNbr();
1127  bool added_stater = false;
1128  age_t AGE = _pop->getCurrentAge();
1129 
1130  if(AGE == (ADULTS | OFFSPRG) || AGE == ALL) {
1131  warning("saving only offspring neutral allele frequencies\n");
1132  AGE = OFFSPRG;
1133  }
1134 
1136 
1137  if(stater == NULL) {
1138  stater = new TTNeutralGenesSH(_FHLinkedTrait);
1139  added_stater = true;
1140  }
1141  //do this here, population may have been resized...
1142  stater->allocateTables(nb_locus, nb_allele);
1143 
1144  stater->setAlleleTables(AGE);
1145  stater->setHeteroTable (AGE);
1146 
1147  DataTable<double> *freqTable = stater->getAlleleFreqTable();
1148  DataTable<double> *hTable = stater->getHeteroTable();
1149 
1150  std::string filename = get_path() +
1152  string file1 = filename + get_extension();
1153  string file2 = filename + ".Ho";
1154  string file3 = filename + ".He";
1155 
1156  unsigned int *sizes = new unsigned int [patchNbr];
1157 
1158  for(unsigned int i = 0; i < patchNbr; ++i)
1159  sizes[i] = _pop->size(AGE, i);
1160 
1161 #ifdef _DEBUG_
1162  message("TTNeutralGenesFH::write_freq (%s)\n",filename.c_str());
1163 #endif
1164 
1165  ofstream FILE1 (file1.c_str(), ios::out);
1166  ofstream FILE2 (file2.c_str(), ios::out);
1167  ofstream FILE3 (file3.c_str(), ios::out);
1168 
1169  if(!FILE1) fatal("could not open neutral allele freq output file!!\n");
1170  if(!FILE2) fatal("could not open neutral allele heterosygosity output file!!\n");
1171  if(!FILE3) fatal("could not open neutral allele heterosygosity output file!!\n");
1172 
1173  for(unsigned int i = 0; i < patchNbr; ++i) {
1174 
1175  for(unsigned int u = 0; u < nb_allele; ++u)
1176  FILE1<<"p"<<i+1<<"a"<<u+1<<" ";
1177 
1178  FILE2<<"p"<<i+1<<"Ho ";
1179  FILE3<<"p"<<i+1<<"He ";
1180  }
1181  FILE1<<endl; FILE2<<endl; FILE3<<endl;
1182 
1183  double p, p2;
1184  for(unsigned int l = 0; l < nb_locus; ++l) {
1185  for(unsigned int i = 0; i < patchNbr; ++i) {
1186  p2 = 0;
1187  for(unsigned int u = 0; u < nb_allele; u++){
1188  p = freqTable->get(i, l, u);
1189  FILE1<< p << " ";
1190  p2 += p*p; //He = 1 - sum of p2
1191  }
1192  FILE2<<hTable->get(i, l, 0)/sizes[i]<<" ";
1193  FILE3<<1-p2<<" ";
1194  }
1195  FILE1<<endl; FILE2<<endl; FILE3<<endl;
1196  }
1197 
1198  FILE1.close(); FILE2.close(); FILE3.close();
1199 
1200  if(added_stater) delete stater;
1201 
1202  delete [] sizes;
1203 }
void setHeteroTable(age_t AGE)
Definition: stats_fstat.cc:166
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
#define ADULTS
Adults age class flag (breeders).
Definition: types.h:54
void allocateTables(unsigned int loci, unsigned int all)
Definition: stats_fstat.cc:43
void setAlleleTables(age_t AGE)
Definition: stats_fstat.cc:74
unsigned int age_t
Age class flags.
Definition: types.h:46
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
unsigned int getPatchNbr()
Definition: metapop.h:270
T get(unsigned int group, unsigned int Class, unsigned int elmnt)
Returns value stored of the element 'elmnt' of the class 'Class' in the group 'group'.
Definition: datatable.h:228
void warning(const char *str,...)
Definition: output.cc:56
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
#define ALL
All ages age class flag.
Definition: types.h:56
#define OFFSPRG
Offspring age class flag.
Definition: types.h:50
unsigned int size()
Get the total number of individuals present in the population, all sex and age classes together...
Definition: metapop.h:302
std::string & get_path()
Definition: filehandler.h:139
DataTable< double > * getAlleleFreqTable()
Accessor to the table of allele frequencies, per patch.
Definition: ttneutralgenes.h:361
DataTable< double > * getHeteroTable()
Definition: ttneutralgenes.h:365
The stat handler for neutral markers.
Definition: ttneutralgenes.h:268
age_t getCurrentAge()
Definition: metapop.h:289
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
TTNeutralGenesSH * get_stater()
Definition: ttneutralgenes.h:191
void TTNeutralGenesFH::write_Fst_i ( )

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, ADULTS, TTNeutralGenesSH::allocateTables(), fatal(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), FileHandler::get_service(), TProtoNeutralGenes::get_stater(), FileServices::getGenerationReplicateFileName(), Metapop::getPatchNbr(), message(), TTNeutralGenesSH::setFst_li(), Metapop::size(), and warning().

Referenced by TProtoNeutralGenes::loadFileServices().

1060 {
1061  if(_pop->size(ADULTS) == 0) {
1062  warning("No adults in pop, not writing the Fst distribution to file.\n");
1063  return;
1064  }
1065 
1066  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1067  unsigned int patchNbr = _pop->getPatchNbr();
1068  double **fst_i;
1069  bool added_stater = false;
1070 
1072 
1073  if(stater == NULL) {
1074  stater = new TTNeutralGenesSH(_FHLinkedTrait);
1076  added_stater = true;
1077  }
1078 
1079  fst_i = new double* [patchNbr];
1080  for(unsigned int i = 0; i < patchNbr; i++)
1081  fst_i[i] = new double [nb_locus];
1082 
1083  stater->setFst_li(patchNbr, nb_locus, fst_i);
1084 
1085  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
1086 
1087 #ifdef _DEBUG_
1088 message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
1089 #endif
1090 
1091 ofstream FILE (filename.c_str(), ios::out);
1092 
1093 if(!FILE) fatal("could not open Fst_i output file!!\n");
1094 
1095 for(unsigned int i = 0; i < patchNbr; ++i)
1096  FILE<<"patch"<<i+1<<" ";
1097 
1098 FILE<<endl;
1099 
1100 for(unsigned int j = 0; j < nb_locus; ++j) {
1101 
1102  for(unsigned int i = 0; i < patchNbr; i++)
1103  FILE<<fst_i[i][j]<<" ";
1104 
1105  FILE<<endl;
1106 }
1107 
1108 FILE.close();
1109 
1110 if(added_stater) delete stater;
1111 
1112 for(unsigned int i = 0; i < patchNbr; i++)
1113  delete [] fst_i[i];
1114 
1115 delete [] fst_i;
1116 
1117 }
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
#define ADULTS
Adults age class flag (breeders).
Definition: types.h:54
void allocateTables(unsigned int loci, unsigned int all)
Definition: stats_fstat.cc:43
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
unsigned int getPatchNbr()
Definition: metapop.h:270
void warning(const char *str,...)
Definition: output.cc:56
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
void setFst_li(unsigned int N, unsigned int L, double **array)
Computes the per-locus per-patch Fst values using Weir&Hill 2002 approach.
Definition: stats_fstat.cc:1057
unsigned int size()
Get the total number of individuals present in the population, all sex and age classes together...
Definition: metapop.h:302
std::string & get_path()
Definition: filehandler.h:139
The stat handler for neutral markers.
Definition: ttneutralgenes.h:268
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
TTNeutralGenesSH * get_stater()
Definition: ttneutralgenes.h:191
void TTNeutralGenesFH::write_FSTAT ( )

The file format is FSTAT-like, with extra info about the individual added (age, sex, pedigree, origin). The file extension is ".dat".

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FileHandler::_pop, ADLTx, fatal(), FEM, Patch::get(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), FileHandler::get_service(), FileServices::getGenerationReplicateFileName(), Individual::getHome(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getPedigreeClass(), Individual::getTrait(), MAL, message(), OFFSx, and Patch::size().

Referenced by TProtoNeutralGenes::loadFileServices().

716 {
719  unsigned int position;
720  unsigned int ploidy = _FHLinkedTrait->get_ploidy();
721  unsigned int nb_all = _FHLinkedTrait->get_allele_num();
722  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
723  unsigned int patchNbr = _pop->getPatchNbr();
724  unsigned char** seq;
725  Patch* current_patch;
726  Individual *ind;
727 
728  position = nb_all > 99 ? 3 : 2; //assumes nb_all not sup. to 999
729 
730  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
731 
732 #ifdef _DEBUG_
733  message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
734 #endif
735 
736  ofstream FILE (filename.c_str(), ios::out);
737 
738  if(!FILE) fatal("could not open FSTAT output file!!\n");
739 
740  FILE<<patchNbr<<" "<< nb_locus + 4 <<" "<<nb_all<<" "<<position<<"\n";
741 
742  for (unsigned int i = 0; i < nb_locus; ++i)
743  FILE<<"loc"<<i+1<<"\n";
744  //add names for the three last fields:
745  FILE<<"age\n"<<"sex\n"<<"ped\n"<<"origin\n";
746 
747  for (unsigned int i = 0; i < patchNbr; ++i) {
748 
749  current_patch = _pop->getPatch(i);
750 
751  for (unsigned int j = 0; j < current_patch->size(FEM, OFFSx); ++j) {
752 
753  FILE<<i+1<<" ";
754  ind = current_patch->get(FEM, OFFSx, j);
755  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
756 
757  for(unsigned int k = 0; k < nb_locus; ++k) {
758  for (unsigned int l = 0; l < ploidy; ++l) {
759  FILE.fill('0');
760  FILE.width(position);
761  FILE<<(unsigned int)(seq[l][k]+1);
762  }
763  FILE<<" ";
764  }
765  FILE << OFFSx << " "<< FEM << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
766  }
767 
768  for (unsigned int j = 0; j < current_patch->size(MAL, OFFSx); ++j) {
769 
770  FILE<<i+1<<" ";
771  ind = current_patch->get(MAL, OFFSx, j);
772  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
773 
774  for(unsigned int k = 0; k < nb_locus; ++k) {
775  for (unsigned int l = 0; l < ploidy; ++l) {
776  FILE.fill('0');
777  FILE.width(position);
778  FILE<<(unsigned int)(seq[l][k]+1);
779  }
780  FILE<<" ";
781  }
782  FILE << OFFSx << " "<< MAL << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
783  }
784 
785  for (unsigned int j = 0; j < current_patch->size(FEM, ADLTx); ++j) {
786 
787  FILE<<i+1<<" ";
788  ind = current_patch->get(FEM, ADLTx, j);
789  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
790 
791  for(unsigned int k = 0; k < nb_locus; ++k) {
792  for (unsigned int l = 0; l < ploidy; ++l) {
793  FILE.fill('0');
794  FILE.width(position);
795  FILE<<(unsigned int)(seq[l][k]+1);
796  }
797  FILE<<" ";
798  }
799  FILE << ADLTx << " "<< FEM << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
800  }
801 
802  for (unsigned int j = 0; j < current_patch->size(MAL, ADLTx); ++j) {
803 
804  FILE<<i+1<<" ";
805  ind = current_patch->get(MAL, ADLTx, j);
806  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
807 
808  for(unsigned int k = 0; k < nb_locus; ++k) {
809  for (unsigned int l = 0; l < ploidy; ++l) {
810  FILE.fill('0');
811  FILE.width(position);
812  FILE<<(unsigned int)(seq[l][k]+1);
813  }
814  FILE<<" ";
815  }
816  FILE << ADLTx << " "<< MAL << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
817  }
818 
819  }
820 
821  FILE.close();
822 
823 }
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
unsigned int size(age_t AGE)
Returns the size of the container of the appropriate age class(es) for both sexes.
Definition: metapop.h:487
TTrait * getTrait(IDX T)
Trait accessor.
Definition: individual.h:277
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
unsigned int get_ploidy()
Definition: ttneutralgenes.h:187
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
Second class in the metapopulation design structure, between the Metapop and Individual classes...
Definition: metapop.h:421
Patch * getPatch(unsigned int i)
Patch accessor, return the ith+1 patch in the metapop.
Definition: metapop.h:251
unsigned int getPatchNbr()
Definition: metapop.h:270
Definition: types.h:42
virtual void ** get_sequence() const =0
sequence accessor.
unsigned int getPedigreeClass()
Returns the pedigree class of the individual, as set during offspring creation.
Definition: individual.h:179
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
unsigned short getHome()
Definition: individual.h:128
Definition: types.h:37
Definition: types.h:42
Individual * get(sex_t SEX, age_idx AGE, unsigned int at)
Returns a pointer to the individual sitting at the index passed.
Definition: metapop.h:523
std::string & get_path()
Definition: filehandler.h:139
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
This class contains traits along with other individual information (sex, pedigree, etc. ).
Definition: individual.h:49
Definition: types.h:37
int _FHLinkedTraitIndex
Definition: filehandler.h:220
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
void TTNeutralGenesFH::write_GENEPOP ( )

The file format is GENEPOP-like, with extra info about the individual added (age, sex, pedigree, origin). The file extension is ".dat".

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FileHandler::_pop, ADLTx, fatal(), FEM, Patch::get(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), FileHandler::get_service(), FileServices::getGenerationReplicateFileName(), Individual::getHome(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getPedigreeClass(), Individual::getTrait(), MAL, message(), OFFSx, and Patch::size().

Referenced by TProtoNeutralGenes::loadFileServices().

828 {
831  unsigned int position;
832  unsigned int ploidy = _FHLinkedTrait->get_ploidy();
833  unsigned int nb_all = _FHLinkedTrait->get_allele_num();
834  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
835  unsigned int patchNbr = _pop->getPatchNbr();
836  unsigned char** seq;
837  Patch* current_patch;
838  Individual *ind;
839 
840  position = nb_all > 99 ? 3 : 2; //assumes nb_all not sup. to 999
841 
842  if(ploidy == 1) position = 1;
843 
844  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
845 
846 #ifdef _DEBUG_
847  message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
848 #endif
849 
850  ofstream FILE (filename.c_str(), ios::out);
851 
852  if(!FILE) fatal("could not open FSTAT output file!!\n");
853 
854  FILE<<"Title line: "<<patchNbr<<" patches, "<<nb_locus<<" loci with "<<nb_all<<" alleles\n";
855 
856  for (unsigned int i = 0; i < nb_locus; ++i)
857  FILE<<"loc"<<i+1<<", ";
858 
859  //add names for the three last fields:
860  FILE<<"age, sex, ped, origin\n";
861 
862  for (unsigned int i = 0; i < patchNbr; ++i) {
863 
864  current_patch = _pop->getPatch(i);
865 
866  FILE<<"POP\n";
867 
868  for (unsigned int j = 0; j < current_patch->size(FEM, OFFSx); ++j) {
869 
870  FILE<<i+1<<", ";
871  ind = current_patch->get(FEM, OFFSx, j);
872  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
873 
874  for(unsigned int k = 0; k < nb_locus; ++k) {
875  for (unsigned int l = 0; l < ploidy; ++l) {
876  FILE.fill('0');
877  FILE.width(position);
878  FILE<<(unsigned int)(seq[l][k]+1);
879  }
880  FILE<<" ";
881  }
882  FILE << OFFSx << " "<< FEM << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
883  }
884 
885  for (unsigned int j = 0; j < current_patch->size(MAL, OFFSx); ++j) {
886 
887  FILE<<i+1<<", ";
888  ind = current_patch->get(MAL, OFFSx, j);
889  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
890 
891  for(unsigned int k = 0; k < nb_locus; ++k) {
892  for (unsigned int l = 0; l < ploidy; ++l) {
893  FILE.fill('0');
894  FILE.width(position);
895  FILE<<(unsigned int)(seq[l][k]+1);
896  }
897  FILE<<" ";
898  }
899  FILE << OFFSx << " "<< MAL << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
900  }
901 
902  for (unsigned int j = 0; j < current_patch->size(FEM, ADLTx); ++j) {
903 
904  FILE<<i+1<<", ";
905  ind = current_patch->get(FEM, ADLTx, j);
906  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
907 
908  for(unsigned int k = 0; k < nb_locus; ++k) {
909  for (unsigned int l = 0; l < ploidy; ++l) {
910  FILE.fill('0');
911  FILE.width(position);
912  FILE<<(unsigned int)(seq[l][k]+1);
913  }
914  FILE<<" ";
915  }
916  FILE << ADLTx << " "<< FEM << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
917  }
918 
919  for (unsigned int j = 0; j < current_patch->size(MAL, ADLTx); ++j) {
920 
921  FILE<<i+1<<", ";
922  ind = current_patch->get(MAL, ADLTx, j);
923  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
924 
925  for(unsigned int k = 0; k < nb_locus; ++k) {
926  for (unsigned int l = 0; l < ploidy; ++l) {
927  FILE.fill('0');
928  FILE.width(position);
929  FILE<<(unsigned int)(seq[l][k]+1);
930  }
931  FILE<<" ";
932  }
933  FILE << ADLTx << " "<< MAL << " " << ind->getPedigreeClass() << " " << ind->getHome()+1<<endl;
934  }
935 
936  }
937 
938  FILE.close();
939 
940 }
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
unsigned int size(age_t AGE)
Returns the size of the container of the appropriate age class(es) for both sexes.
Definition: metapop.h:487
TTrait * getTrait(IDX T)
Trait accessor.
Definition: individual.h:277
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
unsigned int get_ploidy()
Definition: ttneutralgenes.h:187
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
Second class in the metapopulation design structure, between the Metapop and Individual classes...
Definition: metapop.h:421
Patch * getPatch(unsigned int i)
Patch accessor, return the ith+1 patch in the metapop.
Definition: metapop.h:251
unsigned int getPatchNbr()
Definition: metapop.h:270
Definition: types.h:42
virtual void ** get_sequence() const =0
sequence accessor.
unsigned int getPedigreeClass()
Returns the pedigree class of the individual, as set during offspring creation.
Definition: individual.h:179
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
unsigned short getHome()
Definition: individual.h:128
Definition: types.h:37
Definition: types.h:42
Individual * get(sex_t SEX, age_idx AGE, unsigned int at)
Returns a pointer to the individual sitting at the index passed.
Definition: metapop.h:523
std::string & get_path()
Definition: filehandler.h:139
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
This class contains traits along with other individual information (sex, pedigree, etc. ).
Definition: individual.h:49
Definition: types.h:37
int _FHLinkedTraitIndex
Definition: filehandler.h:220
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
void TTNeutralGenesFH::write_TAB ( )

The file format is FSTAT-like, with age class and sex added after the pop id. The file extension is ".dat".

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FileHandler::_pop, ADLTx, fatal(), FEM, Patch::get(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), FileHandler::get_service(), Individual::getFather(), Individual::getFatherID(), FileServices::getGenerationReplicateFileName(), Individual::getHome(), Individual::getID(), Individual::getMother(), Individual::getMotherID(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getPedigreeClass(), Individual::getSex(), Individual::getTrait(), MAL, message(), OFFSx, and Patch::size().

Referenced by TProtoNeutralGenes::loadFileServices().

605 {
608  unsigned int ploidy = _FHLinkedTrait->get_ploidy();
609  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
610  unsigned int patchNbr = _pop->getPatchNbr();
611  unsigned char** seq;
612  Patch* current_patch;
613  Individual *ind;
614 
615  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
616 
617 #ifdef _DEBUG_
618  message("TTNeutralGenesFH::write_TAB (%s)\n",filename.c_str());
619 #endif
620 
621  ofstream FILE (filename.c_str(), ios::out);
622 
623  if(!FILE) fatal("could not open TAB output file!!\n");
624 
625  FILE<<"pop";
626 
627  for (unsigned int i = 0; i < nb_locus; ++i)
628  for (unsigned int j = 0; j < ploidy; ++j)
629  FILE<<" l"<<i+1<<"a"<<j+1;
630  //add names for the three last fields:
631  FILE<<" age sex home ped isMigrant father mother ID\n";
632 
633  for (unsigned int i = 0; i < patchNbr; ++i) {
634 
635  current_patch = _pop->getPatch(i);
636 
637  for (unsigned int j = 0; j < current_patch->size(FEM, OFFSx); ++j) {
638 
639  FILE<<i+1<<" ";
640  ind = current_patch->get(FEM, OFFSx, j);
641  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
642 
643  for(unsigned int k = 0; k < nb_locus; ++k) {
644  for (unsigned int l = 0; l < ploidy; ++l) {
645  FILE<<(unsigned int)(seq[l][k]+1)<<" ";
646  }
647  }
648  FILE << OFFSx <<" "<<ind->getSex()<<" "<<ind->getHome()+1<<" "<<ind->getPedigreeClass()<<" "
649  << (ind->getFather() && ind->getMother() ?
650  (ind->getFather()->getHome()!=i) + (ind->getMother()->getHome()!=i) : 0)
651  <<" "<<ind->getFatherID()<<" "<<ind->getMotherID()<<" "<<ind->getID()<<std::endl;
652 
653  }
654 
655  for (unsigned int j = 0; j < current_patch->size(MAL, OFFSx); ++j) {
656 
657  FILE<<i+1<<" ";
658  ind = current_patch->get(MAL, OFFSx, j);
659  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
660 
661  for(unsigned int k = 0; k < nb_locus; ++k) {
662  for (unsigned int l = 0; l < ploidy; ++l) {
663  FILE<<(unsigned int)(seq[l][k]+1)<<" ";
664  }
665  }
666  FILE << OFFSx <<" "<<ind->getSex()<<" "<<ind->getHome()+1<<" "<<ind->getPedigreeClass()<<" "
667  << (ind->getFather() && ind->getMother() ?
668  (ind->getFather()->getHome() != i) + (ind->getMother()->getHome() != i) : 0)
669  <<" "<<ind->getFatherID()<<" "<<ind->getMotherID()<<" "<<ind->getID()<<std::endl;
670 
671  }
672 
673  for (unsigned int j = 0; j < current_patch->size(FEM, ADLTx); ++j) {
674 
675  FILE<<i+1<<" ";
676  ind = current_patch->get(FEM, ADLTx, j);
677  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
678 
679  for(unsigned int k = 0; k < nb_locus; ++k) {
680  for (unsigned int l = 0; l < ploidy; ++l) {
681  FILE<<(unsigned int)(seq[l][k]+1)<<" ";
682  }
683  }
684  FILE << ADLTx <<" "<<ind->getSex()<<" "<<ind->getHome()+1<<" "<<ind->getPedigreeClass()<<" "
685  << (ind->getFather() && ind->getMother() ?
686  (ind->getFather()->getHome() != i) + (ind->getMother()->getHome() != i) : 0)
687  <<" "<<ind->getFatherID()<<" "<<ind->getMotherID()<<" "<<ind->getID()<<std::endl;
688  }
689 
690  for (unsigned int j = 0; j < current_patch->size(MAL, ADLTx); ++j) {
691 
692  FILE<<i+1<<" ";
693  ind = current_patch->get(MAL, ADLTx, j);
694  seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
695 
696  for(unsigned int k = 0; k < nb_locus; ++k) {
697  for (unsigned int l = 0; l < ploidy; ++l) {
698  FILE<<(unsigned int)(seq[l][k]+1)<<" ";
699  }
700  }
701  FILE << ADLTx <<" "<<ind->getSex()<<" "<<ind->getHome()+1<<" "<<ind->getPedigreeClass()<<" "
702  << (ind->getFather() && ind->getMother() ?
703  (ind->getFather()->getHome() != i) + (ind->getMother()->getHome() != i) : 0)
704  <<" "<<ind->getFatherID()<<" "<<ind->getMotherID()<<" "<<ind->getID()<<std::endl;
705  }
706 
707  }
708 
709  FILE.close();
710 
711 }
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
unsigned long getID()
Definition: individual.h:122
unsigned int size(age_t AGE)
Returns the size of the container of the appropriate age class(es) for both sexes.
Definition: metapop.h:487
TTrait * getTrait(IDX T)
Trait accessor.
Definition: individual.h:277
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
Individual * getFather()
Definition: individual.h:126
unsigned long getMotherID()
Definition: individual.h:125
unsigned int get_ploidy()
Definition: ttneutralgenes.h:187
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
Individual * getMother()
Definition: individual.h:127
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
Second class in the metapopulation design structure, between the Metapop and Individual classes...
Definition: metapop.h:421
Patch * getPatch(unsigned int i)
Patch accessor, return the ith+1 patch in the metapop.
Definition: metapop.h:251
unsigned long getFatherID()
Definition: individual.h:124
unsigned int getPatchNbr()
Definition: metapop.h:270
Definition: types.h:42
virtual void ** get_sequence() const =0
sequence accessor.
unsigned int getPedigreeClass()
Returns the pedigree class of the individual, as set during offspring creation.
Definition: individual.h:179
unsigned short getHome()
Definition: individual.h:128
Definition: types.h:37
Definition: types.h:42
Individual * get(sex_t SEX, age_idx AGE, unsigned int at)
Returns a pointer to the individual sitting at the index passed.
Definition: metapop.h:523
std::string & get_path()
Definition: filehandler.h:139
sex_t getSex()
Definition: individual.h:129
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
This class contains traits along with other individual information (sex, pedigree, etc. ).
Definition: individual.h:49
Definition: types.h:37
int _FHLinkedTraitIndex
Definition: filehandler.h:220
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
void TTNeutralGenesFH::write_varcompWC ( )

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, _output_option, FileHandler::_pop, ADLTx, ADULTS, ALL, TTNeutralGenesSH::allocateTables(), fatal(), TMatrix::get(), DataTable< T >::get(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), FileHandler::get_service(), TProtoNeutralGenes::get_stater(), TTNeutralGenesSH::getAlleleCountTable(), TTNeutralGenesSH::getAlleleFreqTable(), Metapop::getCurrentAge(), FileServices::getGenerationReplicateFileName(), TTNeutralGenesSH::getGlobalFreqs(), TTNeutralGenesSH::getHeteroTable(), Metapop::getPatchNbr(), message(), OFFSPRG, OFFSx, TTNeutralGenesSH::setAlleleTables(), TTNeutralGenesSH::setHeteroTable(), TTNeutralGenesSH::setHo2(), Metapop::size(), and warning().

Referenced by TProtoNeutralGenes::loadFileServices().

1208 {
1209  unsigned int nb_allele = _FHLinkedTrait->get_allele_num();
1210  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1211  unsigned int patchNbr = _pop->getPatchNbr();
1212  bool added_stater = false;
1213  age_t AGE = _pop->getCurrentAge();
1214 
1215  if(AGE == (ADULTS | OFFSPRG) || AGE == ALL) {
1216  warning("saving only offspring neutral allele stats\n");
1217  AGE = OFFSPRG;
1218  }
1219 
1220  age_idx age = (AGE == OFFSPRG ? OFFSx : ADLTx);
1221 
1223 
1224  if(stater == NULL) {
1225  stater = new TTNeutralGenesSH(_FHLinkedTrait);
1227  added_stater = true;
1228  }
1229 
1230  stater->setAlleleTables(AGE);
1231  stater->setHeteroTable(AGE);
1232 
1233  TMatrix *globalFreq = stater->getGlobalFreqs();
1234  DataTable< double > *freqTable = stater->getAlleleFreqTable();
1235  DataTable< double > *heteroTable = stater->getHeteroTable();
1236  DataTable< unsigned int > *alleleCountTable = stater->getAlleleCountTable();
1237 
1238  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
1239 
1240 #ifdef _DEBUG_
1241  message("TTNeutralGenesFH::write_varcompWC (%s)\n",filename.c_str());
1242 #endif
1243 
1244  //init
1245  double* pop_sizes = new double [patchNbr];
1246  double tot_size;
1247  double sum_weights = 0;
1248  double nc;
1249  unsigned int extantPs = 0;
1250 
1251  tot_size = _pop->size(AGE);
1252 
1253  for(unsigned int i = 0; i < patchNbr; i++) {
1254  pop_sizes[i] = _pop->size(AGE, i);
1255  if(pop_sizes[i]) {
1256  extantPs++;
1257  sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
1258  }
1259  }
1260 
1261  nc = (tot_size - sum_weights)/(extantPs-1);
1262 
1263 // unsigned int np = extantPs;
1264  unsigned int npl = extantPs; //all loci typed in all patches
1265 
1266  //p = _alleleFreqTable
1267  //pb = _globalAlleleFreq
1268 
1269  //counting num alleles per locus
1270 
1271  unsigned int *alploc = new unsigned int [nb_locus];
1272 
1273  unsigned int **alploc_table = new unsigned int* [nb_locus];
1274 
1275  for(unsigned int i = 0; i < nb_locus; ++i)
1276  alploc_table[i] = new unsigned int[nb_allele];
1277 
1278  unsigned int tot_num_allele = 0;
1279 
1280  for(unsigned int l = 0; l < nb_locus; ++l){
1281 
1282  alploc[l] = 0;
1283 
1284  for(unsigned int cnt, a = 0; a < nb_allele; ++a) {
1285 
1286  cnt=0;
1287 
1288  for(unsigned int i = 0; i < patchNbr; i++) {
1289 
1290  cnt += alleleCountTable->get(i,l,a);
1291 
1292  }
1293  alploc_table[l][a] = (cnt != 0);
1294  alploc[l] += (cnt != 0);
1295  }
1296  tot_num_allele += alploc[l];
1297  }
1298 
1299  //correspondance with hierfstat implementation:
1300  //n, and nal are given by pop_sizes, same num ind typed at all loci in each patch
1301  //nc is the same for each locus
1302  //nt is given by tot_size, same tot num of ind typed at all loci
1303 
1304  //SSG = het/2 for each allele
1305  double *SSG = new double[tot_num_allele];
1306  double *SSP = new double[tot_num_allele];
1307  double *SSi = new double[tot_num_allele];
1308  double *loc_id = new double[tot_num_allele];
1309  double *al_id = new double [tot_num_allele];
1310 
1311  unsigned int all_cntr = 0;
1312 
1313  double het, freq, var;
1314 
1315  //computing sum of squares
1316 
1317  for(unsigned int l = 0; l < nb_locus; ++l) {
1318 
1319  for(unsigned int a = 0; a < nb_allele & all_cntr < tot_num_allele; ++a) {
1320 
1321  if(alploc_table[l][a] == 0) continue; //do not consider alleles not present in the whole pop
1322 
1323  //store locus and all identifiers for output
1324  loc_id[all_cntr] = l+1;
1325  al_id[all_cntr] = a+1;
1326 
1327  SSG[all_cntr] = 0;
1328  SSi[all_cntr] = 0;
1329  SSP[all_cntr] = 0;
1330 
1331  for(unsigned int p = 0; p < patchNbr; ++p){
1332 
1333  if(!_pop->size(AGE, p)) continue; //skip empty patches
1334 
1335  het = heteroTable->get(p, l, a);
1336 
1337  freq = freqTable->get(p, l, a);
1338 
1339  var = freq - globalFreq->get(l, a); //(p_liu - pbar_u)^2
1340 
1341  var *= var;
1342 
1343  SSG[all_cntr] += het;
1344 
1345  SSi[all_cntr] += 2*pop_sizes[p]*freq*(1-freq) - het/2;
1346 
1347  SSP[all_cntr] += 2*pop_sizes[p]*var;
1348  }
1349 
1350  all_cntr++;
1351  }
1352 
1353  }
1354 
1355  assert(all_cntr == tot_num_allele);
1356 
1357  // open the file
1358 
1359  ofstream FILE (filename.c_str(), ios::out);
1360 
1361  if(!FILE) fatal("could not open neutral vcomp output file!!\n");
1362 
1363  // print column names
1364 
1365  if(_output_option == "allele") {
1366 
1367  FILE<<"locus allele pbar het siga sigb sigw Fst Fis";
1368 
1369  for (unsigned int p = 1; p <= patchNbr; ++p)
1370  FILE<<" het.p"<<p;
1371 
1372  for (unsigned int p = 1; p <= patchNbr; ++p)
1373  FILE<<" freq.p"<<p;
1374 
1375  FILE<<endl;
1376 
1377  } else {
1378 
1379  FILE<<"locus maj.al pbar.maj.al het siga sigb sigw Fst Fis";
1380 
1381  for (unsigned int p = 1; p <= patchNbr; ++p)
1382  FILE<<" het.p"<<p;
1383 
1384  //allele frequencies in each patch
1385  for (unsigned int p = 1; p <= patchNbr; ++p)
1386 // for(unsigned int u = 0; u < nb_allele-1; ++u) //skip last allele, can be deduced...
1387  FILE<<" freq.maj."<<"p"<< p;
1388 
1389  FILE<<endl;
1390  }
1391 
1392  //-----------------------------------------------------------------------------------------
1393  // allele specific stats:
1394 
1395  double *MSG = new double[tot_num_allele];
1396  double *MSP = new double[tot_num_allele];
1397  double *MSI = new double[tot_num_allele];
1398  // double *sigw = new double[tot_num_allele];
1399  double *siga = new double[tot_num_allele];
1400  double *sigb = new double[tot_num_allele];
1401 
1402  for(unsigned int i = 0; i < tot_num_allele; ++i){
1403 
1404  MSG[i] = SSG[i] / (2 * tot_size);
1405  // sigw[i] = MSG[i]; //wasted!
1406 
1407  MSP[i] = SSP[i] / (npl-1);
1408 
1409  MSI[i] = SSi[i]/ (tot_size - npl);
1410 
1411  sigb[i] = 0.5*(MSI[i] - MSG[i]);
1412 
1413  siga[i] = (MSP[i] - MSI[i])/(2*nc);
1414 
1415  if(_output_option == "allele") {
1416  FILE<< loc_id[i] << " " << al_id[i] << " "
1417  //global allele frequency
1418  << globalFreq->get(loc_id[i] - 1, al_id[i] - 1) << " "
1419  //average heterozygosity:
1420  << SSG[i]/tot_size << " "
1421  //variance components:
1422  << siga[i] << " " << sigb[i] << " " << MSG[i] << " "
1423  //Fst
1424  << siga[i]/(siga[i]+sigb[i]+MSG[i]) << " "
1425  //Fis
1426  << sigb[i]/(sigb[i]+MSG[i]) << " ";
1427  //per patch heterozygosity:
1428  for (unsigned int p = 0; p < patchNbr; ++p) {
1429  FILE << heteroTable->get(p, loc_id[i]-1, al_id[i]-1)/pop_sizes[p] << " ";
1430  }
1431  //per patch allele frequency:
1432  for (unsigned int p = 0; p < patchNbr; ++p) {
1433  FILE << freqTable->get(p, loc_id[i]-1, al_id[i]-1) << " ";
1434  }
1435 
1436  FILE<< endl;
1437  }
1438  }
1439 
1440  //-----------------------------------------------------------------------------------------
1441  // locus specific stats:
1442 
1443  if(_output_option == "locus") {
1444 
1445 
1446  double lsiga, lsigb, lsigw, max_all_frq, het;
1447  unsigned int maj_al;
1448 
1449  all_cntr = 0;
1450 
1451  deque <double> loc_het = stater->setHo2(age);
1452 
1453  for(unsigned int i = 0; i < nb_locus; ++i) {
1454 
1455  lsiga = 0;
1456  lsigb = 0;
1457  lsigw = 0;
1458 
1459  max_all_frq = 0;
1460 
1461  for(unsigned int l = 0; l < alploc[i]; ++l) {
1462 
1463  lsiga += siga[all_cntr];
1464  lsigb += sigb[all_cntr];
1465  lsigw += MSG[all_cntr];
1466 
1467  if(max_all_frq < globalFreq->get(i, al_id[all_cntr]-1 ) ) {
1468  max_all_frq = globalFreq->get(i, al_id[all_cntr]-1 );
1469  maj_al = al_id[all_cntr];
1470  }
1471 
1472  all_cntr++;
1473 
1474  }
1475 
1476  FILE << i+1 <<" "<< maj_al <<" "<< max_all_frq <<" ";
1477  FILE<< loc_het[i] <<" "<< lsiga << " " << lsigb <<" ";
1478  FILE<< lsigw << " "<< lsiga /(lsiga + lsigb + lsigw) <<" ";
1479  FILE<< lsigb /(lsigb + lsigw);
1480  // patch heterozygosities:
1481 
1482  for(unsigned int p = 0; p < patchNbr; ++p) {
1483  het = 0;
1484  for(unsigned int a = 0; a < nb_allele; ++a){
1485  het += heteroTable->get(p, i, a);
1486  }
1487  FILE << " " << het/(2.0*pop_sizes[p]);
1488  }
1489  // patch allele freq?
1490  for(unsigned int p = 0; p < patchNbr; ++p) {
1491  FILE<< " " << freqTable->get(p, i, maj_al-1);
1492  }
1493  FILE << endl;
1494  }//END for locus
1495 
1496  } //END per locus output
1497 
1498 
1499 
1500 
1501  FILE.close();
1502 
1503  delete[]pop_sizes;
1504  delete[]alploc;
1505  for(unsigned int i = 0; i < nb_locus; ++i)
1506  delete[]alploc_table[i];
1507  delete[]alploc_table;
1508  delete[]loc_id;
1509  delete[]al_id;
1510  delete[]SSG;
1511  delete[]SSi;
1512  delete[]SSP;
1513  delete[]MSG;
1514  delete[]MSI;
1515  delete[]MSP;
1516  // delete[]sigw;
1517  delete[]siga;
1518  delete[]sigb;
1519 }
void setHeteroTable(age_t AGE)
Definition: stats_fstat.cc:166
void message(const char *message,...)
Definition: output.cc:40
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
#define ADULTS
Adults age class flag (breeders).
Definition: types.h:54
void allocateTables(unsigned int loci, unsigned int all)
Definition: stats_fstat.cc:43
age_idx
Array index of the age classes in the patch sizes and containers arrays.
Definition: types.h:41
void setAlleleTables(age_t AGE)
Definition: stats_fstat.cc:74
unsigned int age_t
Age class flags.
Definition: types.h:46
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:218
TMatrix * getGlobalFreqs()
Accessor to the table of allele frequencies in the whole population.
Definition: ttneutralgenes.h:368
double get(unsigned int i, unsigned int j)
Accessor to element at row i and column j.
Definition: tmatrix.h:120
std::string & get_extension()
Definition: filehandler.h:143
void fatal(const char *str,...)
Definition: output.cc:90
unsigned int getPatchNbr()
Definition: metapop.h:270
deque< double > setHo2(age_idx age_pos)
Definition: stats_fstat.cc:545
T get(unsigned int group, unsigned int Class, unsigned int elmnt)
Returns value stored of the element 'elmnt' of the class 'Class' in the group 'group'.
Definition: datatable.h:228
void warning(const char *str,...)
Definition: output.cc:56
Definition: types.h:42
string _output_option
Definition: ttneutralgenes.h:233
unsigned int get_allele_num()
Definition: ttneutralgenes.h:189
A class to handle matrix in params, coerces matrix into a vector of same total size.
Definition: tmatrix.h:49
#define ALL
All ages age class flag.
Definition: types.h:56
DataTable< unsigned int > * getAlleleCountTable()
Definition: ttneutralgenes.h:363
Definition: types.h:42
#define OFFSPRG
Offspring age class flag.
Definition: types.h:50
unsigned int size()
Get the total number of individuals present in the population, all sex and age classes together...
Definition: metapop.h:302
std::string & get_path()
Definition: filehandler.h:139
DataTable< double > * getAlleleFreqTable()
Accessor to the table of allele frequencies, per patch.
Definition: ttneutralgenes.h:361
DataTable< double > * getHeteroTable()
Definition: ttneutralgenes.h:365
The stat handler for neutral markers.
Definition: ttneutralgenes.h:268
age_t getCurrentAge()
Definition: metapop.h:289
unsigned int get_locus_num()
Definition: ttneutralgenes.h:188
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
TTNeutralGenesSH * get_stater()
Definition: ttneutralgenes.h:191

Member Data Documentation

string TTNeutralGenesFH::_output_option
private

Referenced by setOutputOption(), and write_varcompWC().

void(TTNeutralGenesFH::* TTNeutralGenesFH::write_fct)()
private

Referenced by FHwrite(), and set_write_fct().


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

Generated for Nemo v2.3.0 by  doxygen 1.8.8 --
Catalogued on GSR