diff --git a/doc/Manual/Batchmode/TraitsFile.xlsx b/doc/Manual/Batchmode/TraitsFile.xlsx index 5ebeac2..112181b 100644 Binary files a/doc/Manual/Batchmode/TraitsFile.xlsx and b/doc/Manual/Batchmode/TraitsFile.xlsx differ diff --git a/src/BatchMode.cpp b/src/BatchMode.cpp index 6230642..63412bd 100644 --- a/src/BatchMode.cpp +++ b/src/BatchMode.cpp @@ -2953,7 +2953,7 @@ int CheckTraitsFile(string indir) { string header, colheader; int simNb, nextLineSimNb; - string filename, inTraitType, inSex, inInitDist, inInitParams, + string filename, inTraitType, inSex, inInitDist, inInitParams, inInitDomDist, inInitDomParams, inDominanceDist, inDominanceParams, inIsInherited, inMutationDist, inMutationParams, inPositions, inNbPositions, inExpressionType, inMutationRate, inIsOutput; int nbErrors = 0; @@ -2970,13 +2970,15 @@ int CheckTraitsFile(string indir) bTraitsFile >> header; if (header != "Positions") nbErrors++; bTraitsFile >> header; if (header != "NbrOfPositions") nbErrors++; bTraitsFile >> header; if (header != "ExpressionType") nbErrors++; - bTraitsFile >> header; if (header != "InitialDistribution") nbErrors++; - bTraitsFile >> header; if (header != "InitialParameters") nbErrors++; - bTraitsFile >> header; if (header != "DominanceDistribution") nbErrors++; - bTraitsFile >> header; if (header != "DominanceParameters") nbErrors++; + bTraitsFile >> header; if (header != "InitialAlleleDist") nbErrors++; + bTraitsFile >> header; if (header != "InitialAlleleParams") nbErrors++; + bTraitsFile >> header; if (header != "InitialDomDist") nbErrors++; + bTraitsFile >> header; if (header != "InitialDomParams") nbErrors++; bTraitsFile >> header; if (header != "IsInherited") nbErrors++; bTraitsFile >> header; if (header != "MutationDistribution") nbErrors++; bTraitsFile >> header; if (header != "MutationParameters") nbErrors++; + bTraitsFile >> header; if (header != "DominanceDistribution") nbErrors++; + bTraitsFile >> header; if (header != "DominanceParameters") nbErrors++; bTraitsFile >> header; if (header != "MutationRate") nbErrors++; bTraitsFile >> header; if (header != "OutputValues") nbErrors++; @@ -3005,8 +3007,9 @@ int CheckTraitsFile(string indir) while (!stopReading) { // read and validate columns relating to stage and sex-dependency (NB no IIV here) - bTraitsFile >> inTraitType >> inSex >> inPositions >> inNbPositions >> inExpressionType >> inInitDist >> inInitParams - >> inDominanceDist >> inDominanceParams >> inIsInherited >> inMutationDist >> inMutationParams + bTraitsFile >> inTraitType >> inSex >> inPositions >> inNbPositions + >> inExpressionType >> inInitDist >> inInitParams >> inInitDomDist >> inInitDomParams + >> inIsInherited >> inMutationDist >> inMutationParams >> inDominanceDist >> inDominanceParams >> inMutationRate >> inIsOutput; current = CheckStageSex(whichInputFile, lineNb, simNb, prev, 0, 0, 0, 0, 0, true, false); @@ -3100,24 +3103,19 @@ int CheckTraitsFile(string indir) nbErrors++; } - // Check InitialDistribution + // Check InitialAlleleDist if (tr == NEUTRAL && inInitDist != "uniform") { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "InitialDistribution must be uniform for the neutral trait." << endl; - nbErrors++; - } - if (tr == GENETIC_LOAD && inInitDist != "#") { - BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "InitialDistribution must be blank (#) for genetic load traits." << endl; + batchLog << "InitialAlleleDist must be uniform for the neutral trait." << endl; nbErrors++; } if (isDisp && inInitDist != "normal" && inInitDist != "uniform") { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "InitialDistribution must be either normal or uniform for dispersal traits." << endl; + batchLog << "InitialAlleleDist must be either normal or uniform for dispersal traits." << endl; nbErrors++; } - // Check InitialParameters + // Check InitialAlleleParams const regex patternParamsUnif{ "^\"?min=[-]?([0-9]*[.])?[0-9]+,max=[-]?([0-9]*[.])?[0-9]+\"?$" }; const regex patternParamsNormal{ "^\"?mean=[-]?([0-9]*[.])?[0-9]+,sd=[-]?([0-9]*[.])?[0-9]+\"?$" }; const regex patternParamsGamma{ "^\"?shape=[-]?([0-9]*[.])?[0-9]+,scale=[-]?([0-9]*[.])?[0-9]+\"?$" }; @@ -3129,7 +3127,7 @@ int CheckTraitsFile(string indir) isMatch = regex_search(inInitParams, patternParamsNeutral); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For neutral trait with uniform initialisation, InitialParameters must have form max=int" << endl; + batchLog << "For neutral trait with uniform initialisation, InitialAlleleParams must have form max=int" << endl; nbErrors++; } else { @@ -3144,21 +3142,17 @@ int CheckTraitsFile(string indir) // if not uniform then initDist must be blank, no params else { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For neutral trait with uniform initialisation, InitialParameters must have form max=int" << endl; + batchLog << "For neutral trait with uniform initialisation, InitialAlleleParams must have form max=int" << endl; nbErrors++; } } - if (tr == GENETIC_LOAD && inInitParams != "#") { - BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For genetic load traits, InitialParameters must be blank (#)" << endl; - nbErrors++; - } + if (isDisp) { if (inInitDist == "uniform") { isMatch = regex_search(inInitParams, patternParamsUnif); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For dispersal trait uniform initialisation, InitialParameters must have form min=float,max=float" << endl; + batchLog << "For dispersal trait uniform initialisation, InitialAlleleParams must have form min=float,max=float" << endl; nbErrors++; } } @@ -3166,81 +3160,117 @@ int CheckTraitsFile(string indir) isMatch = regex_search(inInitParams, patternParamsNormal); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For normal initialisation, InitialParameters must have form mean=float,sd=float" << endl; + batchLog << "For normal initialisation, InitialAlleleParams must have form mean=float,sd=float" << endl; nbErrors++; } } } - - // Check DominanceDistribution and DominanceParameters - if (tr == NEUTRAL) { - if (inDominanceDist != "#") { - BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "DominanceDistribution must be left blank (#) for the neutral trait." << endl; - nbErrors++; + if (tr == GENETIC_LOAD) { + if (inInitDist == "uniform") { + isMatch = regex_search(inInitParams, patternParamsUnif); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a uniform distribution, InitialAlleleParams must have form min=float,max=float." << endl; + nbErrors++; + } } - if (inDominanceParams != "#") { - BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "DominanceParameters must be left blank (#) for the neutral trait." << endl; - nbErrors++; + else if (inInitDist == "normal") { + isMatch = regex_search(inInitParams, patternParamsNormal); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a normal distribution, InitialAlleleParams must have form mean=float,sd=float." << endl; + nbErrors++; + } } - } - if (isDisp) { - if (inDominanceDist != "#") { + else if (inInitDist == "gamma") { + isMatch = regex_search(inInitParams, patternParamsGamma); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a Gamma distribution, InitialAlleleParams must have form shape=float,scale=float." << endl; + nbErrors++; + } + } + else if (inInitDist == "negExp") { + isMatch = regex_search(inInitParams, patternParamsMean); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a negative exponential distribution, InitialAlleleParams must have form mean=float." << endl; + nbErrors++; + } + } + else if (inInitDist == "#" && inInitParams != "#") { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "DominanceDistribution must be left blank (#) for dispersal traits." << endl; + batchLog << "If InitialAlleleDist is left blank, InitialAlleleParams must also be blank." << endl; nbErrors++; } - if (inDominanceParams != "#") { + else { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "DominanceParameters must be left blank (#) for dispersal traits." << endl; + batchLog << "For genetic load traits, InitialAlleleDist must be either blank (#), uniform, gamma, negExp or normal" << endl; nbErrors++; } } - if (tr == GENETIC_LOAD) { - if (inDominanceDist == "normal") { - isMatch = regex_search(inDominanceParams, patternParamsNormal); + + // Check InitialDomDist and InitialDomParams + if ((isDisp || tr == NEUTRAL) + && (inInitDomDist != "#" || inInitDomParams != "#")){ + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "InitialDomDist and InitialDomParams must be blank (#) for dispersal and neutral traits." << endl; + nbErrors++; + } + else if (tr == GENETIC_LOAD) { + if (inInitDomDist == "normal") { + isMatch = regex_search(inInitDomParams, patternParamsNormal); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For a normal dominance distribution, DominanceParams must have form mean=float,sd=float" << endl; + batchLog << "For a normal dominance distribution, InitialDomParams must have form mean=float,sd=float" << endl; nbErrors++; } } - else if (inDominanceDist == "gamma") { - isMatch = regex_search(inDominanceParams, patternParamsGamma); + else if (inInitDomDist == "gamma") { + isMatch = regex_search(inInitDomParams, patternParamsGamma); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For a Gamma dominance distribution, DominanceParams must have form shape=float,scale=float" << endl; + batchLog << "For a Gamma dominance distribution, InitialDomParams must have form shape=float,scale=float" << endl; nbErrors++; } } - else if (inDominanceDist == "uniform") { - isMatch = regex_search(inDominanceParams, patternParamsUnif); + else if (inInitDomDist == "uniform") { + isMatch = regex_search(inInitDomParams, patternParamsUnif); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For a uniform dominance distribution, DominanceParams must have form min=float,max=float" << endl; + batchLog << "For a uniform dominance distribution, InitialDomParams must have form min=float,max=float" << endl; nbErrors++; } } - else if (inDominanceDist == "negExp") { - isMatch = regex_search(inDominanceParams, patternParamsMean); + else if (inInitDomDist == "negExp") { + isMatch = regex_search(inInitDomParams, patternParamsMean); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For a negative exponential dominance distribution, DominanceParams must have form mean=float" << endl; + batchLog << "For a negative exponential dominance distribution, InitialDomParams must have form mean=float" << endl; nbErrors++; } } - else if (inDominanceDist == "scaled") { - isMatch = regex_search(inDominanceParams, patternParamsMean); + else if (inInitDomDist == "scaled") { + if (inInitDist == "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "Initial scaled dominance distribution requires InitialAlleleDist to be non-blank." << endl; + nbErrors++; + } + isMatch = regex_search(inInitDomParams, patternParamsMean); if (!isMatch) { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "For a scaled dominance distribution, DominanceParams must have form mean=float" << endl; + batchLog << "For a scaled dominance distribution, InitialDomParams must have form mean=float" << endl; nbErrors++; } } + else if (inInitDomDist == "#" && inInitDomParams != "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "If InitialDomDist is left blank, InitialDomParams must also be blank." << endl; + nbErrors++; + } else { BatchError(whichInputFile, lineNb, 0, " "); - batchLog << "DominanceDistribution must be either normal, gamma, uniform, negExp or scaled for genetic load traits." << endl; + batchLog << "InitialDomDist must be either blank (#), normal, gamma, uniform, negExp or scaled for genetic load traits." << endl; nbErrors++; } } @@ -3365,6 +3395,79 @@ int CheckTraitsFile(string indir) } } + // Check DominanceDistribution and DominanceParameters + if (tr == NEUTRAL) { + if (inDominanceDist != "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "DominanceDistribution must be left blank (#) for the neutral trait." << endl; + nbErrors++; + } + if (inDominanceParams != "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "DominanceParameters must be left blank (#) for the neutral trait." << endl; + nbErrors++; + } + } + if (isDisp) { + if (inDominanceDist != "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "DominanceDistribution must be left blank (#) for dispersal traits." << endl; + nbErrors++; + } + if (inDominanceParams != "#") { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "DominanceParameters must be left blank (#) for dispersal traits." << endl; + nbErrors++; + } + } + if (tr == GENETIC_LOAD) { + if (inDominanceDist == "normal") { + isMatch = regex_search(inDominanceParams, patternParamsNormal); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a normal dominance distribution, DominanceParams must have form mean=float,sd=float" << endl; + nbErrors++; + } + } + else if (inDominanceDist == "gamma") { + isMatch = regex_search(inDominanceParams, patternParamsGamma); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a Gamma dominance distribution, DominanceParams must have form shape=float,scale=float" << endl; + nbErrors++; + } + } + else if (inDominanceDist == "uniform") { + isMatch = regex_search(inDominanceParams, patternParamsUnif); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a uniform dominance distribution, DominanceParams must have form min=float,max=float" << endl; + nbErrors++; + } + } + else if (inDominanceDist == "negExp") { + isMatch = regex_search(inDominanceParams, patternParamsMean); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a negative exponential dominance distribution, DominanceParams must have form mean=float" << endl; + nbErrors++; + } + } + else if (inDominanceDist == "scaled") { + isMatch = regex_search(inDominanceParams, patternParamsMean); + if (!isMatch) { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "For a scaled dominance distribution, DominanceParams must have form mean=float" << endl; + nbErrors++; + } + } + else { + BatchError(whichInputFile, lineNb, 0, " "); + batchLog << "DominanceDistribution must be either normal, gamma, uniform, negExp or scaled for genetic load traits." << endl; + nbErrors++; + } + } + if (inIsOutput != "TRUE" && inIsOutput != "FALSE") { BatchError(whichInputFile, lineNb, 0, " "); batchLog << "OutputValues must be either TRUE or FALSE." << endl; @@ -4918,13 +5021,13 @@ void setUpSpeciesTrait(vector parameters) { const set positions = stringToLoci(parameters[3], parameters[4], genomeSize); const ExpressionType expressionType = stringToExpressionType(parameters[5]); - // Initial distribution parameters + // Initial allele distribution parameters const DistributionType initDist = stringToDistributionType(parameters[6]); const map initParams = stringToParameterMap(parameters[7]); - // Dominance distribution parameters - const DistributionType dominanceDist = stringToDistributionType(parameters[8]); - const map dominanceParams = stringToParameterMap(parameters[9]); + // Initial dominance distribution parameters + const DistributionType initDomDist = stringToDistributionType(parameters[8]); + const map initDomParams = stringToParameterMap(parameters[9]); // Mutation parameters bool isInherited = (parameters[10] == "TRUE"); @@ -4932,27 +5035,33 @@ void setUpSpeciesTrait(vector parameters) { stringToDistributionType(parameters[11]) : DistributionType::NONE; map mutationParameters; - float mutationRate = isInherited ? stof(parameters[13]) : 0.0; if (isInherited) { mutationParameters = stringToParameterMap(parameters[12]); } - int ploidy = gNbSexesDisp; - parameters[14].erase( + // Dominance distribution parameters + const DistributionType dominanceDist = stringToDistributionType(parameters[13]); + const map dominanceParams = stringToParameterMap(parameters[14]); + + float mutationRate = isInherited ? stof(parameters[15]) : 0.0; + + parameters[16].erase( // send windows line endings to hell where they belong - remove(parameters[14].begin(), parameters[14].end(), '\r'), - parameters[14].end() + remove(parameters[16].begin(), parameters[16].end(), '\r'), + parameters[16].end() ); - const bool isOutput = parameters[14] == "TRUE"; + const bool isOutput = parameters[16] == "TRUE"; // Create species trait + int ploidy = gNbSexesDisp; unique_ptr trait(new SpeciesTrait( traitType, sex, positions, expressionType, initDist, initParams, - dominanceDist, dominanceParams, + initDomDist, initDomParams, isInherited, mutationRate, mutationDistribution, mutationParameters, + dominanceDist, dominanceParams, ploidy, isOutput )); diff --git a/src/RScore/GeneticFitnessTrait.cpp b/src/RScore/GeneticFitnessTrait.cpp index 81847de..c20c164 100644 --- a/src/RScore/GeneticFitnessTrait.cpp +++ b/src/RScore/GeneticFitnessTrait.cpp @@ -11,13 +11,106 @@ GeneticFitnessTrait::GeneticFitnessTrait(SpeciesTrait* P) pSpeciesTrait = P; ExpressionType expressionType = pSpeciesTrait->getExpressionType(); - initialise(); - _inherit_func_ptr = (pSpeciesTrait->getPloidy() == 1) ? &GeneticFitnessTrait::inheritHaploid : &GeneticFitnessTrait::inheritDiploid; //this could be changed if we wanted some alternative form of inheritance + // Set initialisation parameters + DistributionType initialDistribution = pSpeciesTrait->getInitialDistribution(); + map initialParameters = pSpeciesTrait->getInitialParameters(); + switch (initialDistribution) { + case UNIFORM: + { + if (initialParameters.count(MAX) != 1) + throw logic_error("Error: initial uniform distribution parameter must contain max value (e.g. max= ) \n"); + if (initialParameters.count(MIN) != 1) + throw logic_error("Error: initial uniform distribution parameter must contain min value (e.g. min= ) \n"); + break; + } + case NORMAL: + { + if (initialParameters.count(MEAN) != 1) + throw logic_error("Error: initial normal distribution parameter must contain mean value (e.g. mean= ) \n"); + if (initialParameters.count(SD) != 1) + throw logic_error("Error: initial normal distribution parameter must contain sdev value (e.g. sdev= ) \n"); + break; + } + case GAMMA: + { + if (initialParameters.count(SHAPE) != 1) + throw logic_error("Error:: genetic load dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n"); + if (initialParameters.count(SCALE) != 1) + throw logic_error("Error:: genetic load dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n"); + break; + } + case NEGEXP: + { + if (initialParameters.count(MEAN) != 1) + throw logic_error("Error:: genetic load dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n"); + break; + } + case NONE: // initialise with default (i.e. zero) values + break; + default: + { + throw logic_error("wrong parameter value for parameter \"initialisation of dispersal traits\", must be uniform/normal \n"); + break; + } + } + + DistributionType initDomDistribution = pSpeciesTrait->getInitDomDistribution(); + map initDomParameters = pSpeciesTrait->getInitDomParameters(); + switch (initDomDistribution) { + case UNIFORM: + { + if (initDomParameters.count(MAX) != 1) + throw logic_error("Error:: genetic load dominance uniform distribution parameter must contain one max value (e.g. max= ) \n"); + if (initDomParameters.count(MIN) != 1) + throw logic_error("Error:: genetic load dominance uniform distribution parameter must contain one min value (e.g. min= ) \n"); + break; + } + case NORMAL: + { + if (initDomParameters.count(MEAN) != 1) + throw logic_error("Error:: genetic load dominance distribution set to normal so parameters must contain one mean value (e.g. mean= ) \n"); + if (initDomParameters.count(SD) != 1) + throw logic_error("Error:: genetic load dominance distribution set to normal so parameters must contain one sdev value (e.g. sdev= ) \n"); + break; + } + case GAMMA: + { + if (initDomParameters.count(SHAPE) != 1) + throw logic_error("Error:: genetic load dominance distribution set to gamma so parameters must contain one shape value (e.g. shape= ) \n"); + if (initDomParameters.count(SCALE) != 1) + throw logic_error("Error:: genetic load dominance distribution set to gamma so parameters must contain one scale value (e.g. scale= ) \n"); + break; + } + case NEGEXP: + { + if (initDomParameters.count(MEAN) != 1) + throw logic_error("Error:: genetic load dominance distribution set to negative exponential (negative decay) so parameters must contain mean value (e.g. mean= ) \n"); + break; + } + case SCALED: + { + if (initDomParameters.count(MEAN) != 1) + throw logic_error("Error:: genetic load dominance distribution set to scaled, so parameters must contain mean dominance value (e.g. mean= ) \n"); + // Set for drawing initial values + setScaledCoeff(initialDistribution, initialParameters); + break; + } + case NONE: // default values, zero-dominance coefficients + break; + default: + { + throw logic_error("Error:: wrong parameter value for genetic load dominance model, must be uniform/normal/gamma/negExp/scaled \n"); + break; + } + } + + // Draw initial values + initialise(); + DistributionType mutationDistribution = pSpeciesTrait->getMutationDistribution(); map mutationParameters = pSpeciesTrait->getMutationParameters(); - switch (mutationDistribution) { case UNIFORM: { @@ -92,24 +185,8 @@ GeneticFitnessTrait::GeneticFitnessTrait(SpeciesTrait* P) if (dominanceParameters.count(MEAN) != 1) throw logic_error("Error:: genetic load dominance distribution set to scaled, so parameters must contain mean dominance value (e.g. mean= ) \n"); - // Calculate mean selection coeff s_d for calculation of k - switch (mutationDistribution) - { - case UNIFORM: - scaledDomMeanSelCoeff = (mutationParameters.find(MIN)->second + mutationParameters.find(MAX)->second) / 2; - break; - case NORMAL: - scaledDomMeanSelCoeff = mutationParameters.find(MEAN)->second; - break; - case GAMMA: - scaledDomMeanSelCoeff = mutationParameters.find(SHAPE)->second * mutationParameters.find(SCALE)->second; - break; - case NEGEXP: - scaledDomMeanSelCoeff = 1 / mutationParameters.find(MEAN)->second; - break; - default: - break; - } + // Set for drawing mutations (overwrite initial value) + setScaledCoeff(mutationDistribution, mutationParameters); break; } default: @@ -120,6 +197,29 @@ GeneticFitnessTrait::GeneticFitnessTrait(SpeciesTrait* P) } } +// Calculate mean selection coeff s_d for calculation of k +void GeneticFitnessTrait::setScaledCoeff(const DistributionType& selCoeffDist, const map& selCoeffParams) +{ + switch (selCoeffDist) + { + case UNIFORM: + scaledDomMeanSelCoeff = (selCoeffParams.find(MIN)->second + selCoeffParams.find(MAX)->second) / 2; + break; + case NORMAL: + scaledDomMeanSelCoeff = selCoeffParams.find(MEAN)->second; + break; + case GAMMA: + scaledDomMeanSelCoeff = selCoeffParams.find(SHAPE)->second * selCoeffParams.find(SCALE)->second; + break; + case NEGEXP: + scaledDomMeanSelCoeff = 1 / selCoeffParams.find(MEAN)->second; + break; + case NONE: + throw logic_error("Scaled dominance distribution cannot be used with default allele distribution."); + default: break; + } +} + // ---------------------------------------------------------------------------------------- // Inheritance constructor // Copies immutable features from a parent trait @@ -134,12 +234,25 @@ GeneticFitnessTrait::GeneticFitnessTrait(const GeneticFitnessTrait& T) : } void GeneticFitnessTrait::initialise() { - // All positions start at wild type, mutations accumulate through simulation - const set genePositions = pSpeciesTrait->getGenePositions(); + float initSelCoeff; + float initDomCoeff; short ploidy = pSpeciesTrait->getPloidy(); - const vector> wildTypeGene(ploidy, wildType); + auto initDist = pSpeciesTrait->getInitialDistribution(); + auto initParams = pSpeciesTrait->getInitialParameters(); + auto initDomDist = pSpeciesTrait->getInitDomDistribution(); + auto initDomParams = pSpeciesTrait->getInitDomParameters(); + + const set genePositions = pSpeciesTrait->getGenePositions(); for (auto position : genePositions) { - genes.insert(make_pair(position, wildTypeGene)); + vector> initialGene(ploidy); + for (int p = 0; p < ploidy; p++) { + initSelCoeff = initDist == NONE ? 0.0 + : drawSelectionCoef(initDist, initParams); + initDomCoeff = initDomDist == NONE ? 0.0 + : drawDominance(initSelCoeff, initDomDist, initDomParams); + initialGene[p] = make_shared(initSelCoeff, initDomCoeff); + } + genes.insert(make_pair(position, initialGene)); } } @@ -171,8 +284,15 @@ void GeneticFitnessTrait::mutate() auto it = genes.find(m); if (it == genes.end()) throw runtime_error("Locus sampled for mutation doesn't exist."); - newSelectionCoef = drawSelectionCoef(); - newDominanceCoef = drawDominance(newSelectionCoef); + newSelectionCoef = drawSelectionCoef( + pSpeciesTrait->getMutationDistribution(), + pSpeciesTrait->getMutationParameters() + ); + newDominanceCoef = drawDominance( + newSelectionCoef, + pSpeciesTrait->getDominanceDistribution(), + pSpeciesTrait->getDominanceParameters() + ); it->second[p] = make_shared(newSelectionCoef, newDominanceCoef); } } @@ -182,24 +302,21 @@ void GeneticFitnessTrait::mutate() // ---------------------------------------------------------------------------------------- // get dominance value for new mutation // ---------------------------------------------------------------------------------------- -float GeneticFitnessTrait::drawDominance(float selCoef) { - - DistributionType dominanceDistribution = pSpeciesTrait->getDominanceDistribution(); - map dominanceParameters = pSpeciesTrait->getDominanceParameters(); +float GeneticFitnessTrait::drawDominance(float selCoef, const DistributionType& domDist, const map& domParams) { float h; - switch (dominanceDistribution) { + switch (domDist) { case UNIFORM: { - float maxD = dominanceParameters.find(MAX)->second; - float minD = dominanceParameters.find(MIN)->second; + float maxD = domParams.find(MAX)->second; + float minD = domParams.find(MIN)->second; h = pRandom->FRandom(minD, maxD); break; } case NORMAL: { - const float mean = dominanceParameters.find(MEAN)->second; - const float sd = dominanceParameters.find(SD)->second; + const float mean = domParams.find(MEAN)->second; + const float sd = domParams.find(SD)->second; do { h = static_cast(pRandom->Normal(mean, sd)); } while (h <= 0.0); @@ -207,20 +324,20 @@ float GeneticFitnessTrait::drawDominance(float selCoef) { } case GAMMA: { - const float shape = dominanceParameters.find(SHAPE)->second; - const float scale = dominanceParameters.find(SCALE)->second; + const float shape = domParams.find(SHAPE)->second; + const float scale = domParams.find(SCALE)->second; h = static_cast(pRandom->Gamma(shape, scale)); break; } case NEGEXP: { - const float mean = dominanceParameters.find(MEAN)->second; + const float mean = domParams.find(MEAN)->second; h = static_cast(pRandom->NegExp(mean)); break; } case SCALED: { - const float h_d = dominanceParameters.find(MEAN)->second; + const float h_d = domParams.find(MEAN)->second; const float k = -log(2 * h_d) / scaledDomMeanSelCoeff; const float max = exp(-k * selCoef); h = pRandom->FRandom(0, max); @@ -243,10 +360,7 @@ float GeneticFitnessTrait::drawDominance(float selCoef) { // if the mutation distribution enable it, take a negative value // down to -1 representing the effect of beneficial mutations // ---------------------------------------------------------------------------------------- -float GeneticFitnessTrait::drawSelectionCoef() { - - DistributionType mutationDistribution = pSpeciesTrait->getMutationDistribution(); - map mutationParameters = pSpeciesTrait->getMutationParameters(); +float GeneticFitnessTrait::drawSelectionCoef(const DistributionType& mutationDistribution, const map& mutationParameters) { float s = 0.0; // default selection coefficient is 0 diff --git a/src/RScore/GeneticFitnessTrait.h b/src/RScore/GeneticFitnessTrait.h index 001e1c2..db9e7ab 100644 --- a/src/RScore/GeneticFitnessTrait.h +++ b/src/RScore/GeneticFitnessTrait.h @@ -59,10 +59,19 @@ class GeneticFitnessTrait : public QuantitativeTrait { // Initialisation void initialise(); + void setScaledCoeff(const DistributionType& selCoeffDist, const map& selCoeffParams); + // Mutation float scaledDomMeanSelCoeff = 0; // s_d, only for scaled dominance distribution - float drawDominance(float); - float drawSelectionCoef(); + float drawDominance( + float selCoef, + const DistributionType& domDist, + const map& domParams + ); + float drawSelectionCoef( + const DistributionType& mutationDistribution, + const map& mutationParameters + ); // Immutable features, set at initialisation // and passed down to every subsequent trait copy diff --git a/src/RScore/SpeciesTrait.cpp b/src/RScore/SpeciesTrait.cpp index 61db3a5..09afa6f 100644 --- a/src/RScore/SpeciesTrait.cpp +++ b/src/RScore/SpeciesTrait.cpp @@ -3,16 +3,20 @@ // Species trait constructor SpeciesTrait::SpeciesTrait( - const TraitType& trType, const sex_t& sx, - const set& pos, + const TraitType& trType, + const sex_t& sx, + const set& pos, const ExpressionType& expr, - const DistributionType& initDist, + const DistributionType& initDist, const map initParams, - const DistributionType& dominanceDist, - const map dominanceParams, - bool isInherited, const float& mutRate, - const DistributionType& mutationDist, + const DistributionType& initDomDist, + const map initDomParams, + bool isInherited, + const float& mutRate, + const DistributionType& mutationDist, const map mutationParams, + const DistributionType& dominanceDist, + const map dominanceParams, const int nPloidy, const bool isOutput) : traitType{ trType }, @@ -21,6 +25,8 @@ SpeciesTrait::SpeciesTrait( expressionType{ expr }, initialDistribution{ initDist }, initialParameters{ initParams }, + initialDomDistribution{ initDomDist }, + initialDomParameters{ initDomParams }, dominanceDistribution{ dominanceDist }, dominanceParameters{ dominanceParams }, inherited{ isInherited }, @@ -31,7 +37,7 @@ SpeciesTrait::SpeciesTrait( traitIsOutput{ isOutput } { // Check distribution parameters - // Initial distribution + // Initial allele distribution for (auto [paramType, paramVal] : initParams) { switch (paramType) { @@ -48,6 +54,23 @@ SpeciesTrait::SpeciesTrait( } } + // Initial dominance distribution + for (auto [paramType, paramVal] : initDomParams) { + switch (paramType) + { + case MIN: case MAX: case MEAN: + if (paramVal < 0.0) + throw logic_error("Invalid parameter value: initial dominance parameter " + to_string(paramType) + " must not be negative."); + break; + case SD: case SHAPE: case SCALE: + if (paramVal <= 0.0) + throw logic_error("Invalid parameter value: initial dominance parameter " + to_string(paramType) + " must be strictly positive"); + break; + default: + break; + } + } + // Mutation distribution for (auto [paramType, paramVal] : mutationParams) { switch (paramType) @@ -205,6 +228,8 @@ SpeciesTrait* createTestEmigSpTrait(const set& genePositions, const bool& i 0.0, // mutation rate DistributionType::UNIFORM, distParams, + DistributionType::NONE, // no dominance + distParams, isDiploid ? 2 : 1, false ); @@ -222,14 +247,16 @@ SpeciesTrait* createTestGenLoadTrait(const set& genePositions, const bool& sex_t::NA, genePositions, ExpressionType::MULTIPLICATIVE, - DistributionType::UNIFORM, + DistributionType::NONE, distParams, - DistributionType::UNIFORM, + DistributionType::NONE, // initialise dominance to zero distParams, true, // isInherited 0.0, // mutation rate DistributionType::UNIFORM, distParams, + DistributionType::UNIFORM, + distParams, isDiploid ? 2 : 1, false ); @@ -249,13 +276,15 @@ SpeciesTrait* createTestNeutralSpTrait(const float& maxAlleleVal, const set ExpressionType::NOTEXPR, // Sample initial values from uniform(0, max) DistributionType::UNIFORM, distParams, - // No dominance - DistributionType::NONE, map{}, + DistributionType::NONE, // No dominance + map{}, true, // isInherited 0.0, // mutation rate // Mutation sampled from a uniform(0, max) DistributionType::KAM, distParams, + DistributionType::NONE, // No dominance + map{}, isDiploid ? 2 : 1, false ); diff --git a/src/RScore/SpeciesTrait.h b/src/RScore/SpeciesTrait.h index 5a323e9..c866dd4 100644 --- a/src/RScore/SpeciesTrait.h +++ b/src/RScore/SpeciesTrait.h @@ -22,12 +22,14 @@ class SpeciesTrait { const ExpressionType& expr, const DistributionType& initDist, const map initParams, - const DistributionType& dominanceDist, - const map dominanceParams, + const DistributionType& initDomDist, + const map initDomParams, bool isInherited, const float& mutationRate, const DistributionType& mutationDist, const map mutationParams, + const DistributionType& dominanceDist, + const map dominanceParams, const int ploidy, const bool isOutput ); @@ -44,12 +46,15 @@ class SpeciesTrait { int getPositionsSize() const { return static_cast(genePositions.size()); } bool isInherited() const { return inherited; } + DistributionType getInitialDistribution() const { return initialDistribution; }; + map getInitialParameters() const { return initialParameters; }; + DistributionType getInitDomDistribution() const { return initialDomDistribution; }; + map getInitDomParameters() const { return initialDomParameters; }; DistributionType getMutationDistribution() const { return mutationDistribution; }; map getMutationParameters() const { return mutationParameters; }; DistributionType getDominanceDistribution() const { return dominanceDistribution; }; map getDominanceParameters() const { return dominanceParameters; }; - DistributionType getInitialDistribution() const { return initialDistribution; }; - map getInitialParameters() const { return initialParameters; }; + ExpressionType getExpressionType() const { return expressionType; }; int getNbNeutralAlleles() const { @@ -78,6 +83,8 @@ class SpeciesTrait { ExpressionType expressionType; DistributionType initialDistribution; map initialParameters; + DistributionType initialDomDistribution; + map initialDomParameters; DistributionType dominanceDistribution; map dominanceParameters; bool inherited;