...
 
Commits (10)
......@@ -2,6 +2,9 @@
This file reflects the most important changes between releases. For further details, please chech the commits in the releases branches in gitlab.
## v2.4 (2020-05-31)
- Important improvements in performance
## v2.3 (2020-03-06)
- Fixes using restriction matrix for management
- Fixes using relationship matrix
......
FROM alpine:3.4
LABEL version="v2.3"
LABEL version="v2.4"
LABEL description="Metapop2 provides an analysis of gene and allelic diversity in subdivided populations from molecular genotype or coancestry data as well as a tool for the management of genetic diversity in conservation programs."
LABEL maintainer="[email protected]"
......
# Metapop2 (v2.3) [![build status](https://gitlab.com/elcortegano/metapop2/badges/master/pipeline.svg)](https://gitlab.com/elcortegano/metapop2/commits/master)
# Metapop2 (v2.4) [![build status](https://gitlab.com/elcortegano/metapop2/badges/master/pipeline.svg)](https://gitlab.com/elcortegano/metapop2/commits/master)
## Overview
......
......@@ -15,7 +15,7 @@
* --------------------------------------------------------------------------------------------
* Convert Metapop files into Genepop files
* ********************************************************************************************/
void convert2gnp ( std::string filename, std::vector<population>& deme ) {
void convert2gnp (std::string filename, std::vector<population>& deme, const allele2D& alleles) {
std::size_t loci (deme[0].getIndividual(0).getLociNumber());
std::ofstream genepop_file;
......@@ -33,31 +33,26 @@ void convert2gnp ( std::string filename, std::vector<population>& deme ) {
std::size_t digits (1);
allele_t max (0);
for (std::size_t p(0); p < deme.size(); ++p) {
if (DEBUG) std::cout << "\r Reading population " << p+1 << " / " << deme.size();
for (std::size_t i(0); i < deme[p].getPopulationSize(); ++i) {
for (std::size_t l(0); l < loci; ++l) {
if ( deme[p].getIndividual(i).getLocus(l).getA(0) > max ) {
max = deme[p].getIndividual(i).getLocus(l).getA(0);
}
if ( deme[p].getIndividual(i).getLocus(l).getA(1) > max ) {
max = deme[p].getIndividual(i).getLocus(l).getA(1);
}
for (std::size_t l(0);l < alleles.size(); ++l) {
if (DEBUG) std::cout << "\r Reading locus " << l+1 << " / " << alleles.size();
for (std::size_t a(0); a < alleles[l].size(); ++a) {
if (alleles[l][a] > max ) {
max = alleles[l][a];
}
}
}
while (max /= 10) ++digits;
if (DEBUG) std::cout << std::endl;
for (std::size_t p (0); p < deme.size(); p++){
for (std::size_t p (0); p < deme.size(); ++p) {
if (DEBUG) std::cout << "\r Writting population " << p+1 << " / " << deme.size();
genepop_file << "Pop" << std::endl;
for (std::size_t i(0); i<deme[p].getPopulationSize(); i++){
for (std::size_t i(0); i<deme[p].getPopulationSize(); ++i) {
individual ind (deme[p].getIndividual(i));
genepop_file << deme[p].getLabel() << p << i << ", ";
for (std::size_t l(0); l< loci; l++){
allele_t a (deme[p].getIndividual(i).getLocus(l).getA(0));
allele_t b (deme[p].getIndividual(i).getLocus(l).getA(1));
allele_t a (ind.getLocus(l).getA(0));
allele_t b (ind.getLocus(l).getA(1));
genepop_file << std::setfill('0') << setw(digits) << a << setw(digits) << b << "\t";
}
genepop_file << std::endl;
......@@ -163,6 +158,7 @@ void convert (const Metapop& input, const MTPConfigFile& config) {
std::string mode (config.get_convert_mode());
bool use_nucleotides (config.get_use_nucleotides());
int ploidy (config.get_ploidy());
allele2D alleles (input.get_alleles_per_locus());
if (mode == "mtp2gp") {
if (use_nucleotides) {
......@@ -171,7 +167,7 @@ void convert (const Metapop& input, const MTPConfigFile& config) {
}
std::string genepop_filename (caso + ".gen");
std::cout << "Converting " << input.get_filename() << " (Metapop format) into " << genepop_filename << " (Genepop format)" << std::endl;
convert2gnp (genepop_filename, deme);
convert2gnp (genepop_filename, deme, alleles);
} else if (mode == "mtp2ped") {
std::string ped_filename (caso + ".ped");
std::cout << "Converting " << input.get_filename() << " (Metapop format) into " << ped_filename << " (Plink PED format)" << std::endl;
......
......@@ -6,10 +6,10 @@
#include "output.h"
// ===== FUNCTION'S PROTOTYPES =====
void convert2gnp ( std::string , const std::vector<population>& );
void convert2ped ( std::string , const std::vector<population>& , bool , int );
void convert_deme2mtp ( std::string , const std::vector<population>& , bool );
void convert ( const Metapop& , const MTPConfigFile& );
void convert2gnp (std::string, const std::vector<population>&, const allele2D&);
void convert2ped (std::string, const std::vector<population>&, bool, int);
void convert_deme2mtp (std::string, const std::vector<population>&, bool);
void convert (const Metapop&, const MTPConfigFile&);
#endif
......@@ -33,7 +33,7 @@ individual::individual (std::string fname, unshort ID, unshort sex, const std::v
* --------------------------------------------------------------------------------------------
* List of methods to return protected attributes
* ********************************************************************************************/
locus individual::getLocus ( std::size_t index) const {
locus individual::getLocus (const std::size_t& index) const {
return loci[index];
}
......
......@@ -20,7 +20,7 @@ public:
individual ( std::string , unshort , unshort , const std::vector<locus>& );
locus getLocus ( std::size_t ) const;
locus getLocus ( const std::size_t& ) const;
std::size_t getLociNumber ();
std::string getLabel ();
int getSex ();
......
......@@ -23,7 +23,7 @@ locus::locus(const allele1D& as) {
this->check();
}
locus::locus(allele_t a1, allele_t a2) {
locus::locus(const allele_t& a1, const allele_t& a2) {
alleles.push_back(a1);
alleles.push_back(a2);
this->check();
......@@ -34,11 +34,11 @@ locus::locus(allele_t a1, allele_t a2) {
* --------------------------------------------------------------------------------------------
* Return information of alleles in the locus.
* ********************************************************************************************/
allele_t locus::getA (std::size_t n) const {
allele_t locus::getA (const std::size_t& n) const {
return alleles[n];
}
allele_t locus::getNAllels (allele_t a) const {
allele_t locus::getNAllels (const allele_t& a) const {
unshort nalel (0);
for (auto ai: alleles) if (a==ai) ++nalel;
return nalel;
......
......@@ -13,10 +13,10 @@ class locus {
public:
locus ( const std::vector<allele_t>& );
locus ( allele_t , allele_t );
allele_t getA ( std::size_t ) const;
locus ( const allele_t& , const allele_t& );
allele_t getA ( const std::size_t& ) const;
allele1D getAllels () const { return alleles; }
allele_t getNAllels ( allele_t ) const;
allele_t getNAllels ( const allele_t& ) const;
bool status () const;
protected:
......
......@@ -10,10 +10,10 @@
* files, and how the program works. *
* *
* *
* Version: 2.3 *
* Version: 2.4 *
* Created: 19/12/2011 *
* Updated: 06/03/2020 *
* Compiler: g++ (9.2.1) *
* Updated: 31/05/2020 *
* Compiler: g++ (9.3.0) *
* *
* Author: Eugenio López-Cortegano *
* Email: [email protected] *
......
......@@ -32,7 +32,7 @@ using namespace std;
/* GLOBAL VARIABLES */
const std::string version = "v2.3";
const std::string version = "v2.4";
std::mt19937 random_generator; // Mersenne Twister pseudorandom number generator engine
std::uniform_real_distribution<double> random_dist(0.0, 1.0);
bool DEBUG (false);
......
......@@ -50,6 +50,8 @@ Metapop::Metapop (std::string file_name, const MTPConfigFile& config) {
if (config.get_datafile_sformat()=="pedigree") read_pedigree = true;
// Read and format input file
bool compute_coancestry (true);
if (config.get_convert_mode() != "none") compute_coancestry = false;
use_nucleotides = config.get_use_nucleotides();
std::string data ((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
file.close();
......@@ -58,7 +60,7 @@ Metapop::Metapop (std::string file_name, const MTPConfigFile& config) {
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); // input_metapop
else this->update (data, compute_coancestry); // input_metapop
if (for_simulation) use_nucleotides = false;
// Read management file
......@@ -168,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 ();
}
......@@ -178,7 +180,7 @@ Metapop::Metapop ( const Metapop& original) { // Copy constructor for bootstrap
* --------------------------------------------------------------------------------------------
* Read a string with the content of the datafile, and update its attributes.
* ********************************************************************************************/
void Metapop::update (std::string file) {
void Metapop::update (std::string file, bool compute_coancestry) {
// Initialize variables
n_subpops = 0;
......@@ -199,11 +201,12 @@ void Metapop::update (std::string file) {
// Read metapop file
this->input_metapop(file);
//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
......@@ -239,6 +242,7 @@ void Metapop::compute_diversity ( ) {
/* CALCULATE INDIVIDUAL FREQUENCIES () */
if (use_molecular_markers) {
calculate_ind_freq ();
coancestry_from_freq();
}
/* CALCULATE GENETIC PARAMETERS */
......@@ -422,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();
......@@ -468,6 +473,42 @@ 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 tmpf += (*ind_freq[i][ind_i])[k][a] * (*ind_freq[j][ind_j])[k][a];
}
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);
}
}
}
/* ********************************************************************************************
......@@ -1161,9 +1202,6 @@ void Metapop::input_genepop (std::string str) {
std::cout << std::endl;
}
// Calculate coancestry
this->coancestry_matrix();
}
/* ********************************************************************************************
......@@ -1303,8 +1341,7 @@ void Metapop::input_pedfile (std::string str, bool calc_pedigree, std::string rm
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 ();
......
......@@ -148,7 +148,7 @@ public:
Metapop ( std::string , const MTPConfigFile& );
Metapop ( std::string , const MTPConfigFile& , std::string ); // user for conversion to metapop format
Metapop ( const Metapop& ); // copy constructor for bootstrap
void update ( std::string );
void update ( std::string , bool );
void check ( const MTPConfigFile& );
~Metapop ();
......@@ -210,7 +210,7 @@ protected:
// Input methods
void input_metapop ( std::string );
void input_genepop ( std::string );
void input_pedfile ( std::string, bool, 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 );
......
......@@ -778,7 +778,7 @@ void printSolution ( const MTPConfigFile& config, Metapop& metapop, const Calcul
}
//Establish the input for the new generation
if (simMode) metapop.update (dat.str());
if (simMode) metapop.update (dat.str(), true);
fr.close();
/* SAVE SIMULATION RESULTS */
......