plugins/base1/se.lu.onk/trunk/ZTest/src/ztest/ZTest.java

Code
Comments
Other
Rev Date Author Line
187 13 Oct 06 enell 1 /*
187 13 Oct 06 enell 2  $Id$
187 13 Oct 06 enell 3
187 13 Oct 06 enell 4  Copyright (C) 2006 Johan Enell
187 13 Oct 06 enell 5
187 13 Oct 06 enell 6  This file is part of BASE - BioArray Software Environment.
187 13 Oct 06 enell 7  Available at http://base.thep.lu.se/
187 13 Oct 06 enell 8
187 13 Oct 06 enell 9  BASE is free software; you can redistribute it and/or modify it
187 13 Oct 06 enell 10  under the terms of the GNU General Public License as published by
187 13 Oct 06 enell 11  the Free Software Foundation; either version 2 of the License, or
187 13 Oct 06 enell 12  (at your option) any later version.
187 13 Oct 06 enell 13
187 13 Oct 06 enell 14  BASE is distributed in the hope that it will be useful, but
187 13 Oct 06 enell 15  WITHOUT ANY WARRANTY; without even the implied warranty of
187 13 Oct 06 enell 16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
187 13 Oct 06 enell 17  General Public License for more details.
187 13 Oct 06 enell 18
187 13 Oct 06 enell 19  You should have received a copy of the GNU General Public License
187 13 Oct 06 enell 20  along with this program; if not, write to the Free Software
187 13 Oct 06 enell 21  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
187 13 Oct 06 enell 22  02111-1307, USA.
187 13 Oct 06 enell 23  */
187 13 Oct 06 enell 24 package ztest;
187 13 Oct 06 enell 25
203 24 Nov 06 enell 26 import org.apache.commons.math.MathException;
203 24 Nov 06 enell 27 import org.apache.commons.math.special.Erf;
203 24 Nov 06 enell 28
187 13 Oct 06 enell 29 import java.util.HashMap;
187 13 Oct 06 enell 30 import java.util.List;
187 13 Oct 06 enell 31
187 13 Oct 06 enell 32 import basefile.BASEFileSpotSection;
187 13 Oct 06 enell 33
187 13 Oct 06 enell 34 public class ZTest
187 13 Oct 06 enell 35 {
217 08 Dec 06 enell 36   private HashMap<Integer, Double> means;
217 08 Dec 06 enell 37   
217 08 Dec 06 enell 38   private HashMap<Integer, Double> sds;
217 08 Dec 06 enell 39   
187 13 Oct 06 enell 40   public void calculate(BASEFileSpotSection<Reporter, Spot> bfss)
187 13 Oct 06 enell 41   {
187 13 Oct 06 enell 42     List<Integer> assays = bfss.findFieldIntList("assays");
217 08 Dec 06 enell 43     means = getAssayMeans(bfss, assays);
217 08 Dec 06 enell 44     sds = getAssaySDS(bfss, assays, means);
197 21 Nov 06 enell 45     
187 13 Oct 06 enell 46     for (int i = 0; i < bfss.getReporterSize(); i++)
187 13 Oct 06 enell 47     {
187 13 Oct 06 enell 48       int popSize = 0;
187 13 Oct 06 enell 49       double popSD = 0;
187 13 Oct 06 enell 50       double popMean = 0;
187 13 Oct 06 enell 51       double meanM = 0;
187 13 Oct 06 enell 52
187 13 Oct 06 enell 53       for (int j = 0; j < assays.size(); j++)
187 13 Oct 06 enell 54       {
197 21 Nov 06 enell 55         double m = bfss.getSpot((i*assays.size()) + j).getM();
187 13 Oct 06 enell 56         if (!Double.isNaN(m))
187 13 Oct 06 enell 57         {
187 13 Oct 06 enell 58           meanM += m;
187 13 Oct 06 enell 59           popMean += means.get(assays.get(j));
187 13 Oct 06 enell 60           popSD += Math.pow(sds.get(assays.get(j)), 2);
187 13 Oct 06 enell 61           popSize++;
187 13 Oct 06 enell 62         }
187 13 Oct 06 enell 63       }
197 21 Nov 06 enell 64       popSD = Math.sqrt(popSD) / popSize;
187 13 Oct 06 enell 65       popMean = popMean / popSize;
187 13 Oct 06 enell 66       meanM = meanM / popSize;
187 13 Oct 06 enell 67
187 13 Oct 06 enell 68       double z = z(meanM, popMean, popSD);
197 21 Nov 06 enell 69       
187 13 Oct 06 enell 70       double normsdist = normsdist(Math.abs(z));
187 13 Oct 06 enell 71       double p = p(normsdist);
197 21 Nov 06 enell 72
197 21 Nov 06 enell 73       Reporter r = bfss.getReporter(i); 
197 21 Nov 06 enell 74       r.setP(p);
197 21 Nov 06 enell 75       r.setM(meanM);
197 21 Nov 06 enell 76       r.setSize(popSize);
187 13 Oct 06 enell 77     }
187 13 Oct 06 enell 78   }
217 08 Dec 06 enell 79   
217 08 Dec 06 enell 80   public double getTotalSD()
217 08 Dec 06 enell 81   {
217 08 Dec 06 enell 82     double sd = 0;
217 08 Dec 06 enell 83     for(Double d : sds.values())
217 08 Dec 06 enell 84     {
217 08 Dec 06 enell 85       sd += Math.pow(d, 2);
217 08 Dec 06 enell 86     }
217 08 Dec 06 enell 87     return Math.sqrt(sd) / sds.size();
217 08 Dec 06 enell 88   }
217 08 Dec 06 enell 89   
217 08 Dec 06 enell 90   public double getTotalMean()
217 08 Dec 06 enell 91   {
217 08 Dec 06 enell 92     double mean = 0;
217 08 Dec 06 enell 93     for(Double d : means.values())
217 08 Dec 06 enell 94     {
217 08 Dec 06 enell 95       mean += d;
217 08 Dec 06 enell 96     }
217 08 Dec 06 enell 97     return mean/means.size();
217 08 Dec 06 enell 98   }
187 13 Oct 06 enell 99
187 13 Oct 06 enell 100   private double p(double normsdist)
187 13 Oct 06 enell 101   {
206 24 Nov 06 enell 102     return 2.0 * (1.0 - normsdist);
187 13 Oct 06 enell 103   }
187 13 Oct 06 enell 104
187 13 Oct 06 enell 105   private double normsdist(double z)
187 13 Oct 06 enell 106   {
203 24 Nov 06 enell 107     double value = 0;
203 24 Nov 06 enell 108     try
197 21 Nov 06 enell 109     {
203 24 Nov 06 enell 110       if (z < -26) value = 0;
203 24 Nov 06 enell 111       else if (z > 26) value = 1;
203 24 Nov 06 enell 112       else value = 0.5 + Erf.erf(z/Math.sqrt(2))/2;
197 21 Nov 06 enell 113     }
203 24 Nov 06 enell 114     catch (MathException e)
203 24 Nov 06 enell 115     {
203 24 Nov 06 enell 116       System.err.println("This should never happen but if it do... The plug-in is crap.");
203 24 Nov 06 enell 117       System.err.println(z);
203 24 Nov 06 enell 118       e.printStackTrace();
203 24 Nov 06 enell 119       System.exit(0);
203 24 Nov 06 enell 120     }
197 21 Nov 06 enell 121     return value;
187 13 Oct 06 enell 122   }
197 21 Nov 06 enell 123   
187 13 Oct 06 enell 124   private double z(double x, double popMean, double popSD)
187 13 Oct 06 enell 125   {
187 13 Oct 06 enell 126     return (x - popMean) / popSD;
187 13 Oct 06 enell 127   }
187 13 Oct 06 enell 128
197 21 Nov 06 enell 129   private HashMap<Integer, Double> getAssaySDS(BASEFileSpotSection<Reporter, Spot> bfss, List<Integer> assays, HashMap<Integer, Double> means)
187 13 Oct 06 enell 130   {
197 21 Nov 06 enell 131     HashMap<Integer, Double> sds = new HashMap<Integer, Double>();
197 21 Nov 06 enell 132     for (int i = 0; i < assays.size(); i++)
187 13 Oct 06 enell 133     {
197 21 Nov 06 enell 134       int size = 0;
197 21 Nov 06 enell 135       double sqrSum = 0;
197 21 Nov 06 enell 136       for (int j = 0; j < bfss.getSpotSize(); j += assays.size())
187 13 Oct 06 enell 137       {
197 21 Nov 06 enell 138         double m = bfss.getSpot(i + j).getM();
197 21 Nov 06 enell 139         if (!Double.isNaN(m))
187 13 Oct 06 enell 140         {
197 21 Nov 06 enell 141           sqrSum += Math.pow(m - means.get(assays.get(i)), 2);
197 21 Nov 06 enell 142           size++;
187 13 Oct 06 enell 143         }
187 13 Oct 06 enell 144       }
197 21 Nov 06 enell 145       sds.put(assays.get(i), Math.sqrt(sqrSum / size));
187 13 Oct 06 enell 146     }
187 13 Oct 06 enell 147
187 13 Oct 06 enell 148     return sds;
187 13 Oct 06 enell 149   }
187 13 Oct 06 enell 150
197 21 Nov 06 enell 151   private HashMap<Integer, Double> getAssayMeans(BASEFileSpotSection<Reporter, Spot> bfss, List<Integer> assays)
187 13 Oct 06 enell 152   {
197 21 Nov 06 enell 153     HashMap<Integer, Double> means = new HashMap<Integer, Double>();
197 21 Nov 06 enell 154     for (int i = 0; i < assays.size(); i++)
187 13 Oct 06 enell 155     {
197 21 Nov 06 enell 156       int size = 0;
197 21 Nov 06 enell 157       double sum = 0;
197 21 Nov 06 enell 158       for (int j = 0; j < bfss.getSpotSize(); j += assays.size())
187 13 Oct 06 enell 159       {
197 21 Nov 06 enell 160         double m = bfss.getSpot(i + j).getM();
197 21 Nov 06 enell 161         if (!Double.isNaN(m))
187 13 Oct 06 enell 162         {
197 21 Nov 06 enell 163           sum += m;
197 21 Nov 06 enell 164           size++;
187 13 Oct 06 enell 165         }
187 13 Oct 06 enell 166       }
197 21 Nov 06 enell 167       means.put(assays.get(i), sum / size);
187 13 Oct 06 enell 168     }
187 13 Oct 06 enell 169     return means;
187 13 Oct 06 enell 170   }
187 13 Oct 06 enell 171 }

Generated on Tue Nov 21 03:41:11 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)