Nemo  2.3.56
Simulate forward-in-time genetic evolution in a spatially explicit, individual-based stochastic simulator
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_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)())
 
- Public Member Functions inherited from TraitFileHandler< TProtoNeutralGenes >
 TraitFileHandler (TProtoNeutralGenes *trait_proto, const char *ext)
 
virtual ~TraitFileHandler ()
 
virtual void FHwrite ()=0
 
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)
 
- 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 FHwrite ()=0
 Default behavior of the class, called by Handler::update(). More...
 
virtual void FHread (string &filename)=0
 Default input function. More...
 
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 void init ()=0
 Inits state. More...
 
virtual void update ()=0
 Updates the handler state. More...
 
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::TTNeutralGenesFH ( TProtoNeutralGenes TP)
inline
244 { }
void(TTNeutralGenesFH::* write_fct)()
Definition: ttneutralgenes.h:238
Template class for the trait's FileHandler.
Definition: filehandler.h:217

◆ ~TTNeutralGenesFH()

virtual TTNeutralGenesFH::~TTNeutralGenesFH ( )
inlinevirtual
246{ }

Member Function Documentation

◆ FHread()

void TTNeutralGenesFH::FHread ( string &  filename)
virtual

reads in only FSTAT format

Implements FileHandler.

1200{
1202 // make sure we are using the whole pop, not a sub-sampled one:
1203 Metapop *pop = get_service()->get_pop_ptr();
1204
1205 unsigned int digit, nloci_infile, nall_infile, npatch_infile;
1206 // unsigned int ploidy = _FHLinkedTrait->get_ploidy();
1207 unsigned int nb_all = _FHLinkedTrait->get_allele_num();
1208 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1209 unsigned int patchNbr = pop->getPatchNbr();
1210
1211 bool is_extended = false;
1212
1213 ifstream FILE(filename.c_str(),ios::in);
1214
1215 if(!FILE) fatal("could not open FSTAT input file \"%s\"\n", filename.c_str());
1216
1217 FILE >> npatch_infile >> nloci_infile >> nall_infile >> digit;
1218
1219 if(npatch_infile != patchNbr) fatal("# of patch in FSTAT file differs from simulation settings\n");
1220
1221 if(nloci_infile != nb_locus && nloci_infile != nb_locus +4) {
1222
1223 fatal("# of loci in FSTAT file (%i loci) differs from simulation settings (%i loci)\n",nloci_infile, nb_locus);
1224
1225 } else if (nloci_infile == nb_locus + 4) {
1226 is_extended = true;
1227 } else if (nloci_infile == nb_locus){
1228 is_extended = false;
1229 } else {
1230
1231 error("when reading the first line of the input FSTAT file \"%s\" \n", filename.c_str());
1232 fatal("please specify either (num locus) or (num locus + 5) as the second number on the first line of the input file.\n");
1233
1234 }
1235
1236 if(nall_infile != nb_all) fatal("# of alleles in FSTAT file differs from simulation settings\n");
1237
1238 digit = (int) pow(10.0,(double)digit);
1239 string str;
1240
1241 for(unsigned int i = 0; i < nloci_infile; ++i) {
1242 FILE >> str; //we don't do anything with this here
1243 }
1244
1245
1246 unsigned int patch, genot, all0, all1, age, sex, ped, origin;
1247 age_idx agex;
1248 Individual *ind;
1249 unsigned char* seq[2];
1250 int lnbr = nloci_infile +2; //line counter
1251 seq[0] = new unsigned char[nb_locus];
1252 seq[1] = new unsigned char[nb_locus];
1253
1254 while(FILE>>patch) {
1255
1256 if(FILE.bad()) {
1257 error("Reading input genotypes from \"%s\" failed at line %i.\n",filename.c_str(), lnbr);
1258 FILE.clear();
1259 FILE >> str;
1260 fatal("Expecting population number as first element on the line, but received this: %s \n", str.c_str());
1261 }
1262
1263 //check patch value
1264 if(patch > patchNbr) fatal("found an illegal patch identifier in FSTAT file (%i) at line %i",patch, lnbr);
1265
1266 for(unsigned int i = 0; i < nb_locus; ++i) {
1267 FILE>>genot;
1268
1269 all0 = genot/digit;
1270 all1 = genot%digit;
1271
1272 if( all0 > nb_all ) {
1273
1274 error("in FSTAT input file at line %i, locus %i : \
1275first allele value %d is greater than the max value specified (%i)!\n", lnbr, i+1, all0, nb_all);
1276
1277 fatal("Please check the input file.\n");
1278
1279 } else if ( all0 > 0 ) {
1280
1281 seq[0][i] = all0 - 1;
1282
1283 } else {
1284
1285 fatal("in FSTAT input file at line %i, locus %i first allele value: \
1286allele value 0 is not allowed!\n*** Please check the input file.\n", lnbr, i+1, all0, nb_all);
1287
1288 }
1289
1290 if( all1 > nb_all ){
1291
1292 error("in FSTAT input file at line %i, locus %i : \
1293second allele value %i is greater than the max value specified (%i)!\n", lnbr, i+1, all1, nb_all);
1294
1295 fatal("Please check the input file.\n");
1296
1297 } else if (all1 > 0) {
1298
1299 seq[1][i] = all1 - 1;
1300
1301 } else {
1302
1303 fatal("in FSTAT input file at line %i, locus %i second allele value: \
1304allele value 0 is not allowed!\n*** Please check the input file.\n", lnbr, i+1, all0, nb_all);
1305
1306 }
1307 }
1308
1309 if(is_extended) {
1310
1311 FILE >> age >> sex >> ped >> origin;
1312 //age index now saved in the FSTAT file
1313 agex = (static_cast<age_idx> (age) == ADLTx ? ADLTx : OFFSx);
1314
1315 } else {
1316 // by default, individuals are adult females
1317 agex = ADLTx;
1318 sex = FEM;
1319 ped = 0;
1320 origin = 0;
1321
1322 }
1323
1324
1325 ind = pop->makeNewIndividual(0, 0, sex_t(sex), origin - 1);
1326 ind->setPedigreeClass((unsigned char)ped);
1327 ind->getTrait(_FHLinkedTraitIndex)->set_sequence((void**)seq);
1328
1329 pop->getPatch(patch-1)->add(sex_t(sex), agex, ind);
1330
1331 lnbr++;
1332
1333 }
1334
1335 FILE.close();
1336
1337 delete [] seq[0];
1338 delete [] seq[1];
1339}
FileServices * get_service()
Returns pointer to the FileServices.
Definition: filehandler.h:135
virtual Metapop * get_pop_ptr()
Accessor to the pointer to the main population.
Definition: fileservices.h:111
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
This class contains traits along with other individual information (sex, pedigree,...
Definition: individual.h:49
void setPedigreeClass(Individual *mother, Individual *father)
Definition: individual.h:115
TTrait * getTrait(IDX T)
Trait accessor.
Definition: individual.h:277
Top class of the metapopulation structure, contains the patches.
Definition: metapop.h:80
unsigned int getPatchNbr()
Definition: metapop.h:276
Patch * getPatch(unsigned int i)
Patch accessor, return the ith+1 patch in the metapop.
Definition: metapop.h:257
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:549
unsigned int get_locus_num()
Definition: ttneutralgenes.h:190
unsigned int get_allele_num()
Definition: ttneutralgenes.h:191
virtual void set_sequence(void **seq)=0
Called to set the sequence pointer to an existing trait.
int _FHLinkedTraitIndex
Definition: filehandler.h:220
TProtoNeutralGenes * _FHLinkedTrait
Definition: filehandler.h:219
void fatal(const char *str,...)
Definition: output.cc:96
int error(const char *str,...)
Definition: output.cc:77
sex_t
Sex types, males are always 0 and females 1!!
Definition: types.h:36
@ FEM
Definition: types.h:37
age_idx
Array index of the age classes in the patch sizes and containers arrays.
Definition: types.h:41
@ OFFSx
Definition: types.h:42
@ ADLTx
Definition: types.h:42

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

◆ FHwrite()

void TTNeutralGenesFH::FHwrite ( )
virtual

Implements TraitFileHandler< TProtoNeutralGenes >.

715{
716 if(!_pop->isAlive()) return;
717
718 // reset the pop ptr from the file services, will be the main metapop without sub sampling
719 // or a sub sampled pop otherwise:
720
722
723 (this->*write_fct)();
724
725 // reset the pop ptr to the main pop
727
728}
Metapop * _pop
Pointer to the current metapop, set during initialization within the init function.
Definition: filehandler.h:103
Metapop * getSampledPop()
Sets the down-sampled population and provides accessor to file handlers.
Definition: fileservices.cc:197
bool isAlive()
Checks if the population still contains at least one individual in any sex or age class.
Definition: metapop.h:307

References FileHandler::_pop, FileServices::get_pop_ptr(), FileHandler::get_service(), FileServices::getSampledPop(), Metapop::isAlive(), and write_fct.

◆ print_PLINK_PED()

void TTNeutralGenesFH::print_PLINK_PED ( ofstream &  FH,
age_idx  Ax,
Patch patch 
)
906{
907 Individual *ind;
908 double** seq;
909 char BASE[2] = {'A','G'};
910 unsigned int loc;
911 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
912 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
913
914 // SPECIFICATION FOR THE .fam FILE = 6 first values in .ped files:
915// .fam (PLINK sample information file)
916//
917// Sample information file accompanying a .bed binary genotype table.
918// Also generated by "--recode lgen" and "--recode rlist".
919//
920// A text file with no header line, and one line per sample with the following six fields:
921//
922// 1. Family ID ('FID')
923// 2. Within-family ID ('IID'; cannot be '0')
924// 3. Within-family ID of father ('0' if father isn't in dataset)
925// 4. Within-family ID of mother ('0' if mother isn't in dataset)
926// 5. Sex code ('1' = male, '2' = female, '0' = unknown)
927// 6. Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)
928//
929// If there are any numeric phenotype values other than {-9, 0, 1, 2}, the phenotype is interpreted as a quantitative trait instead of case/control status. In this case, -9 normally still designates a missing phenotype; use --missing-phenotype if this is problematic.
930
931
932 for (unsigned int j = 0; j < patch->size(FEM, Ax); ++j) {
933
934 ind = patch->get(FEM, Ax, j);
935
936 FH<<"fam"<<ind->getHome()+1<<" "<<ind->getID()<<" 0 0" //<<ind->getFatherID()<<" "<<ind->getMotherID()
937 <<" 2 -9";
938
939 seq = (double**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
940
941 for(unsigned int k = 0; k < nb_locus; ++k) {
942 for (unsigned int l = 0; l < ploidy; ++l) {
943 FH<< BASE[ (unsigned int)(seq[l][k]) ]<<" ";
944 }
945 }
946
947 FH <<std::endl;
948
949 }
950
951 for (unsigned int j = 0; j < patch->size(MAL, Ax); ++j) {
952
953 ind = patch->get(MAL, Ax, j);
954
955 FH<<"fam"<<ind->getHome()+1<<" 0 0" //<<ind->getID()+1<<ind->getFatherID()<<" "<<ind->getMotherID()
956 <<" 1 -9";
957
958 seq = (double**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
959
960 for(unsigned int k = 0; k < nb_locus; ++k) {
961 for (unsigned int l = 0; l < ploidy; ++l) {
962 FH<<" "<< BASE[ (unsigned int)(seq[l][k]) ];
963 }
964 }
965
966 FH<<std::endl;
967
968 }
969}
unsigned long getID()
Definition: individual.h:122
unsigned short getHome()
Definition: individual.h:128
unsigned int size(age_t AGE)
Returns the size of the container of the appropriate age class(es) for both sexes.
Definition: metapop.h:496
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:532
unsigned int get_ploidy()
Definition: ttneutralgenes.h:189
virtual void ** get_sequence() const =0
sequence accessor.
@ MAL
Definition: types.h:37

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FEM, Patch::get(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), Individual::getHome(), Individual::getID(), Individual::getTrait(), MAL, and Patch::size().

Referenced by write_PLINK().

+ Here is the caller graph for this function:

◆ set_write_fct()

void TTNeutralGenesFH::set_write_fct ( void(TTNeutralGenesFH::*)()  fct_ptr)
inline
266{write_fct = fct_ptr;}

References write_fct.

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ setOutputOption()

void TTNeutralGenesFH::setOutputOption ( string  opt)
704{
705
706 if(opt != "1") //it means that the param received an argument value
707 _output_option = opt;
708 else
709 _output_option = "locus"; //default
710}
string _output_option
Definition: ttneutralgenes.h:237

References _output_option.

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_Fst_i()

void TTNeutralGenesFH::write_Fst_i ( )
1344{
1345 // make sure we are using the whole pop, not a sub-sampled one:
1346 Metapop *pop = get_service()->get_pop_ptr();
1347
1348 if(pop->size(ADULTS) == 0) {
1349 warning("No adults in pop, not writing the Fst distribution to file.\n");
1350 return;
1351 }
1352
1353 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1354 unsigned int patchNbr = pop->getPatchNbr();
1355 double **fst_i;
1356 bool added_stater = false;
1357
1359
1360 if(stater == NULL) {
1361 stater = new TTNeutralGenesSH(_FHLinkedTrait);
1363 added_stater = true;
1364 }
1365
1366 fst_i = new double* [patchNbr];
1367 for(unsigned int i = 0; i < patchNbr; i++)
1368 fst_i[i] = new double [nb_locus];
1369
1370 stater->setFst_li(patchNbr, nb_locus, fst_i);
1371
1372 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
1373
1374#ifdef _DEBUG_
1375 message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
1376#endif
1377
1378 ofstream FILE (filename.c_str(), ios::out);
1379
1380 if(!FILE) fatal("could not open Fst_i output file!!\n");
1381
1382 for(unsigned int i = 0; i < patchNbr; ++i)
1383 FILE<<"patch"<<i+1<<" ";
1384
1385 FILE<<endl;
1386
1387 for(unsigned int j = 0; j < nb_locus; ++j) {
1388
1389 for(unsigned int i = 0; i < patchNbr; i++)
1390 FILE<<fst_i[i][j]<<" ";
1391
1392 FILE<<endl;
1393 }
1394
1395 FILE.close();
1396
1397 if(added_stater) delete stater;
1398
1399 for(unsigned int i = 0; i < patchNbr; i++)
1400 delete [] fst_i[i];
1401
1402 delete [] fst_i;
1403
1404}
std::string & get_extension()
Definition: filehandler.h:143
std::string & get_path()
Definition: filehandler.h:139
string getGenerationReplicateFileName()
Accessor to the current file name with generation and replicate counters added.
Definition: fileservices.cc:390
unsigned int size()
Get the total number of individuals present in the population, all sex and age classes together.
Definition: metapop.h:310
TTNeutralGenesSH * get_stater()
Definition: ttneutralgenes.h:193
The stat handler for neutral markers.
Definition: ttneutralgenes.h:277
void allocateTables(unsigned int loci, unsigned int all)
Definition: stats_fstat.cc:43
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
void warning(const char *str,...)
Definition: output.cc:58
void message(const char *message,...)
Definition: output.cc:40
#define ADULTS
Adults age class flag (breeders).
Definition: types.h:54

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

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_FSTAT()

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".

1033{
1036 unsigned int position;
1037 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
1038 unsigned int nb_all = _FHLinkedTrait->get_allele_num();
1039 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1040 unsigned int patchNbr = _pop->getPatchNbr();
1041 unsigned char** seq;
1042 Patch* current_patch;
1043 Individual *ind;
1044
1045 position = nb_all > 99 ? 3 : nb_all > 9 ? 2 : 1; //assumes nb_all not sup. to 999
1046
1047 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
1048 + get_extension();
1049
1050#ifdef _DEBUG_
1051 message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
1052#endif
1053
1054 ofstream FILE (filename.c_str(), ios::out);
1055
1056 if(!FILE) fatal("could not open FSTAT output file!!\n");
1057
1058 FILE<<patchNbr<<" "<< nb_locus + 4 <<" "<<nb_all<<" "<<position<<"\n";
1059
1060 for (unsigned int i = 0; i < nb_locus; ++i)
1061 FILE<<"loc"<<i+1<<"\n";
1062
1063 //add names for the three last fields:
1064 FILE<<"age\n"<<"sex\n"<<"ped\n"<<"origin\n";
1065
1066 for (unsigned int i = 0; i < patchNbr; ++i) {
1067
1068 current_patch = _pop->getPatch(i);
1069
1070 write_patch_FSTAT(current_patch, FEM, OFFSx, FILE, position);
1071 write_patch_FSTAT(current_patch, MAL, OFFSx, FILE, position);
1072 write_patch_FSTAT(current_patch, FEM, ADLTx, FILE, position);
1073 write_patch_FSTAT(current_patch, MAL, ADLTx, FILE, position);
1074
1075 }
1076
1077 FILE.close();
1078
1079}
Second class in the metapopulation design structure, between the Metapop and Individual classes.
Definition: metapop.h:430
void write_patch_FSTAT(Patch *patch, sex_t SEX, age_idx AGE, ofstream &FH, unsigned int digits)
Definition: ttneutralgenes.cc:1083

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, ADLTx, fatal(), FEM, TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), FileHandler::get_service(), FileServices::getGenerationReplicateFileName(), Metapop::getPatch(), Metapop::getPatchNbr(), MAL, message(), OFFSx, and write_patch_FSTAT().

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_GENEPOP()

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 ".txt".

1114{
1117 unsigned int position;
1118 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
1119 unsigned int nb_all = _FHLinkedTrait->get_allele_num();
1120 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1121 unsigned int patchNbr = _pop->getPatchNbr();
1122 unsigned char** seq;
1123 Patch* current_patch;
1124 Individual *ind;
1125
1126 position = nb_all > 99 ? 3 : 2; //assumes nb_all not sup. to 999
1127
1128 if(ploidy == 1) position = 1;
1129
1130 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
1131 + get_extension();
1132
1133#ifdef _DEBUG_
1134 message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
1135#endif
1136
1137 ofstream FILE (filename.c_str(), ios::out);
1138
1139 if(!FILE) fatal("could not open FSTAT output file!!\n");
1140
1141 FILE<<"Title line: "<<patchNbr<<" patches, "<<nb_locus<<" loci with "<<nb_all<<" alleles\n";
1142
1143 for (unsigned int i = 0; i < nb_locus; ++i)
1144 FILE<<"loc"<<i+1<<", ";
1145
1146 //add names for the three last fields:
1147 FILE<<"age, sex, ped, origin\n";
1148
1149 for (unsigned int i = 0; i < patchNbr; ++i) {
1150
1151 current_patch = _pop->getPatch(i);
1152
1153 FILE<<"POP\n";
1154
1155 write_patch_GENEPOP(current_patch, FEM, OFFSx, FILE, position);
1156 write_patch_GENEPOP(current_patch, MAL, OFFSx, FILE, position);
1157 write_patch_GENEPOP(current_patch, FEM, ADLTx, FILE, position);
1158 write_patch_GENEPOP(current_patch, MAL, ADLTx, FILE, position);
1159
1160
1161 }
1162
1163 FILE.close();
1164
1165}
void write_patch_GENEPOP(Patch *patch, sex_t SEX, age_idx AGE, ofstream &FH, unsigned int digits)
Definition: ttneutralgenes.cc:1169

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, ADLTx, fatal(), FEM, TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), FileHandler::get_service(), FileServices::getGenerationReplicateFileName(), Metapop::getPatch(), Metapop::getPatchNbr(), MAL, message(), OFFSx, and write_patch_GENEPOP().

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_patch_FSTAT()

void TTNeutralGenesFH::write_patch_FSTAT ( Patch patch,
sex_t  SEX,
age_idx  AGE,
ofstream &  FH,
unsigned int  digits 
)
1084{
1085 unsigned char** seq;
1086 Individual *ind;
1087
1088 for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
1089
1090 FH<<patch->getID() + 1<<" ";
1091 ind = patch->get(SEX, AGE, j);
1092 seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
1093
1094 for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
1095 for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
1096 FH.fill('0');
1097 FH.width(digits);
1098 FH<<(unsigned int)(seq[l][k]+1);
1099 }
1100 FH<<" ";
1101 }
1102
1103 FH << AGE <<" "
1104 <<SEX<<" "
1105 <<ind->getPedigreeClass()<<" "
1106 <<ind->getHome()+1<<std::endl;
1107 }
1108
1109}
unsigned int getPedigreeClass()
Returns the pedigree class of the individual, as set during offspring creation.
Definition: individual.h:179
unsigned int getID()
Definition: metapop.h:478

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, Patch::get(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), Individual::getHome(), Patch::getID(), Individual::getPedigreeClass(), Individual::getTrait(), and Patch::size().

Referenced by write_FSTAT().

+ Here is the caller graph for this function:

◆ write_patch_GENEPOP()

void TTNeutralGenesFH::write_patch_GENEPOP ( Patch patch,
sex_t  SEX,
age_idx  AGE,
ofstream &  FH,
unsigned int  digits 
)
1170{
1171 unsigned char** seq;
1172 Individual *ind;
1173
1174 for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
1175
1176 FH<<patch->getID() + 1<<", ";
1177 ind = patch->get(SEX, AGE, j);
1178 seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
1179
1180 for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
1181 for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
1182 FH.fill('0');
1183 FH.width(digits);
1184 FH<<(unsigned int)(seq[l][k]+1);
1185 }
1186 FH<<" ";
1187 }
1188
1189 FH << AGE <<" "
1190 <<SEX<<" "
1191 <<ind->getPedigreeClass()<<" "
1192 <<ind->getHome()+1<<std::endl;
1193 }
1194
1195}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, Patch::get(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), Individual::getHome(), Patch::getID(), Individual::getPedigreeClass(), Individual::getTrait(), and Patch::size().

Referenced by write_GENEPOP().

+ Here is the caller graph for this function:

◆ write_patch_TAB()

void TTNeutralGenesFH::write_patch_TAB ( Patch patch,
sex_t  SEX,
age_idx  AGE,
ofstream &  FH 
)
788{
789 unsigned char** seq;
790 Individual *ind;
791
792 for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
793
794 FH<<patch->getID() + 1<<" ";
795
796 ind = patch->get(SEX, AGE, j);
797 seq = (unsigned char**)ind->getTrait(_FHLinkedTraitIndex)->get_sequence();
798
799 if(_output_option == "snp") {
800
801 for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k)
802 FH<<(unsigned int)(seq[0][k])+(unsigned int)(seq[1][k])<<" ";
803 }
804 else {
805 for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
806 for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
807 FH<<(unsigned int)(seq[l][k]+1)<<" ";
808 }
809 }
810 }
811 FH << AGE <<" "
812 <<SEX<<" "
813 <<ind->getHome()+1<<" "
814 <<ind->getPedigreeClass()<<" "
815 << (ind->getFather() && ind->getMother() ?
816 (ind->getFather()->getHome()!=patch->getID()) + (ind->getMother()->getHome()!= patch->getID()) : 0)
817 <<" "
818 <<ind->getFatherID()<<" "
819 <<ind->getMotherID()<<" "
820 <<ind->getID()<<std::endl;
821
822 }
823
824}
Individual * getFather()
Definition: individual.h:126
Individual * getMother()
Definition: individual.h:127
unsigned long getMotherID()
Definition: individual.h:125
unsigned long getFatherID()
Definition: individual.h:124

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, _output_option, Patch::get(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), TTrait::get_sequence(), Individual::getFather(), Individual::getFatherID(), Individual::getHome(), Individual::getID(), Patch::getID(), Individual::getMother(), Individual::getMotherID(), Individual::getPedigreeClass(), Individual::getTrait(), and Patch::size().

Referenced by write_TAB().

+ Here is the caller graph for this function:

◆ write_PLINK()

void TTNeutralGenesFH::write_PLINK ( )
829{
830
832 fatal("PLINK file output for the neutral loci is only possible for di-allelic loci (for now)\n");
833
834
835 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
836 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
837 unsigned int patchNbr = _pop->getPatchNbr();
838 Patch* current_patch;
839 age_t pop_age = _pop->getCurrentAge(); //flag telling which age class should contain individuals
840
841 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + ".ped";
842
843#ifdef _DEBUG_
844 message("TTNeutralGenesFH::write_PLINK (%s)\n",filename.c_str());
845#endif
846
847 // the PED file -------------------------------------------------------------------------
848 ofstream PED (filename.c_str(), ios::out);
849
850 if(!PED) fatal("could not open plink .ped output file!!\n");
851
852 for (unsigned int i = 0; i < patchNbr; ++i) {
853
854 current_patch = _pop->getPatch(i);
855
856 if( pop_age & ADULTS )
857 print_PLINK_PED(PED, ADLTx, current_patch);
858
859 if( pop_age & OFFSPRG )
860 print_PLINK_PED(PED, OFFSx, current_patch);
861
862 }
863
864
865 PED.close();
866
867 // the MAP file -------------------------------------------------------------------------
868 filename = get_path() + this->get_service()->getGenerationReplicateFileName() + ".map";
869
870 ofstream MAP (filename.c_str(), ios::out);
871
872 if(!MAP) fatal("could not open plink .map output file!!\n");
873
874 double **map;
875 map = new double* [nb_locus];
876 for(unsigned int k = 0; k < nb_locus; ++k) map[k] = new double[2];
877
878 bool found = _FHLinkedTrait->_map.getGeneticMap(_FHLinkedTrait->get_type(), map, nb_locus);
879
880 if( found ) {
881
882 // MAP FORMAT (PLINK1.9): chrmsm ID; Loc ID; position (cM); bp ID
883 for(unsigned int k = 0; k < nb_locus; ++k) {
884 MAP<<map[k][0]+1<<" "<<"loc"<<k+1<<" "<<map[k][1]<<" "<<k+1<<endl;
885 }
886 } else { // trait didn't register a genetic map, loci are unlinked (free recombination)
887
888 warning("PLINK .map file: we assume ntrl loci are unlinked, separated by 50M and on a single chromosome\n");
889
890 // we're gonna set all loci on a single chrmsme, but 50M apart
891 for(unsigned int k = 0; k < nb_locus; ++k) {
892 MAP<<"1 "<<"loc"<<k+1<<" "<< k*5000.0 + 1.0<<" "<<k+1<<endl;
893 }
894 }
895
896 MAP.close();
897
898 for(unsigned int k = 0; k < nb_locus; ++k) delete [] map[k];
899 delete [] map;
900}
bool getGeneticMap(trait_t trait, double **table, unsigned int table_length)
Definition: ttrait_with_map.cc:874
age_t getCurrentAge()
Definition: metapop.h:297
virtual trait_t get_type() const
Definition: ttneutralgenes.h:207
void print_PLINK_PED(ofstream &FH, age_idx Ax, Patch *patch)
Definition: ttneutralgenes.cc:905
static GeneticMap _map
Definition: ttrait_with_map.h:201
unsigned int age_t
Age class flags.
Definition: types.h:46
#define OFFSPRG
Offspring age class flag.
Definition: types.h:50

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TTProtoWithMap::_map, FileHandler::_pop, ADLTx, ADULTS, fatal(), TProtoNeutralGenes::get_allele_num(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), FileHandler::get_service(), TProtoNeutralGenes::get_type(), Metapop::getCurrentAge(), FileServices::getGenerationReplicateFileName(), GeneticMap::getGeneticMap(), Metapop::getPatch(), Metapop::getPatchNbr(), message(), OFFSPRG, OFFSx, print_PLINK_PED(), and warning().

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_PLINK_BED()

void TTNeutralGenesFH::write_PLINK_BED ( ofstream &  BED)
975{
976 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
977 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
978 unsigned int patchNbr = _pop->getPatchNbr();
979 unsigned int popsize = _pop->size();
980 unsigned int n_blocs = (unsigned int)ceil(popsize / 4.0);
981 Patch* current_patch;
982
983 unsigned char BIN_GENO[4] = {0,1,2,3}; //00: homoz A/0; 01: missing; 10: heteroz; 11: homoz G/1
984
985 unsigned char BIT_MASK = 0x3; //number 3 or '11'
986
987 unsigned char MAGIC[4] = {0x6c,0x1b,0x01,'\0'};
988
989 BED<<MAGIC;
990
991 unsigned char BYTE = 0;
992
993 bitstring DATA(n_blocs * 8); //one bloc is one byte
994
995 for(unsigned int k = 0, pos; k < nb_locus; ++k) {
996
997 DATA.reset();
998
999 pos = 0;
1000
1001 for (unsigned int i = 0; i < patchNbr; ++i) {
1002
1003 current_patch = _pop->getPatch(i);
1004
1005 // store_BED(FEM, OFFSx, BED, current_patch, DATA, &pos, k);
1006 // store_BED(MAL, OFFSx, BED, current_patch, DATA, &pos, k);
1007 // store_BED(FEM, ADLTx, BED, current_patch, DATA, &pos, k);
1008 // store_BED(MAL, ADLTx, BED, current_patch, DATA, &pos, k);
1009 }
1010
1011 }
1012}
Non-template and faster implementation of std::bitset.
Definition: bitstring.h:56

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), Metapop::getPatch(), Metapop::getPatchNbr(), bitstring::reset(), and Metapop::size().

◆ write_TAB()

void TTNeutralGenesFH::write_TAB ( )

The file extension is ".txt".

734{
736 unsigned int ploidy = _FHLinkedTrait->get_ploidy();
737 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
738 unsigned int patchNbr = _pop->getPatchNbr();
739 unsigned char** seq;
740 Patch* current_patch;
741 Individual *ind;
742
743 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
744 + get_extension();
745
746#ifdef _DEBUG_
747 message("TTNeutralGenesFH::write_TAB (%s)\n",filename.c_str());
748#endif
749
750 ofstream FILE (filename.c_str(), ios::out);
751
752 if(!FILE) fatal("could not open TAB output file!!\n");
753
754 FILE<<"pop";
755
756 for (unsigned int i = 0; i < nb_locus; ++i) {
757
758 if( _output_option == "snp") {
759
760 FILE<<" l"<<i+1;
761
762 } else {
763 for (unsigned int j = 0; j < ploidy; ++j)
764 FILE<<" l"<<i+1<<"a"<<j+1;
765 }
766 }
767 //add names for the three last fields:
768 FILE<<" age sex home ped isMigrant father mother ID\n";
769
770 for (unsigned int i = 0; i < patchNbr; ++i) {
771
772 current_patch = _pop->getPatch(i);
773
774 write_patch_TAB(current_patch, FEM, OFFSx, FILE);
775 write_patch_TAB(current_patch, MAL, OFFSx, FILE);
776 write_patch_TAB(current_patch, FEM, ADLTx, FILE);
777 write_patch_TAB(current_patch, MAL, ADLTx, FILE);
778
779 }
780
781 FILE.close();
782
783}
void write_patch_TAB(Patch *patch, sex_t SEX, age_idx AGE, ofstream &FH)
Definition: ttneutralgenes.cc:787

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, _output_option, FileHandler::_pop, ADLTx, fatal(), FEM, FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), TProtoNeutralGenes::get_ploidy(), FileHandler::get_service(), FileServices::getGenerationReplicateFileName(), Metapop::getPatch(), Metapop::getPatchNbr(), MAL, message(), OFFSx, and write_patch_TAB().

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

◆ write_varcompWC()

void TTNeutralGenesFH::write_varcompWC ( )
1410{
1411 // make sure we are using the whole pop, not a sub-sampled one:
1412 Metapop *pop = get_service()->get_pop_ptr();
1413
1414 unsigned int nb_allele = _FHLinkedTrait->get_allele_num();
1415 unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
1416 unsigned int patchNbr = pop->getPatchNbr();
1417 bool added_stater = false;
1418 age_t AGE = pop->getCurrentAge();
1419
1420 if(AGE == (ADULTS | OFFSPRG) || AGE == ALL) {
1421 // warning("saving only offspring neutral allele stats\n");
1422 AGE = OFFSPRG;
1423 }
1424
1425 age_idx age = (AGE == OFFSPRG ? OFFSx : ADLTx);
1426
1428
1429 if(stater == NULL) {
1430 stater = new TTNeutralGenesSH(_FHLinkedTrait);
1432 added_stater = true;
1433 }
1434
1435 stater->setAlleleTables(AGE);
1436 stater->setHeteroTable(AGE);
1437
1438 TMatrix *globalFreq = stater->getGlobalFreqs();
1439 DataTable< double > *freqTable = stater->getAlleleFreqTable();
1440 DataTable< double > *heteroTable = stater->getHeteroTable();
1441 DataTable< unsigned int > *alleleCountTable = stater->getAlleleCountTable();
1442
1443 std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
1444
1445#ifdef _DEBUG_
1446 message("TTNeutralGenesFH::write_varcompWC (%s)\n",filename.c_str());
1447#endif
1448
1449 //init
1450 double* pop_sizes = new double [patchNbr];
1451 double tot_size;
1452 double sum_weights = 0;
1453 double nc;
1454 unsigned int extantPs = 0;
1455
1456 tot_size = pop->size(AGE);
1457
1458 for(unsigned int i = 0; i < patchNbr; i++) {
1459 pop_sizes[i] = pop->size(AGE, i);
1460 if(pop_sizes[i]) {
1461 extantPs++;
1462 sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
1463 }
1464 }
1465
1466 nc = (tot_size - sum_weights)/(extantPs-1);
1467
1468 // unsigned int np = extantPs;
1469 unsigned int npl = extantPs; //all loci typed in all patches
1470
1471 //p = _alleleFreqTable
1472 //pb = _globalAlleleFreq
1473
1474 //counting num alleles per locus
1475
1476 unsigned int *alploc = new unsigned int [nb_locus];
1477
1478 unsigned int **alploc_table = new unsigned int* [nb_locus];
1479
1480 for(unsigned int i = 0; i < nb_locus; ++i)
1481 alploc_table[i] = new unsigned int[nb_allele];
1482
1483 unsigned int tot_num_allele = 0;
1484
1485 for(unsigned int l = 0; l < nb_locus; ++l){
1486
1487 alploc[l] = 0;
1488
1489 for(unsigned int cnt, a = 0; a < nb_allele; ++a) {
1490
1491 cnt=0;
1492
1493 for(unsigned int i = 0; i < patchNbr; i++) {
1494
1495 cnt += alleleCountTable->get(i,l,a);
1496
1497 }
1498 alploc_table[l][a] = (cnt != 0);
1499 alploc[l] += (cnt != 0);
1500 }
1501 tot_num_allele += alploc[l];
1502 }
1503
1504 //correspondance with hierfstat implementation:
1505 //n, and nal are given by pop_sizes, same num ind typed at all loci in each patch
1506 //nc is the same for each locus
1507 //nt is given by tot_size, same tot num of ind typed at all loci
1508
1509 //SSG = het/2 for each allele
1510 double *SSG = new double[tot_num_allele];
1511 double *SSP = new double[tot_num_allele];
1512 double *SSi = new double[tot_num_allele];
1513 double *loc_id = new double[tot_num_allele];
1514 double *al_id = new double [tot_num_allele];
1515
1516 unsigned int all_cntr = 0;
1517
1518 double het, freq, var;
1519
1520 //computing sum of squares
1521
1522 for(unsigned int l = 0; l < nb_locus; ++l) {
1523
1524 for(unsigned int a = 0; a < nb_allele & all_cntr < tot_num_allele; ++a) {
1525
1526 if(alploc_table[l][a] == 0) continue; //do not consider alleles not present in the whole pop
1527
1528 //store locus and all identifiers for output
1529 loc_id[all_cntr] = l+1;
1530 al_id[all_cntr] = a+1;
1531
1532 SSG[all_cntr] = 0;
1533 SSi[all_cntr] = 0;
1534 SSP[all_cntr] = 0;
1535
1536 for(unsigned int p = 0; p < patchNbr; ++p){
1537
1538 if(!pop->size(AGE, p)) continue; //skip empty patches
1539
1540 het = heteroTable->get(p, l, a);
1541
1542 freq = freqTable->get(p, l, a);
1543
1544 var = freq - globalFreq->get(l, a); //(p_liu - pbar_u)^2
1545
1546 var *= var;
1547
1548 SSG[all_cntr] += het;
1549
1550 SSi[all_cntr] += 2*pop_sizes[p]*freq*(1-freq) - het/2;
1551
1552 SSP[all_cntr] += 2*pop_sizes[p]*var;
1553 }
1554
1555 all_cntr++;
1556 }
1557
1558 }
1559
1560 assert(all_cntr == tot_num_allele);
1561
1562 // open the file
1563
1564 ofstream FILE (filename.c_str(), ios::out);
1565
1566 if(!FILE) fatal("could not open neutral vcomp output file!!\n");
1567
1568 // print column names
1569
1570 if(_output_option == "allele") {
1571
1572 FILE<<"locus allele pbar het siga sigb sigw Fst Fis";
1573
1574 for (unsigned int p = 1; p <= patchNbr; ++p)
1575 FILE<<" het.p"<<p;
1576
1577 for (unsigned int p = 1; p <= patchNbr; ++p)
1578 FILE<<" freq.p"<<p;
1579
1580 FILE<<endl;
1581
1582 } else {
1583
1584 FILE<<"locus maj.al pbar.maj.al het siga sigb sigw Fst Fis";
1585
1586 for (unsigned int p = 1; p <= patchNbr; ++p)
1587 FILE<<" het.p"<<p;
1588
1589 //allele frequencies in each patch
1590 for (unsigned int p = 1; p <= patchNbr; ++p)
1591 // for(unsigned int u = 0; u < nb_allele-1; ++u) //skip last allele, can be deduced...
1592 FILE<<" freq.maj."<<"p"<< p;
1593
1594 FILE<<endl;
1595 }
1596
1597 //-----------------------------------------------------------------------------------------
1598 // allele specific stats:
1599
1600 double *MSG = new double[tot_num_allele];
1601 double *MSP = new double[tot_num_allele];
1602 double *MSI = new double[tot_num_allele];
1603 // double *sigw = new double[tot_num_allele];
1604 double *siga = new double[tot_num_allele];
1605 double *sigb = new double[tot_num_allele];
1606
1607 for(unsigned int i = 0; i < tot_num_allele; ++i){
1608
1609 MSG[i] = SSG[i] / (2 * tot_size);
1610 // sigw[i] = MSG[i]; //wasted!
1611
1612 MSP[i] = SSP[i] / (npl-1);
1613
1614 MSI[i] = SSi[i]/ (tot_size - npl);
1615
1616 sigb[i] = 0.5*(MSI[i] - MSG[i]);
1617
1618 siga[i] = (MSP[i] - MSI[i])/(2*nc);
1619
1620 if(_output_option == "allele") {
1621 FILE<< loc_id[i] << " " << al_id[i] << " "
1622 //global allele frequency
1623 << globalFreq->get(loc_id[i] - 1, al_id[i] - 1) << " "
1624 //average heterozygosity:
1625 << SSG[i]/tot_size << " "
1626 //variance components:
1627 << siga[i] << " " << sigb[i] << " " << MSG[i] << " "
1628 //Fst
1629 << siga[i]/(siga[i]+sigb[i]+MSG[i]) << " "
1630 //Fis
1631 << sigb[i]/(sigb[i]+MSG[i]) << " ";
1632 //per patch heterozygosity:
1633 for (unsigned int p = 0; p < patchNbr; ++p) {
1634 FILE << heteroTable->get(p, loc_id[i]-1, al_id[i]-1)/pop_sizes[p] << " ";
1635 }
1636 //per patch allele frequency:
1637 for (unsigned int p = 0; p < patchNbr; ++p) {
1638 FILE << freqTable->get(p, loc_id[i]-1, al_id[i]-1) << " ";
1639 }
1640
1641 FILE<< endl;
1642 }
1643 }
1644
1645 //-----------------------------------------------------------------------------------------
1646 // locus specific stats:
1647
1648 if(_output_option == "locus") {
1649
1650
1651 double lsiga, lsigb, lsigw, max_all_frq;
1652 unsigned int maj_al;
1653
1654 all_cntr = 0;
1655
1656 deque <double> loc_het = stater->setHo2(age);
1657
1658 for(unsigned int i = 0; i < nb_locus; ++i) {
1659
1660 lsiga = 0;
1661 lsigb = 0;
1662 lsigw = 0;
1663
1664 max_all_frq = 0;
1665
1666 for(unsigned int l = 0; l < alploc[i]; ++l) {
1667
1668 lsiga += siga[all_cntr];
1669 lsigb += sigb[all_cntr];
1670 lsigw += MSG[all_cntr];
1671
1672 if(max_all_frq < globalFreq->get(i, al_id[all_cntr]-1 ) ) {
1673 max_all_frq = globalFreq->get(i, al_id[all_cntr]-1 );
1674 maj_al = al_id[all_cntr];
1675 }
1676
1677 all_cntr++;
1678
1679 }
1680
1681 FILE << i+1 <<" "<< maj_al <<" "<< max_all_frq <<" ";
1682 FILE<< loc_het[i] <<" "<< lsiga << " " << lsigb <<" ";
1683 FILE<< lsigw << " "<< lsiga /(lsiga + lsigb + lsigw) <<" ";
1684 FILE<< lsigb /(lsigb + lsigw);
1685 // patch heterozygosities:
1686
1687 for(unsigned int p = 0; p < patchNbr; ++p) {
1688 het = 0;
1689 for(unsigned int a = 0; a < nb_allele; ++a){
1690 het += heteroTable->get(p, i, a);
1691 }
1692 FILE << " " << het/(2.0*pop_sizes[p]);
1693 }
1694 // patch allele freq?
1695 for(unsigned int p = 0; p < patchNbr; ++p) {
1696 FILE<< " " << freqTable->get(p, i, maj_al-1);
1697 }
1698 FILE << endl;
1699 }//END for locus
1700
1701 } //END per locus output
1702
1703
1704
1705
1706 FILE.close();
1707
1708 delete[]pop_sizes;
1709 delete[]alploc;
1710 for(unsigned int i = 0; i < nb_locus; ++i)
1711 delete[]alploc_table[i];
1712 delete[]alploc_table;
1713 delete[]loc_id;
1714 delete[]al_id;
1715 delete[]SSG;
1716 delete[]SSi;
1717 delete[]SSP;
1718 delete[]MSG;
1719 delete[]MSI;
1720 delete[]MSP;
1721 // delete[]sigw;
1722 delete[]siga;
1723 delete[]sigb;
1724}
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
A class to handle matrix in params, coerces matrix into a vector of same total size.
Definition: tmatrix.h:49
double get(unsigned int i, unsigned int j)
Accessor to element at row i and column j.
Definition: tmatrix.h:147
DataTable< unsigned int > * getAlleleCountTable()
Definition: ttneutralgenes.h:372
deque< double > setHo2(age_idx age_pos)
Definition: stats_fstat.cc:545
void setAlleleTables(age_t AGE)
Definition: stats_fstat.cc:74
TMatrix * getGlobalFreqs()
Accessor to the table of allele frequencies in the whole population.
Definition: ttneutralgenes.h:377
DataTable< double > * getHeteroTable()
Definition: ttneutralgenes.h:374
void setHeteroTable(age_t AGE)
Definition: stats_fstat.cc:166
DataTable< double > * getAlleleFreqTable()
Accessor to the table of allele frequencies, per patch.
Definition: ttneutralgenes.h:370
#define ALL
All ages age class flag.
Definition: types.h:56

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, _output_option, ADLTx, ADULTS, ALL, TTNeutralGenesSH::allocateTables(), fatal(), DataTable< T >::get(), TMatrix::get(), TProtoNeutralGenes::get_allele_num(), FileHandler::get_extension(), TProtoNeutralGenes::get_locus_num(), FileHandler::get_path(), FileServices::get_pop_ptr(), 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(), and Metapop::size().

Referenced by TProtoNeutralGenes::loadFileServices().

+ Here is the caller graph for this function:

Member Data Documentation

◆ _output_option

string TTNeutralGenesFH::_output_option
private

◆ write_fct

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.56 by  doxygen 1.9.0 -- Nemo is hosted on  Download Nemo

Locations of visitors to this page
Catalogued on GSR