output dominance coefficients, solve #60

RangeShifter · Aug 20, 2024 · 02d86a0 · 02d86a0
1 parent e5d112a
commit 02d86a0
Show file tree

Hide file tree

Showing 8 changed files with 38 additions and 12 deletions.
diff --git 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
@@ -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
@@ -42,6 +42,7 @@ class DispersalTrait : public QuantitativeTrait {
 	void inheritGenes(const bool& fromMother, QuantitativeTrait* parent, set<unsigned int> 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
@@ -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
diff --git 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<unsigned int> 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
@@ -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
@@ -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
@@ -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;