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

Private Attributes
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

    : TraitFileHandler<TProtoNeutralGenes> (TP, ".txt"), write_fct(0)
  { }

~TTNeutralGenesFH()

virtual TTNeutralGenesFH::~TTNeutralGenesFH ( )

inlinevirtual

{ }

Member Function Documentation

FHread()

void TTNeutralGenesFH::FHread ( string &  filename )

virtual

reads in only FSTAT format

Implements FileHandler.

{ 
  // make sure we are using the whole pop, not a sub-sampled one:
//  Metapop *pop = get_service()->get_pop_ptr();
 
  unsigned int digit, nloci_infile, nall_infile, npatch_infile;
  //  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_all   = _FHLinkedTrait->get_allele_num();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  unsigned int extended_num_colums = 0;
 
  bool is_extended = false;
 
  ifstream FILE(filename.c_str(),ios::in);
 
  if(!FILE) fatal("could not open FSTAT input file \"%s\"\n", filename.c_str());
 
  FILE >> npatch_infile >> nloci_infile >> nall_infile >> digit;
 
  if(npatch_infile != patchNbr) fatal("number of patch in FSTAT file (%i) differs from simulation settings (%i)\n", npatch_infile, patchNbr);
 
  if (nloci_infile > nb_locus) {
 
    is_extended = true;
 
    extended_num_colums = nloci_infile - nb_locus;
 
    if(extended_num_colums != FSTAT_EXTRA_INFO_LOCI) {
      error("FSTAT input genotype file \"%s\" has a non-regular number of extra info loci of %i\n", filename.c_str(), extended_num_colums);
      error("Nemo expects %i extra info loci in the extented FSTAT format (age, sex, ped, origin)\n", FSTAT_EXTRA_INFO_LOCI);
      error("Extra info loci will be ignored\n");
    }
 
  } else if (nloci_infile == nb_locus){
 
    is_extended = false;
 
  } else {
 
    error("when reading the first line of the input FSTAT file \"%s\" \n", filename.c_str());
    fatal("please specify either (num locus) or (num locus + %i) as the second number on the first line of the input file.\n", FSTAT_EXTRA_INFO_LOCI);
 
  }
 
 
 
  if(nall_infile != nb_all)  fatal("number of alleles in FSTAT file differs from simulation settings\n");
 
  digit = (int) pow(10.0,(double)digit);
 
  string str;
 
 
#ifdef _DEBUG_
  message(">>>>  -- FSTAT%s file: \n", (is_extended ? " extended" : ""));
  message(">>>>     %i patch, %i loci, %i alleles, %i digits \n", npatch_infile, nloci_infile, nall_infile, digit);
#endif
 
  // read the locus names
  for(unsigned int i = 0; i < nloci_infile; ++i) {
 
    FILE >> str;
 
  }
 
  unsigned int patch, genot, all0, all1, age, sex, ped, origin;
  age_idx agex;
  Individual *ind;
  TTrait* trait;
  int lnbr = nloci_infile +2; //line counter
 
  unsigned int *loc_all0 = new unsigned int[nb_locus];
  unsigned int *loc_all1 = new unsigned int[nb_locus];
 
  while(FILE>>patch) {
 
    if(FILE.bad()) {
      error("Reading input genotypes from \"%s\" failed at line %i.\n",filename.c_str(), lnbr);
      FILE.clear();
      FILE >> str;
      fatal("Expecting population number as first element on the line, but received this: %s \n", str.c_str());
    }
 
    //check patch value
    if(patch > patchNbr) fatal("found an illegal patch identifier in FSTAT file (%i) at line %i",patch, lnbr);
 
    for(unsigned int i = 0; i < nb_locus; ++i) {
      FILE>>genot;
 
      all0 = genot/digit;
      all1 = genot%digit;
 
      if( all0 > nb_all ) {
 
        error("in FSTAT input file at line %i, locus %i : \
first allele value %d is greater than the max value specified (%i)!\n", lnbr, i+1, all0, nb_all);
 
        fatal("Please check the input file.\n");
 
      } else if ( all0 > 0 ) {
 
        loc_all0[i] = all0 - 1;
 
      } else {
 
        fatal("in FSTAT input file at line %i, locus %i first allele value: \
allele value 0 is not allowed!\n*** Please check the input file.\n", lnbr, i+1, all0, nb_all);
 
      }
 
      if( all1 > nb_all ){
 
        error("in FSTAT input file at line %i, locus %i : \
second allele value %i is greater than the max value specified (%i)!\n", lnbr, i+1, all1, nb_all);
 
        fatal("Please check the input file.\n");
 
      } else if (all1 > 0) {
 
        loc_all1[i] = all1 - 1;
 
      } else {
 
        fatal("in FSTAT input file at line %i, locus %i second allele value: \
allele value 0 is not allowed!\n*** Please check the input file.\n", lnbr, i+1, all0, nb_all);
 
      }
    }
 
    if(is_extended) {
 
      if(extended_num_colums == FSTAT_EXTRA_INFO_LOCI) {
 
        FILE >> age >> sex >> ped >> origin;
 
        //age index now saved in the FSTAT file
        agex = (static_cast<age_idx> (age) == ADLTx ? ADLTx : OFFSx);
 
        if(sex > 1)
          fatal("in FSTAT input file at line %i, extra column \"sex\" has irregular value %i, must be 0 or 1\n", lnbr, sex);
 
      } else { //irregular file
       getline(FILE, str); // we ignore the extra fields, and use default values
       agex = ADLTx;
       sex = FEM;
       ped = 0;
       origin = 1;
      }
 
    } else {
      // by default, individuals are adult females
      agex = ADLTx;
      sex = FEM;
      ped = 0;
      origin = 1;
 
    }
 
 
    ind = _pop->makeNewIndividual(0, 0, sex_t(sex), origin - 1);
    ind->setPedigreeClass((unsigned char)ped);
    trait = ind->getTrait(_FHLinkedTraitIndex);
    for(unsigned int i = 0; i < nb_locus; ++i) {
      trait->set_allele_value(i, 0, (double)loc_all0[i]);
      trait->set_allele_value(i, 1, (double)loc_all1[i]);
    }
 
    _pop->getPatch(patch-1)->add(sex_t(sex), agex, ind);
 
    lnbr++;
 
  }
 
  delete [] loc_all0;
  delete [] loc_all1;
 
  FILE.close();
 
#ifdef _DEBUG_
  message(">>>>     read FSTAT file with %i lines, %i individuals in pop \n", lnbr, _pop->size());
#endif
}

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

FHwrite()

void TTNeutralGenesFH::FHwrite ( )

virtual

Implements TraitFileHandler< TProtoNeutralGenes >.

{
  if(!_pop->isAlive()) return;
 
  // reset the pop ptr from the file services, will be the main metapop without sub sampling
  // or a sub sampled pop otherwise:
 
  _pop = get_service()->getSampledPop();
 
  (this->*write_fct)();
 
  // reset the pop ptr to the main pop
  _pop = get_service()->get_pop_ptr();
 
}

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
	)

{
  Individual *ind;
  TTrait* trait;
  char BASE[2] = {'A','G'};
  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
 
  // SPECIFICATION FOR THE .fam FILE = 6 first values in .ped files:
//  .fam (PLINK sample information file)
//
//  Sample information file accompanying a .bed binary genotype table.
//  Also generated by "--recode lgen" and "--recode rlist".
//
//  A text file with no header line, and one line per sample with the following six fields:
//
//      1. Family ID ('FID')
//      2. Within-family ID ('IID'; cannot be '0')
//      3. Within-family ID of father ('0' if father isn't in dataset)
//      4. Within-family ID of mother ('0' if mother isn't in dataset)
//      5. Sex code ('1' = male, '2' = female, '0' = unknown)
//      6. Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)
//
//  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.
 
 
  for (unsigned int j = 0; j < patch->size(FEM, Ax); ++j) {
 
    ind = patch->get(FEM, Ax, j);
 
    FH<<"fam"<<ind->getHome()+1
        <<" "<<ind->getID();         // don't add +1 to keep IDs consistent across files
 
    if(Ax == OFFSx)
      FH<<" "<<ind->getFatherID()<<" "<<ind->getMotherID(); //parents may be in file for offspring, although not guaranteed
    else
      FH<<" 0 0"; //parents not in file for adults
 
    FH<<" 2 -9";
 
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    for(unsigned int k = 0; k < nb_locus; ++k) {
        FH<<" "<< BASE[ trait->get_allele(k, FEM) ]<<" "<< BASE[ trait->get_allele(k, MAL) ];  // the maternally inherited allele comes first
    }
 
    FH <<std::endl;
 
  }
 
  for (unsigned int j = 0; j < patch->size(MAL, Ax); ++j) {
 
    ind = patch->get(MAL, Ax, j);
 
    FH<<"fam"<<ind->getHome()+1
        <<" "<<ind->getID();
 
    if(Ax == OFFSx)
      FH<<" "<<ind->getFatherID()<<" "<<ind->getMotherID(); //parents may be in file for offspring, although not guaranteed
    else
      FH<<" 0 0"; //parents not in file for adults
 
    FH<<" 1 -9";
 
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    for(unsigned int k = 0; k < nb_locus; ++k) {
      FH<<" "<< BASE[ trait->get_allele(k, FEM) ]<<" "<< BASE[ trait->get_allele(k, MAL) ];  // the maternally inherited allele comes first
    }
 
    FH<<std::endl;
 
  }
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, FEM, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), Individual::getFatherID(), Individual::getHome(), Individual::getID(), Individual::getMotherID(), Individual::getTrait(), MAL, OFFSx, 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

{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

)

{
 
  if(opt != "1")              //it means that the param received an argument value
    _output_option = opt;
  else
    _output_option = "locus"; //default
}

References _output_option.

Referenced by TProtoNeutralGenes::loadFileServices().

Here is the caller graph for this function:

write_Fst_i()

void TTNeutralGenesFH::write_Fst_i

(

)

{  
  // make sure we are using the whole pop, not a sub-sampled one:
  Metapop *pop = get_service()->get_pop_ptr();
 
  if(pop->size(ADULTS) == 0) {
    warning("No adults in pop, not writing the Fst distribution to file.\n");
    return;
  }
 
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = pop->getPatchNbr();
  double **fst_i;
  bool added_stater = false;
 
  TTNeutralGenesSH* stater = _FHLinkedTrait->get_stater();
 
  if(stater == NULL) {
    stater = new TTNeutralGenesSH(_FHLinkedTrait);
    stater->allocateTables(_FHLinkedTrait->get_locus_num() , _FHLinkedTrait->get_allele_num());
    added_stater = true;
  }
 
  fst_i = new double* [patchNbr];
  for(unsigned int i = 0; i < patchNbr; i++)
    fst_i[i] = new double [nb_locus];
 
  stater->setFst_li(patchNbr, nb_locus, fst_i);
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
#endif
 
  ofstream FILE (filename.c_str(), ios::out);
 
  if(!FILE) fatal("could not open Fst_i output file!!\n");
 
  for(unsigned int i = 0; i < patchNbr; ++i)
    FILE<<"patch"<<i+1<<" ";
 
  FILE<<endl;
 
  for(unsigned int j = 0; j < nb_locus; ++j) {
 
    for(unsigned int i = 0; i < patchNbr; i++)
      FILE<<fst_i[i][j]<<" ";
 
    FILE<<endl;
  }
 
  FILE.close();
 
  if(added_stater) delete stater;
 
  for(unsigned int i = 0; i < patchNbr; i++)
    delete [] fst_i[i];
 
  delete [] fst_i;
 
}

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

{
  unsigned int position;
  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_all   = _FHLinkedTrait->get_allele_num();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  Patch* current_patch;
 
  position = nb_all > 99 ? 3 : nb_all > 9 ? 2 : 1; //assumes nb_all not sup. to 999
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
                                        + get_extension();
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
#endif
 
  ofstream FILE (filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
  if(!FILE) fatal("could not open FSTAT output file!!\n");
 
  FILE<<patchNbr<<" "<< nb_locus + FSTAT_EXTRA_INFO_LOCI <<" "<<nb_all<<" "<<position<<"\n";
 
  for (unsigned int i = 0; i < nb_locus; ++i)
    FILE<<"loc"<<i+1<<"\n";
 
  //add names for the three last fields:
  FILE<<"age\n"<<"sex\n"<<"ped\n"<<"origin\n";
 
  for (unsigned int i = 0; i < patchNbr; ++i) {
 
    current_patch = _pop->getPatch(i);
 
    write_patch_FSTAT(current_patch, FEM, OFFSx, FILE, position);
    write_patch_FSTAT(current_patch, MAL, OFFSx, FILE, position);
    write_patch_FSTAT(current_patch, FEM, ADLTx, FILE, position);
    write_patch_FSTAT(current_patch, MAL, ADLTx, FILE, position);
 
  }
 
  FILE.close();
  std::ios_base::sync_with_stdio(true); // reset
 
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, FileHandler::_pop, ADLTx, fatal(), FEM, FSTAT_EXTRA_INFO_LOCI, 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".

{
  unsigned int position;
  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_all   = _FHLinkedTrait->get_allele_num();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  Patch* current_patch;
 
  position = nb_all > 99 ? 3 : 2; //assumes nb_all not sup. to 999
 
  if(ploidy == 1) position = 1;
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
                                        + get_extension();
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::FHwrite (%s)\n",filename.c_str());
#endif
 
  ofstream FILE (filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
  if(!FILE) fatal("could not open FSTAT output file!!\n");
 
  FILE<<"Title line: "<<patchNbr<<" patches, "<<nb_locus<<" loci with "<<nb_all<<" alleles\n";
 
  for (unsigned int i = 0; i < nb_locus; ++i)
    FILE<<"loc"<<i+1<<", ";
 
  //add names for the three last fields:
  FILE<<"age, sex, ped, origin\n";
 
  for (unsigned int i = 0; i < patchNbr; ++i) {
 
    current_patch = _pop->getPatch(i);
 
    FILE<<"POP\n";
 
    write_patch_GENEPOP(current_patch, FEM, OFFSx, FILE, position);
    write_patch_GENEPOP(current_patch, MAL, OFFSx, FILE, position);
    write_patch_GENEPOP(current_patch, FEM, ADLTx, FILE, position);
    write_patch_GENEPOP(current_patch, MAL, ADLTx, FILE, position);
 
 
  }
 
  FILE.close();
  std::ios_base::sync_with_stdio(true); // reset
 
}

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
	)

{
  TTrait* trait;
  Individual *ind;
 
  for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
 
    FH<<patch->getID() + 1<<" ";
    ind = patch->get(SEX, AGE, j);
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
      for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
        FH.fill('0');
        FH.width(digits);
        FH<<(trait->get_allele(k, l)+1);
      }
      FH<<" ";
    }
 
    FH << AGE <<" "
        <<SEX<<" "
        <<ind->getPedigreeClass()<<" "
        <<ind->getHome()+1<<std::endl;
  }
 
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), 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
	)

{
  TTrait* trait;
  Individual *ind;
 
  for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
 
    FH<<patch->getID() + 1<<", ";
    ind = patch->get(SEX, AGE, j);
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
      for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
        FH.fill('0');
        FH.width(digits);
        FH<<(trait->get_allele(k, l)+1);
      }
      FH<<" ";
    }
 
    FH << AGE <<" "
        <<SEX<<" "
        <<ind->getPedigreeClass()<<" "
        <<ind->getHome()+1<<std::endl;
  }
 
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), 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
	)

{
  TTrait* trait;
  Individual *ind;
 
  for (unsigned int j = 0; j < patch->size(SEX, AGE); ++j) {
 
    FH<<patch->getID() + 1<<" ";
 
    ind = patch->get(SEX, AGE, j);
    trait = ind->getTrait(_FHLinkedTraitIndex);
 
    if(_output_option == "snp") {
 
      for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k)
        FH<<trait->get_allele(k, 0)+trait->get_allele(k, 1)<<" ";
 
    } else {
 
      for(unsigned int k = 0; k < _FHLinkedTrait->get_locus_num(); ++k) {
        for (unsigned int l = 0; l < _FHLinkedTrait->get_ploidy(); ++l) {
          FH<<(trait->get_allele(k, l)+1)<<" ";
        }
      }
    }
    FH << AGE <<" "
        <<SEX<<" "
        <<ind->getHome()+1<<" "
        <<ind->getPedigreeClass()<<" "
        << (ind->getFather() && ind->getMother() ?
            (ind->getFather()->getHome()!=patch->getID()) + (ind->getMother()->getHome()!= patch->getID()) : 0)
            <<" "
        <<ind->getFatherID()<<" "
        <<ind->getMotherID()<<" "
        <<ind->getID()<<std::endl;
 
  }
 
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTraitIndex, _output_option, Patch::get(), TTrait::get_allele(), TProtoNeutralGenes::get_locus_num(), TProtoNeutralGenes::get_ploidy(), 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

(

)

{
 
  if(_FHLinkedTrait->get_allele_num() > 2)
    fatal("PLINK file output for the neutral loci is only possible for di-allelic loci\n");
 
 
  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  Patch* current_patch;
  age_t pop_age = _pop->getCurrentAge(); //flag telling which age class should contain individuals
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + "-ntrl.ped";
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::write_PLINK (%s)\n",filename.c_str());
#endif
 
  // the PED file -------------------------------------------------------------------------
  ofstream PED (filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
  if(!PED) fatal("could not open plink .ped output file!!\n");
 
  for (unsigned int i = 0; i < patchNbr; ++i) {
 
    current_patch = _pop->getPatch(i);
 
    if( pop_age & OFFSPRG )
      print_PLINK_PED(PED, OFFSx, current_patch);
 
    if( pop_age & ADULTS )
      print_PLINK_PED(PED, ADLTx, current_patch);
 
  }
 
 
  PED.close();
 
  // the MAP file -------------------------------------------------------------------------
  filename = get_path() + this->get_service()->getGenerationReplicateFileName() + "-ntrl.map";
 
  ofstream MAP (filename.c_str(), ios::out);
 
  if(!MAP) fatal("could not open plink .map output file!!\n");
 
  double **map;
  map = new double* [nb_locus];
  for(unsigned int k = 0; k < nb_locus; ++k) map[k] = new double[2];
 
  bool found = _FHLinkedTrait->_map.getGeneticMap(_FHLinkedTrait->get_type(), map, nb_locus);
 
  if( found ) {
 
    // MAP FORMAT (PLINK1.9): chrmsm ID; Loc ID; position (cM); bp ID
    for(unsigned int k = 0; k < nb_locus; ++k) {
      MAP<<map[k][0]+1<<" "<<"loc"<<k+1<<" "<<map[k][1]<<" "<<k+1<<endl;
    }
  } else { // trait didn't register a genetic map, loci are unlinked (free recombination)
 
    warning("PLINK .map file: we assume ntrl loci are unlinked, separated by 50M and on a single chromosome\n");
 
    // we're gonna set all loci on a single chrmsme, but 50M apart
    for(unsigned int k = 0; k < nb_locus; ++k) {
      MAP<<"1 "<<"loc"<<k+1<<" "<< k*5000.0 + 1.0<<" "<<k+1<<endl;
    }
  }
 
  MAP.close();
 
  std::ios_base::sync_with_stdio(true); // reset
 
  for(unsigned int k = 0; k < nb_locus; ++k) delete [] map[k];
  delete [] map;
}

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

)

write_TAB()

void TTNeutralGenesFH::write_TAB

(

)

The file extension is ".txt".

{
  unsigned int ploidy   = _FHLinkedTrait->get_ploidy();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = _pop->getPatchNbr();
  Patch* current_patch;
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName()
                                        + get_extension();
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::write_TAB (%s)\n",filename.c_str());
#endif
 
  ofstream FILE (filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
  if(!FILE) fatal("could not open TAB output file!!\n");
 
  FILE<<"pop";
 
  for (unsigned int i = 0; i < nb_locus; ++i) {
 
    if( _output_option == "snp") {
 
      FILE<<" l"<<i+1;
 
    } else {
      for (unsigned int j = 0; j < ploidy; ++j)
        FILE<<" l"<<i+1<<1<<j+1;
    }
  }
  //add names for the three last fields:
  FILE<<" age sex home ped isMigrant father mother ID\n";
 
  for (unsigned int i = 0; i < patchNbr; ++i) {
 
    current_patch = _pop->getPatch(i);
 
    write_patch_TAB(current_patch, FEM, OFFSx, FILE);
    write_patch_TAB(current_patch, MAL, OFFSx, FILE);
    write_patch_TAB(current_patch, FEM, ADLTx, FILE);
    write_patch_TAB(current_patch, MAL, ADLTx, FILE);
 
  }
 
  FILE.close();
 
  std::ios_base::sync_with_stdio(true); // reset
}

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

(

)

{  
 
  // make sure we are using the whole pop, not a sub-sampled one:
  Metapop *pop = get_service()->get_pop_ptr();
 
  unsigned int nb_allele = _FHLinkedTrait->get_allele_num();
  unsigned int nb_locus = _FHLinkedTrait->get_locus_num();
  unsigned int patchNbr = pop->getPatchNbr();
  bool added_stater = false;
 
  age_t AGE = pop->getCurrentAge();
 
  if(AGE == ALL) {
//    warning("saving only offspring neutral allele stats\n");
    AGE = OFFSPRG;
  }
 
  age_idx age = (AGE == OFFSPRG ? OFFSx : ADLTx);
 
//  cout << "----------TTNeutralGenesFH::write_varcompWC----------\n"<<endl;
//
//  cout << "   num alleles: "<<nb_allele<<endl;
//  cout << "   num loci   : "<<nb_locus <<endl;
//  cout << "   pop AGE    : "<<AGE<<endl;
//  cout << "   age pos    : "<<age<<endl;
 
//  string pad = "    ";
  TTNeutralGenesSH* stater = _FHLinkedTrait->get_stater();
 
  if(stater == NULL) {
    stater = new TTNeutralGenesSH(_FHLinkedTrait);
    stater->allocateTables(_FHLinkedTrait->get_locus_num() , _FHLinkedTrait->get_allele_num());
    added_stater = true;
  }
 
  stater->setAlleleTables(AGE);
  stater->setHeteroTable(AGE);
 
  TMatrix *globalFreq = stater->getGlobalFreqs();
  DataTable< double > *freqTable = stater->getAlleleFreqTable();
  DataTable< double > *heteroTable = stater->getHeteroTable();
  DataTable< unsigned int > *alleleCountTable = stater->getAlleleCountTable();
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + get_extension();
 
#ifdef _DEBUG_
  message("TTNeutralGenesFH::write_varcompWC (%s)\n",filename.c_str());
#endif
 
  //init
  double* pop_sizes = new double [patchNbr];
  double tot_size;
  double sum_weights = 0;
  double nc;
  unsigned int extantPs = 0;
 
  tot_size = pop->size(AGE);
 
  for(unsigned int i = 0; i < patchNbr; i++) {
    pop_sizes[i] = pop->size(AGE, i);
    if(pop_sizes[i]) {
      extantPs++;
      sum_weights += (pop_sizes[i] * pop_sizes[i] / tot_size);
    }
  }
 
  nc = (tot_size - sum_weights)/(extantPs-1);
 
  //      unsigned int np = extantPs;
  unsigned int npl = extantPs; //all loci typed in all patches
 
  //p = _alleleFreqTable
  //pb = _globalAlleleFreq
 
  unsigned int *alploc = new unsigned int [nb_locus];
 
  unsigned int **alploc_table = new unsigned int* [nb_locus];
 
  for(unsigned int i = 0; i < nb_locus; ++i)
    alploc_table[i] = new unsigned int[nb_allele];
 
  unsigned int tot_num_allele = 0;
 
  for(unsigned int l = 0; l < nb_locus; ++l){
 
    alploc[l] = 0;
 
    for(unsigned int cnt,  a = 0; a < nb_allele; ++a) {
 
      cnt=0;
 
      for(unsigned int i = 0; i < patchNbr; i++) {
 
        cnt += alleleCountTable->get(i,l,a);
 
      }
      alploc_table[l][a] = (cnt != 0);
      alploc[l] += (cnt != 0);
    }
 
    tot_num_allele += alploc[l];
  }
 
 
  //correspondance with hierfstat implementation:
  //n, and nal are given by pop_sizes, same num ind typed at all loci in each patch
  //nc is the same for each locus
  //nt is given by tot_size, same tot num of ind typed at all loci
 
  //SSG = het/2 for each allele
  double *SSG = new double[tot_num_allele];
  double *SSP = new double[tot_num_allele];
  double *SSi = new double[tot_num_allele];
  double *loc_id = new double[tot_num_allele];
  double *al_id = new double [tot_num_allele];
 
  unsigned int all_cntr = 0;
 
  double het, freq, var;
 
  for(unsigned int l = 0; l < nb_locus; ++l) {
 
    for(unsigned int a = 0; a < nb_allele & all_cntr < tot_num_allele; ++a) {
 
      if(alploc_table[l][a] == 0) continue; //do not consider alleles not present in the whole pop
 
      //store locus and all identifiers for output
      loc_id[all_cntr] = l+1;
      al_id[all_cntr] = a+1;
 
      SSG[all_cntr] = 0;
      SSi[all_cntr] = 0;
      SSP[all_cntr] = 0;
 
      for(unsigned int p = 0; p < patchNbr; ++p){
 
        if(!pop->size(AGE, p)) continue; //skip empty patches
 
        het = heteroTable->get(p, l, a);
 
        freq = freqTable->get(p, l, a);
 
        var = freq - globalFreq->get(l, a); //(p_liu - pbar_u)^2
 
        var *= var;
 
        SSG[all_cntr] += het;
 
        SSi[all_cntr] += 2*pop_sizes[p]*freq*(1-freq) - het/2;
 
        SSP[all_cntr] += 2*pop_sizes[p]*var;
      }
 
      all_cntr++;
    }
 
  }
 
  assert(all_cntr == tot_num_allele);
 
  // open the file
  ofstream FILE (filename.c_str(), ios::out);
 
  if(!FILE) fatal("could not open neutral vcomp output file!!\n");
 
  // print column names
 
 
  if(_output_option == "allele") {
 
    FILE<<"locus allele pbar het siga sigb sigw Fst Fis";
 
    for (unsigned int p = 1; p <= patchNbr; ++p)
      FILE<<" het.p"<<p;
 
    for (unsigned int p = 1; p <= patchNbr; ++p)
      FILE<<" freq.p"<<p;
 
    FILE<<endl;
 
  } else {
 
    FILE<<"locus maj.al pbar.maj.al het siga sigb sigw Fst Fis";
 
    for (unsigned int p = 1; p <= patchNbr; ++p)
      FILE<<" het.p"<<p;
 
    //allele frequencies in each patch
    for (unsigned int p = 1; p <= patchNbr; ++p)
      //                  for(unsigned int u = 0; u < nb_allele-1; ++u) //skip last allele, can be deduced...
      FILE<<" freq.maj."<<"p"<< p;
 
    FILE<<endl;
  }
 
  //-----------------------------------------------------------------------------------------
  // allele specific stats:
 
  double *MSG = new double[tot_num_allele];
  double *MSP = new double[tot_num_allele];
  double *MSI = new double[tot_num_allele];
  //      double *sigw = new double[tot_num_allele];
  double *siga = new double[tot_num_allele];
  double *sigb = new double[tot_num_allele];
 
  for(unsigned int i = 0; i < tot_num_allele; ++i){
 
    MSG[i] = SSG[i] / (2 * tot_size);
    //            sigw[i] = MSG[i]; //wasted!
 
    MSP[i] = SSP[i] / (npl-1);
 
    MSI[i] = SSi[i]/ (tot_size - npl);
 
    sigb[i] = 0.5*(MSI[i] - MSG[i]);
 
    siga[i] = (MSP[i] - MSI[i])/(2*nc);
 
    if(_output_option == "allele") {
 
      FILE<< loc_id[i] << " " << al_id[i] << " "
          //global allele frequency
          << globalFreq->get(loc_id[i] - 1, al_id[i] - 1) << " "
          //average heterozygosity:
          << SSG[i]/tot_size << " "
          //variance components:
          << siga[i] << " " << sigb[i] << " " << MSG[i] << " "
          //Fst
          << siga[i]/(siga[i]+sigb[i]+MSG[i]) << " "
          //Fis
          << sigb[i]/(sigb[i]+MSG[i]) << " ";
      //per patch heterozygosity:
      for (unsigned int p = 0; p < patchNbr; ++p) {
        FILE << heteroTable->get(p, loc_id[i]-1, al_id[i]-1)/pop_sizes[p] << " ";
      }
      //per patch allele frequency:
      for (unsigned int p = 0; p < patchNbr; ++p) {
        FILE << freqTable->get(p, loc_id[i]-1, al_id[i]-1) << " ";
      }
 
      FILE<< endl;
    }
  }
 
  //-----------------------------------------------------------------------------------------
  if(_output_option == "locus") {
 
    double lsiga, lsigb, lsigw, max_all_frq;
    unsigned int maj_al = 0;
 
    all_cntr = 0;
 
    deque <double> loc_het = stater->setHo2(age);
 
    for(unsigned int i = 0; i < nb_locus; ++i) {
 
      lsiga = 0;
      lsigb = 0;
      lsigw = 0;
 
      max_all_frq = 0;
 
      for(unsigned int l = 0; l < alploc[i]; ++l) {
 
        lsiga += siga[all_cntr];
        lsigb += sigb[all_cntr];
        lsigw += MSG[all_cntr];
 
        if(max_all_frq < globalFreq->get(i, al_id[all_cntr]-1 ) ) {
          max_all_frq = globalFreq->get(i, al_id[all_cntr]-1 );
          maj_al = al_id[all_cntr];
        }
 
        all_cntr++;
 
      }
      FILE << i+1 <<" "<< maj_al <<" "<< max_all_frq <<" ";
      FILE<< loc_het[i] <<" "<< lsiga << " " << lsigb <<" ";
      FILE<< lsigw << " "<< lsiga /(lsiga + lsigb + lsigw) <<" ";
      FILE<< lsigb /(lsigb + lsigw)<<fflush;
 
      for(unsigned int p = 0; p < patchNbr; ++p) {
        het = 0;
        for(unsigned int a = 0; a < nb_allele; ++a){
          het += heteroTable->get(p, i, a);
        }
        FILE << " " << het/(2.0*pop_sizes[p]);
      }
      for(unsigned int p = 0; p < patchNbr; ++p) {
        FILE<< " " << freqTable->get(p, i, maj_al-1);
      }
      FILE << endl;
    }//END for locus
 
  } //END per locus output
 
 
  FILE.close();
 
  delete[]pop_sizes;
  delete[]alploc;
  for(unsigned int i = 0; i < nb_locus; ++i)
    delete[]alploc_table[i];
  delete[]alploc_table;
  delete[]loc_id;
  delete[]al_id;
  delete[]SSG;
  delete[]SSi;
  delete[]SSP;
  delete[]MSG;
  delete[]MSI;
  delete[]MSP;
  //      delete[]sigw;
  delete[]siga;
  delete[]sigb;
 
  if(added_stater) delete stater;
}

References TraitFileHandler< TProtoNeutralGenes >::_FHLinkedTrait, _output_option, ADLTx, 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

Referenced by setOutputOption(), write_patch_TAB(), write_TAB(), and write_varcompWC().

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:

• ttneutralgenes.h
• ttneutralgenes.cc