    001    package calhoun.analysis.crf.features.supporting.phylogenetic;
    002
    003    import org.apache.commons.logging.Log;
    004    import org.apache.commons.logging.LogFactory;
    005
    006    import calhoun.util.Assert;
    007
    008    public class HKY85Model extends NucleotideSubstitutionModel {
    009            private static final long serialVersionUID = -1555350959092722736L;
    010            private static final Log log = LogFactory.getLog(HKY85Model.class);
    011            // Model for nucleotide substitution proposed by Kimura in 1980.
    012            // Models the overall rate of nucleotide substitution and the difference
    013            // between rate of transitions and transversions.
    014
    015            double ts;
    016            double tv;
    017            double piA;
    018            double piC;
    019            double piG;
    020            double piT;
    021
    022            /* HKY85 Nucleotide substitution model
    023             * 1. Website:         http://www.mun.ca/biology/scarr/Transitions_vs_Transversions.htm
    024             * 2. Reviw article:   Pietro Lio and Nick Goldman, "Models of Molecular Evolution and Phylogeny", Genome Research, 1998.
    025             * 3. Primary article: Hasegawa, Kishino, Yano, 1985, "Dating of the human-ape splitting by a molecular clock of mitochondrial DNA", J. Molecular Evolution. 22, 160-174.
    026             */
    027
    028            public HKY85Model(double[] parms) {
    029                    setParameters(parms);
    030            }
    031
    032            public HKY85Model() {
    033                    // Provide some reasonable seed values
    034                    ts  = 0.6;
    035                    tv  = 0.2;
    036                    piA = 0.25;
    037                    piC = 0.25;
    038                    piG = 0.25;
    039                    piT = 0.25;
    040            }
    041
    042            @Override
    043            public String getEvolutionaryModelName() {
    044                    return "HKY85_nucleotide_evolution_model";
    045            }
    046
    047            @Override
    048            public double[] getParameters() {
    049                    double[] ret = new double[6];
    050                    ret[0] = ts;
    051                    ret[1] = tv;
    052                    ret[2] = piA;
    053                    ret[3] = piC;
    054                    ret[4] = piG;
    055                    ret[5] = piT;
    056                    return ret;
    057            }
    058
    059
    060            @Override
    061            public void setParameters(double[] parms) {
    062                    Assert.a(parms.length == 5);
    063                    for (int j=0; j<5; j++) {
    064                            Assert.a(parms[j] > 0);
    065                    }
    066
    067                    ts  = parms[0];
    068                    tv  = parms[1];
    069                    piA = parms[2];
    070                    piC = parms[3];
    071                    piG = parms[4];
    072                    piT = 1 - piA - piC - piG;
    073                    Assert.a(piT>0);
    074
    075                    // Maybe here I should be using a hashing function, but let's first just write it out;
    076                    double[][] X = new double[4][4];
    077                    int A=0, C=1, G=2, T=3;
    078                    // Given an i base, X[i][j] is rate of changing to base j.
    079                    X[A][A]=0;         X[A][C]=tv*piC;   X[A][G]=ts*piG;   X[A][T]=tv*piT;
    080                    X[C][A]=tv*piA;    X[C][C]=0;        X[C][G]=tv*piG;   X[C][T]=ts*piT;
    081                    X[G][A]=ts*piA;    X[G][C]=tv*piC;   X[G][G]=0;        X[G][T]=tv*piT;
    082                    X[T][A]=tv*piA;    X[T][C]=ts*piC;   X[T][G]=tv*piG;   X[T][T]=0;
    083                    for (int i=0; i<4; i++) {
    084                            for (int j=0; j<4; j++) {
    085                                    if (j!=i) {
    086                                            X[i][i] -= X[i][j];
    087                                    }
    088                            }
    089                    }
    090
    091                    R.assign(X);
    092            }
    093
    094            @Override
    095            public void summarize() {
    096                    log.debug("HKY85 Nucleotide substitution rate matrix; ts=" + ts + "  tv=" + tv + "  piA="+piA+"  piC="+piC+"  piG="+piG+"  piT="+piT);
    097            }
    098    }