Commit 287f1ab7 authored by Etienne Klein's avatar Etienne Klein
Browse files

t push -u origin master

Codes associated to the Gauzere et al. 2016 Molecular Ecology Resources paper
parents
File added
##### Codes-Script Article MolEcolRes ######
- "code C++" directory : contains all the C++ codes (.cpp, .o, .hpp). NB: the parameters defining the mating system properties of the ancestral and current populations (alpha, Npc) are declared in the "Simul.cpp" code.
A Code::Blocks project (.cbp) is also provided as a possibility to generate the exe file on different platforms
Note that several functions from the Boost library are used in the C++ code
- "wrapper6.py" : python script that defines almost all the input parameters (Nfam, Noff, Nrepet, Nloci...). The script also implements the random draw of the allelic frequencies. The compiled C++ code is called from the Python application.
- "Simul.exe" : exemple of compiled C++ code
- "Script_AsREML_Method_PhiAmatrix.as" : AsREML script to analyse the simulated maternal progenies using an animal model
\ No newline at end of file
!LOGFILE !QUIET !WORKSPACE 200
simulations to estimate heritability from an estimated relatedness PhiA_matrix
Ind !A # 1st colomn of phenotype data = individual label
Phenotype !M0 # 2nd colomn of phenotype data = phenotypic value
PhiA_matrix.grm !NSD !DENSE # relationship matrix in dense format and non inversible
Pheno_repet1.asd !Skip1 # Data file with phenotypic values and fixed/random effects - first line is headings
Phenotype ~ mu !r giv(Ind,1) # animal model with no fixed effects
!PIN !define # define the contents of the output file - here calculates the total phenotypic variance and narrow-sense heritability
F phenvar 1 + 2
F addvar 1 * 4
H herit 1 3
\ No newline at end of file
File added
This diff is collapsed.
#ifndef CURRENTPOP
#define CURRENTPOP
#include "Indiv.hpp"
#include "Locus.hpp"
#include "Random.hpp"
#include "Params.hpp"
#include<vector> // SDM: je vire tous les vector contenus dans cette classe
#include <boost/math/distributions/gamma.hpp>
#include <boost/random/lagged_fibonacci.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/math/special_functions/round.hpp>
using namespace std;
class CurrentPop {
// SDM: je mets le default constructeur en prive pour verifier qu'il ne sert pas... si le compilateur se plaint de CurrentPop( ) is private ... in this context a viendra de l
CurrentPop();
public:
CurrentPop(int,int,int,int,float,int,vector<vector< float > > ,Random*,float,float ,float,int/*,float*/,vector<float> ); //le dernier float est pour FreqQuant
~CurrentPop();
Indiv** pop;
Indiv** popF;
Locus*** allelequant;
void Generation(boost::variate_generator<boost::lagged_fibonacci19937&, boost::uniform_real<> >& tirage,
boost::math::gamma_distribution<>& gamma, int Ngroup);
void Ech(boost::variate_generator<boost::lagged_fibonacci19937&, boost::uniform_real<> >& tirage,
boost::math::gamma_distribution<>& gamma, int Ngroup);
void PrintBreedVal();
private:
int length; // devient membre de la classe pour connaitre la taille des tableaux
int lengthQuant;
inline int round(float);
float GlobalSum; // SDM sum est sauvegarde pour pas avoir a recalculer les frequences
float* MargSum;
float Loiselle(int,int);
float Ritland96(int,int);
float APSbis(int,int);
float GlobalSumQuant;
float** FreqMatrixQuant;
float* MargSumQuant;
int Npop,Nfam,Noff,Nloc,acc,nullall,Ngeneration;
int NlocQuant,NallQuant;
float s,sigmaA,sigmaE,sigmaM;
float** FreqMatrix;
float** KinshipMatrix;
float** tpKinshipMatrix;
float** KinFinal;
Indiv** popCache;
int g;
ofstream genofileQuant, genofile, pedigree, phenodata, AddValues, loiselfile, loiselpatfile, freqall, phiAfile, phiAdense;
Random* random;
};
#endif
/***************************************************
*Classe Indiv:
***************************************************/
#include<iostream>
#include "Indiv.hpp"
Indiv::Indiv() {
Mere=0;
Pere=0;
// this->Nloc=Nloc;
// this->NlocQuant=NlocQuant;
// this->NallQuant=NallQuant;
breedVal=0; //?
phenVal=0; //?
Genotype[0]= new int[Nloc];
Genotype[1]= new int[Nloc];
GenotypeQuant[0]= new int[NlocQuant];
GenotypeQuant[1]= new int[NlocQuant];
// AllelicEffect = new double [NlocQuant][NallQuant];
}
Indiv::~Indiv(){//on libre la mmoire alloue pour les objets une fois qu'ils sont dtruits
delete[] Genotype[0];
delete[] Genotype[1];
delete[] GenotypeQuant[0];
delete[] GenotypeQuant[1];
}
int Indiv::Nloc=0; //on dfinit le membre statique Nloc pour allouer la mmoire ncessaire
// son stockage (bien que le membre statique ne dpende en ralit pas de ce scope)
int Indiv::NlocQuant=0;
void Indiv::SetGenotypes(int allele,int chromosome,int locus){
Genotype[chromosome][locus]=allele;
}
void Indiv::SetGenotypesQuant(int allele,int chromosome,int locus){
GenotypeQuant[chromosome][locus]=allele;
}
/*int Indiv::getGenotypes(int chromosome, int locus){
return Genotype[chromosome][locus];
}*/
void Indiv::videtonsac(){//on pourra virer tout a quand a marchera
for(int i=0;i<Nloc;i++) cout<<this<<" "<<i<<" "<<Genotype[0][i]<<" "<<Genotype[1][i]<<endl;
for(int i=0;i<NlocQuant;i++) cout<<this<<" "<<i<<" "<<GenotypeQuant[0][i]<<" "<<GenotypeQuant[1][i]<<endl;
}
/* void Indiv::SetAllelicEffect(int locus,int allele){
AllelicEffect[locus][allele]=random->nrand(0.,sigmaA);
} */
void Indiv::SetbreedVal(double value){
breedVal=value;
}
double Indiv::getbreedVal(){
return breedVal;
}
void Indiv::SetphenVal(double value){
phenVal=value;
}
double Indiv::getphenVal(){
return phenVal;
}
#ifndef INDIV
#define INDIV
using namespace std;
class Indiv {
public:
Indiv();
~Indiv();
void SetGenotypes(int,int,int);
void SetGenotypesQuant(int,int,int);
// void SetAllelicEffect(int,int);
void SetbreedVal(double);
void SetphenVal(double);
double getbreedVal();//accesseurs
double getphenVal();//accesseurs
//int getGenotypes(int,int);
int Indexfather;
int Indexmother;
int* Genotype[2];// mis public pour eviter de trs nombreux appels un accesseur
int* GenotypeQuant[2];
// double* AllelicEffect;
Indiv* Mere;
Indiv* Pere;
void videtonsac();
static int Nloc;
static int NlocQuant;
// static int NallQuant;
private:
double breedVal;//la breeding value de chaque individu
double phenVal;//valeur phno de l'individu
};
#endif
/***************************************************
*Classe Locus:
***************************************************/
#include<iostream>
#include "Locus.hpp"
using namespace std;
Locus::Locus() {
// this->NlocQuant=NlocQuant;
AllelicEffect = 0;
}
Locus::~Locus() {
// delete m_param;
}
void Locus::SetAllelicEffect(int allele){
AllelicEffect=random->nrand(0.,(sigmaA));
}
float Locus::GetAllelicEffect() const {
return AllelicEffect;
}
float Locus::sigmaA = 0.;
//int Locus::NlocQuant = 0.;
Random* Locus::random = NULL;
#ifndef LOCUS
#define LOCUS
#include "Random.hpp"
//#include "Params.hpp"
using namespace std;
class Locus {
public:
Locus(); // Constructeur
~Locus();
// static int NlocQuant;
// float AllelicEffect;
// static int NallQuant;
static float sigmaA;
// static int NlocQuant;
static Random* random;
void SetAllelicEffect(int allele);
float GetAllelicEffect() const;
private:
// int NlocQuant, NallQuant;
// Params *NallQuant;
float AllelicEffect;
};
#endif
/***************************************************
*Classe Init permet d'initialiser le jeu de donnes*
* partir duquel on va simuler. *
***************************************************/
#include<iostream>
#include<fstream>
#include <vector>
using namespace std;
#include "Params.hpp"
#include <cstdlib>
Params::Params(int argn, char** argv) {
switch (argn) {
case 1: // pas d'argument: affichage de l'aide
usage();
//inutile de breaker
case 2:
arguments=argv;
readfile(argv[1]);
checkparameters();
break;
/* case 10: // bon nombre d'arguments: lecture et verification
arguments = argv;
readparameters();
checkparameters();
break; */
default: // nombre incorrect d'arguments
cout<<"Ligne de commande incorrecte\n"<<endl;
exit(-1);
}
}
int Params::round(float f) {
return (int)(f + 0.5);
}
void Params::usage() {
printf("*********\n");
printf("[Simul] : usage\n");
printf("Type [executable name] sigmaA sigmaE Npop Nloc s Ngenerations Ngeno NlocQuanti\n");
printf("or [executable name] Parameter file\n");
printf("\tsigmaA = standard deviation of allelic effects\n");
printf("\tsigmaE = error standard deviation\n");
printf("\tNlocQuant = number of quantitative loci (positive integer)\n");
// printf("\tFreqAllQuant = frequency of allele A at quantitative locus\n");
printf("\tNallQuant = number of alleles by locus quanti\n");
printf("\tNpop = population size (positive integer)\n");
printf("\tNfam = nb familles sampled (positive integer)\n");
printf("\tNoff = nb offspring per familles sampled (positive integer)\n");
printf("\tNloc = number of loci (positive integer)\n");
printf("\ts = selfing rate (positive real or null)\n");
printf("\tNgenerations = integer larger or equal to 1\n");
printf("\tNall = number of alleles by locus (type N if all individuals have\n\tdifferent alleles in the base population)\n");
// printf("**********\n");
// printf("\t***************\n"); //allele freq for quant loci
system("PAUSE");
exit(1);
}
/*
// litparametres(): convertit les arguments de la ligne de commande
void Params::readparameters() {
int i;
float f,p;
sscanf(arguments[1],"%f",&sigmaA);
sscanf(arguments[2],"%f",&sigmaE);
NlocQuant = atoi(arguments[3]);
NallQuant = atoi(arguments[4]);
Npop = atoi(arguments[5]);
Nloc = atoi(arguments[6]);
sscanf(arguments[7],"%f",&s);
Ngeneration=atoi(arguments[8]);
Ngeno=atoi(arguments[9]);
// sscanf(arguments[10],"%f",&FreqAllQuant);
for (i=0;i<Ngeno;i++) {
cin>>f;
Freq.push_back(f);
//tant que le bon nombre de frequence n'est pas saisi, le programme attend
}
for (i=0;i<NallQuant;i++) { //idem pour freq locus quanti
cin>>p;
FreqQuant.push_back(p);
}
}
*/
void Params::readfile(char* infile) {
ifstream file(infile);
if (!file.is_open()) {cerr<<"cannot open infile"<<endl;exit(-1);}
next(file); file >> sigmaA;
next(file); file >> sigmaE;
next(file); file >> sigmaM;
next(file); file >> NlocQuant;
next(file); file >> NallQuant;
next(file); file >> Npop;
next(file); file >> Nfam;
next(file); file >> Noff;
next(file); file >> Nloc;
next(file); file >> s;
next(file); file >> Ngeneration;
next(file); file >> Ngeno;
// next(file); file >> FreqAllQuant;
// Freqloc.resize(Ngeno);//il faut faire a sinon par dfaut il est vide!
Freq.resize(Nloc,vector<float>(Ngeno));
for(int i=0;i<Nloc;i++) {
for(int j=0;j<Ngeno;j++) {
next(file);
file >> Freq[i][j];
cout<<"locus "<<i<<"allele"<<j<<"freq"<<Freq[i][j]<<endl;
}
}
// file.ignore(8192, '\n');
FreqQuant.resize(NallQuant);
for(int i=0;i<NallQuant;i++) {
next(file);
file >> FreqQuant[i];
}
// for(int i=0;i<NallQuant;i++) cout<<FreqQuant[i]<<"\t";
file.close();
}
void Params::next(ifstream& infile) {
while (infile.get()!='#') if(infile.peek()==EOF){
cerr<<"problem with parameter file format"<<endl;exit(-1);
}
}
// verifie la validite des parametres
void Params::checkparameters() {
int i,j;
float somme=0;
if (Npop>10000) {cerr<<"Error: population too big\n";exit(-1);}
if (Npop<Nfam) {cerr<<"Error: too many families sampled\n";exit(-1);}
// if (FreqAllQuant>1) {cerr<<"Error: Frequencies must be within[0-1]\n";exit(-1);}
if (s<0 || s>1) {cerr<<"Error: incorrect selfing rate\n";exit(-1);}
if (Ngeneration==0) {cerr<<"Error: Number of generations must be at least 1\n";exit(-1);}
for(i=0;i<Nloc;i++) {
for(j=0;j<Ngeno;j++) {
if (Freq[i][j]<0||Freq[i][j]>1) {cerr<<"Error: Frequencies must be within[0-1]\n";exit(1);}
// if ((round(Freq[i][j]*Npop)-Freq[i][j]*Npop)>0.000001||(round(Freq[i][j]*Npop)-Freq[i][j]*Npop)<-0.000001) {cerr<<"erreur: frequences ne doivent pas exceder precision (1/Npop)"<<endl;exit(-1);}
somme+=Freq[i][j];
}
}
if ((Nloc-somme)>0.0001||(Nloc-somme)<-0.0001) {cerr<<"Error: Frequencies must sum to 1 \n"<<endl;exit(-1);}
}
void Params::save(ofstream& file) const {
int i,j;
file << "sigmaA;" << sigmaA << endl;
file << "sigmaE;" << sigmaE << endl;
file << "sigmaM;" << sigmaM << endl;
file << "Number of loci quanti;" << NlocQuant << endl;
file << "Number of alleles by locus quanti;" << NallQuant << endl;
file << "Number of individuals;" << Npop << endl;
file << "Number of families sampled;" << Nfam << endl;
file << "Number of offspring per families sampled;" << Noff << endl;
file << "Number of loci;" << Nloc << endl;
file << "Selfing rate;" << s << endl;
file << "Number of generations;" << Ngeneration << endl;
file << "Number of alleles by locus;" << Ngeno << endl;
// file << "freq A quanti #" << FreqAllQuant << endl;
file << "Allele frequencies;";
for (i=0;i<Nloc;i++) {
file <<"(";
for (j=0;j<Ngeno;j++) file << Freq[i][j]<<"_";
file << ")"<<"\t";
}
file << endl;
file << "Allele frequencies quanti;";
for (unsigned int i=0;i<FreqQuant.size();i++) file << FreqQuant[i]<<"_";
}
#ifndef PARAMS
#define PARAMS
#include<vector>
#include <iostream>
class Params {
public:
Params(int, char**);
int Npop,Nfam,Noff,Nloc,Ngeno,Ngeneration,NlocQuant,NallQuant;
float s,sigmaA,sigmaE, sigmaM/*,FreqAllQuant*/;
// vector<float> Freqloc;
vector<vector< float > > Freq;
vector<float> FreqQuant;
void save(ofstream&) const;
private:
char** arguments;
void usage();
void readparameters();
void readfile(char*);
void next(ifstream&);
void checkparameters();
int round(float);
};
#endif
#include<iostream>
using namespace std;
#include<fstream>
#include "Random.hpp"
#include<cmath>
Random::Random(){
default_seed_1 = 123456789;
default_seed_2 = 987654321;
filename ="seeds.txt";
FILE* file;
if (!(file = fopen(filename,"r")) || fscanf(file,"%lf %lf",&S1,&S2) !=2) {
S1 = default_seed_1;
S2 = default_seed_2;
}
if (file) fclose(file);
}
/* ifstream file7("seeds.txt");
int S1=-1; S2=-1;
if (file7.eof()) { cerr << "no seeds" << endl; exit(1); }
file7 >>S1;
if (file7.eof()) { cerr << "miss one seed" << endl; exit(1); }
file7 >>S2;
if (S1<0 || S2<0) { cerr << "invalid values for seeds" << endl; exit(1); }
//stockage des graines dans l'objet
this->S1=S1;
this->S2=S2;
}*/
Random::~Random() {
FILE* file=fopen(filename,"w");
if (!file) cerr<<"cannot open seed file\n";
printSeeds(file);
fclose(file);
}
double Random::uniform (){
double z, r;
int i;
//printf("seeds %f %f\n",S1,S2); //si besoin de savoir les seeds chaque boucle
r = S1 / 53668;
i = (int) r;//prend la partie entire de r
r = (double) i;//on remet le nombre en double prcision
S1 = 40014 * (S1 - r * 53668) - r * 12211; // redfinition de la premire graine
if (S1 < 0)
S1 += 2147483563;
r = S2 / 52774;
i = (int) r;
r = (double) i;
S2 = 40692 * (S2 - r * 52774) - r * 3791;//redfinition de la deuxime graine
if (S2 < 0)
S2 += 2147483399;
z = S1 - S2;
if (z < 1)
z += 2147483562;
return (z * 4.656613e-10);
}
int Random::irand(int x) {
return (int) (uniform()*x);
}
double Random::erand(double esperance) {
double tp = 0.0;
while (tp==0.0) tp = uniform();
return ( -(esperance)*log(tp));
}
double Random::nrand(double mean,double std){
const int NbTirages = 12; //augmenter pour une meilleure precision. 12 est bien.
//increase for better precision. 12 works fine.
double valeur = 0;
for(int i=0 ; i < NbTirages ; ++i) valeur += uniform();
//on recentre la somme / centering the sum
valeur -= NbTirages/2;
//on etale suivant l'ecartype / spread with standard deviation
//le 12 n'a rien a voir avec NbTirages, mais explique pourquoi justement, on prend souvent
//NbTirages = 12
//the 12 is not related to NbTirages, but it explains why it is often chosen that NbTirages=12
valeur *= (NbTirages == 12) ? std
: sqrt(12/static_cast<double>(NbTirages))*std;
//on centre sur la moyenne / debias
valeur += mean;
return valeur;
}
void Random::printSeeds(FILE* file) {
fprintf(file,"%ld %ld\n\n",(long)S1,(long)S2);
}
#ifndef RANDOM
#define RANDOM
class Random {
public:
Random();
~Random();
double uniform();
double erand(double);
int irand(int);//pas besoin de drand ici on utilisera uniform
double nrand(double,double);
void printSeeds(FILE*);
private:
double S1, S2;//les seeds qui servent l'algorithme
long default_seed_1;
long default_seed_2;
const char* filename;
};
#endif
#include <iostream>
#include <fstream>
#include <math.h>
#include <vector>
#include <boost/math/distributions/gamma.hpp>
#include <boost/random/lagged_fibonacci.hpp>