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...
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)
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
 
 
    //---- 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(), and setLookupTable().

Referenced by TTProtoWithMap::registerGeneticMap().

Here is the caller graph for this function:

clear()

void GeneticMap::clear

(

)

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

References _nTrait, _traits, and reset_tables().

Referenced by IndFactory::clearPrototype().

Here is the caller graph for this function:

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.

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

Here is the caller graph for this function:

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("failed to find trait \"%s\" in genetic map"), 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]*_resolution;
 
  }
 
  return true;
}

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

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

Here is the caller graph for this function:

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.

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

Here is the caller graph for this function:

getResolution()

double GeneticMap::getResolution ( )

inline

{return _resolution;}

References _resolution.

Referenced by TTProtoWithMap::recordRandomMap().

Here is the caller graph for this function:

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
 
//  cout << "\n@@@ GeneticMap::recombine( SEX = "<<SEX<<" nbRec = "<<nbRec<<")"<<endl;
 
  //adjust
  if(nbRec > _totalNumLoci - 1) nbRec =  _totalNumLoci - 1; 
  
  //reset the list of recombination events.
  _junctions.resize(nbRec);
 
  //now determine where the x-over will take place
  if(nbRec > 0) {
 
    auto drawRecPosition = [this](unsigned int& n){n = (unsigned int)(RAND::Uniform()*_recombLength);};
//    auto print = [this](const unsigned int& n){cout<<" "<<n;};
 
    //draw the positions
    for_each(_junctions.begin(), _junctions.end(), drawRecPosition);
 
    //sort them
    sort(_junctions.begin(), _junctions.end());
 
//    cout<< " --- sorted junctions ("<< _junctions.size()<<"): ";
//    for_each(_junctions.cbegin(), _junctions.cend(), print);
 
    //remove the duplicates (they cancel each other):
//    unique(_junctions.begin(), _junctions.end());
 
//    cout<<"\n --- after removing duplicates ("<< _junctions.size()<<"): ";
//    for_each(_junctions.cbegin(), _junctions.cend(), print);
 
 
    //the (automatically sorted) table of < locus number, num of x-over at locus position >
    map< unsigned int, unsigned int > cumul_rec_positions;
    
    //the iterator to read within the map
    map< unsigned int, unsigned int>::const_iterator recIter;
    
//    cout<<"\n --- computing rec map position: \n";
    //set up trait recombination spots, more than one x-over may happen b/n two loci of a given trait
    //two consecutive x-overs b/n two loci cancel each other
    //we need to count the number of consecutive x-overs for each trait separately
    //the use of a map to store the x-over positions insures that they will be sorted
 
    nbRec = _junctions.size();
 
    for (unsigned int t = 0; t < _nTrait; t++) {
      
      cumul_rec_positions.clear();
      
//      cout<<"   trait "<< t << endl;
 
      for(unsigned int i = 0, hit, num_interleaved_xover; i < nbRec; ++i) {
      
        
        // RAND::Uniform() : draw position of x-over
        // _totalLength+1 : the +1 is necessary because the lookupTable (the map) starts with position 0
        //      position 0 should never be drawn as it would suppress recombination b/n two first loci
        //      it is safe so because the lookupTable has _totalLength +1 elements
        // _lociLookupTable : gives the locus to the right of that 'hit' position
        
        size_t pos = _junctions[i];
        
        hit = _lociLookupTable[t][pos];
        
//        cout<< "   @"<<pos<<" -> loc "<< hit <<"\n";
        
        num_interleaved_xover = cumul_rec_positions[hit];
        //returns 0 if position hit was not yet recorded
        
        cumul_rec_positions[hit] = ++num_interleaved_xover; 
        //indicates that one more x-over happened b/n two loci
        
       }
//      cout<<"   recording rec locus position: ";
      // record only the positions at which an odd number of x-over happened:
      for (recIter = cumul_rec_positions.begin(); recIter != cumul_rec_positions.end(); recIter++) {
        
        if(recIter->second % 2) {
 
          _recPositions[SEX][t].push_back(recIter->first); //add locus position
        
//        cout<<recIter->first<<", ";
        }
      }
//      cout<<endl;
    }// 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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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 !!!!

!!! this means mumLocus may be returned, adjust _recPositions accordingly !!!!

{  
  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"<<"{ ";
    if(c > 0) 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];
    
    while (pos < _totalLength && loc >= _numLociPerTrait[idx]) 
    {
    
      _lociLookupTable[idx][pos++] = _numLociPerTrait[idx]; 
      //(i.e. must contain num loci + 1 elements)!!!!
 
 
//      cout << loc<<":"<<pos-1<<"(tab["<<pos-1<<"]="<<_lociLookupTable[idx][pos-1]<<") ";
      
    }
//    cout << "}\n";
  }
   
  
#ifdef _DEBUG_
  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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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(), 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(), 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(), 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(), and reset_tables().

_resolution

double GeneticMap::_resolution

private

Referenced by addTrait(), getGeneticMap(), 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(), clear(), getGeneticMap(), getLocusPositionTable(), and unregisterTrait().

---

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

• ttrait_with_map.h
• ttrait_with_map.cc