Commit 5ccd9327 authored by E.L.Cortegano's avatar E.L.Cortegano

Improved performance of coancestry computation

parent a4f0f72e
Pipeline #151322048 passed with stages
in 3 minutes and 38 seconds
......@@ -58,8 +58,8 @@ Metapop::Metapop (std::string file_name, const MTPConfigFile& config) {
clean_datafile(data);
filename = file_name;
if (input_format == GNP) input_genepop (data, compute_coancestry);
else if (input_format == PED) input_pedfile (data, read_pedigree, config.get_rm_file(), compute_coancestry);
if (input_format == GNP) input_genepop (data);
else if (input_format == PED) input_pedfile (data, read_pedigree, config.get_rm_file());
else this->update (data, compute_coancestry); // input_metapop
if (for_simulation) use_nucleotides = false;
......@@ -119,8 +119,8 @@ Metapop::Metapop ( std::string file_name, const MTPConfigFile& config, std::stri
filename = file_name;
input_format = config.get_datafile_format();
if (mode == "gp2mtp") input_genepop (data, false);
else if (mode == "ped2mtp") input_pedfile (data, false,"", false);
if (mode == "gp2mtp") input_genepop (data);
else if (mode == "ped2mtp") input_pedfile (data, false,"");
else {
std::cerr << "Unknown conversion mode. Select one of 'none', 'mtp2gp' or 'mtp2ped'." << std::endl;
exit(-1);
......@@ -170,7 +170,7 @@ Metapop::Metapop ( const Metapop& original) { // Copy constructor for bootstrap
if (n_loci != tmp_loci) print_error (ERROR_MTP_BOOTL);
n_loci = original.n_loci;
clean_alleles ();
this->coancestry_matrix();
//this->coancestry_matrix();
set_struc ();
}
......@@ -201,12 +201,12 @@ void Metapop::update (std::string file, bool compute_coancestry) {
// Read metapop file
this->input_metapop(file);
if (!compute_coancestry) return;
//if (!compute_coancestry) return;
// Calculate coancestry matrix from molecular markers (if it hasn't been inputed)
// should this be part of metapop class?
if (use_molecular_markers) {
this->coancestry_matrix();
//this->coancestry_matrix();
}
// Set struc for management
......@@ -242,6 +242,7 @@ void Metapop::compute_diversity ( ) {
/* CALCULATE INDIVIDUAL FREQUENCIES () */
if (use_molecular_markers) {
calculate_ind_freq ();
coancestry_from_freq();
}
/* CALCULATE GENETIC PARAMETERS */
......@@ -425,9 +426,10 @@ void Metapop::read_fmatrix ( std::stringstream& str) {
/* ********************************************************************************************
* MTPInput:: coancestry ()
* --------------------------------------------------------------------------------------------
* Calculates the matrix of molecular coancestry between individuals
* Calculates the matrix of molecular coancestry between individuals from individual genotypes
* NOTE: This function is disabled, and coancestries are now computed from allele frequencies
* ********************************************************************************************/
void Metapop::coancestry_matrix () {
/*void Metapop::coancestry_matrix () {
if (DEBUG) std::cout << "Calculating coancestry matrix..." << endl;
f_mol.clear();
......@@ -471,6 +473,46 @@ void Metapop::coancestry_matrix () {
std::cout << "Done" << std::endl;
std::cout << std::endl;
}
}*/
/* ********************************************************************************************
* MTPInput:: coancestry_from_freq ()
* --------------------------------------------------------------------------------------------
* Calculates the matrix of molecular coancestry from allele frequencies
* ********************************************************************************************/
void Metapop::coancestry_from_freq () {
f_mol.clear();
for (std::size_t i(0); i < n_subpops; ++i) {
for (std::size_t ind_i (0), N_i (deme[i].getPopulationSize()); ind_i < N_i; ++ind_i) {
vdou1D f_vector;
for (std::size_t j(0); j<i+1; ++j) {
for (std::size_t ind_j(0), N_j(deme[j].getPopulationSize()); ind_j < N_j; ++ind_j) {
double f (0.0);
size_t eval_loci (n_loci);
for (std::size_t k(0); k<n_loci; ++k) {
bool all_loci (true);
double tmpf (0.0);
for (std::size_t a(0); a<alleles_per_locus[k].size(); ++a) {
if (!(*eval_ind_locus[i][ind_i])[k]) { all_loci = false; break; }
else if (!(*eval_ind_locus[j][ind_j])[k]) { all_loci = false; break; }
else if ((*ind_freq[i][ind_i])[k][a]== 0.0 || (*ind_freq[j][ind_j])[k][a] == 0.0) continue;
else if ((*ind_freq[i][ind_i])[k][a] == 1.0 && (*ind_freq[j][ind_j])[k][a] == 1.0) tmpf += 1.0;
else if ((*ind_freq[i][ind_i])[k][a] == 0.5 && (*ind_freq[j][ind_j])[k][a] == 0.5) tmpf += 0.25;
else if ((*ind_freq[i][ind_i])[k][a] == 0.5 || (*ind_freq[j][ind_j])[k][a] == 0.5) {
tmpf += 0.5;
}
}
if (!all_loci) --eval_loci; // If no allele, no calculation for that loci
else f += tmpf;
}
if (eval_loci) f /= eval_loci;
f_vector.push_back(f);
}
}
f_mol.push_back(f_vector);
}
}
}
/* ********************************************************************************************
......@@ -1057,7 +1099,7 @@ void Metapop::input_metapop (std::string file) {
* --------------------------------------------------------------------------------------------
* This function creates an intern and simplified metapop input file from a genepop input file
* ********************************************************************************************/
void Metapop::input_genepop (std::string str, bool compute_coancestry) {
void Metapop::input_genepop (std::string str) {
use_molecular_markers = true;
......@@ -1164,10 +1206,6 @@ void Metapop::input_genepop (std::string str, bool compute_coancestry) {
std::cout << std::endl;
}
// Calculate coancestry
if (!compute_coancestry) return;
this->coancestry_matrix();
}
/* ********************************************************************************************
......@@ -1175,7 +1213,7 @@ void Metapop::input_genepop (std::string str, bool compute_coancestry) {
* --------------------------------------------------------------------------------------------
* This function creates an intern and simplified metapop input file from a plink PED file
* ********************************************************************************************/
void Metapop::input_pedfile (std::string str, bool calc_pedigree, std::string rmatrix_filename, bool compute_coancestry) {
void Metapop::input_pedfile (std::string str, bool calc_pedigree, std::string rmatrix_filename) {
if (!calc_pedigree) use_molecular_markers = true;
else {
......@@ -1304,12 +1342,10 @@ void Metapop::input_pedfile (std::string str, bool calc_pedigree, std::string rm
}
// Test pedigree
if (!compute_coancestry) return;
if (calc_pedigree) test_pedigree(pedigree);
// Calculate coancestry
if (!calc_pedigree) this->coancestry_matrix();
else this->genealogical_matrix(pedigree);
if (calc_pedigree) this->genealogical_matrix(pedigree);
// Set struc and prgeny for management
set_struc ();
......
......@@ -209,8 +209,8 @@ protected:
// Input methods
void input_metapop ( std::string );
void input_genepop ( std::string , bool );
void input_pedfile ( std::string, bool, std::string , bool );
void input_genepop ( std::string );
void input_pedfile ( std::string, bool, std::string );
void default_mng ();
void input_management ( std::string );
void end_body ( std::size_t& , std::size_t& , std::size_t& , std::string& , std::vector<individual>& );
......@@ -247,7 +247,8 @@ protected:
allele_t get_allele ( std::string& );
// Coancestry methods
void coancestry_matrix ();
//void coancestry_matrix ();
void coancestry_from_freq ();
void genealogical_matrix ( const vint2D& );
double MeanCoancestryInd ( int , int , int , int );
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment