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

void	placeRecombinationEvents (sex_t SEX, vector< unsigned int > junctions, const vector< bool > &firstRecPos)
	Deterministic counterpart of recombine(sex_t): instead of drawing the x-over map positions and per-chromosome start sides at random, they are supplied by the caller. More...

void	placeRecombinationEvents (sex_t SEX)
	Shared no-allocation core of recombine()/placeRecombinationEvents(): maps the x-over positions already stored in the member _junctions to per-trait locus indices (dedup + lookup-table mapping + trailing-gap/boundary skip + sentinel), filling _recPositions[SEX]. 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 (const 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

101 {reset_tables();}

GeneticMap::reset_tables

void reset_tables()

Definition: ttrait_with_map.cc:1404

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;
     
   }
   
 #if defined(_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;
  
  
 #if defined(_DEBUG_) or defined(_REC_DEBUG_)
      cout << "GeneticMap::addTrait: "<<endl;
      cout << "    map stats for trait \""<<trait<<"\" at "<<traitIdx<<"\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);
   }
  
 #if defined(_DEBUG_) or defined(_REC_DEBUG_)
      cout << "GeneticMap::addTrait: "<<endl;
      cout << "  traits in table:\n";
      map<trait_t, unsigned int>::const_iterator tt_it = _traits.begin();
      while(tt_it != _traits.end()) {
        cout << "   ["<<tt_it->second<<"] "<<tt_it->first<<endl;
        tt_it++;
      }
 #endif
   
   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().

Here is the caller graph for this function:

◆ 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 TTProtoWithMap::is_mapped(), and TTQuantiFH::setOutputOption().

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.

191 { return _chrsmFirstRecombPosition[SEX]; }

GeneticMap::_chrsmFirstRecombPosition

vector< bool > _chrsmFirstRecombPosition[2]

Two vectors holding the starting copy of each chromosome to use when creating the two gametes that ar...

Definition: ttrait_with_map.h:78

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[0][i] = chr;
       table[1][i] = _locPositionsPerTrait[idx][i];
  
   }
  
   return true;
 }

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

Referenced by GenotyperFH::FHwrite(), GenotyperFH::print_PLINK(), TTNeutralGenesFH::write_PLINK(), and TTQuantiFH::write_PLINK().

Here is the caller graph for this function:

◆ getLocusPositionTable()

unsigned int* GeneticMap::getLocusPositionTable ( const trait_t trait )

inline

194 {return _locPositionsPerTrait[ _traits[ trait ] ];}

References _locPositionsPerTrait, and _traits.

Referenced by TTProtoWithMap::get_locus_map_positions().

Here is the caller graph for this function:

◆ 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

105 {return _resolution;}

References _resolution.

Referenced by TTProtoWithMap::recordRandomMap().

Here is the caller graph for this function:

◆ placeRecombinationEvents() [1/2]

void GeneticMap::placeRecombinationEvents ( sex_t SEX )

Shared no-allocation core of recombine()/placeRecombinationEvents(): maps the x-over positions already stored in the member _junctions to per-trait locus indices (dedup + lookup-table mapping + trailing-gap/boundary skip + sentinel), filling _recPositions[SEX].

It does NOT set _chrsmFirstRecombPosition: each caller sets the per-chromosome start copies itself (recombine() draws them, the deterministic overload copies them from firstRecPos).

 {
   //empty the array of x-over positions
   for (unsigned int t = 0; t < _nTrait; t++) {
     _recPositions[SEX][t].clear();
   }
  
   //_junctions are assumed in range [0,_totalLength]: the random path (recombine) draws them
   //that way by construction and only calls this when nbRec>0; the deterministic entry point
   //clamps untrusted input before calling. So no per-gamete bounds scan happens here.
  
   unsigned int nbRec = _junctions.size();
  
   //now determine where the x-over will take place
   if(nbRec > 0) {
  
     //sort the positions
     sort(_junctions.begin(), _junctions.end());
  
 #ifdef _REC_DEBUG_
     cout<< " --- sorted junctions ("<< _junctions.size()<<"): ";
     auto print = [](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;
  
       // chromosome of a contiguous map position. A boundary coordinate B_c is shared by
       // chromosome c-1's max-position locus and chromosome c's min-position locus;
       // setLookupTable fills B_c during chromosome c-1's loop (its bound is inclusive:
       // pos <= _perChrsmLength[c-1]), so B_c belongs to the LOWER chromosome. We must use
       // the SAME convention here ("<=", not "<"): otherwise a crossover landing exactly on
       // B_c is classified as chromosome c here but mapped to chromosome c-1's last locus by
       // the lookup, and traits disagree on which chromosome it belongs to -> desync at the
       // sparse trait's first locus of the next chromosome.
       auto chromOfPos = [this](size_t p) -> unsigned int {
         unsigned int end = 0;
         for (unsigned int c = 0; c < _numChromosome; ++c) {
           end += _perChrsmLength[c];
           if (p <= end) return c;
         }
         return (_numChromosome ? _numChromosome - 1 : 0);
       };
  
       for (unsigned int t = 0; t < _nTrait; t++) {
  
         // chromosome of a locus index for trait t; the sentinel region (loc past the
         // last locus of the last chromosome) returns _numChromosome.
         auto chromOfLocus = [this, t](unsigned int loc) -> unsigned int {
           unsigned int end = 0;
           for (unsigned int c = 0; c < _numChromosome; ++c) {
             end += _numLociPerChrsmPerTrait[t][c];
             if (loc < end) return c;
           }
           return _numChromosome;
         };
  
 #ifdef _REC_DEBUG_
         cout<<"   trait "<< t << endl;
 #endif
  
         for(unsigned int i = 0, hit; i < nbRec; ++i) {
  
           pos = _junctions[i];
  
           hit = _lociLookupTable[t][pos];
  
           // A crossover that falls on chromosome c but past this trait's LAST locus on c
           // (its "trailing gap") maps to cStart(c+1), i.e. a locus on a LATER chromosome
           // than the crossover itself. Such a crossover is a no-op for this trait and must
           // NOT be recorded: otherwise reduceJunctions' zip-forward would mis-account it as
           // a flip of chromosome c+1 and decouple this trait from the others across that
           // whole chromosome. (Mirror of the leading-gap case handled in reduceJunctions.)
           if ( chromOfLocus(hit) > chromOfPos(pos) ) {
 #ifdef _REC_DEBUG_
             cout<< "  ["<<i<<"] @"<<pos<<" -> loc "<< hit <<" SKIPPED (trailing gap on chrm "
                 << chromOfPos(pos) <<")\n";
 #endif
             continue;
           }
           _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]);
   }
   
  
 }

References _junctions, _lociLookupTable, _nTrait, _numChromosome, _numLociPerChrsmPerTrait, _numLociPerTrait, _perChrsmLength, and _recPositions.

◆ placeRecombinationEvents() [2/2]

void GeneticMap::placeRecombinationEvents	(	sex_t	SEX,
		vector< unsigned int >	junctions,
		const vector< bool > &	firstRecPos
	)

Deterministic counterpart of recombine(sex_t): instead of drawing the x-over map positions and per-chromosome start sides at random, they are supplied by the caller.

The SAME dedup + per-trait lookup-table mapping + sentinel logic is then applied, so the resulting GeneticMap::_recPositions are identical to what recombine() would have produced for those same draws. Used by tests/tools to reproduce an exact gamete.

Parameters

SEX	the origin of the gamete
junctions	crossover positions in contiguous map coordinates, in 0,_totalLength
firstRecPos	which parental chromosome copy each chromosome starts with (size must be _numChromosome)

 {
   //adopt the caller-supplied positions into the reusable member.
   _junctions.swap(junctions);
  
   //sanitise untrusted input: clamp out-of-range positions (the lookup table has
   //_totalLength+1 elements, so the largest valid index is _totalLength). The random
   //hot path skips this because its draws are in range by construction.
   for (unsigned int& n : _junctions) {
     if (n > _totalLength) {
       warning("GeneticMap::placeRecombinationEvents: junction position %i out of range "
               "(map length %i), clamping.\n", n, _totalLength);
       n = _totalLength;
     }
   }
  
   //run the shared core (same dedup + lookup + sentinel mapping as the random hot path).
   placeRecombinationEvents(SEX);
  
   //set the starting copy for each chromosome from the supplied values
   _chrsmFirstRecombPosition[SEX].assign(_numChromosome, 0); //re-initialise the array
   for(unsigned int c = 0; c < _numChromosome && c < firstRecPos.size(); ++c) {
     _chrsmFirstRecombPosition[SEX][c] = firstRecPos[c];
   }
 }

References _chrsmFirstRecombPosition, _junctions, _numChromosome, _totalLength, and warning().

Referenced by recombine().

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;
  
  
   
   //create a gamete:
   
   //draw the number of x-over depending on the total map size
   nbRec = (unsigned int)RAND::Poisson( _totRecombEventsMean );
  
  
 #ifdef _REC_DEBUG_
   cout << "\n@@@ GeneticMap::recombine( SEX = "<<SEX<<" nbRec = "<<nbRec<<")"<<endl;
 #endif
  
   //adjust
   if(nbRec > _totalNumLoci - 1) nbRec =  _totalNumLoci - 1; 
  
   if(nbRec > 0) {
  
     //draw where the x-over will take place (contiguous map positions) into the reusable
     //_junctions buffer; resize() keeps the previously allocated capacity, so no per-gamete
     //heap allocation happens after warm-up. The draw is in range [0,_totalLength] by
     //construction (Uniform() < 1, _recombLength == _totalLength + 1), so no clamp is needed.
     _junctions.resize(nbRec);
     for (unsigned int& n : _junctions)
       n = (unsigned int)(RAND::Uniform() * _recombLength);
  
     //map the drawn positions to per-trait loci (clears _recPositions, dedup + lookup, sentinel)
     placeRecombinationEvents(SEX);
  
   } else {
  
     //common case (no crossover, the majority of gametes on a short map): just clear and
     //append the end sentinel for each trait, skipping the buffer/dedup/lookup machinery and
     //the extra function call entirely.
     for (unsigned int t = 0; t < _nTrait; t++) {
       _recPositions[SEX][t].clear();
       _recPositions[SEX][t].push_back(_numLociPerTrait[t]);
     }
   }
  
   //draw the starting copy for each chromosome straight into the destination (no temp vector)
   _chrsmFirstRecombPosition[SEX].assign(_numChromosome, 0);
   for(unsigned int c = 0; c < _numChromosome; ++c)
     _chrsmFirstRecombPosition[SEX][c] = RAND::RandBool();
 }

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

Referenced by TTProtoWithMap::recombine().

Here is the caller graph for this function:

◆ 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++; flipper = !flipper;}
  
     // 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(), TTNeutralGenes_bitstring::inherit_low(), TProtoQuanti::inherit_low(), and TProtoNeutralGenes::inherit_low().

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

 {  
   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<<" num loci: "<<_numLociPerChrsmPerTrait[idx][c]
          <<" ("<<_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: (from first pos:"<<_chrsmFirstLocusPosition[c]+offset<<") [ ";
     
     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)\nlookup table: {position:locus} (first 100 + 100th) (offset = "<<offset<<", stride = "<<stride<<")\n   {";
     
     for(unsigned int i = 0; i < _numLociPerChrsmPerTrait[idx][c]; ++i) {
       if( i < 100 || !(i % 100) )
         cout<<_locPositionsPerTrait[idx][i+stride]+offset<<":"<<_lociLookupTable[idx][_locPositionsPerTrait[idx][i+stride]+offset]<<", ";
       if( i == 100) cout << "|>100) ";
     }
     // last locus
     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(), placeRecombinationEvents(), 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 placeRecombinationEvents(), and 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(), placeRecombinationEvents(), 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(), placeRecombinationEvents(), recombine(), rescaleMap(), and unregisterTrait().

◆ _numChromosome

unsigned int GeneticMap::_numChromosome

private

Referenced by addTrait(), placeRecombinationEvents(), 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(), placeRecombinationEvents(), reduceJunctions(), reset_tables(), and setLookupTable().

◆ _numLociPerTrait

vector< unsigned int > GeneticMap::_numLociPerTrait

private

Vector of number of loci for each trait.

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

◆ _perChrsmLength

unsigned int* GeneticMap::_perChrsmLength

private

Referenced by addTrait(), placeRecombinationEvents(), 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(), placeRecombinationEvents(), 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(), placeRecombinationEvents(), 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:

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ GeneticMap()

◆ ~GeneticMap()

Member Function Documentation

◆ addTrait()

◆ checkRegisteredTrait()

◆ clear()

◆ getFirstRecPosition()

◆ getGeneticMap()

◆ getLocusPositionTable()

◆ getRecLoci()

◆ getResolution()

◆ placeRecombinationEvents() [1/2]

◆ placeRecombinationEvents() [2/2]

◆ recombine()

◆ reduceJunctions()

◆ registerIndForRecombine()

◆ rescaleMap()

◆ reset_tables()

◆ setLookupTable()

◆ setResolution()

◆ unregisterTrait()

Member Data Documentation

◆ _chrsmFirstLocusPosition

◆ _chrsmFirstRecombPosition

◆ _currentIndividual

◆ _junctions

◆ _lociLookupTable

◆ _locPositionsPerTrait

◆ _nTrait

◆ _numChromosome

◆ _numChrsmPerTrait

◆ _numLociPerChrsmPerTrait

◆ _numLociPerTrait

◆ _perChrsmLength

◆ _recombinationRate

◆ _recombLength

◆ _recPositions

◆ _resolution

◆ _totalLength

◆ _totalNumLoci

◆ _totRecombEventsMean

◆ _traits