343 lines
11 KiB
C++
343 lines
11 KiB
C++
#include "DerevyankoReport2014.h"
|
|
#include "individual/Trait.h"
|
|
#include "individual/Phenotype.h"
|
|
//#include "population/BreedingStrategies/PopulationBreedingStrategy.h"
|
|
#include "population/BreedingStrategies/NeutralEvolutionBreedStrat.h"
|
|
#include "processor/Settings.h"
|
|
|
|
#include <list>
|
|
#include <fstream>
|
|
#include <cstdio>
|
|
#include <ctime>
|
|
#include <cstdlib>
|
|
#include <boost/lexical_cast.hpp>
|
|
#include <boost/algorithm/string.hpp>
|
|
#include <boost/tokenizer.hpp>
|
|
|
|
using std::string;
|
|
|
|
std::string DerevyankoReport2014::getInitGene(int type, int num){
|
|
|
|
string seq, initSeq = "";
|
|
srand((unsigned int)time(NULL));
|
|
for(unsigned int i = 0; i < Settings::locusLength; i++){
|
|
int nucl = rand()%4;
|
|
if(nucl == 0){initSeq += 'A';}
|
|
else if(nucl == 1){initSeq += 'T';}
|
|
else if(nucl == 2){initSeq += 'G';}
|
|
else if(nucl == 3){initSeq += 'C';}
|
|
}
|
|
//initSeq = "TTCTTTCATGGGGAAGCAGATTTGGGTACCACCCAAGTATTGACTCACCCATCAACAACC";
|
|
|
|
if(type == 0){
|
|
// Ìèòîõîíäðèàëüíàÿ ÄÍÊ (ïîïóëÿöèÿ "ëþäåé")
|
|
seq = initSeq;
|
|
}
|
|
else{
|
|
// Ìóòàíòíàÿ ìòÄÍÊ (ïîïóëÿöèÿ íåàíäåíèñîâöåâ)
|
|
seq = initSeq;
|
|
seq.replace(10,1,"T");
|
|
seq.replace(20,1,"A");
|
|
seq.replace(30,1,"C");
|
|
seq.replace(40,1,"G");
|
|
seq.replace(50,1,"G");
|
|
}
|
|
|
|
return seq;
|
|
}
|
|
|
|
BisexualPopulation* DerevyankoReport2014::populationFactory(int size, const std::vector<std::string>& genome){
|
|
BisexualPopulation* ans;
|
|
|
|
std::list<Individual*> males;
|
|
std::list<Individual*> females;
|
|
|
|
// Ïîäãîòîâêà èíèöèàòîðíîãî ãåíîìà
|
|
Chromosome chrom1("Mitochondrial chromosome");
|
|
for(unsigned int i = 0; i < genome.size(); i++){
|
|
string gName = "mtDna_" + boost::lexical_cast<std::string>(i+1);
|
|
Gene gene(Gene::Sequence, gName, genome.at(i));
|
|
chrom1.insertGeneToEnd(gene);
|
|
}
|
|
|
|
std::vector<Chromosome> fGenome;
|
|
std::vector<Chromosome> mGenome;
|
|
fGenome.push_back(chrom1);
|
|
mGenome.push_back(chrom1);
|
|
|
|
// Ôåíîòèï
|
|
Trait trait1(Trait::Discrete, "Age", 20);
|
|
// Ïóë ñóáñòðàòîâ
|
|
InnerSubstratesPool* subPool = 0;
|
|
// (END) Ïóë ñóáñòðàòîâ
|
|
|
|
for(int i = 0; i < size/2; i++){
|
|
Genotype* genotype = new Genotype(fGenome, mGenome);
|
|
Phenotype* phenotype = new Phenotype(trait1);
|
|
males.push_back(new Individual(genotype, phenotype, subPool, Individual::male, 0));
|
|
|
|
genotype = new Genotype(fGenome, mGenome);
|
|
phenotype = new Phenotype(trait1);
|
|
females.push_back(new Individual(genotype, phenotype, subPool, Individual::female, 0));
|
|
}
|
|
|
|
ans = new BisexualPopulation(males,females);
|
|
return ans;
|
|
}
|
|
|
|
void DerevyankoReport2014::report2014_01(int argc, char** argv){
|
|
std::vector<std::string> calcScript;
|
|
if(argc > 1){
|
|
calcScript = parseScript(argv[1]);
|
|
}
|
|
else{
|
|
calcScript = parseScript("script.txt");
|
|
}
|
|
// Ãåíåðàöèÿ ñòàðòîâûõ ïîïóëÿöèé
|
|
int numLoci = Settings::numLoci;
|
|
std::vector<string> initMtDNA(numLoci);
|
|
std::vector<string> initMtDNA_ancient(numLoci);
|
|
for(unsigned int i = 0; i < initMtDNA.size(); i++){
|
|
initMtDNA.at(i) = getInitGene(0,i);
|
|
if(i<= numLoci*Settings::percentDiffLoci){
|
|
string seq = initMtDNA.at(i);
|
|
for(int j = 0; j < 5; j++){
|
|
int pos = (int)(j/5.*Settings::locusLength);
|
|
if(seq.at(pos)=='T')
|
|
seq.replace(pos,1,"G");
|
|
else
|
|
seq.replace(pos,1,"T");
|
|
}
|
|
initMtDNA_ancient.at(i) = seq;
|
|
}
|
|
else{
|
|
initMtDNA_ancient.at(i) = initMtDNA.at(i);
|
|
}
|
|
}
|
|
|
|
int popSizeH = Settings::PopulationHomoInitSize;
|
|
int popSizeA = Settings::PopulationAncientInitSize;
|
|
float initBirthRate = Settings::BirthRate;
|
|
float deathRate = Settings::NaturalDeathRate;
|
|
|
|
BisexualPopulation* pop_Sapiens = populationFactory(popSizeH, initMtDNA);
|
|
NeutralEvolutionBreedingStrategy* breedStrat_Sapiens =
|
|
dynamic_cast<NeutralEvolutionBreedingStrategy*>(PopulationBreedingStrategy::getInstance("neutral"));
|
|
breedStrat_Sapiens->setBirthRate(initBirthRate);
|
|
breedStrat_Sapiens->setDeathRate(deathRate);
|
|
pop_Sapiens->setBreedingStrategy(breedStrat_Sapiens);
|
|
|
|
BisexualPopulation* pop_Ancient = populationFactory(popSizeA, initMtDNA_ancient);
|
|
NeutralEvolutionBreedingStrategy* breedStrat_Ancient =
|
|
dynamic_cast<NeutralEvolutionBreedingStrategy*>(PopulationBreedingStrategy::getInstance("neutral"));
|
|
breedStrat_Ancient->setBirthRate(initBirthRate);
|
|
breedStrat_Ancient->setDeathRate(deathRate);
|
|
pop_Ancient->setBreedingStrategy(breedStrat_Ancient);
|
|
|
|
///////////////////////////////////////////
|
|
//
|
|
// Ýâîëþöèîííûé ïðîöåññ
|
|
//
|
|
///////////////////////////////////////////
|
|
//Settings::ProbMtDNARecomb = 1e-4;
|
|
|
|
// Øàã 1.
|
|
// Âàðèàíò 1
|
|
int genPerMigr_H_A = Settings::generationsPerMigration;
|
|
float part_H_to_A = Settings::migrationPartHomoToAncient;
|
|
float part_A_to_H = Settings::migrationPartAncientToHomo;
|
|
|
|
unsigned int generation = 0;
|
|
for(unsigned int curLine = 0; curLine < calcScript.size(); curLine++){
|
|
if(calcScript.at(curLine).at(0) == '/' || calcScript.at(curLine).at(0) == '#')
|
|
continue;
|
|
|
|
std::vector<std::string > tokensVector; // Search for tokens
|
|
boost::split(tokensVector, calcScript.at(curLine), boost::is_any_of("=") );
|
|
std::string instr = tokensVector.at(0);
|
|
boost::to_lower(instr);
|
|
boost::trim(instr);
|
|
|
|
// Èíñòðóêöèè
|
|
|
|
if(instr == "evolution"){
|
|
int iter = 0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &iter);
|
|
for(unsigned int i = 0; i < iter; i++, generation++){
|
|
std::cout<<"gen\t"<<generation<<std::endl;
|
|
pop_Sapiens->breedAll();
|
|
pop_Sapiens->mutationAll();
|
|
|
|
pop_Ancient->breedAll();
|
|
pop_Ancient->mutationAll();
|
|
}
|
|
}
|
|
else if(instr == "evolution_m"){
|
|
int iter = 0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &iter);
|
|
for(unsigned int i = 0; i < iter; i++, generation++){
|
|
std::cout<<"gen\t"<<generation<<std::endl;
|
|
pop_Sapiens->breedAll();
|
|
pop_Sapiens->mutationAll();
|
|
|
|
pop_Ancient->breedAll();
|
|
pop_Ancient->mutationAll();
|
|
if(i != 0 && ((i % genPerMigr_H_A) == 0)){
|
|
std::string code = BisexualPopulation::mutualMigration(pop_Sapiens, pop_Ancient, part_H_to_A, part_A_to_H);
|
|
std::cout<<"Migration Sapienti<->Ancient\t"<<code<<"\tIter\t"<<generation<<std::endl;
|
|
}
|
|
}
|
|
}
|
|
else if(instr == "migration"){
|
|
double part_H_to_A = 0.0;
|
|
double part_A_to_H = 0.0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg %lg", &part_H_to_A, &part_A_to_H);
|
|
std::string code = BisexualPopulation::mutualMigration(pop_Sapiens, pop_Ancient, part_H_to_A, part_A_to_H);
|
|
std::cout<<"Migration Sapienti<->Ancient\t"<<code<<"\tIter\t"<<generation<<std::endl;
|
|
}
|
|
} // (END)for(unsigned int curLine = 0; curLine < calcScript.size(); curLine++)
|
|
|
|
|
|
/* for(generation = 0; generation < 5000; generation++){
|
|
std::cout<<"gen\t"<<generation<<std::endl;
|
|
pop_Sapiens->breedAll();
|
|
pop_Sapiens->mutationAll();
|
|
|
|
pop_Ancient->breedAll();
|
|
pop_Ancient->mutationAll();
|
|
|
|
if(generation != 0 && ((generation % genPerMigr_H_A) == 0)){
|
|
std::string code = BisexualPopulation::mutualMigration(pop_Sapiens, pop_Ancient, part_H_to_A, part_A_to_H);
|
|
std::cout<<"Migration Sapienti<->Ancient\t"<<code<<std::endl;
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
// Êîíåö ýâîëþöèè - âûâîä ñòàòèñòèêè â ôàéëû
|
|
std::ofstream file("GenStatistics.pop_Sapiens.fasta");
|
|
pop_Sapiens->putGeneticStatisticsToStream(file);
|
|
file.close();
|
|
|
|
file.open("GenStatistics.pop_Sapiens.txt");
|
|
pop_Sapiens->putGeneticSimpleStatisticsToStream(file);
|
|
file.close();
|
|
|
|
////////////////////////////////////////
|
|
file.open("GenStatistics.pop_Ancient.fasta");
|
|
pop_Ancient->putGeneticStatisticsToStream(file);
|
|
file.close();
|
|
|
|
file.open("GenStatistics.pop_Ancient.txt");
|
|
pop_Ancient->putGeneticSimpleStatisticsToStream(file);
|
|
file.close();
|
|
////////////////////////////////////////
|
|
|
|
}
|
|
|
|
std::vector<std::string> DerevyankoReport2014::parseScript(std::string fileName){
|
|
std::vector<std::string> script;
|
|
std::ifstream file(fileName.c_str());
|
|
std::string str;
|
|
|
|
while(getline(file,str)){
|
|
if(str.length() > 2)
|
|
script.push_back(str);
|
|
}
|
|
file.close();
|
|
|
|
// process Script
|
|
unsigned int curLine;
|
|
for(curLine = 0; curLine < script.size(); curLine++){
|
|
if(script.at(curLine).at(0) == '/' || script.at(curLine).at(0) == '#')
|
|
continue;
|
|
|
|
std::vector<std::string > tokensVector; // Search for tokens
|
|
boost::split(tokensVector, script.at(curLine), boost::is_any_of("=") );
|
|
|
|
std::string instr = tokensVector.at(0);
|
|
boost::to_lower(instr);
|
|
boost::trim(instr);
|
|
|
|
// Èíñòðóêöèè
|
|
|
|
if(instr == "model_start"){
|
|
vector<std::string>::const_iterator first = script.begin() + curLine;
|
|
vector<std::string>::const_iterator last = script.end();
|
|
std::vector<std::string> ans(first, last);
|
|
return ans;
|
|
}
|
|
else if(instr == "pop_homo_init"){
|
|
int flag = 0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &flag);
|
|
if(flag != 0)
|
|
Settings::PopulationHomoInitSize = flag;
|
|
}
|
|
else if(instr == "pop_ancient_init"){
|
|
int flag = 0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &flag);
|
|
if(flag != 0)
|
|
Settings::PopulationAncientInitSize = flag;
|
|
}
|
|
else if(instr == "generations_per_migration"){
|
|
int flag = Settings::generationsPerMigration;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &flag);
|
|
if(flag != 0)
|
|
Settings::generationsPerMigration = flag;
|
|
}
|
|
else if(instr == "migr_homo_ancient_ratio"){
|
|
double ratio = Settings::migrationPartHomoToAncient;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &ratio);
|
|
Settings::migrationPartHomoToAncient = ratio;
|
|
}
|
|
else if(instr == "migr_ancient_homo_ratio"){
|
|
double ratio = Settings::migrationPartAncientToHomo;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &ratio);
|
|
Settings::migrationPartAncientToHomo = ratio;
|
|
}
|
|
else if(instr == "birth" || instr == "birthrate" || instr == "birth_rate" ){
|
|
double rate = Settings::BirthRate;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &rate);
|
|
Settings::BirthRate = rate;
|
|
}
|
|
else if(instr == "death" || instr == "deathrate" || instr == "death_rate" ){
|
|
double rate = Settings::NaturalDeathRate;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &rate);
|
|
Settings::NaturalDeathRate = rate;
|
|
}
|
|
else if(instr == "percent_diff_loci" ){
|
|
double value = Settings::percentDiffLoci;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &value);
|
|
Settings::percentDiffLoci = value;
|
|
}
|
|
else if(instr == "num_loci"){
|
|
int flag = 0;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%d", &flag);
|
|
if(flag != 0)
|
|
Settings::numLoci = flag;
|
|
}
|
|
else if(instr == "prob_recomb" ){
|
|
double value = Settings::ProbMtDNARecomb;
|
|
sscanf_s(tokensVector.at(1).c_str(),"%lg", &value);
|
|
Settings::ProbMtDNARecomb = value;
|
|
}
|
|
|
|
|
|
|
|
//*****************************
|
|
else if(instr == "breedings"){
|
|
//long int kMax = Settings::KMaxParameter;
|
|
//sscanf_s(tokensVector.at(1).c_str(),"%ld", &kMax);
|
|
std::string strategy = tokensVector.at(1);
|
|
boost::to_lower(strategy);
|
|
boost::trim(strategy);
|
|
Settings::BreedingStrategy = strategy;
|
|
}
|
|
} // (END) for(curLine = 0; curLine < script.size(); curLine++)
|
|
std::vector<std::string> a;
|
|
return a;
|
|
}
|
|
|
|
//void DerevyankoReport2014::parseCalc(std::vector<std::string> script){
|
|
//
|
|
//}
|