mev-4.0.01/source/org/tigr/microarray/mev/cluster/algorithm/impl/PavlidisTemplateMatching.java

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2 Copyright @ 1999-2003, The Institute for Genomic Research (TIGR).
2 26 Feb 07 jari 3 All rights reserved.
2 26 Feb 07 jari 4 */
2 26 Feb 07 jari 5 /*
2 26 Feb 07 jari 6  * $RCSfile: PavlidisTemplateMatching.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.4 $
2 26 Feb 07 jari 8  * $Date: 2005/03/10 15:45:20 $
2 26 Feb 07 jari 9  * $Author: braistedj $
2 26 Feb 07 jari 10  * $State: Exp $
2 26 Feb 07 jari 11  */
2 26 Feb 07 jari 12 package org.tigr.microarray.mev.cluster.algorithm.impl;
2 26 Feb 07 jari 13
2 26 Feb 07 jari 14 import java.util.Vector;
2 26 Feb 07 jari 15
2 26 Feb 07 jari 16 import javax.swing.JButton;
2 26 Feb 07 jari 17
2 26 Feb 07 jari 18 import org.tigr.microarray.mev.cluster.Cluster;
2 26 Feb 07 jari 19 import org.tigr.microarray.mev.cluster.Node;
2 26 Feb 07 jari 20 import org.tigr.microarray.mev.cluster.NodeList;
2 26 Feb 07 jari 21 import org.tigr.microarray.mev.cluster.NodeValue;
2 26 Feb 07 jari 22 import org.tigr.microarray.mev.cluster.NodeValueList;
2 26 Feb 07 jari 23 import org.tigr.microarray.mev.cluster.algorithm.AbortException;
2 26 Feb 07 jari 24 import org.tigr.microarray.mev.cluster.algorithm.AbstractAlgorithm;
2 26 Feb 07 jari 25 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmData;
2 26 Feb 07 jari 26 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmEvent;
2 26 Feb 07 jari 27 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmException;
2 26 Feb 07 jari 28 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmParameters;
2 26 Feb 07 jari 29 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 30
2 26 Feb 07 jari 31 import JSci.maths.statistics.TDistribution;
2 26 Feb 07 jari 32
2 26 Feb 07 jari 33 public class PavlidisTemplateMatching extends AbstractAlgorithm {
2 26 Feb 07 jari 34     private boolean stop = false;
2 26 Feb 07 jari 35     
2 26 Feb 07 jari 36     private int function;
2 26 Feb 07 jari 37     private float factor;
2 26 Feb 07 jari 38     private boolean absolute;
2 26 Feb 07 jari 39     
2 26 Feb 07 jari 40     private int number_of_genes;
2 26 Feb 07 jari 41     private int number_of_samples;
2 26 Feb 07 jari 42     
2 26 Feb 07 jari 43     private FloatMatrix expMatrix;
2 26 Feb 07 jari 44     private FloatMatrix newMatrix;
2 26 Feb 07 jari 45     
2 26 Feb 07 jari 46     private long StartTime;
2 26 Feb 07 jari 47     private long CalculationTime;
2 26 Feb 07 jari 48     private boolean ptmGenes;
2 26 Feb 07 jari 49     
2 26 Feb 07 jari 50     private double[] pValues, rValues;
2 26 Feb 07 jari 51     
2 26 Feb 07 jari 52     
2 26 Feb 07 jari 53     boolean useAbsolute; // true = use absolute value of correlation; false = use the signed value of the correlation
2 26 Feb 07 jari 54     boolean useR, drawSigTreesOnly;
2 26 Feb 07 jari 55     JButton abortButton;
2 26 Feb 07 jari 56     
2 26 Feb 07 jari 57     Vector templateVector;
2 26 Feb 07 jari 58     FloatMatrix templateVectorMatrix; //NEED TO IMPORT TEMPLATE VECTOR AS A MATRIX BECAUSE "ALGORITHMDATA" CLASS ONLY HAS getMatrix() METHOD
2 26 Feb 07 jari 59     //NO WAY TO IMPORT A VECTOR DIRECTLY, EITHER USING CONFMAP OR ALGORITHMDATA CLASS
2 26 Feb 07 jari 60     float[] geneTemplate;    
2 26 Feb 07 jari 61     float threshold;    
2 26 Feb 07 jari 62     int origNumGenes, numSamples;
2 26 Feb 07 jari 63     
2 26 Feb 07 jari 64     private int hcl_function;
2 26 Feb 07 jari 65     private boolean hcl_absolute;
2 26 Feb 07 jari 66     
2 26 Feb 07 jari 67     public AlgorithmData execute(AlgorithmData data) throws AlgorithmException {
2 26 Feb 07 jari 68   
2 26 Feb 07 jari 69   
2 26 Feb 07 jari 70   AlgorithmParameters map = data.getParams();
2 26 Feb 07 jari 71   
2 26 Feb 07 jari 72   function = map.getInt("distance-function", EUCLIDEAN);
2 26 Feb 07 jari 73   factor   = map.getFloat("distance-factor", 1.0f);
2 26 Feb 07 jari 74   absolute = map.getBoolean("distance-absolute", false);
2 26 Feb 07 jari 75                 
2 26 Feb 07 jari 76         hcl_function = map.getInt("hcl-distance-function", EUCLIDEAN);
2 26 Feb 07 jari 77         hcl_absolute = map.getBoolean("hcl-distance-absolute", false);
2 26 Feb 07 jari 78         
2 26 Feb 07 jari 79         ptmGenes = map.getBoolean("ptm-cluster-genes", true);
2 26 Feb 07 jari 80   threshold = map.getFloat("threshold", 0.8f);
2 26 Feb 07 jari 81   useAbsolute = map.getBoolean("use-absolute", false);
2 26 Feb 07 jari 82   useR = map.getBoolean("useR", false);
2 26 Feb 07 jari 83   templateVectorMatrix = data.getMatrix("templateVectorMatrix");
2 26 Feb 07 jari 84   templateVector = convertToVector(templateVectorMatrix);
2 26 Feb 07 jari 85   
2 26 Feb 07 jari 86   boolean hierarchical_tree = map.getBoolean("hierarchical-tree", false);
2 26 Feb 07 jari 87         if (hierarchical_tree) {
2 26 Feb 07 jari 88             drawSigTreesOnly = map.getBoolean("draw-sig-trees-only");
2 26 Feb 07 jari 89         }         
2 26 Feb 07 jari 90   int method_linkage = map.getInt("method-linkage", 0);
2 26 Feb 07 jari 91   boolean calculate_genes = map.getBoolean("calculate-genes", false);
2 26 Feb 07 jari 92   boolean calculate_experiments = map.getBoolean("calculate-experiments", false);
2 26 Feb 07 jari 93   
2 26 Feb 07 jari 94   this.expMatrix = data.getMatrix("experiment");
2 26 Feb 07 jari 95   
2 26 Feb 07 jari 96   number_of_genes   = this.expMatrix.getRowDimension();
2 26 Feb 07 jari 97   number_of_samples = this.expMatrix.getColumnDimension();
2 26 Feb 07 jari 98   
2 26 Feb 07 jari 99   origNumGenes = this.expMatrix.getRowDimension();
2 26 Feb 07 jari 100   numSamples = this.expMatrix.getColumnDimension();
2 26 Feb 07 jari 101   
2 26 Feb 07 jari 102   geneTemplate = new float[templateVector.size()];
2 26 Feb 07 jari 103         
2 26 Feb 07 jari 104         pValues = new double[number_of_genes];
2 26 Feb 07 jari 105         rValues = new double[number_of_genes];
2 26 Feb 07 jari 106   
2 26 Feb 07 jari 107   for (int i = 0; i < templateVector.size(); i++) {
2 26 Feb 07 jari 108       geneTemplate[i] = ((Float)templateVector.get(i)).floatValue();
2 26 Feb 07 jari 109   }
2 26 Feb 07 jari 110   
2 26 Feb 07 jari 111   
2 26 Feb 07 jari 112   Vector clusters[] = calculate();
2 26 Feb 07 jari 113   
2 26 Feb 07 jari 114   FloatMatrix means = getMeans(clusters);
2 26 Feb 07 jari 115   FloatMatrix variances = getVariances(clusters, means);
2 26 Feb 07 jari 116   
2 26 Feb 07 jari 117   AlgorithmEvent event = null;
2 26 Feb 07 jari 118   if (hierarchical_tree) {
2 26 Feb 07 jari 119       event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, clusters.length, "Calculate Hierarchical Trees");
2 26 Feb 07 jari 120       fireValueChanged(event);
2 26 Feb 07 jari 121       event.setIntValue(0);
2 26 Feb 07 jari 122       event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 123       fireValueChanged(event);
2 26 Feb 07 jari 124   }
2 26 Feb 07 jari 125   
2 26 Feb 07 jari 126   Cluster result_cluster = new Cluster();
2 26 Feb 07 jari 127   NodeList nodeList = result_cluster.getNodeList();
2 26 Feb 07 jari 128   int[] features;
2 26 Feb 07 jari 129   for (int i=0; i<clusters.length; i++) {
2 26 Feb 07 jari 130       if (stop) {
2 26 Feb 07 jari 131     throw new AbortException();
2 26 Feb 07 jari 132       }
2 26 Feb 07 jari 133       features = convert2int(clusters[i]);
2 26 Feb 07 jari 134       Node node = new Node(features);
2 26 Feb 07 jari 135       nodeList.addNode(node);
2 26 Feb 07 jari 136       if (hierarchical_tree) {
2 26 Feb 07 jari 137                 if (drawSigTreesOnly) {
2 26 Feb 07 jari 138                     if (i == 0) {
2 26 Feb 07 jari 139                         node.setValues(calculateHierarchicalTree(features, method_linkage, calculate_genes, calculate_experiments));
2 26 Feb 07 jari 140                         event.setIntValue(i+1);
2 26 Feb 07 jari 141                         fireValueChanged(event);                       
2 26 Feb 07 jari 142                     }
2 26 Feb 07 jari 143                 } else {
2 26 Feb 07 jari 144                     node.setValues(calculateHierarchicalTree(features, method_linkage, calculate_genes, calculate_experiments));
2 26 Feb 07 jari 145                     event.setIntValue(i+1);
2 26 Feb 07 jari 146                     fireValueChanged(event);
2 26 Feb 07 jari 147                 }
2 26 Feb 07 jari 148       }
2 26 Feb 07 jari 149   }
2 26 Feb 07 jari 150         
2 26 Feb 07 jari 151         FloatMatrix rValuesMatrix = new FloatMatrix(number_of_genes, 1);
2 26 Feb 07 jari 152         FloatMatrix pValuesMatrix = new FloatMatrix(number_of_genes, 1);     
2 26 Feb 07 jari 153         
2 26 Feb 07 jari 154         for (int i = 0; i < pValues.length; i++) {
2 26 Feb 07 jari 155             rValuesMatrix.A[i][0] = (float)(rValues[i]);
2 26 Feb 07 jari 156             pValuesMatrix.A[i][0] = (float)(pValues[i]);
2 26 Feb 07 jari 157         }
2 26 Feb 07 jari 158   
2 26 Feb 07 jari 159   // prepare the result
2 26 Feb 07 jari 160   AlgorithmData result = new AlgorithmData();
2 26 Feb 07 jari 161   result.addCluster("cluster", result_cluster);
2 26 Feb 07 jari 162   result.addMatrix("clusters_means", means);
2 26 Feb 07 jari 163   result.addMatrix("clusters_variances", variances);
2 26 Feb 07 jari 164         result.addMatrix("rValuesMatrix", rValuesMatrix);
2 26 Feb 07 jari 165         result.addMatrix("pValuesMatrix", pValuesMatrix);
2 26 Feb 07 jari 166   return result;
2 26 Feb 07 jari 167     }
2 26 Feb 07 jari 168     
2 26 Feb 07 jari 169     
2 26 Feb 07 jari 170     private NodeValueList calculateHierarchicalTree(int[] features, int method, boolean genes, boolean experiments) throws AlgorithmException {
2 26 Feb 07 jari 171   NodeValueList nodeList = new NodeValueList();
2 26 Feb 07 jari 172   AlgorithmData data = new AlgorithmData();
2 26 Feb 07 jari 173   FloatMatrix experiment;
2 26 Feb 07 jari 174         
2 26 Feb 07 jari 175         if(ptmGenes)
2 26 Feb 07 jari 176             experiment = getSubExperiment(this.expMatrix, features);
2 26 Feb 07 jari 177         else
2 26 Feb 07 jari 178             experiment = this.getSubExperimentReducedCols(this.expMatrix, features);
2 26 Feb 07 jari 179         
2 26 Feb 07 jari 180         data.addMatrix("experiment", experiment);
2 26 Feb 07 jari 181         data.addParam("hcl-distance-function", String.valueOf(this.hcl_function));
2 26 Feb 07 jari 182         data.addParam("hcl-distance-absolute", String.valueOf(this.hcl_absolute));
2 26 Feb 07 jari 183   data.addParam("method-linkage", String.valueOf(method));
2 26 Feb 07 jari 184         
2 26 Feb 07 jari 185   HCL hcl = new HCL();
2 26 Feb 07 jari 186   AlgorithmData result;
2 26 Feb 07 jari 187   if (genes) {
2 26 Feb 07 jari 188       data.addParam("calculate-genes", String.valueOf(true));
2 26 Feb 07 jari 189       result = hcl.execute(data);
2 26 Feb 07 jari 190       validate(result);
2 26 Feb 07 jari 191       addNodeValues(nodeList, result);
2 26 Feb 07 jari 192   }
2 26 Feb 07 jari 193   if (experiments) {
2 26 Feb 07 jari 194       data.addParam("calculate-genes", String.valueOf(false));
2 26 Feb 07 jari 195       result = hcl.execute(data);
2 26 Feb 07 jari 196       validate(result);
2 26 Feb 07 jari 197       addNodeValues(nodeList, result);
2 26 Feb 07 jari 198   }
2 26 Feb 07 jari 199   return nodeList;
2 26 Feb 07 jari 200     }
2 26 Feb 07 jari 201     
2 26 Feb 07 jari 202     private void addNodeValues(NodeValueList target_list, AlgorithmData source_result) {
2 26 Feb 07 jari 203   target_list.addNodeValue(new NodeValue("child-1-array", source_result.getIntArray("child-1-array")));
2 26 Feb 07 jari 204   target_list.addNodeValue(new NodeValue("child-2-array", source_result.getIntArray("child-2-array")));
2 26 Feb 07 jari 205   target_list.addNodeValue(new NodeValue("node-order", source_result.getIntArray("node-order")));
2 26 Feb 07 jari 206   target_list.addNodeValue(new NodeValue("height", source_result.getMatrix("height").getRowPackedCopy()));
2 26 Feb 07 jari 207     }
2 26 Feb 07 jari 208     
2 26 Feb 07 jari 209     private FloatMatrix getSubExperiment(FloatMatrix experiment, int[] features) {
2 26 Feb 07 jari 210   FloatMatrix subExperiment = new FloatMatrix(features.length, experiment.getColumnDimension());
2 26 Feb 07 jari 211   for (int i=0; i<features.length; i++) {
2 26 Feb 07 jari 212       subExperiment.A[i] = experiment.A[features[i]];
2 26 Feb 07 jari 213   }
2 26 Feb 07 jari 214   return subExperiment;
2 26 Feb 07 jari 215     }
2 26 Feb 07 jari 216     
2 26 Feb 07 jari 217         /**
2 26 Feb 07 jari 218      *  Creates a matrix with reduced columns (samples) as during experiment clustering
2 26 Feb 07 jari 219      */
2 26 Feb 07 jari 220     private FloatMatrix getSubExperimentReducedCols(FloatMatrix experiment, int[] features) {
2 26 Feb 07 jari 221         FloatMatrix copyMatrix = experiment.copy();
2 26 Feb 07 jari 222         FloatMatrix subExperiment = new FloatMatrix(features.length, copyMatrix.getColumnDimension());
2 26 Feb 07 jari 223         for (int i=0; i<features.length; i++) {
2 26 Feb 07 jari 224             subExperiment.A[i] = copyMatrix.A[features[i]];
2 26 Feb 07 jari 225         }
2 26 Feb 07 jari 226         subExperiment = subExperiment.transpose();
2 26 Feb 07 jari 227         return subExperiment;
2 26 Feb 07 jari 228     }
2 26 Feb 07 jari 229     
2 26 Feb 07 jari 230     /**
2 26 Feb 07 jari 231      * Checking the result of hcl algorithm calculation.
2 26 Feb 07 jari 232      * @throws AlgorithmException, if the result is incorrect.
2 26 Feb 07 jari 233      */
2 26 Feb 07 jari 234     private void validate(AlgorithmData result) throws AlgorithmException {
2 26 Feb 07 jari 235   if (result.getIntArray("child-1-array") == null) {
2 26 Feb 07 jari 236       throw new AlgorithmException("parameter 'child-1-array' is null");
2 26 Feb 07 jari 237   }
2 26 Feb 07 jari 238   if (result.getIntArray("child-2-array") == null) {
2 26 Feb 07 jari 239       throw new AlgorithmException("parameter 'child-2-array' is null");
2 26 Feb 07 jari 240   }
2 26 Feb 07 jari 241   if (result.getIntArray("node-order") == null) {
2 26 Feb 07 jari 242       throw new AlgorithmException("parameter 'node-order' is null");
2 26 Feb 07 jari 243   }
2 26 Feb 07 jari 244   if (result.getMatrix("height") == null) {
2 26 Feb 07 jari 245       throw new AlgorithmException("parameter 'height' is null");
2 26 Feb 07 jari 246   }
2 26 Feb 07 jari 247     }
2 26 Feb 07 jari 248     
2 26 Feb 07 jari 249     private int[] convert2int(Vector source) {
2 26 Feb 07 jari 250   int[] int_matrix = new int[source.size()];
2 26 Feb 07 jari 251   for (int i=0; i<int_matrix.length; i++) {
2 26 Feb 07 jari 252       int_matrix[i] = ((Integer) source.get(i)).intValue();
2 26 Feb 07 jari 253   }
2 26 Feb 07 jari 254   return int_matrix;
2 26 Feb 07 jari 255     }
2 26 Feb 07 jari 256     
2 26 Feb 07 jari 257     public void abort() {
2 26 Feb 07 jari 258   stop = true;
2 26 Feb 07 jari 259     }
2 26 Feb 07 jari 260     
2 26 Feb 07 jari 261     private Vector convertToVector(FloatMatrix tempMatrix) {
2 26 Feb 07 jari 262   Vector temp = new Vector();
2 26 Feb 07 jari 263   for (int i = 0; i < tempMatrix.A[0].length; i++) {
2 26 Feb 07 jari 264       temp.add(new Float(tempMatrix.A[0][i]));
2 26 Feb 07 jari 265   }
2 26 Feb 07 jari 266   return temp;
2 26 Feb 07 jari 267     }
2 26 Feb 07 jari 268     
2 26 Feb 07 jari 269     
2 26 Feb 07 jari 270     public Vector itf(Vector integerVector) {
2 26 Feb 07 jari 271   Vector floatVector = new Vector();
2 26 Feb 07 jari 272   
2 26 Feb 07 jari 273   for (int i = 0; i < integerVector.size(); i++) {
2 26 Feb 07 jari 274       floatVector.addElement(new Float(((Integer) integerVector.elementAt(i)).intValue()));
2 26 Feb 07 jari 275   }
2 26 Feb 07 jari 276   
2 26 Feb 07 jari 277   return floatVector;
2 26 Feb 07 jari 278     }
2 26 Feb 07 jari 279     
2 26 Feb 07 jari 280     
2 26 Feb 07 jari 281     Vector[] calculate() {
2 26 Feb 07 jari 282   
2 26 Feb 07 jari 283   Vector allUniqueIDIndices = new Vector();
2 26 Feb 07 jari 284   for(int i = 0; i < origNumGenes; i++) {
2 26 Feb 07 jari 285       allUniqueIDIndices.add(new Integer(i));
2 26 Feb 07 jari 286   }
2 26 Feb 07 jari 287   
2 26 Feb 07 jari 288   Vector remainingGenes = (Vector)(allUniqueIDIndices).clone();
2 26 Feb 07 jari 289   StartTime=System.currentTimeMillis();
2 26 Feb 07 jari 290   Vector similarGenes = findSimilarGenes();
2 26 Feb 07 jari 291   remainingGenes.removeAll(similarGenes);
2 26 Feb 07 jari 292   CalculationTime=System.currentTimeMillis()-StartTime;
2 26 Feb 07 jari 293   
2 26 Feb 07 jari 294   Vector[] clusters = new Vector[2];
2 26 Feb 07 jari 295   clusters[0] = similarGenes;
2 26 Feb 07 jari 296   clusters[1] = remainingGenes;
2 26 Feb 07 jari 297   
2 26 Feb 07 jari 298   return clusters;
2 26 Feb 07 jari 299     }
2 26 Feb 07 jari 300     
2 26 Feb 07 jari 301     
2 26 Feb 07 jari 302     FloatMatrix addTemplateToexpMatrix() {
2 26 Feb 07 jari 303   FloatMatrix newMatrix = new FloatMatrix(origNumGenes + 1, numSamples);
2 26 Feb 07 jari 304   for (int i = 0; i < origNumGenes; i++) { //copy all the elements of expMatrix into newMatrix
2 26 Feb 07 jari 305       for (int j = 0; j < numSamples; j++) {
2 26 Feb 07 jari 306     newMatrix.A[i][j] = expMatrix.A[i][j];
2 26 Feb 07 jari 307       }
2 26 Feb 07 jari 308   }
2 26 Feb 07 jari 309   
2 26 Feb 07 jari 310   for (int k = 0; k < geneTemplate.length; k++) {
2 26 Feb 07 jari 311       newMatrix.A[origNumGenes][k] = geneTemplate[k];
2 26 Feb 07 jari 312   }
2 26 Feb 07 jari 313   /*
2 26 Feb 07 jari 314   for (int i = 0; i < origNumGenes + 1; i++) {
2 26 Feb 07 jari 315       for (int j = 0; j < numSamples; j++) {
2 26 Feb 07 jari 316     
2 26 Feb 07 jari 317       }
2 26 Feb 07 jari 318   }
2 26 Feb 07 jari 319          */
2 26 Feb 07 jari 320   
2 26 Feb 07 jari 321   return newMatrix;
2 26 Feb 07 jari 322     }
2 26 Feb 07 jari 323     
2 26 Feb 07 jari 324     
2 26 Feb 07 jari 325     Vector findSimilarGenes() {
2 26 Feb 07 jari 326   Vector similarGenes = new Vector();
2 26 Feb 07 jari 327   double pearsonR;
2 26 Feb 07 jari 328   newMatrix = addTemplateToexpMatrix();
2 26 Feb 07 jari 329   
2 26 Feb 07 jari 330   
2 26 Feb 07 jari 331   Vector allUniqueIDIndices = new Vector();
2 26 Feb 07 jari 332   for(int i = 0; i < origNumGenes; i++) {
2 26 Feb 07 jari 333       allUniqueIDIndices.add(new Integer(i));
2 26 Feb 07 jari 334   }
2 26 Feb 07 jari 335   
2 26 Feb 07 jari 336   for (int i = 0; i < origNumGenes; i++){
2 26 Feb 07 jari 337       pearsonR = ExperimentUtil.genePearson(newMatrix, null, i, origNumGenes, factor);
2 26 Feb 07 jari 338             rValues[i] = pearsonR;
2 26 Feb 07 jari 339             pValues[i] = getProb(pearsonR);
2 26 Feb 07 jari 340       if(useR) {
2 26 Feb 07 jari 341     
2 26 Feb 07 jari 342     if (useAbsolute == true) {
2 26 Feb 07 jari 343         pearsonR = Math.abs(pearsonR);
2 26 Feb 07 jari 344         if (pearsonR >= threshold) {
2 26 Feb 07 jari 345       similarGenes.add(allUniqueIDIndices.elementAt(i));
2 26 Feb 07 jari 346         }
2 26 Feb 07 jari 347     } else {
2 26 Feb 07 jari 348         if (pearsonR >= threshold) {
2 26 Feb 07 jari 349       similarGenes.add(allUniqueIDIndices.elementAt(i));
2 26 Feb 07 jari 350         }
2 26 Feb 07 jari 351     }
2 26 Feb 07 jari 352     
2 26 Feb 07 jari 353       } else { // if (!useR)
2 26 Feb 07 jari 354     
2 26 Feb 07 jari 355     if (useAbsolute == true) {
2 26 Feb 07 jari 356         if (getProb(pearsonR)<= threshold) {
2 26 Feb 07 jari 357       similarGenes.add(allUniqueIDIndices.elementAt(i));
2 26 Feb 07 jari 358         }
2 26 Feb 07 jari 359         
2 26 Feb 07 jari 360     } else {
2 26 Feb 07 jari 361         if ((pearsonR > 0)&&(getProb(pearsonR)<= threshold)) {
2 26 Feb 07 jari 362       similarGenes.add(allUniqueIDIndices.elementAt(i));
2 26 Feb 07 jari 363         }
2 26 Feb 07 jari 364     }
2 26 Feb 07 jari 365       }
2 26 Feb 07 jari 366       
2 26 Feb 07 jari 367   }
2 26 Feb 07 jari 368   
2 26 Feb 07 jari 369   return similarGenes;
2 26 Feb 07 jari 370     }
2 26 Feb 07 jari 371     
2 26 Feb 07 jari 372     
2 26 Feb 07 jari 373     double getProb(double pearsonR){
2 26 Feb 07 jari 374   
2 26 Feb 07 jari 375   double prob;
2 26 Feb 07 jari 376   double tValue = getTValue(Math.abs(pearsonR));
2 26 Feb 07 jari 377   int df = geneTemplate.length - 2;
2 26 Feb 07 jari 378   
2 26 Feb 07 jari 379   TDistribution tDist = new TDistribution(df);
2 26 Feb 07 jari 380   //prob = tDist.probability(tValue);
2 26 Feb 07 jari 381         prob = 2*(1-tDist.cumulative(tValue));
2 26 Feb 07 jari 382         
2 26 Feb 07 jari 383         if (prob > 1.0d) {
2 26 Feb 07 jari 384             prob = 1.0d;
2 26 Feb 07 jari 385         }
2 26 Feb 07 jari 386   
2 26 Feb 07 jari 387   return prob;
2 26 Feb 07 jari 388     }
2 26 Feb 07 jari 389     //
2 26 Feb 07 jari 390     
2 26 Feb 07 jari 391     double getTValue(double pearsonR){ //as explained on pg 333, Jaccard and Becker
2 26 Feb 07 jari 392   double tValue;
2 26 Feb 07 jari 393   double stdError;
2 26 Feb 07 jari 394   int n = geneTemplate.length;
2 26 Feb 07 jari 395   
2 26 Feb 07 jari 396   stdError = Math.sqrt((1-pearsonR*pearsonR)/(n - 2));
2 26 Feb 07 jari 397   
2 26 Feb 07 jari 398   tValue = pearsonR/stdError;
2 26 Feb 07 jari 399   
2 26 Feb 07 jari 400   return tValue;
2 26 Feb 07 jari 401     }
2 26 Feb 07 jari 402     
2 26 Feb 07 jari 403     
2 26 Feb 07 jari 404     private FloatMatrix getMeans(Vector[] clusters) {
2 26 Feb 07 jari 405   FloatMatrix means = new FloatMatrix(clusters.length, number_of_samples);
2 26 Feb 07 jari 406   FloatMatrix mean;
2 26 Feb 07 jari 407   for (int i=0; i<clusters.length; i++) {
2 26 Feb 07 jari 408       mean = getMean(clusters[i]);
2 26 Feb 07 jari 409       means.A[i] = mean.A[0];
2 26 Feb 07 jari 410   }
2 26 Feb 07 jari 411   return means;
2 26 Feb 07 jari 412     }
2 26 Feb 07 jari 413     
2 26 Feb 07 jari 414     private FloatMatrix getMean(Vector cluster) {
2 26 Feb 07 jari 415   FloatMatrix mean = new FloatMatrix(1, number_of_samples);
2 26 Feb 07 jari 416   float currentMean;
2 26 Feb 07 jari 417   int n = cluster.size();
2 26 Feb 07 jari 418   int denom = 0;
2 26 Feb 07 jari 419   float value;
2 26 Feb 07 jari 420   for (int i=0; i<number_of_samples; i++) {
2 26 Feb 07 jari 421       currentMean = 0f;
2 26 Feb 07 jari 422       denom  = 0;
2 26 Feb 07 jari 423       for (int j=0; j<n; j++) {
2 26 Feb 07 jari 424     value = expMatrix.get(((Integer) cluster.get(j)).intValue(), i);
2 26 Feb 07 jari 425     if (!Float.isNaN(value)) {
2 26 Feb 07 jari 426         currentMean += value;
2 26 Feb 07 jari 427         denom++;
2 26 Feb 07 jari 428     }
2 26 Feb 07 jari 429       }
2 26 Feb 07 jari 430       mean.set(0, i, currentMean/(float)denom);
2 26 Feb 07 jari 431   }
2 26 Feb 07 jari 432   
2 26 Feb 07 jari 433   return mean;
2 26 Feb 07 jari 434     }
2 26 Feb 07 jari 435     
2 26 Feb 07 jari 436     private FloatMatrix getVariances(Vector[] clusters, FloatMatrix means) {
2 26 Feb 07 jari 437   final int rows = means.getRowDimension();
2 26 Feb 07 jari 438   final int columns = means.getColumnDimension();
2 26 Feb 07 jari 439   FloatMatrix variances = new FloatMatrix(rows, columns);
2 26 Feb 07 jari 440   for (int row=0; row<rows; row++) {
2 26 Feb 07 jari 441       for (int column=0; column<columns; column++) {
2 26 Feb 07 jari 442     variances.set(row, column, getSampleVariance(clusters[row], column, means.get(row, column)));
2 26 Feb 07 jari 443       }
2 26 Feb 07 jari 444   }
2 26 Feb 07 jari 445   return variances;
2 26 Feb 07 jari 446     }
2 26 Feb 07 jari 447     
2 26 Feb 07 jari 448     int validN;
2 26 Feb 07 jari 449     
2 26 Feb 07 jari 450     private float getSampleNormalizedSum(Vector cluster, int column, float mean) {
2 26 Feb 07 jari 451   final int size = cluster.size();
2 26 Feb 07 jari 452   float sum = 0f;
2 26 Feb 07 jari 453   validN = 0;
2 26 Feb 07 jari 454   float value;
2 26 Feb 07 jari 455   for (int i=0; i<size; i++) {
2 26 Feb 07 jari 456       value = expMatrix.get(((Integer) cluster.get(i)).intValue(), column);
2 26 Feb 07 jari 457       if (!Float.isNaN(value)) {
2 26 Feb 07 jari 458     sum += Math.pow(value-mean, 2);
2 26 Feb 07 jari 459     validN++;
2 26 Feb 07 jari 460       }
2 26 Feb 07 jari 461   }
2 26 Feb 07 jari 462   return sum;
2 26 Feb 07 jari 463     }
2 26 Feb 07 jari 464     
2 26 Feb 07 jari 465     private float getSampleVariance(Vector cluster, int column, float mean) {
2 26 Feb 07 jari 466   return(float)Math.sqrt(getSampleNormalizedSum(cluster, column, mean)/(float)(validN-1));
2 26 Feb 07 jari 467     }
2 26 Feb 07 jari 468     
2 26 Feb 07 jari 469 }