GeneticMap Class Reference

GeneticMap. More...

#include <ttrait_with_map.h>

Collaboration diagram for GeneticMap:

Public Member Functions
	GeneticMap ()
	~GeneticMap ()
bool	getGeneticMap (trait_t trait, double **table, unsigned int table_length)
double	getResolution ()
double	setResolution (double val)
void	rescaleMap (double val)
void	reset_tables ()
void	clear ()
void	setLookupTable (unsigned int idx)
	Bbuilds the lookup table for each trait. More...
void	recombine (sex_t SEX)
	Called by TTProtoWithMap::recombine twice to create the two gametes necessary for the creation of a new individual. More...
vector< pair< unsigned int, unsigned int > >	reduceJunctions (sex_t SEX, unsigned int trait_idx)
	Remove multiple x-over at the same locus when traits differ in number of loci. More...
bool	registerIndForRecombine (unsigned long ID)
	Called by TTProtoWithMap::recombine with individual ID passed down from Individual::recombine. More...
unsigned int	addTrait (trait_t trait, unsigned int nChrm, unsigned int nLoc, unsigned int *nLocChrm, double resolution, unsigned int *locPositions)
	Returns the table index for the registered trait. More...
void	unregisterTrait (trait_t trait)
bool	checkRegisteredTrait (trait_t trait)
	Returns true if trait 'trait' has registered a genetic map, false otherwise. More...
vector< unsigned int > &	getRecLoci (sex_t SEX, unsigned int trait)
	Returns a vector of the loci where crossing-overs take place. More...
vector< bool > &	getFirstRecPosition (sex_t SEX)
	Returns the vector of the first chromosome position for recombination, used for all traits. More...
unsigned int *	getLocusPositionTable (trait_t trait)

Private Attributes
unsigned long	_currentIndividual
map< trait_t, unsigned int >	_traits
	Table mapping trait type to its position index in the following tables. More...
unsigned int	_nTrait
	Number of traits registered in the map. More...
vector< unsigned int * >	_lociLookupTable
	A list of tables that map the map position (cM) to a locus, for each trait. More...
vector< unsigned int >	_numChrsmPerTrait
	Vector of number of chromosomes for each trait. More...
vector< unsigned int >	_numLociPerTrait
	Vector of number of loci for each trait. More...
vector< unsigned int * >	_numLociPerChrsmPerTrait
	Vector containing a table of number of loci per chromosome for each trait. More...
vector< unsigned int * >	_locPositionsPerTrait
	Vector containing the table of map position for the loci of each trait. More...
vector< vector< unsigned int > >	_recPositions [2]
	Vector of tables containing, for each trait, the locus number at which x-overs happen. More...
vector< bool >	_chrsmFirstRecombPosition [2]
	Two vectors holding the starting copy of each chromosome to use when creating the two gametes that are used to create a new individual. More...
vector< unsigned int >	_junctions
	A vector to store the position of the recombination events. More...
unsigned int	_numChromosome
unsigned int *	_perChrsmLength
unsigned int *	_chrsmFirstLocusPosition
unsigned int	_totalLength
unsigned int	_recombLength
unsigned int	_totalNumLoci
double	_resolution
double	_totRecombEventsMean
double	_recombinationRate

Detailed Description

GeneticMap.

Constructor & Destructor Documentation

GeneticMap()

GeneticMap::GeneticMap ( )

inline

               :_currentIndividual(0), _nTrait(0),
  _numChromosome(0), _perChrsmLength(0), 
  _chrsmFirstLocusPosition(0), _totalLength(0), _recombLength(0), _totalNumLoci(0), _resolution(1),
  _totRecombEventsMean(0), _recombinationRate(0.5) {}

~GeneticMap()

GeneticMap::~GeneticMap ( )

inline

{reset_tables();}

References reset_tables().

Member Function Documentation

addTrait()

unsigned int GeneticMap::addTrait	(	trait_t	trait,
		unsigned int	nChrm,
		unsigned int	nLoc,
		unsigned int *	nLocChrm,
		double	resolution,
		unsigned int *	locPositions
	)

Returns the table index for the registered trait.

{  
  map<trait_t, unsigned int>::iterator tIter;
  
  unsigned int traitIdx = 0;
  
  bool do_add = true;
    
  tIter = _traits.find(trait);
  
  if ( tIter != _traits.end() ) {
    //this shouldn't happen... unless we are loading a pop from source!!
    
    traitIdx = tIter->second;
    
    if(_numChrsmPerTrait[traitIdx] != nChrm) 
      fatal("mismatch while loading source population and resetting genetic map for trait \"%s\" \n\
>>>>> number of chromosomes differ (%i != %i)\n\
>>>>> please check match between init file and source population\n",trait.c_str(),_numChrsmPerTrait[traitIdx],nChrm);
    
    if(_numLociPerTrait[traitIdx] != nLoc) 
      fatal("mismatch while loading source population and resetting genetic map for trait \"%s\"\n\
>>>>> number of loci differ (%i != %i)\n\
>>>>> please check match between init file and source population\n",trait.c_str(),_numLociPerTrait[traitIdx], nLoc);
    
    unsigned int *locTable = _numLociPerChrsmPerTrait[traitIdx];
    
    for (unsigned int i=0; i<nChrm; i++)
      if(locTable[i] != nLocChrm[i]) 
        fatal("mismatch while loading source population and resetting genetic map for trait \"%s\"\n\
>>>>> number of loci differ (%i != %i) on chromosome %i\n\
>>>>> please check match between init file and source population\n",trait.c_str(),locTable[i], nLocChrm[i], i+1);
    
    do_add = false;
    
  }
  
#ifdef _DEBUG_
  cout << "GeneticMap::addTrait: \n";
  cout << "   _numChromosome = "<<nChrm<<endl;
  cout << "   _totalNumLoci  = "<<_totalNumLoci<<endl;
  cout << "   _resolution    = "<<_resolution<<endl;
  cout << "   _totalLength   = "<<_totalLength<<endl;
  cout << ":::"<<trait<<"::: params:\n";
  cout << "   num chromosome = "<<nChrm<<endl;
  cout << "   num loci       = "<<nLoc<<endl;
  cout << "   map resolution = "<<resolution<<endl;
#endif
 
  if (_nTrait == 0) { //this is the first trait to register its map
 
    _numChromosome = nChrm;
    _totalNumLoci = nLoc;
    _resolution = resolution;
    _totalLength = locPositions[nLoc-1];
    
    //destroy all previously allocated tables
    reset_tables();
    
    _perChrsmLength = new unsigned int[_numChromosome];
    _chrsmFirstLocusPosition = new unsigned int[_numChromosome];
    
  } else {  
    //we have already registered some trait's genetic map, we append this trait's info and update the map
    
    if (nChrm != _numChromosome) //check chromosome numbers
      fatal("Traits in init file don't have same number of chromosomes, cannot set the genetic map!\n");
 
    //number of loci
    _totalNumLoci = max(_totalNumLoci, nLoc);
    
    //setting the resolution:
    if (resolution < _resolution) {
      
//      cout << "GeneticMap::addTrait: rescaling map resolution (old "<< _resolution << "; new "<< resolution<< ")\n";
 
      rescaleMap(resolution);
      
    } else if (resolution > _resolution) {
      
//      cout << "GeneticMap::addTrait: rescaling locus position (old "<< _resolution << "; new "<< resolution<< "; ratio "<< resolution/_resolution <<")\n";
 
      double ratio = resolution/_resolution;
      
      for(unsigned int i = 0; i < nLoc; ++i)
        locPositions[i] = (unsigned int)((double)locPositions[i]*ratio);
    }
    
  }
  
  if(do_add) {  //only at the start of a simulation, not when loading pop from binary source
    
    traitIdx = _nTrait;
        
    _traits[trait] = _nTrait++;
  
    
    //---- record num chrmsm
    if ( _numChrsmPerTrait.size() != _nTrait -1 ) {
      fatal("Genetic map::wrong size of table of num chrms per trait (%i != %i)", 
            _numChrsmPerTrait.size(), _nTrait -1);
    } else
      _numChrsmPerTrait.push_back(nChrm);
    
    //---- record num loci
    if ( _numLociPerTrait.size() != _nTrait -1 ) {
      fatal("Genetic map::wrong size of table of num loci per trait (%i != %i)", 
            _numLociPerTrait.size(), _nTrait -1);
    } else
      _numLociPerTrait.push_back(nLoc);
    
    //---- record table of num loci per chrmsm
    if ( _numLociPerChrsmPerTrait.size() != _nTrait -1 )
      fatal("Genetic map::wrong size of table of num loci per chrms per trait (%i != %i)", 
            _numLociPerChrsmPerTrait.size(), _nTrait-1);
    
    unsigned int* locTable = new unsigned int[nChrm];
    
    for (unsigned int i = 0; i < nChrm; i++) locTable[i] = nLocChrm[i];
    
    _numLociPerChrsmPerTrait.push_back(locTable);
        
    //---- record the locus positions 
    if ( _locPositionsPerTrait.size() != _nTrait -1 ) {
      fatal("Genetic map::wrong size of table of loc position per trait (%i != %i)", 
            _locPositionsPerTrait.size(), _nTrait-1);
    }
    
    //copy the locus positions, they are at the correct resolution already
    unsigned int* posTable = new unsigned int[nLoc];
    
    for (unsigned int i = 0; i < nLoc; ++i) {
      posTable[i] = locPositions[i];
    }
    
    _locPositionsPerTrait.push_back(posTable);
    
    
    //---- create tables to record positions of x-over, for each gamete (male and female)
    
    for(unsigned int s = 0; s < 2; ++s){
      if ( _recPositions[s].size() != _nTrait -1 )
        fatal("Genetic map::wrong size of rec positions table (%i != %i)", 
              _recPositions[s].size(), _nTrait-1);
      
      _recPositions[s].push_back(vector<unsigned int>()); //add an empty array for each trait
    }
    
    
    //---- set the map positions according to other trait's
    //---- need to correctly set starting position and length of each chromosome
    
    if (_nTrait == 1) { //first trait to register
 
      for (unsigned int c = 0, stride = 0; c < _numChromosome; ++c) {
        _chrsmFirstLocusPosition[c] = locPositions[stride];
        _perChrsmLength[c] = locPositions[stride + nLocChrm[c] - 1] - locPositions[stride]; //position of last locus - position of first locus on that chromosome.
        stride += nLocChrm[c];
      }
      
    } else {
      
      for (unsigned int c = 0, lastPos, stride = 0; c < _numChromosome; ++c) {
        
        lastPos = max(locPositions[stride + nLocChrm[c] - 1], 
                      _chrsmFirstLocusPosition[c] + _perChrsmLength[c]);
        
        _chrsmFirstLocusPosition[c] = min(_chrsmFirstLocusPosition[c], locPositions[stride]);
        
        _perChrsmLength[c] = lastPos - _chrsmFirstLocusPosition[c];
        
        stride += nLocChrm[c];
      }
    }
    
    //---- compute total chromosome length
    _totalLength = 0;
    for (unsigned int c = 0; c < _numChromosome; ++c) _totalLength += _perChrsmLength[c];
    
    //---- mean num recombination events, map length of 1M = 1 x-over on average
    _totRecombEventsMean = (double)_totalLength * 0.01 * _resolution;
    
    //---- set the map size used to draw recombination spots; one less operation in recombine()...
    _recombLength = _totalLength + 1;
 
    //---- set the recombination rate between adjacent positions on the map;
    _recombinationRate = 0.01 * _resolution;
 
 
#ifdef _DEBUG_
     cout << "GeneticMap::addTrait: "<<endl;
     cout << "    map stats for trait \""<<trait<<"\"\n";
     cout << "   _totalLength         = "<<_totalLength<<endl;
     cout << "   _totRecombEventsMean = "<<_totRecombEventsMean<<endl;
#endif
 
     if(_totalLength > 9e+8)
       warning("Your genetic map is larger than %.2fcM with %i elements, you may experience memory allocation problems\
 as well as problems at runtime.\n",
           _totalLength*_resolution, _totalLength);
 
    //---- create new slot for the lookup table of that trait
    //     will be allocated and filled in setLookupTable
    if ( _lociLookupTable.size() != _nTrait -1 )
      fatal("Genetic map::wrong size of loci lookup tables (%i != %i)", 
            _lociLookupTable.size(), _nTrait-1);
    
    _lociLookupTable.push_back(NULL);
    
    //---- set all lookup tables (all genetic maps):
    //     we reset all tables each time a trait is added because _totalLength changes
    for (unsigned int t = 0; t < _nTrait; ++t) setLookupTable(t);
  }
  
  return traitIdx;
}

References _chrsmFirstLocusPosition, _lociLookupTable, _locPositionsPerTrait, _nTrait, _numChromosome, _numChrsmPerTrait, _numLociPerChrsmPerTrait, _numLociPerTrait, _perChrsmLength, _recombinationRate, _recombLength, _recPositions, _resolution, _totalLength, _totalNumLoci, _totRecombEventsMean, _traits, fatal(), rescaleMap(), reset_tables(), setLookupTable(), and warning().

Referenced by TTProtoWithMap::registerGeneticMap().

checkRegisteredTrait()

bool GeneticMap::checkRegisteredTrait

(

trait_t

trait

)

Returns true if trait 'trait' has registered a genetic map, false otherwise.

Parameters

trait

The trait name to check for in the trait table.

{
  map<trait_t, unsigned int>::iterator tIter = _traits.find(trait);
 
  return ( tIter != _traits.end() );
}

References _traits.

Referenced by TTQuantiFH::setOutputOption().

clear()

void GeneticMap::clear

(

)

{
  reset_tables();
  _traits.clear();
  _nTrait = 0;
}

References _nTrait, _traits, and reset_tables().

Referenced by IndFactory::clearPrototype().

getFirstRecPosition()

vector< bool >& GeneticMap::getFirstRecPosition ( sex_t  SEX )

inline

Returns the vector of the first chromosome position for recombination, used for all traits.

  { return _chrsmFirstRecombPosition[SEX]; }

References _chrsmFirstRecombPosition.

getGeneticMap()

bool GeneticMap::getGeneticMap	(	trait_t	trait,
		double **	table,
		unsigned int	table_length
	)

{
  map< trait_t, unsigned int>::const_iterator ITER = _traits.find(trait);
 
  if(ITER == _traits.end())
    return(error("GeneticMap::getGeneticMap::failed to find trait \"%s\" in genetic map\n", trait.c_str()));
 
  unsigned int idx = ITER->second;
 
  assert(table_length == _numLociPerTrait[idx]);
 
  unsigned int chr = 0, loc_cnt = _numLociPerChrsmPerTrait[idx][0];
 
  for(unsigned int i = 0; i < table_length; ++i) {
 
      if(i >= loc_cnt) {
        chr++;
        loc_cnt += _numLociPerChrsmPerTrait[idx][chr];
      }
 
      table[i][0] = chr;
      table[i][1] = _locPositionsPerTrait[idx][i];
 
  }
 
  return true;
}

References _locPositionsPerTrait, _numLociPerChrsmPerTrait, _numLociPerTrait, _traits, and error().

Referenced by TTNeutralGenesFH::write_PLINK(), and TTQuantiFH::write_PLINK().

getLocusPositionTable()

unsigned int* GeneticMap::getLocusPositionTable ( trait_t  trait )

inline

  {return _locPositionsPerTrait[ _traits[ trait ] ];}

References _locPositionsPerTrait, and _traits.

getRecLoci()

vector< unsigned int>& GeneticMap::getRecLoci ( sex_t  SEX, unsigned int  trait  )

inline

Returns a vector of the loci where crossing-overs take place.

  {
    return _recPositions[SEX][trait];
  }

References _recPositions.

getResolution()

double GeneticMap::getResolution ( )

inline

{return _resolution;}

References _resolution.

Referenced by TTProtoWithMap::recordRandomMap().

recombine()

void GeneticMap::recombine

(

sex_t

SEX

)

Called by TTProtoWithMap::recombine twice to create the two gametes necessary for the creation of a new individual.

This function sets the location of the x-over on the chromosome map and records the loci positions at which x-over happen for each trait separately. The positions are stored in GeneticMap::_recPositions, which are then passed to the trait-specific inheritance functions when setting up the genetics of the traits in new individuals.

Parameters

SEX	the origin of the gamete, whether from the father or the mother

{
  unsigned int nbRec = 0;
 
 
  //empty the array of x-over positions
  for (unsigned int t = 0; t < _nTrait; t++) {
    _recPositions[SEX][t].clear();
  }
  
  //create a gamete:
  
  //draw the number of x-over depending on the total map size
  nbRec = (unsigned int)RAND::Poisson( _totRecombEventsMean );
 
//  nbRec = (unsigned int)RAND::Binomial( _recombinationRate , _totalLength );
//  call to Poisson() is faster
 
#ifdef _REC_DEBUG_
  cout << "\n@@@ GeneticMap::recombine( SEX = "<<SEX<<" nbRec = "<<nbRec<<")"<<endl;
#endif
 
  //adjust
  if(nbRec > _totalNumLoci - 1) nbRec =  _totalNumLoci - 1; 
 
  //now determine where the x-over will take place
  if(nbRec > 0) {
 
    //reset the list of recombination events.
    _junctions.resize(nbRec);
 
    auto drawRecPosition = [this](unsigned int& n){n = (unsigned int)(RAND::Uniform()*_recombLength);};
 
    //draw the positions
    for_each(_junctions.begin(), _junctions.end(), drawRecPosition);
 
    //sort them
    sort(_junctions.begin(), _junctions.end());
 
#ifdef _REC_DEBUG_
    cout<< " --- sorted junctions ("<< _junctions.size()<<"): ";
    auto print = [this](unsigned int& n){cout<<" "<<n;};
    for_each(_junctions.begin(), _junctions.end(), print);
    cout<<endl;
#endif
    
    //remove the duplicates, we keep only one when that happens (in principle very rarely):
    auto last = unique(_junctions.begin(), _junctions.end());
 
 
    // UPDATE NUMBER OF RECOMBINATION EVENTS
 
    nbRec = last - _junctions.begin(); // these are the unique x-over map positions
 
#ifdef _REC_DEBUG_
      cout<<" --- after removing duplicates ("<< last - _junctions.begin()<<"): ";
      for_each(_junctions.begin(), last, print);
      cout<<endl;
 
      cout<<" --- computing rec map position: \n";
#endif
 
      size_t pos;
 
      for (unsigned int t = 0; t < _nTrait; t++) {
 
#ifdef _REC_DEBUG_
        cout<<"   trait "<< t << endl;
#endif
 
        for(unsigned int i = 0, hit; i < nbRec; ++i) {
 
          pos = _junctions[i];
 
          hit = _lociLookupTable[t][pos];
 
          _recPositions[SEX][t].push_back( hit );
 
 
#ifdef _REC_DEBUG_
          cout<< "  ["<<i<<"] @"<<pos<<" -> loc "<< hit <<"\n";
#endif
 
        }// end_for_bRec
      }// end_for_traits
  } // end_if_nbRec 
  
  // we always add the last position possible, to avoid returning an empty vector
  for (unsigned int t = 0; t < _nTrait; t++) {
    _recPositions[SEX][t].push_back(_numLociPerTrait[t]);
  }
  
  //set the starting point for each chromosome randomly (which chromosome copy to use first)
  
  _chrsmFirstRecombPosition[SEX].assign(_numChromosome, 0); //re-initialise the array
  
  for(unsigned int c = 0; c < _numChromosome; ++c) {
    _chrsmFirstRecombPosition[SEX][c] = RAND::RandBool();
  }
 
}

References _chrsmFirstRecombPosition, _junctions, _lociLookupTable, _nTrait, _numChromosome, _numLociPerTrait, _recombLength, _recPositions, _totalNumLoci, _totRecombEventsMean, RAND::Poisson(), RAND::RandBool(), and RAND::Uniform().

Referenced by TTProtoWithMap::recombine().

reduceJunctions()

vector< pair< unsigned int, unsigned int > > GeneticMap::reduceJunctions	(	sex_t	SEX,
		unsigned int	trait_idx
	)

Remove multiple x-over at the same locus when traits differ in number of loci.

{
  unsigned int prevLoc = 0, cEnd = 0;
  bool flipper;
  //number of x-overs
  vector< unsigned int >& recTable = _recPositions[SEX][trait_idx];
  unsigned int nbRec = recTable.size();
 
  vector< pair<unsigned int, unsigned int> > junctions;
 
  unsigned int num_copied = 0;
  
#ifdef _REC_DEBUG_
  cout << "GeneticMap::reduceJunctions; trait "<<trait_idx<< "; nb Rec = "<<nbRec - 1<<endl;
  cout << "   --sex = "<<SEX<<"\n";
#endif
 
  // c is the chromosome number
  // stride is the number of loci considered so-far, sum of chromosome loci
  // rec is the number of x-over done so-far
  for(unsigned int c = 0, cStart = 0, rec = 0; c < _numChromosome; ++c) {
 
    //the copy of the chromosome with which with start
    flipper = _chrsmFirstRecombPosition[SEX][c];
 
    //number of loci to consider, including previous chromosomes
    cEnd = cStart + _numLociPerChrsmPerTrait[trait_idx][c] ;
 
    //last locus at which a x-over happened, will be first position on current chromosome
    prevLoc = cStart;
 
    num_copied = 0;
 
#ifdef _REC_DEBUG_
    cout<<"   --chrm "<<c<<", start="<<cStart<<", end="<<cEnd<<" (rec="<<rec<<"; rec pos = "
        <<recTable[rec]<<")"<<", side="<<flipper<<endl;
#endif
 
    // zip forward if x-over on previous chrm were past the last locus on previous chrm
    // we need to catch up with position in the current chrm without changing the flipper
    while(recTable[rec] == cStart) rec++;
 
    // copy blocs of loci between x-over points on current chromosome
    // skip x-over if locus is not on this chromosome but a latter one (junction pos > cEnd)
    for(; recTable[rec] < cEnd && rec < nbRec; rec++) {
 
      //skip recombination events that would lead to copying 0 byte because even num x-over at same locus
      if(recTable[rec] == prevLoc){
        //but switch side for next bloc to copy on this chromosome
        //absolutely necessary to be in sync with other traits on the map
        flipper = !flipper;
        continue;
      }
 
 
#ifdef _REC_DEBUG_
      cout<< "   ["<<rec<<"] junction "<<recTable[rec]<< " on side "<<flipper<<std::endl;
#endif
 
      junctions.push_back(pair<unsigned int, unsigned int>(recTable[rec], flipper));
 
      num_copied += recTable[rec] - prevLoc;
 
      //update the starting locus to the next recombination point
      prevLoc = recTable[rec];
 
      //switch side for next bloc to copy on this chromosome
      flipper = !flipper;
    }
 
#ifdef _REC_DEBUG_
    cout << "   --copy til end of chrmsm from "<<prevLoc<<" to "<<cEnd<<" on side "<<flipper<<endl;
#endif
 
    num_copied += cEnd - prevLoc;
 
    junctions.push_back(pair<unsigned int, unsigned int>(cEnd, flipper));
 
    if(num_copied != _numLociPerChrsmPerTrait[trait_idx][c])
      error("number of loci copied on chromosome %i is %i != %i\n", c, num_copied, _numLociPerChrsmPerTrait[trait_idx][c]);
 
    cStart += _numLociPerChrsmPerTrait[trait_idx][c];
  }
 
  return junctions;
}

References _chrsmFirstRecombPosition, _numChromosome, _numLociPerChrsmPerTrait, _recPositions, and error().

Referenced by TProtoBDMI::inherit_low(), TProtoDeletMutations_bitstring::inherit_low(), TProtoQuanti::inherit_low(), and TProtoNeutralGenes::inherit_low().

registerIndForRecombine()

bool GeneticMap::registerIndForRecombine ( unsigned long  ID )

inline

Called by TTProtoWithMap::recombine with individual ID passed down from Individual::recombine.

Returns false when already called by the same individual. The function is called for each trait separately but recombination must be computed once per individual only. It is so because some functions create new individual's trait separately, for instance in breed_selection when traits under viability selection are created before neutral traits.

Parameters

ID	the ID number of the individual calling the recombination function

  {
    if (ID == _currentIndividual) return false;
    else  _currentIndividual = ID;
    return true;
  }

References _currentIndividual.

Referenced by TTProtoWithMap::recombine().

rescaleMap()

void GeneticMap::rescaleMap

(

double

val

)

{
  assert(val < _resolution);
  
  double ratio = _resolution/val; //it is assumed that val < _resolution
  
  _resolution = val;
  
  //we have to update all positions 
  for (unsigned int i = 0; i < _nTrait; ++i) {
    for (unsigned int j = 0; j < _numLociPerTrait[i]; ++j) {
      _locPositionsPerTrait[i][j] *= ratio;
    }
  }
  
  for (unsigned int i = 0; i < _numChromosome; ++i) {
    _perChrsmLength[i] *= ratio;
    _chrsmFirstLocusPosition[i] *= ratio;
  }
  
  _totalLength *= ratio;
  _totRecombEventsMean = (double)_totalLength * 0.01 * _resolution;
 
  //reset all lookup tables, the total length has changed!
  for (unsigned int i = 0; i < _nTrait; ++i) setLookupTable(i);
}

References _chrsmFirstLocusPosition, _locPositionsPerTrait, _nTrait, _numChromosome, _numLociPerTrait, _perChrsmLength, _resolution, _totalLength, _totRecombEventsMean, and setLookupTable().

Referenced by addTrait().

reset_tables()

void GeneticMap::reset_tables

(

)

{
  if(_perChrsmLength != NULL) delete [] _perChrsmLength;
  _perChrsmLength = NULL;
  
  if(_chrsmFirstLocusPosition != NULL) delete [] _chrsmFirstLocusPosition;
  _chrsmFirstLocusPosition = NULL;
  
  _numChrsmPerTrait.clear();
  _numLociPerTrait.clear();
  
  for (unsigned int i = 0; i < _lociLookupTable.size(); ++i) {
    if(_lociLookupTable[i] != NULL) delete [] _lociLookupTable[i];
    else error("GeneticMap::reset_tables::found null pointer in _lociLookupTable\n");
  }
  _lociLookupTable.clear();
  
  for (unsigned int i = 0; i < _numLociPerChrsmPerTrait.size(); ++i) {
    if(_numLociPerChrsmPerTrait[i] != NULL) delete [] _numLociPerChrsmPerTrait[i];
    else error("GeneticMap::reset_tables::found null pointer in _numLociPerChrsmPerTrait\n");
  }
  _numLociPerChrsmPerTrait.clear();
  
  for (unsigned int i = 0; i < _locPositionsPerTrait.size(); ++i) {
    if(_locPositionsPerTrait[i] != NULL) delete [] _locPositionsPerTrait[i];
    else error("GeneticMap::reset_tables::found null pointer in _locPositionsPerTrait\n");
  }
  _locPositionsPerTrait.clear();
  
  _recPositions[0].clear();
  _recPositions[1].clear();
  
  _chrsmFirstRecombPosition[0].clear();
  _chrsmFirstRecombPosition[1].clear();
  
}

References _chrsmFirstLocusPosition, _chrsmFirstRecombPosition, _lociLookupTable, _locPositionsPerTrait, _numChrsmPerTrait, _numLociPerChrsmPerTrait, _numLociPerTrait, _perChrsmLength, _recPositions, and error().

Referenced by addTrait(), clear(), and ~GeneticMap().

setLookupTable()

void GeneticMap::setLookupTable

(

unsigned int

idx

)

Bbuilds the lookup table for each trait.

A lookup table maps a chromosomal position to a locus so that when a x-over is placed at a given map position, the corresponding locus can be directly found. The size of the lookup table depends on the _totalLength and _resolution of the map. The lookup table for a 1M map at the 0.01 cM scale will have 10,000 elements.

Parameters

idx	the index of the trait in the table (stored in the trait prototype)

!!! +1 needed here because map must start with 0 !!!!

{  
  if(_lociLookupTable[idx] != NULL)  delete [] _lociLookupTable[idx];
 
  _lociLookupTable[idx] = new unsigned int [_totalLength +1];
  
  if(_lociLookupTable[idx] == NULL)
    fatal("GeneticMap::set lookup table: memory exhausted?\n");
  
  //------------------------------------------------------------------------
  int pre_offset = 0, post_offset = 0, offset;
 
//  cout<<"\n-- Genetic map #"<<idx<<":"<<endl;
  
  for (unsigned int loc = 0, pos = 0, stride = 0, 
       c = 0; c < _numChromosome; ++c) {
 
//    cout<<"Chrm "<<c+1<<" ("<<_perChrsmLength[c]<<" ["<<_resolution<<"cM]):\n"<<"{ ";
 
    loc = stride;
    
    stride += _numLociPerChrsmPerTrait[idx][c];
 
    //the map must start with position zero for the first locus with lowest position
    pre_offset = -_chrsmFirstLocusPosition[c];
 
    //chromosomes are contiguous 
    offset = pre_offset + post_offset;
          
    for(; loc < stride && loc < _numLociPerTrait[idx]; ++loc) 
    {
      
      while (pos <= ( _locPositionsPerTrait[idx][loc] + offset) &&
             ((int)pos - post_offset) <= (int)_perChrsmLength[c]) 
             //table has 1 more elmnt, hence "<="
      
            _lociLookupTable[idx][pos++] = loc;
 
//      cout << loc<<":"<<pos-1<<"("<<(int)pos - post_offset<<"; tab["<<pos-1<<"]="<<_lociLookupTable[idx][pos-1]<<") ";
    
    }
    
    post_offset += _perChrsmLength[c];
    
    // need to pad to the right when we reached last locus on a chromosome but not at the last map position of that chrm
    // this must apply only to all chrm but the last, hence loc < max num loci
    while( pos < post_offset && loc == stride && loc < _numLociPerTrait[idx])
    {
      _lociLookupTable[idx][pos++] = loc;
//      cout << loc<<":"<<pos-1<<"(tab["<<pos-1<<"]="<<_lociLookupTable[idx][pos-1]<<") ";
 
    }
 
    // we need to pad to the right when position of last locus is within last chromosome (chrom length set by other traits too)
    // pad with last locus number+1 (tot num loci) on that chromosome
    while (pos < _totalLength+1 && loc >= _numLociPerTrait[idx])
    {
    
      _lociLookupTable[idx][pos++] = _numLociPerTrait[idx]; 
 
//      cout << loc<<":"<<pos-1<<"(tab["<<pos-1<<"]="<<_lociLookupTable[idx][pos-1]<<") ";
      
    }
//    cout << "}\n";
  }
   
  
#ifdef _DEBUG_
//  map< trait_t, unsigned int>::iterator ITER = _traits.begin() + idx;
 
  cout << "\nGeneticMap::setLookupTable ("<<idx<<")\n";
  
  cout << "lookup table: [";
  for(unsigned int i = 0; i < _totalLength+1 && i < 101; ++i) {
    cout << _lociLookupTable[idx][i] << " ";
    if( i == 100) cout << "|>100) ";
  }
  if(_totalLength >= 100) cout << "... last pos ("<<_totalLength<<") "<<_lociLookupTable[idx][_totalLength];
  cout << "] (truncated if > 100 positions)\n";
  
  pre_offset = 0; post_offset = 0;
  for(unsigned int c = 0, stride = 0; c < _numChromosome; c++) {
    
    cout<<"+++Chrm "<<c+1<<" ("<<_perChrsmLength[c]<<" ["<<_resolution<<"cM]):\n";
    cout << "Loc positions: [ ";
    
    for(unsigned int i = 0; i < _numLociPerChrsmPerTrait[idx][c] && i < 101; ++i){
      if( i < 100 || !(i % 100) ) 
        cout << _locPositionsPerTrait[idx][i+stride] << " " ;
      if( i == 100) cout << "|>100) ";
    }
    cout << "]\n";
   
    pre_offset = -_chrsmFirstLocusPosition[c];
    offset = pre_offset + post_offset;
 
    cout << "lookup table: [ ";
    
    for(unsigned int cnt=0, i = _chrsmFirstLocusPosition[c]+offset; i < _chrsmFirstLocusPosition[c]+_perChrsmLength[c]+offset && cnt < 100; ++i, ++cnt){
      cout << _lociLookupTable[idx][i] << " " ;
    }
    cout << "] (first 100 positions only)\n{position:locus} {";
    
    for(unsigned int i = 0; i < _numLociPerChrsmPerTrait[idx][c]-1; ++i) {
      if( i < 100 || !(i % 100) ) cout<<_locPositionsPerTrait[idx][i+stride]+offset<<":"<<_lociLookupTable[idx][_locPositionsPerTrait[idx][i+stride]+offset]<<", ";
      if( i == 100) cout << "|>100) ";
    }
    cout<<_locPositionsPerTrait[idx][stride + _numLociPerChrsmPerTrait[idx][c]-1]+offset
        <<":"<<_lociLookupTable[idx][_locPositionsPerTrait[idx][stride + _numLociPerChrsmPerTrait[idx][c]-1]+offset]<<"}\n";
    
    post_offset += _perChrsmLength[c];
    stride += _numLociPerChrsmPerTrait[idx][c];
  }
  cout<<endl;
//    cout<<"Loc pos lookup table:"<<endl<<"{";
//    for(unsigned int i = 0; i < _superChrsmLength; ++i)
//      cout<<_recLociPositionTable[i]<<", ";
//    cout<<"}"<<endl;  
 
  cout<<"   Tot map length: "<<_totalLength<<endl;
  cout<<"   Mean num recombination events: "<<_totRecombEventsMean<<endl;
#endif
}

References _chrsmFirstLocusPosition, _lociLookupTable, _locPositionsPerTrait, _numChromosome, _numLociPerChrsmPerTrait, _numLociPerTrait, _perChrsmLength, _resolution, _totalLength, _totRecombEventsMean, and fatal().

Referenced by addTrait(), and rescaleMap().

setResolution()

double GeneticMap::setResolution ( double  val )

inline

  {
    _resolution = (val < _resolution ? val : _resolution); 
    return _resolution;
  }

References _resolution.

unregisterTrait()

void GeneticMap::unregisterTrait

(

trait_t

trait

)

{
  map<trait_t, unsigned int>::iterator tIter;
  
  tIter = _traits.find(trait);
    
  if ( tIter != _traits.end() )  {
    
//    cout<<"GeneticMap::unregisterTrait::"<<tIter->first<<" (idx: "<<tIter->second<<")"<<endl;
    
    _traits.erase(tIter);
    
    //@TODO trait lookup tables must be erased as well
 
//      cout<<"GeneticMap::unregisterTrait::done\n";
  }
  else fatal("Genetic map::unregisterTrait: trait \"%s\" is not registered\n", trait.c_str());
  
  _nTrait--;
}

References _nTrait, _traits, and fatal().

Referenced by TTProtoWithMap::reset(), and TTProtoWithMap::unregisterFromGeneticMap().

Member Data Documentation

_chrsmFirstLocusPosition

unsigned int* GeneticMap::_chrsmFirstLocusPosition

private

Referenced by addTrait(), rescaleMap(), reset_tables(), and setLookupTable().

_chrsmFirstRecombPosition

vector< bool > GeneticMap::_chrsmFirstRecombPosition[2]

private

Two vectors holding the starting copy of each chromosome to use when creating the two gametes that are used to create a new individual.

Referenced by getFirstRecPosition(), recombine(), reduceJunctions(), and reset_tables().

_currentIndividual

unsigned long GeneticMap::_currentIndividual

private

Referenced by registerIndForRecombine().

_junctions

vector< unsigned int > GeneticMap::_junctions

private

A vector to store the position of the recombination events.

Referenced by recombine().

_lociLookupTable

vector< unsigned int* > GeneticMap::_lociLookupTable

private

A list of tables that map the map position (cM) to a locus, for each trait.

The length of the table is the length of the genetic map, which depends on the map resolution (cM by default).

Referenced by addTrait(), recombine(), reset_tables(), and setLookupTable().

_locPositionsPerTrait

vector< unsigned int* > GeneticMap::_locPositionsPerTrait

private

Vector containing the table of map position for the loci of each trait.

Positions are recorded according to the minimum map resolution as used in the lookup table.

Referenced by addTrait(), getGeneticMap(), getLocusPositionTable(), rescaleMap(), reset_tables(), and setLookupTable().

_nTrait

unsigned int GeneticMap::_nTrait

private

Number of traits registered in the map.

Length of following arrays.

Referenced by addTrait(), clear(), recombine(), rescaleMap(), and unregisterTrait().

_numChromosome

unsigned int GeneticMap::_numChromosome

private

Referenced by addTrait(), recombine(), reduceJunctions(), rescaleMap(), and setLookupTable().

_numChrsmPerTrait

vector< unsigned int > GeneticMap::_numChrsmPerTrait

private

Vector of number of chromosomes for each trait.

Referenced by addTrait(), and reset_tables().

_numLociPerChrsmPerTrait

vector< unsigned int* > GeneticMap::_numLociPerChrsmPerTrait

private

Vector containing a table of number of loci per chromosome for each trait.

Referenced by addTrait(), getGeneticMap(), reduceJunctions(), reset_tables(), and setLookupTable().

_numLociPerTrait

vector< unsigned int > GeneticMap::_numLociPerTrait

private

Vector of number of loci for each trait.

Referenced by addTrait(), getGeneticMap(), recombine(), rescaleMap(), reset_tables(), and setLookupTable().

_perChrsmLength

unsigned int* GeneticMap::_perChrsmLength

private

Referenced by addTrait(), rescaleMap(), reset_tables(), and setLookupTable().

_recombinationRate

double GeneticMap::_recombinationRate

private

Referenced by addTrait().

_recombLength

unsigned int GeneticMap::_recombLength

private

Referenced by addTrait(), and recombine().

_recPositions

vector< vector < unsigned int > > GeneticMap::_recPositions[2]

private

Vector of tables containing, for each trait, the locus number at which x-overs happen.

Updated at each generation. There is matrix per sex/gamete with one trait per row and a variable number of elements in each row.

Referenced by addTrait(), getRecLoci(), recombine(), reduceJunctions(), and reset_tables().

_resolution

double GeneticMap::_resolution

private

Referenced by addTrait(), getResolution(), rescaleMap(), setLookupTable(), and setResolution().

_totalLength

unsigned int GeneticMap::_totalLength

private

Referenced by addTrait(), rescaleMap(), and setLookupTable().

_totalNumLoci

unsigned int GeneticMap::_totalNumLoci

private

Referenced by addTrait(), and recombine().

_totRecombEventsMean

double GeneticMap::_totRecombEventsMean

private

Referenced by addTrait(), recombine(), rescaleMap(), and setLookupTable().

_traits

map< trait_t, unsigned int > GeneticMap::_traits

private

Table mapping trait type to its position index in the following tables.

Referenced by addTrait(), checkRegisteredTrait(), clear(), getGeneticMap(), getLocusPositionTable(), and unregisterTrait().

---

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

• ttrait_with_map.h
• ttrait_with_map.cc