From 02d86a01ef0fd14c58c0e8ef2d324be3cd108d69 Mon Sep 17 00:00:00 2001 From: Theo Pannetier Date: Tue, 20 Aug 2024 15:04:42 +0100 Subject: [PATCH] output dominance coefficients, solve #60 --- src/RScore/Community.cpp | 4 ++-- src/RScore/DispersalTrait.cpp | 7 +++++++ src/RScore/DispersalTrait.h | 1 + src/RScore/GeneticFitnessTrait.cpp | 14 +++++++------- src/RScore/GeneticFitnessTrait.h | 2 +- src/RScore/NeutralTrait.h | 4 ++++ src/RScore/Population.cpp | 17 +++++++++++++++-- src/RScore/QuantitativeTrait.h | 1 + 8 files changed, 38 insertions(+), 12 deletions(-) diff --git a/src/RScore/Community.cpp b/src/RScore/Community.cpp index c940316..ccf3662 100644 --- a/src/RScore/Community.cpp +++ b/src/RScore/Community.cpp @@ -1584,10 +1584,10 @@ bool Community::openOutGenesFile(const bool& isDiploid, const int landNr, const ofsGenes.open(name.c_str()); ofsGenes << "Year\tGeneration\tIndID\ttraitType\tlocusPosition"; if (isDiploid) { - ofsGenes << "\talleleValueA\talleleValueB"; + ofsGenes << "\talleleValueA\tdomCoefA\talleleValueBA\tdomCoefB"; } else { - ofsGenes << "\talleleValue"; + ofsGenes << "\talleleValueA\tdomCoefA"; } ofsGenes << endl; return ofsGenes.is_open(); diff --git a/src/RScore/DispersalTrait.cpp b/src/RScore/DispersalTrait.cpp index 1cf59df..8413775 100644 --- a/src/RScore/DispersalTrait.cpp +++ b/src/RScore/DispersalTrait.cpp @@ -455,6 +455,13 @@ float DispersalTrait::getAlleleValueAtLocus(short whichChromosome, int position) return it->second[whichChromosome].get()->getAlleleValue(); } +float DispersalTrait::getDomCoefAtLocus(short whichChromosome, int position) const { + auto it = genes.find(position); + if (it == genes.end()) + throw runtime_error("The genetic load locus queried for its dominance coefficient does not exist."); + return it->second[whichChromosome]->getDominanceCoef(); +} + #if RSDEBUG // Testing only // Get allele ID at locus diff --git a/src/RScore/DispersalTrait.h b/src/RScore/DispersalTrait.h index 3f9815c..76a76ab 100644 --- a/src/RScore/DispersalTrait.h +++ b/src/RScore/DispersalTrait.h @@ -42,6 +42,7 @@ class DispersalTrait : public QuantitativeTrait { void inheritGenes(const bool& fromMother, QuantitativeTrait* parent, set const& recomPositions, int startingChromosome) override; float getAlleleValueAtLocus(short chromosome, int i) const override; + float getDomCoefAtLocus(short chromosome, int position) const override; int countHeterozygoteLoci() const; bool isHeterozygoteAtLocus(int locus) const override; diff --git a/src/RScore/GeneticFitnessTrait.cpp b/src/RScore/GeneticFitnessTrait.cpp index 771a676..cd266c5 100644 --- a/src/RScore/GeneticFitnessTrait.cpp +++ b/src/RScore/GeneticFitnessTrait.cpp @@ -426,6 +426,13 @@ float GeneticFitnessTrait::getAlleleValueAtLocus(short whichChromosome, int posi return it->second[whichChromosome] == 0 ? wildType->getAlleleValue() : it->second[whichChromosome]->getAlleleValue(); } +float GeneticFitnessTrait::getDomCoefAtLocus(short whichChromosome, int position) const { + auto it = genes.find(position); + if (it == genes.end()) + throw runtime_error("The genetic load locus queried for its dominance coefficient does not exist."); + return it->second[whichChromosome] == 0 ? wildType->getDominanceCoef() : it->second[whichChromosome]->getDominanceCoef(); +} + #if RSDEBUG // Testing only // Get allele ID at locus int GeneticFitnessTrait::getAlleleIDAtLocus(short whichChromosome, int position) const { @@ -435,11 +442,4 @@ int GeneticFitnessTrait::getAlleleIDAtLocus(short whichChromosome, int position) return it->second[whichChromosome].get()->getId(); } -float GeneticFitnessTrait::getDomCoefAtLocus(short whichChromosome, int position) const { - auto it = genes.find(position); - if (it == genes.end()) - throw runtime_error("The genetic load locus queried for its dominance coefficient does not exist."); - return it->second[whichChromosome] == 0 ? wildType->getDominanceCoef() : it->second[whichChromosome]->getDominanceCoef(); -} - #endif // RSDEBUG \ No newline at end of file diff --git a/src/RScore/GeneticFitnessTrait.h b/src/RScore/GeneticFitnessTrait.h index 8bd60ac..ebb6e76 100644 --- a/src/RScore/GeneticFitnessTrait.h +++ b/src/RScore/GeneticFitnessTrait.h @@ -46,11 +46,11 @@ class GeneticFitnessTrait : public QuantitativeTrait { virtual void inheritGenes(const bool& fromMother, QuantitativeTrait* parent, set const& recomPositions, int startingChromosome) override; virtual float getAlleleValueAtLocus(short chromosome, int position) const override; + virtual float getDomCoefAtLocus(short chromosome, int position) const override; virtual int countHeterozygoteLoci() const; virtual bool isHeterozygoteAtLocus(int locus) const override; #if RSDEBUG // for testing only int getAlleleIDAtLocus(short whichChromosome, int position) const; - float getDomCoefAtLocus(short chromosome, int position) const; #endif private: diff --git a/src/RScore/NeutralTrait.h b/src/RScore/NeutralTrait.h index e2a0e6d..ed12339 100644 --- a/src/RScore/NeutralTrait.h +++ b/src/RScore/NeutralTrait.h @@ -47,6 +47,10 @@ class NeutralTrait : public QuantitativeTrait { } virtual float getAlleleValueAtLocus(short chromosome, int position) const override; + virtual float getDomCoefAtLocus(short chromosome, int position) const override { + return 0.0; + } + virtual int countHeterozygoteLoci() const; virtual bool isHeterozygoteAtLocus(int locus) const override; #if RSDEBUG // for testing only diff --git a/src/RScore/Population.cpp b/src/RScore/Population.cpp index 97ee943..59b302f 100644 --- a/src/RScore/Population.cpp +++ b/src/RScore/Population.cpp @@ -1938,6 +1938,7 @@ void Population::outputGeneValues(ofstream& ofsGenes, const int& yr, const int& auto traitTypes = pSpecies->getTraitTypes(); int indID; float alleleOnChromA, alleleOnChromB; + float domCoefA, domCoefB; // Fetch map to positions for each trait // Presumably faster than fetching for every individual @@ -1955,10 +1956,22 @@ void Population::outputGeneValues(ofstream& ofsGenes, const int& yr, const int& auto indTrait = ind->getTrait(trType); for (auto pos : positions) { alleleOnChromA = indTrait->getAlleleValueAtLocus(0, pos); - ofsGenes << yr << '\t' << gen << '\t' << indID << '\t' << trType << '\t' << pos << '\t' << alleleOnChromA; + if (trType == GENETIC_LOAD1 || trType == GENETIC_LOAD2 || trType == GENETIC_LOAD3 || trType == GENETIC_LOAD4 || trType == GENETIC_LOAD5) { + domCoefA = indTrait->getDomCoefAtLocus(0, pos); + } + else { + domCoefA = 0.0; + } + ofsGenes << yr << '\t' << gen << '\t' << indID << '\t' << to_string(trType) << '\t' << pos << '\t' << alleleOnChromA << '\t' << domCoefA; if (isDiploid) { alleleOnChromB = indTrait->getAlleleValueAtLocus(1, pos); - ofsGenes << '\t' << alleleOnChromB; + if (trType == GENETIC_LOAD1 || trType == GENETIC_LOAD2 || trType == GENETIC_LOAD3 || trType == GENETIC_LOAD4 || trType == GENETIC_LOAD5) { + domCoefB = indTrait->getDomCoefAtLocus(1, pos); + } + else { + domCoefB = 0.0; + } + ofsGenes << '\t' << alleleOnChromB << '\t' << domCoefB; } ofsGenes << endl; } diff --git a/src/RScore/QuantitativeTrait.h b/src/RScore/QuantitativeTrait.h index e19a9bf..7d0a0e4 100644 --- a/src/RScore/QuantitativeTrait.h +++ b/src/RScore/QuantitativeTrait.h @@ -21,6 +21,7 @@ class QuantitativeTrait { virtual float getMutationRate() const = 0; virtual bool isInherited() const = 0; virtual float getAlleleValueAtLocus(short chromosome, int i) const = 0; + virtual float getDomCoefAtLocus(short whichChromosome, int position) const = 0; virtual int countHeterozygoteLoci() const = 0; virtual bool isHeterozygoteAtLocus(int loci) const = 0; virtual float express() = 0;