#include "DerevyankoReport2015.h" #include "individual/Trait.h" #include "individual/Phenotype.h" //#include "population/BreedingStrategies/PopulationBreedingStrategy.h" #include "population/BreedingStrategies/NeutralEvolutionBreedStrat.h" #include "processor/Settings.h" #include #include #include #include #include #include #include #include using std::string; std::string DerevyankoReport2015::getInitGene(int type, int num){ // type - тип человека (0 - кроманьонец, 1 - условный неандерталец) // num - номер гаплогруппы (0-4 - для кроманьонцев, 5-9 - для неандертальцев) string initMtDNA = "TTCTTTCATGGGGAAGCAGATTTGGGTACCACCCAAGTATTGACTCACCCATCAACAACC"; initMtDNA += "GCTATGTATTTCGTACATTACTGCCAGCCACCATGAATATTGTACGGTACCATAAATACT"; initMtDNA += "TGACCACCTGTAGTACATAAAAACCCAATCCACATCAAAACCCCCCCCTCATGCTTACAA"; initMtDNA += "GCAAGTACAGCAATCAACCTTCAACTATCACACATCAACTGCAACTCCAAAGCCACCCCT"; initMtDNA += "CACCCACTAGGATATCAACAAACCTACCCATCCTTAACAGTACATGGTACATAAAGCCAT"; initMtDNA += "TTACCGTACATAGCACATTACAGTCAAATCCCTTCTCGTCCCCATGGATGACCCCCCTCA"; string subSeqs[] = { "AAAAAAAAAA", "TTTTTTTTTT", "GGGGGGGGGG", "CCCCCCCCCC", "AAAAAAAAAA", "AAAAATTTTT", "TTTTTAAAAA", "GGGGGTTTTT", "CCCCCGGGGG", "GGGGGCCCCC"}; // Митохондриальная ДНК int pos = num*10; num = num > 9 ? 9 : num; initMtDNA.replace(pos, subSeqs[num].length(), subSeqs[num]); return initMtDNA; } BisexualPopulation* DerevyankoReport2015::populationFactory(int type, int size){ BisexualPopulation* ans; std::list males; std::list females; // Фенотип Trait trait1(Trait::Discrete, "Age", 20); // Пул субстратов InnerSubstratesPool* subPool = 0; // (END) Пул субстратов for(int j = 0; j < size/2; j++){ // Подготовка инициаторного генома Chromosome chrom1("Mitochondrial chromosome"); string gName = type == 0 ? "mtDna_KR" : "mtDna_Neand"; int num = 5*type + (j%5); Gene gene(Gene::Sequence, gName, getInitGene(type,num)); chrom1.insertGeneToEnd(gene); Chromosome chrom2("Numt genome"); gName = "Numt"; Gene gene2(Gene::Sequence, gName, ""); chrom2.insertGeneToEnd(gene2); std::vector fGenome; std::vector mGenome; fGenome.push_back(chrom1); fGenome.push_back(chrom2); mGenome.push_back(chrom1); mGenome.push_back(chrom2); 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 DerevyankoReport2015::report2015_01(int argc, char** argv){ std::vector calcScript; if(argc > 1){ calcScript = parseScript(argv[1]); } else{ calcScript = parseScript("script.txt"); } // Генерация стартовых популяций const int krom = 0; const int neand = 1; int popSizeH = Settings::PopulationHomoInitSize; int popSizeA = Settings::PopulationAncientInitSize; float initBirthRate = Settings::BirthRate; float deathRate = Settings::NaturalDeathRate; BisexualPopulation* pop_Sapiens = populationFactory(krom, popSizeH); NeutralEvolutionBreedingStrategy* breedStrat_Sapiens = dynamic_cast(PopulationBreedingStrategy::getInstance("neutral")); breedStrat_Sapiens->setBirthRate(initBirthRate); breedStrat_Sapiens->setDeathRate(deathRate); pop_Sapiens->setBreedingStrategy(breedStrat_Sapiens); BisexualPopulation* pop_Ancient = populationFactory(neand, popSizeA); NeutralEvolutionBreedingStrategy* breedStrat_Ancient = dynamic_cast(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 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"<breedAll(); //pop_Sapiens->mutationAll(); // TODO Сделать мутационные стратегии pop_Sapiens->mutationDerevyanko2015(1e-6, 1e-8); pop_Ancient->breedAll(); //pop_Ancient->mutationAll(); pop_Ancient->mutationDerevyanko2015(1e-6, 1e-8); } } 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"<breedAll(); //pop_Sapiens->mutationAll(); pop_Sapiens->mutationDerevyanko2015(1e-6, 1e-8); pop_Ancient->breedAll(); //pop_Ancient->mutationAll(); pop_Ancient->mutationDerevyanko2015(1e-6, 1e-8); 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"<Ancient\t"<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"<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 DerevyankoReport2015::parseScript(std::string fileName){ std::vector 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 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::const_iterator first = script.begin() + curLine; vector::const_iterator last = script.end(); std::vector 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 a; return a; } //void DerevyankoReport2015::parseCalc(std::vector script){ // //}