package presentation.epiphyte;
import java.util.ArrayList;
import java.util.Calendar;
import java.util.Iterator;
import java.util.TreeSet;
import presentation.dataOutput.C_OutputData;
import repast.simphony.essentials.RepastEssentials;
import thing.C_Rodent;
import thing.dna.C_ChromosomePair;
import thing.dna.C_GenePair;
import thing.dna.C_GenomeEucaryote;
import thing.dna.C_XsomePairMicrosat;
import com.vividsolutions.jts.geom.Coordinate;
import data.I_sim_constants;
/** Retrieves information on genes and alleles of the whole population.
* The rodents' list (listRodents) is obtained from the population inspector and updated when needed. This
* implementation to avoid permanent calls to C_InspectorPopulation.listRodents static field
* @see C_InspectorPopulation#listRodents
* @author A Realini 05.2011 rev. J. Le Fur 03.2013 - SimMasto project (http://simmasto.org)*/
public class C_InspectorGenetic extends A_Inspector {
protected TreeSet listRodents = C_InspectorPopulation.listRodents;
protected int genePopSaveStep = 0;// steps interval for writing genePop files
protected C_OutputData genesFile;
protected C_OutputData genePopFile;
protected double observedHeterozygosityHO = 0, expectedHeterozygosityHE = 0, meanAllelicRichness = 0;
/** List of list of allelic frequencies for each locus */
protected ArrayList[] allelicFrequencies = new ArrayList[I_sim_constants.NUMBER_GENES];
protected double[] expectedHeterozygosityByLocus = new double[I_sim_constants.NUMBER_GENES];
protected int[] richnessByLocus = new int[I_sim_constants.NUMBER_GENES];
public C_InspectorGenetic() {
super();
indicatorsHeader = "Tick;ObservedHetero;ExpectedHetero;FIS;MeanAllelicRichness";
genePopFile = new C_OutputData("GenePop");
initTabGenePopFile();
genesFile = new C_OutputData("Genes.csv");
initTabGenesFile();
}
@Override
public void step() {
setListRodents(); // update the list before computing
super.step(); // // compute and store indicators values
// save private files
recordGenesInFile();
// write genePop file only at an interval defined in I_sim_constants
if (RepastEssentials.GetTickCount() == I_sim_constants.INTERVAL_ECRITURE_GENE_POP * genePopSaveStep) {
recordGenePopInFile();
genePopSaveStep++;
}
}
/** Update the private rodent list */
private void setListRodents() {
listRodents = C_InspectorPopulation.listRodents;
}
/** Initialize the list of different alleles, for the locus passed in parameter, with the firs rodent of the
* listRodents list (list in C_InspectorPopulation).
* @param locus : locus number for which allelic richness is requested */
private ArrayList initDifferingAllelesList(int locus) {
ArrayList listDifferingAllelesByLocus = new ArrayList();
C_XsomePairMicrosat microsat = ((C_GenomeEucaryote) listRodents.first().getGenome())
.getMicrosatXsome();
C_GenePair genesPair = (C_GenePair) microsat.getLocusAllele(locus);
listDifferingAllelesByLocus.add(genesPair.getGene(C_ChromosomePair.PARENT_1).getAllele());
if (isHeterozygous(genesPair)) listDifferingAllelesByLocus.add(genesPair.getGene(
C_ChromosomePair.PARENT_2).getAllele());
return listDifferingAllelesByLocus;
}
/** Compute allelic richness for each locus (nb differing alleles for this locus).
* Results are put in the table RichnessByLocus */
private void computeAllelicRichnessEveryLocus() {
for (int locus = 0; locus < I_sim_constants.NUMBER_GENES; locus++) {
ArrayList listDifferingAllelesByLocus = initDifferingAllelesList(locus);
Iterator rodents = listRodents.iterator();
while (rodents.hasNext()) {
C_XsomePairMicrosat microsat = ((C_GenomeEucaryote) rodents.next().getGenome())
.getMicrosatXsome();
C_GenePair genePair = (C_GenePair) microsat.getLocusAllele(locus);
Object allele = genePair.getGene(C_ChromosomePair.PARENT_1).getAllele();
if (!listDifferingAllelesByLocus.contains(allele)) listDifferingAllelesByLocus.add(allele);
if (isHeterozygous(genePair)) {
allele = genePair.getGene(C_ChromosomePair.PARENT_2).getAllele();
if (!listDifferingAllelesByLocus.contains(allele)) listDifferingAllelesByLocus
.add(allele);
}
}
richnessByLocus[locus] = listDifferingAllelesByLocus.size();
}
}
/** Compute the mean allelic richness for the whole population
* Formula: sum(allelic richness by locus)/nb genes */
private void computeAllelicRichnessPopMean() {
double somme = 0;
for (int i = 0; i < richnessByLocus.length; i++)
somme += richnessByLocus[i];
meanAllelicRichness = convertNumber((somme / I_sim_constants.NUMBER_GENES), 3);
}
/** Compute allelic frequencies for the locus passed in parameter and writes in the allelicFrequencies list
* Formula: nb occurences of one allele at locus i / (2 * nb indiv.) */
private void computeAllelicRichnessOneLocus(int locus) {
ArrayList freq = new ArrayList();
if (!listRodents.isEmpty()) {
ArrayList listeAlleles = getListeAllAllelesOneLocus(locus);
while (!listeAlleles.isEmpty()) {
Object allele = listeAlleles.get(0);// retrieve the first allele of the list to compare with the others
listeAlleles.remove(0); // remove the reference allele
// Retrieves the index where the first occurence of the allele is found. If not found, index=-1
int cpt = 1;
int index = listeAlleles.indexOf(allele);
while (index != -1) {
listeAlleles.remove(index);
cpt++;
index = listeAlleles.indexOf(allele);
}
freq.add(convertNumber((double) cpt / (2 * listRodents.size()), 5));
}
}
allelicFrequencies[locus] = freq;
}
/** Return the list of every alleles existing for the locus passed as a parameter (used to compute the allelic
* frequency */
private ArrayList getListeAllAllelesOneLocus(int locus) {
ArrayList liste = new ArrayList();
// add the chromosomes left and right to the list and not the pair per se (for easier computation)
Iterator rodents = listRodents.iterator();
while (rodents.hasNext()) {
C_Rodent rodent = rodents.next();
C_XsomePairMicrosat microsat = ((C_GenomeEucaryote) rodent.getGenome()).getMicrosatXsome();
liste.add(((C_GenePair) microsat.getLocusAllele(locus)).getGene(C_ChromosomePair.PARENT_1)
.getAllele());
liste.add(((C_GenePair) microsat.getLocusAllele(locus)).getGene(C_ChromosomePair.PARENT_2)
.getAllele());
}
return liste;
}
/** Compute global observed heterozygosity
* Formula : sum(nb heterozygotes at locus i/nb indiv) / nb locus */
private void computeHeterozygosityObservedHO() {
// for each locus, count the number of heterozygotes then divide by the number of individuals
double sum = 0;
for (int locus = 0; locus < I_sim_constants.NUMBER_GENES; locus++) {
int cpt = 0;
Iterator rodents = listRodents.iterator();
while (rodents.hasNext()) {
C_Rodent rodent = rodents.next();
C_XsomePairMicrosat microsat = ((C_GenomeEucaryote) rodent.getGenome()).getMicrosatXsome();
if (isHeterozygous((C_GenePair) microsat.getLocusAllele(locus))) cpt++;
}
sum += (double) cpt / (double) listRodents.size();
}
// divide the sum by the number of loci //
observedHeterozygosityHO = convertNumber((sum / I_sim_constants.NUMBER_GENES), 5);
}
/** Compute expected heterozygosity for the locus passed as parameter
* Formula : 1 - sum((allelic frequency at locus i)^2) */
private void computeHeterozygosityExpectedHE(int locus) {
double sum = 0.0;
computeAllelicRichnessOneLocus(locus);
ArrayList freq = allelicFrequencies[locus];
for (int i = 0; i < freq.size(); i++)
sum += Math.pow(freq.get(i), 2);
expectedHeterozygosityByLocus[locus] = convertNumber(1 - sum, 3);
}
/** Compute expected heterozygosity over the whole population
* Formula: mean of expected heterozygosities */
private void computeHeterozygosityExpectedGlobal() {
double sum = 0;
for (int locus = 0; locus < I_sim_constants.NUMBER_GENES; locus++) {
computeHeterozygosityExpectedHE(locus);
sum += expectedHeterozygosityByLocus[locus];
}
expectedHeterozygosityHE = convertNumber((sum / I_sim_constants.NUMBER_GENES), 5);
}
/** Compute and return the fixation index */
public double getFixationIndex() {
setListRodents(); // update the list before computing
double fisValue = 1 - (observedHeterozygosityHO / expectedHeterozygosityHE);
// patch: at the end of simulation, FIS drops to -1 which ruins the graph.
if ((convertNumber(fisValue, 5) == -1.) || (listRodents.size() == 1)) return 0.;
else return convertNumber(fisValue, 5);
}
/** Return true is the gene pair passed as parameter is heterozygous, else return false */
private boolean isHeterozygous(C_GenePair paire) {
Object left = paire.getGene(C_ChromosomePair.PARENT_1).getAllele();
Object right = paire.getGene(C_ChromosomePair.PARENT_2).getAllele();
if (left.equals(right)) return false;
else return true;
}
@Override
public void computeIndicators() {
setListRodents(); // update the list before computing
super.computeIndicators();
computeAllelicRichnessEveryLocus();
computeHeterozygosityExpectedGlobal();
computeHeterozygosityObservedHO();
computeAllelicRichnessPopMean();
}
@Override
/** Store the current state of indicators in the field including the super ones / JLF 01.2013*/
public String storeIndicatorsValues() {
indicatorsValues = RepastEssentials.GetTickCount() + ";" + observedHeterozygosityHO + ";"
+ expectedHeterozygosityHE + ";" + getFixationIndex() + ";" + meanAllelicRichness;
return indicatorsValues;
}
// OUTPUT FILES MANAGEMENT //
/** Initialize the header of the genetic indicators file */
private void initTabGenesFile() {
genesFile.write("Tick;");
loopOnGenesForGenesFileHeader("locus", ';');
loopOnGenesForGenesFileHeader("LocusFrequency", ';');
loopOnGenesForGenesFileHeader("expectedHetero", ';');
genesFile.writeln("ObservedHetero;ExpectedHetero;FIS;MeanAllelicRichness");
}
/** Write the title with an index at the end as many times as there are requested genes to output
* @param title : column title to rewrite NUMBER_GENES times
* @param separator : separator between two strings of the headers */
private void loopOnGenesForGenesFileHeader(String title, char separator) {
for (int i = 0; i < I_sim_constants.NUMBER_GENES; i++)
genesFile.write(title + i + separator);
}
/** Initialize the header for genePop output file */
private void initTabGenePopFile() {
genePopFile.writeln(Calendar.getInstance().getTime() + " Run n°" + genePopFile.numRun
+ " | Species: ");
for (int i = 0; i < I_sim_constants.NUMBER_GENES; i++)
genePopFile.writeln("locus" + i);
}
/** Output data in genePop file */
private void recordGenePopInFile() {
C_GenePair[] paire = new C_GenePair[I_sim_constants.NUMBER_GENES];
genePopFile.writeln("Pop " + RepastEssentials.GetTickCount());
Iterator rodents = listRodents.iterator();
while (rodents.hasNext()) {
C_Rodent microtus = rodents.next();
Coordinate point = microtus.getCoord_Umeters();
genePopFile.write(point.x + " " + point.y + " , ");
for (int locus = 0; locus < I_sim_constants.NUMBER_GENES; locus++) {
C_XsomePairMicrosat microsat = ((C_GenomeEucaryote) microtus.getGenome()).getMicrosatXsome();
paire[locus] = (C_GenePair) microsat.getLocusAllele(locus);
}
genePopFile.writeln(paire[0] + " " + paire[1] + " " + paire[2]);
}
}
/** Output data in genetic indicators file */
private void recordGenesInFile() {
genesFile.write(RepastEssentials.GetTickCount() + ";");
for (int i = 0; i < I_sim_constants.NUMBER_GENES; i++)
genesFile.write(richnessByLocus[i] + ";");
for (int i = 0; i < I_sim_constants.NUMBER_GENES; i++)
genesFile.write(allelicFrequencies[i] + ";");
for (int i = 0; i < I_sim_constants.NUMBER_GENES; i++)
genesFile.write(expectedHeterozygosityByLocus[i] + ";");
genesFile.writeln(observedHeterozygosityHO + ";" + expectedHeterozygosityHE + ";"
+ getFixationIndex() + ";" + meanAllelicRichness);
}
/** close private files */
public void closeSimulation() {
meanAllelicRichness = 0;
expectedHeterozygosityHE = 0;
observedHeterozygosityHO = 0;
genesFile.writeln("0;0;0;0;0;0;0;0;0;0;0;" + RepastEssentials.GetTickCount() + ";0;0;0;"
+ genesFile.getNumRun());
genesFile.closeFile();
genePopFile.closeFile();
}
/** Utility: reduce to the requested number of digit after the decimal point(above 5 the double is automatically
* transformed in scientific format).
* @param val : the double to transform
* @param n : the number of digits after the decimal point */
private double convertNumber(double val, int n) {
final int MAX = n + 2; // n digits after the decimal point + decimal point + 0 before the point
final int facteur = (int) Math.pow(10, n); // to shift the decimal point of n digits
// Convert only if the number of digit after the dec. point exceeds n
if (((Double) val).toString().length() > MAX) {
int i = (int) (val * facteur); // On supprime les chiffres après n
return (double) i / facteur; // On replace la virgule
}
else return val;
}
/** getter
* Warning, the field is not updated at the call, this is done in step() procedure JLF 02.2013
* @see #step() */
public double getMeanAllelicRichness() {
return meanAllelicRichness;
}
}