TTQuantiFH Class Reference

TTQuantiFH. More...

#include <ttquanti.h>

Inheritance diagram for TTQuantiFH:

Collaboration diagram for TTQuantiFH:

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

Private Attributes
string	_output_option
bool	_has_genetic_map

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

Detailed Description

TTQuantiFH.

Constructor & Destructor Documentation

TTQuantiFH()

TTQuantiFH::TTQuantiFH ( TProtoQuanti *  T )

inline

                              : TraitFileHandler<TProtoQuanti>(T,".quanti"),
    _has_genetic_map(0) {}

~TTQuantiFH()

virtual TTQuantiFH::~TTQuantiFH ( )

inlinevirtual

{}

Member Function Documentation

FHread()

void TTQuantiFH::FHread ( string &  filename )

virtual

Implements FileHandler.

{
//  unsigned int nb_locus = _FHLinkedTrait->get_num_locus();
  unsigned int nb_trait = _FHLinkedTrait->get_num_traits();
  unsigned int seqLength = _FHLinkedTrait->get_seq_length();
  bool has_genotype = (!_FHLinkedTrait->get_env_var().empty());
  bool has_additive_genotype = (_FHLinkedTrait->get_dominance_model() != 0);
 
  unsigned int patchNbr = _pop->getPatchNbr();
 
  TTQuanti_continuous* trait;
  TTQuanti* DA_trait;
 
  ifstream FILE(filename.c_str(),ios::in);
 
  if(!FILE) fatal("could not open QUANTI input file \"%s\"\n",filename.c_str());
 
  unsigned int age, sex, ped, origin, dummy;
  unsigned int loc;
  age_idx stagex;
  unsigned int pop;
//  age_idx agex;
  Individual *ind;
  double all0, all1, pheno, geno;
  double *seq[2];
  seq[0] = new double [seqLength];
  seq[1] = new double [seqLength];
 
  unsigned char *binseq[2];
  binseq[0] = new unsigned char [seqLength];
  binseq[1] = new unsigned char [seqLength];
 
  int lnbr = 1;
 
  // swallow the header
  FILE.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
 
  // swallow the second line if it stores the dominance effects:
  if(has_additive_genotype)
    FILE.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
 
  string str;
 
  while(FILE>>pop) {
 
    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());
    }
 
    if(pop > patchNbr)
      fatal("Patch number found in file exceeds number of patches in the population.\n");
 
    for(unsigned int k = 0; k < nb_trait; ++k) {
      for(unsigned int i = 0; i < _FHLinkedTrait->get_num_locus(k); ++i) {
 
        loc = _FHLinkedTrait->get_locus_seq_pos(i, k);
 
        FILE>>all0>>all1;
 
        if(FILE.bad())
          fatal("Reading input genotypes from \"%s\" failed at line %i.\n",filename.c_str(), lnbr);
 
        seq[0][loc] = all0;
        seq[1][loc] = all1;
      }
    }
 
 
    for(unsigned int k = 0; k < nb_trait; k++) {
      FILE >> pheno;                   //read phenotypic value
      if(has_genotype) FILE >> geno;   // genotypic values, only present if trait is influenced by env. var.
      if(has_additive_genotype) FILE >> geno; // geno unused;
    }
 
    // read extra fields
    FILE >> age >> sex >> origin >> ped >> dummy >> dummy >> dummy >> dummy;
 
    stagex = age_idx(age);
 
    ind = _pop->makeNewIndividual(0, 0, sex_t(sex), origin - 1);
    ind->setPedigreeClass((unsigned char)ped);
    ind->setAge(age);
 
    if(_FHLinkedTrait->get_allele_model() > 2) { // non-diallelic loci
 
      trait = dynamic_cast<TTQuanti_continuous*> (ind->getTrait(_FHLinkedTraitIndex));
      trait->set_sequence((void**)seq);
      trait->set_value();
 
    } else { // diallelic loci
 
      DA_trait = dynamic_cast<TTQuanti*> (ind->getTrait(_FHLinkedTraitIndex));
 
      //have to make a binary sequence out of the data read:
      for (unsigned int l = 0; l < seqLength; ++l) {
 
        binseq[0][l] = (seq[0][l] != _FHLinkedTrait->get_seq_diallele_value(l, 0));
        binseq[1][l] = (seq[1][l] != _FHLinkedTrait->get_seq_diallele_value(l, 0));
 
      }
 
      DA_trait->set_sequence((void**)binseq);
      DA_trait->set_value();
    }
 
 
    //      ind->show_up();
 
    _pop->getPatch(pop-1)->add(sex_t(sex), stagex, ind);
 
    lnbr++;
  }
 
 
  FILE.close();
 
  delete seq[0];
  delete seq[1];
  delete binseq[0];
  delete binseq[1];
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, TraitFileHandler< TProtoQuanti >::_FHLinkedTraitIndex, FileHandler::_pop, Patch::add(), error(), fatal(), TProtoQuanti::get_allele_model(), TProtoQuanti::get_dominance_model(), TProtoQuanti::get_env_var(), TProtoQuanti::get_locus_seq_pos(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), TProtoQuanti::get_seq_diallele_value(), TProtoQuanti::get_seq_length(), Metapop::getPatch(), Metapop::getPatchNbr(), Individual::getTrait(), IndFactory::makeNewIndividual(), TTQuanti_continuous::set_sequence(), TTrait::set_sequence(), TTQuanti::set_value(), Individual::setAge(), and Individual::setPedigreeClass().

FHwrite()

void TTQuantiFH::FHwrite ( )

virtual

Implements TraitFileHandler< TProtoQuanti >.

{
  if (!_pop->isAlive()) return;
 
  if(!get_service()) fatal("TTQuantiFH: link to file services not set!!\n");
 
  // sets the pop ptr to a sub sampled pop, if sub sampling happened
  _pop = get_service()->getSampledPop();
 
  if(_output_option == "plink")
    write_PLINK();
  else
    write_TABLE();
 
 
  // reset to pop ptr to main pop:
  _pop = get_service()->get_pop_ptr();
 
}

References _output_option, FileHandler::_pop, fatal(), FileServices::get_pop_ptr(), FileHandler::get_service(), FileServices::getSampledPop(), Metapop::isAlive(), write_PLINK(), and write_TABLE().

print()

void TTQuantiFH::print	(	ofstream &	FH,
		Patch *	patch,
		sex_t	SEX,
		age_idx	Ax,
		unsigned int	print_gene,
		bool	print_genotype,
		bool	print_additive_genotype
	)

{
  Individual* ind;
  TTQuanti* trait;
 
  double* Tval;
 
  unsigned int pos, nb_trait=_FHLinkedTrait->get_num_traits();
 
  for(unsigned int j = 0, size = patch->size(SEX, Ax); j < size; j++) {
 
    ind = patch->get(SEX, Ax, j);
    trait = dynamic_cast<TTQuanti*> (ind->getTrait(_FHLinkedTraitIndex));
 
    FH<<patch->getID() +1<<" ";
 
 
    if(print_gene) { // print the genotype values at each locus
 
      FH.precision(6);
 
      for(unsigned int k = 0; k < nb_trait; k++) {
        for(unsigned int l = 0; l < _FHLinkedTrait->get_num_locus(k); l++) {
 
          pos = _FHLinkedTrait->get_locus_seq_pos(l, k);
 
          if(print_gene == 1) {  // print the allelic values, not the allelic state
 
            FH << trait->get_allele_value(pos, FEM) <<" "
               << trait->get_allele_value(pos, MAL) <<" ";
 
          } else { // must be option 2 or more; print the genotype state, 0/1/2
 
            //SNP genotypes, for di-allelic loci only (checked already)
 
            if(print_gene == 2) {
              // SNP genotypes are 0/1/2 for 00/01/11 genotypes
              FH << trait->get_allele_bit(pos, FEM) + trait->get_allele_bit(pos, MAL) << " ";
            } else {
              // print the allele ID 0/1
              FH << trait->get_allele_bit(pos, FEM)<< " "
                 << trait->get_allele_bit(pos, MAL) << " ";
            }
          }
        }
      }
    }
 
    FH.precision(4);
 
    Tval = (double*)trait->getValue();
 
    for(unsigned int k = 0; k < nb_trait; k++) {
      FH<<Tval[k]<<" ";
      if(print_genotype) FH << trait->get_full_genotype(k) << " ";
      if(print_additive_genotype) FH << trait->get_additive_genotype(k) << " ";
    }
 
    FH<<Ax<<" "<<ind->getSex()<<" "<<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< TProtoQuanti >::_FHLinkedTrait, TraitFileHandler< TProtoQuanti >::_FHLinkedTraitIndex, FEM, Patch::get(), TTQuanti::get_additive_genotype(), TTQuanti::get_allele_bit(), TTrait::get_allele_value(), TTQuanti::get_full_genotype(), TProtoQuanti::get_locus_seq_pos(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), Individual::getFather(), Individual::getFatherID(), Individual::getHome(), Individual::getID(), Patch::getID(), Individual::getMother(), Individual::getMotherID(), Individual::getPedigreeClass(), Individual::getSex(), Individual::getTrait(), TTQuanti::getValue(), MAL, and Patch::size().

Referenced by write_TABLE().

Here is the caller graph for this function:

print_PLINK_FAM()

void TTQuantiFH::print_PLINK_FAM	(	ofstream &	FH,
		age_idx	Ax,
		Patch *	patch
	)

{
  Individual *ind;
  unsigned int nb_trait = _FHLinkedTrait->get_num_traits();
 
  // 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.
 
 
  // here, we save the phenotypic value of the first trait by default
  // the complete set of phenotypic values for all traits is saved in the .fam file
 
  for (unsigned int j = 0; j < patch->size(FEM, Ax); ++j) {
 
    ind = patch->get(FEM, 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<<" 2";
 
    for(unsigned int k = 0; k < nb_trait; ++k) {
      FH << " " << ((double*)ind->getTraitValue(_FHLinkedTraitIndex))[k];
    }
 
    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";
 
    for(unsigned int k = 0; k < nb_trait; ++k) {
      FH << " " << ((double*)ind->getTraitValue(_FHLinkedTraitIndex))[k];
    }
 
    FH<<std::endl;
 
  }
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, TraitFileHandler< TProtoQuanti >::_FHLinkedTraitIndex, FEM, Patch::get(), TProtoQuanti::get_num_traits(), Individual::getFatherID(), Individual::getHome(), Individual::getID(), Individual::getMotherID(), Individual::getTraitValue(), MAL, OFFSx, and Patch::size().

Referenced by write_PLINK().

Here is the caller graph for this function:

print_PLINK_PED()

void TTQuantiFH::print_PLINK_PED	(	ofstream &	FH,
		age_idx	Ax,
		Patch *	patch
	)

{
  Individual *ind;
//  double** seq;
  TTQuanti* trait;
  char BASE[2] = {'A','G'};
//  const TMatrix& ref_all = _FHLinkedTrait->get_diallele_values();
  unsigned int nb_trait = _FHLinkedTrait->get_num_traits(), pos;
 
  // 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.
 
 
  // here, we save the phenotypic value of the first trait by default
  // the complete set of phenotypic values for all traits is saved in the .fam file
 
  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 "
      <<((double*)ind->getTraitValue(_FHLinkedTraitIndex))[0]<<" ";
 
    trait = dynamic_cast<TTQuanti*>(ind->getTrait(_FHLinkedTraitIndex));
 
    for(unsigned int k = 0; k < nb_trait; ++k) {
      for (unsigned int l = 0; l < _FHLinkedTrait->get_num_locus(k); ++l) {
 
        pos = _FHLinkedTrait->get_locus_seq_pos(l, k);
 
 
        FH << BASE[ trait->get_allele_bit(pos, FEM)] << " "
           << BASE[ trait->get_allele_bit(pos, MAL)] << " ";
 
      }
    }
 
    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 "
      <<((double*)ind->getTraitValue(_FHLinkedTraitIndex))[0]<<" ";
 
    trait = dynamic_cast<TTQuanti*>(ind->getTrait(_FHLinkedTraitIndex));
 
    for(unsigned int k = 0; k < nb_trait; ++k) {
      for (unsigned int l = 0; l < _FHLinkedTrait->get_num_locus(k); ++l) {
 
        pos = _FHLinkedTrait->get_locus_seq_pos(l, k);
 
        FH << BASE[ trait->get_allele_bit(pos, FEM)] << " "
           << BASE[ trait->get_allele_bit(pos, MAL)] << " ";
 
      }
    }
    FH<<std::endl;
 
  }
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, TraitFileHandler< TProtoQuanti >::_FHLinkedTraitIndex, FEM, Patch::get(), TProtoQuanti::get_locus_seq_pos(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), Individual::getFatherID(), Individual::getHome(), Individual::getID(), Individual::getMotherID(), Individual::getTrait(), Individual::getTraitValue(), MAL, OFFSx, and Patch::size().

Referenced by write_PLINK().

Here is the caller graph for this function:

setOutputOption()

void TTQuantiFH::setOutputOption

(

string

opt

)

{
  _output_option = opt;
 
  std::for_each(_output_option.begin(), _output_option.end(), [](char& c){c = ::tolower(c);});
 
  if( !(_output_option == "genotypes"   || _output_option == "genotype"     ||
      _output_option == "snp_genotypes" || _output_option == "snp_genotype" ||
      _output_option == "snp_id"        || _output_option == "plink"        ||
      _output_option == "1"             || _output_option == "0")) {
    fatal("option \"%s\" to parameter \"quanti_output\" is not recognized\n", opt.c_str());
  }
 
  _has_genetic_map = _FHLinkedTrait->_map.checkRegisteredTrait(trait_t(_FHLinkedTrait->get_type()));
 
  if(_output_option == "plink" && !_has_genetic_map)
    warning("PLINK .map file for quant trait: no genetic map found for QTL, we assume they are unlinked, separated by 50M and on a single chromosome.\n");
 
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, _has_genetic_map, TTProtoWithMap::_map, _output_option, GeneticMap::checkRegisteredTrait(), fatal(), TProtoQuanti::get_type(), and warning().

Referenced by TProtoQuanti::loadFileServices().

Here is the caller graph for this function:

write_PLINK()

void TTQuantiFH::write_PLINK

(

)

{
  if(_FHLinkedTrait->get_allele_model() > 2)
    fatal("PLINK file output for the quanti trait is only possible for di-allelic loci.\n");
 
  unsigned int patchNbr = _pop->getPatchNbr();
  Patch* current_patch;
 
  std::string filename = get_path() + this->get_service()->getGenerationReplicateFileName() + "-quanti.ped";
 
  // 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");
 
#ifdef _DEBUG_
  message("TTQuantiFH::write_PLINK (%s)\n",filename.c_str());
#endif
 
  for (unsigned int i = 0; i < patchNbr; ++i) {
 
    current_patch = _pop->getPatch(i);
 
    if( _pop->getCurrentAge() & OFFSPRG)
      print_PLINK_PED(PED, OFFSx, current_patch);
 
    if( _pop->getCurrentAge() & ADULTS) {
      for(unsigned int a = 1; a < _pop->getNumAgeClasses(); ++a)
        print_PLINK_PED(PED, age_idx(a), current_patch);
    }
 
  }
 
  PED.close();
 
  // the FAM file -------------------------------------------------------------------------
  // we write the .fam file only when more than one trait; phenotype columns are added here
 
  if(_FHLinkedTrait->get_num_traits() > 1) {
 
    filename = get_path() + this->get_service()->getGenerationReplicateFileName() + "-quanti.fam";
 
    ofstream FAM (filename.c_str(), ios::out);
    std::ios_base::sync_with_stdio(false); // better for writing performances
 
    if(!FAM) fatal("could not open plink .fam output file!!\n");
 
#ifdef _DEBUG_
    message("TTQuantiFH::write_PLINK (%s)\n",filename.c_str());
#endif
 
    for (unsigned int i = 0; i < patchNbr; ++i) {
 
      current_patch = _pop->getPatch(i);
 
      if( _pop->getCurrentAge() & OFFSPRG)
        print_PLINK_FAM(FAM, OFFSx, current_patch);
 
      if( _pop->getCurrentAge() & ADULTS){
        for(unsigned int a = 1; a < _pop->getNumAgeClasses(); ++a)
          print_PLINK_FAM(FAM, age_idx(a), current_patch);
      }
 
    }
 
    FAM.close();
  }
 
  // the MAP file -------------------------------------------------------------------------
  // !! careful here with pleiotropy, same locus will be registered for each trait it affects
  // output for pleiotropic locis has not be tested with PLINK
 
  filename = get_path() + this->get_service()->getGenerationReplicateFileName() + "-quanti.map";
 
  ofstream MAP (filename.c_str(), ios::out);
 
  if(!MAP) fatal("could not open plink .map output file!!\n");
 
#ifdef _DEBUG_
  message("TTQuantiFH::write_PLINK (%s)\n",filename.c_str());
#endif
 
 
  if( _has_genetic_map ) {
 
    unsigned int nb_locus = _FHLinkedTrait->get_num_locus();
 
    double *map[2];
    map[0] = new double [nb_locus]; // chromosome ID
    map[1] = new double [nb_locus]; // position
 
    _FHLinkedTrait->_map.getGeneticMap(trait_t(_FHLinkedTrait->get_type()), map, nb_locus);
 
    for(unsigned int k = 0; k < _FHLinkedTrait->get_num_traits() ; ++k) {
      for (unsigned int l = 0, LOCUS; l < _FHLinkedTrait->get_num_locus(k); ++l) { // this will give the locus ID
 
        LOCUS = _FHLinkedTrait->get_locus_ID(l, k);
 
        // MAP FORMAT (PLINK1.9): chrmsm ID; Loc ID; position (cM); bp ID
        MAP <<map[0][ LOCUS ]+1
            <<" T"<< k + 1 <<"loc"<< LOCUS + 1 <<" "
            <<map[1][ LOCUS ]<<" "<< LOCUS + 1 <<endl;
 
      }
    }
 
 
    delete [] map[0];
    delete [] map[1];
 
  } else { // trait didn't register a genetic map, loci are unlinked (free recombination)
 
    // we're gonna set all loci on a single chrmsme, but 50M apart
    for(unsigned int k = 0; k < _FHLinkedTrait->get_num_traits() ; ++k) {  //_FHLinkedTrait->get_num_traits()
      for (unsigned int l = 0, LOCUS; l < _FHLinkedTrait->get_num_locus(k); ++l) { // this will give the locus ID
 
        LOCUS = _FHLinkedTrait->get_locus_ID(l, k) + 1;
 
        // write: chromosome, locus name, map position, locus number
        MAP<< 1 <<" T"<<k+1<<"loc"<< LOCUS <<" "<< l*5000.0 + 1.0 <<" "<< LOCUS <<endl;
 
      }
    }
 
  }
  MAP.close();
 
  std::ios_base::sync_with_stdio(true); // reset
 
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, _has_genetic_map, TTProtoWithMap::_map, FileHandler::_pop, ADULTS, fatal(), TProtoQuanti::get_allele_model(), TProtoQuanti::get_locus_ID(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), FileHandler::get_path(), FileHandler::get_service(), TProtoQuanti::get_type(), Metapop::getCurrentAge(), FileServices::getGenerationReplicateFileName(), GeneticMap::getGeneticMap(), Metapop::getNumAgeClasses(), Metapop::getPatch(), Metapop::getPatchNbr(), message(), OFFSPRG, OFFSx, print_PLINK_FAM(), and print_PLINK_PED().

Referenced by FHwrite().

Here is the caller graph for this function:

write_TABLE()

void TTQuantiFH::write_TABLE

(

)

{
  std::string filename = get_filename();
 
  std::ofstream FILE (filename.c_str(), ios::out);
  std::ios_base::sync_with_stdio(false); // better for writing performances
 
 
  if(!FILE) fatal("could not open \"%s\" output file!!\n",filename.c_str());
 
#ifdef _DEBUG_
  message("TTQuantiFH::write_TABLE (%s)\n", filename.c_str());
#endif
 
  //print the genotype at each locus
  unsigned int print_gene = 0;
 
  if(_output_option == "genotypes" || _output_option == "genotype")
    print_gene = 1;
  else if(_output_option == "snp_genotypes" || _output_option == "snp_genotype")
    print_gene = 2;
  else if(_output_option == "snp_id")
    print_gene = 3;
 
 
  //print the genotypic value of the individual in sus of the phenotype
  bool print_genotype = (!_FHLinkedTrait->get_env_var().empty()) || print_gene != 0;
 
  bool print_additive_genotype = (_FHLinkedTrait->get_dominance_model() != 0);
 
  // NOTE: it is not possible to separate total genotypic value from pure additive value in the epistatic case
 
  // First column is the patch ID:
  FILE<<"pop ";
 
  // print the allelic values?
  if(print_gene) {
 
    // we print allelic values ordered trait:locus, alleles at loci affecting one trait are grouped
    for(unsigned int k = 0; k < _FHLinkedTrait->get_num_traits() ; k++) {
      for(unsigned int l = 0, LOCUS; l < _FHLinkedTrait->get_num_locus(k) ; l++) {
 
        LOCUS = _FHLinkedTrait->get_locus_ID(l, k) + 1;
 
        if(print_gene == 1 || print_gene == 3) {  // print two values per locus
          // "x" maternally inherited, "y" paternally inherited
          FILE<<"t"<<k+1<<"l"<< LOCUS <<"x "<<"t"<<k+1<<"l"<< LOCUS <<"y ";
 
        } else { // means > 1; print snp_genotype 0/1/2, instead of alleles state/value
 
          //check allele model:
          if(_FHLinkedTrait->get_allele_model() > 2)
            fatal("trait quant output::SNP genotypes can only be printed for di-allelic loci.\n");
 
          // snp-encoded genotypes: 0, 1, 2
          FILE<<"t"<<k+1<<"l"<< LOCUS <<" ";
        }
      }
    }
  }
 
  for(unsigned int k = 0; k < _FHLinkedTrait->get_num_traits(); k++) {
    FILE<<"P"<<k+1<< " ";
    if(print_genotype) FILE<<"G"<<k+1<< " ";
    if(print_additive_genotype) FILE<<"A"<<k+1<< " ";
  }
 
  FILE<<"age sex home ped isMigrant father mother ID\n";
 
 
  Patch* current_patch;
 
  for(unsigned int i = 0; i < _pop->getPatchNbr(); ++i) {
 
    current_patch = _pop->getPatch(i);
 
    if( current_patch->size(OFFSx) != 0 ) {
      print(FILE, current_patch, FEM, OFFSx, print_gene, print_genotype, print_additive_genotype);
      print(FILE, current_patch, MAL, OFFSx, print_gene, print_genotype, print_additive_genotype);
    }
 
    if( current_patch->size(ADLTx) != 0 ) {
      print(FILE, current_patch, FEM, ADLTx, print_gene, print_genotype, print_additive_genotype);
      print(FILE, current_patch, MAL, ADLTx, print_gene, print_genotype, print_additive_genotype);
    }
 
  }
 
  FILE.close();
  std::ios_base::sync_with_stdio(true); // reset
 
}

References TraitFileHandler< TProtoQuanti >::_FHLinkedTrait, _output_option, FileHandler::_pop, ADLTx, fatal(), FEM, TProtoQuanti::get_allele_model(), TProtoQuanti::get_dominance_model(), TProtoQuanti::get_env_var(), FileHandler::get_filename(), TProtoQuanti::get_locus_ID(), TProtoQuanti::get_num_locus(), TProtoQuanti::get_num_traits(), Metapop::getPatch(), Metapop::getPatchNbr(), MAL, message(), OFFSx, print(), and Patch::size().

Referenced by FHwrite().

Here is the caller graph for this function:

Member Data Documentation

_has_genetic_map

bool TTQuantiFH::_has_genetic_map

private

Referenced by setOutputOption(), and write_PLINK().

_output_option

string TTQuantiFH::_output_option

private

Referenced by FHwrite(), setOutputOption(), and write_TABLE().

---

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

• ttquanti.h
• ttquanti.cc