View Javadoc
1   /*
2    * The baseCode project
3    *
4    * Copyright (c) 2011 University of British Columbia
5    *
6    * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with
7    * the License. You may obtain a copy of the License at
8    *
9    * http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on
12   * an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the
13   * specific language governing permissions and limitations under the License.
14   */
15  package ubic.basecode.math.linearmodels;
16  
17  import org.junit.Test;
18  
19  import ubic.basecode.dataStructure.matrix.DenseDoubleMatrix;
20  import ubic.basecode.dataStructure.matrix.DoubleMatrix;
21  import ubic.basecode.dataStructure.matrix.StringMatrix;
22  import ubic.basecode.io.reader.DoubleMatrixReader;
23  import ubic.basecode.io.reader.StringMatrixReader;
24  import ubic.basecode.math.MatrixStats;
25  import ubic.basecode.math.Smooth;
26  import ubic.basecode.math.linearmodels.DesignMatrix;
27  import cern.colt.matrix.DoubleMatrix1D;
28  import cern.colt.matrix.DoubleMatrix2D;
29  import cern.colt.matrix.impl.DenseDoubleMatrix2D;
30  import ubic.basecode.util.RegressionTesting;
31  
32  import java.io.File;
33  import java.io.OutputStream;
34  import java.io.PrintStream;
35  
36  import static org.junit.Assert.*;
37  
38  /**
39   * @author ptan
40   */
41  public class MeanVarianceEstimatorTest {
42  
43      /**
44       * Tests two things: 1. Input data's X values are not sorted 2. Interpolated X values are outside the domain.
45       *
46       * @throws Exception R code:
47       *                   x = c(9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0)
48       *                   y = c(1.0, 0.25, 0.1111111111111111, 0.0625, 0.04, 0.027777777777777776, 0.02040816326530612, 0.015625, 0.012345679012345678, 0.01)
49       *                   interpolate = c(9.5, 8.5, 7.5, 6.5, 5.5, 4.5, 3.5, 2.5, 1.5, 0.5)
50       *                   f<-approxfun(x,y,rule=2)
51       *                   f(c(interpolate))
52       */
53      @Test
54      public void testInterpolate() throws Exception {
55          double x[] = new double[10];
56          double y[] = new double[x.length];
57          double xInterpolate[] = new double[x.length];
58          for (int i = 0; i < x.length; i++) {
59              x[i] = 10 - (i + 1);
60              y[i] = 1.0 / Math.pow((i + 1), 2);
61              xInterpolate[i] = 10 - (i + 1) + 0.5;
62          }
63  
64          double expected[] = new double[]{1.000000000000000000000, 0.625000000000000000000, 0.180555555555555552472,
65                  0.086805555555555552472, 0.051250000000000003886
66                  , 0.033888888888888885065, 0.024092970521541946793, 0.018016581632653058676, 0.013985339506172839164, 0.011172839506172840135};
67  
68          double yOut[] = Smooth.interpolate(x, y, xInterpolate);
69          for (int i = 0; i < expected.length; i++) {
70              assertEquals(expected[i], yOut[i], 1e-8);
71          }
72      }
73  
74      /**
75       * Duplicate row
76       *
77       * @throws Exception
78       */
79      @Test
80      public void testDuplicateRows() throws Exception {
81          DoubleMatrixReader f = new DoubleMatrixReader();
82          DoubleMatrix<String, String> testMatrix = f.read(this.getClass()
83                  .getResourceAsStream("/data/lmtest11.dat.txt"));
84  
85          StringMatrixReader of = new StringMatrixReader();
86          StringMatrix<String, String> sampleInfo = of.read(this.getClass().getResourceAsStream(
87                  "/data/lmtest11.des.txt"));
88  
89          // add a duplicate entry with a different row id but similar expression levels
90          DoubleMatrix<String, String> duplMatrix = new DenseDoubleMatrix<>(testMatrix.asDoubles());
91          duplMatrix.viewRow(1).assign(duplMatrix.viewRow(0));
92          DoubleMatrix1D libSize = MatrixStats.colSums(duplMatrix);
93          duplMatrix = MatrixStats.convertToLog2Cpm(duplMatrix, libSize);
94  
95          DesignMatrix d = new DesignMatrix(sampleInfo, true);
96          MeanVarianceEstimator est = new MeanVarianceEstimator(d, duplMatrix, libSize);
97          DoubleMatrix2D actuals = est.getWeights();
98  
99          int[] expectedIndices = new int[]{0, 1, 149};
100         double[][] expected = new double[][]{
101                 {0.755036, 0.812033, 0.938803, 0.973691, 0.970172, 0.968853, 0.755036},
102                 {0.755036, 0.812033, 0.938803, 0.973691, 0.970172, 0.968853, 0.755036},
103                 {29.462294, 32.386751, 36.690275, 37.844764, 40.823251, 40.776166, 23.402288}};
104         for (int i = 0; i < expectedIndices.length; i++) {
105             assertArrayEquals(expected[i], actuals.viewRow(expectedIndices[i]).toArray(), 0.1);
106         }
107     }
108 
109     /**
110      * Test calculation of weights for LeastSquaresFitting. See test-squeezevar.R
111      *
112      * @throws Exception
113      */
114     @Test
115     public void testVoom() throws Exception {
116         DoubleMatrixReader reader = new DoubleMatrixReader();
117         DoubleMatrix<String, String> testMatrix = reader.read(this.getClass()
118                 .getResourceAsStream("/data/lmtest11.dat.txt"));
119         DoubleMatrix1D libSize = MatrixStats.colSums(testMatrix);
120 
121         assertArrayEquals(new double[]{6792, 8101, 10133, 10672, 11346, 11329, 4005}, libSize.toArray(), 0.000001);
122 
123         testMatrix = MatrixStats.convertToLog2Cpm(testMatrix, libSize);
124         DoubleMatrix expectedlog2cpm = reader.read(this.getClass()
125                 .getResourceAsStream("/data/lmtest11.log2cpm.txt"));
126         assertTrue(RegressionTesting.closeEnough(expectedlog2cpm, testMatrix, 1e-4));
127 
128         StringMatrixReader of = new StringMatrixReader();
129         StringMatrix<String, String> sampleInfo = of.read(this.getClass().getResourceAsStream(
130                 "/data/lmtest11.des.txt"));
131 
132         DesignMatrix d = new DesignMatrix(sampleInfo, true);
133         MeanVarianceEstimator est = new MeanVarianceEstimator(d, testMatrix, libSize);
134 
135         DoubleMatrix2D actuals = est.getWeights();
136         // Our lowess result is quite close
137         DoubleMatrix2D expectedLowess = new DenseDoubleMatrix2D(reader.read(this.getClass().getResourceAsStream(
138                 "/data/lmtest11.lowess.txt")).asArray());
139         DoubleMatrix2D actualLoess = est.getLoess();
140         assertTrue(RegressionTesting.closeEnough(expectedLowess, actualLoess, 0.01));
141 
142         DoubleMatrix2D expectedWeights = new DenseDoubleMatrix2D(reader.read(this.getClass()
143                 .getResourceAsStream("/data/lmtest11.voomweights.txt")).asArray());
144         // The tolerance for this test is low because our weights are slightly different from the limma:voom ones probably because of the interpolation.
145         assertTrue(RegressionTesting.closeEnough(expectedWeights, actuals, 0.12));
146 
147     }
148 
149     /**
150      * Data has missing values, no Design matrix provided so plot a generic mean-variance plot
151      *
152      * @throws Exception
153      */
154     @Test
155     public void testMeanVarianceNoDesignWithColsRowAllMissing() throws Exception {
156 
157         double[][] dataPrimitive = new double[][]{{Double.NaN, -0.2, -0.3, -4, -5},
158                 {Double.NaN, 0, 7, Double.NaN, 8}, {Double.NaN, 9, -0.3, 0.5, Double.NaN},
159                 {Double.NaN, 3, -0.3, -0.1, Double.NaN},
160                 {Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN}};
161         DoubleMatrix<String, String> data = new DenseDoubleMatrix<>(dataPrimitive);
162         DoubleMatrix1D libSize = MatrixStats.colSums(data);
163         data = MatrixStats.convertToLog2Cpm(data, libSize);
164         MeanVarianceEstimator est = new MeanVarianceEstimator(new DenseDoubleMatrix2D(data.asArray()));
165         DoubleMatrix2D actuals = null;
166 
167         actuals = est.getMeanVariance();
168         assertEquals(5, actuals.rows());
169         assertEquals(16.55, actuals.get(0, 0), 0.01);
170         assertEquals(7.260, actuals.get(0, 1), 0.001);
171         assertEquals(16.42, actuals.get(3, 0), 0.01);
172         assertEquals(2.688, actuals.get(3, 1), 0.001);
173         assertEquals(Double.NaN, actuals.get(4, 0), 0.01);
174         assertEquals(Double.NaN, actuals.get(4, 1), 0.001);
175 
176         actuals = est.getLoess();
177         assertEquals(4, actuals.rows());
178         assertEquals(16.42, actuals.get(0, 0), 0.01);
179         assertEquals(2.688, actuals.get(0, 1), 0.001);
180         assertEquals(18.76, actuals.get(3, 0), 0.01);
181         assertEquals(6.321, actuals.get(3, 1), 0.001);
182     }
183 
184     /**
185      * Data has missing values, no Design matrix provided so plot a generic mean-variance plot (no voom)
186      *
187      * @throws Exception
188      */
189     @Test
190     public void testMeanVarianceNoDesignWithMissing() throws Exception {
191         DoubleMatrixReader f = new DoubleMatrixReader();
192         DoubleMatrix<String, String> testMatrix = f.read(this.getClass().getResourceAsStream(
193                 "/data/example.madata.withmissing.small.txt"));
194         DoubleMatrix<String, String> data = new DenseDoubleMatrix<>(testMatrix.asDoubles());
195         DoubleMatrix1D libSize = MatrixStats.colSums(data);
196         data = MatrixStats.convertToLog2Cpm(data, libSize);
197         MeanVarianceEstimator est = new MeanVarianceEstimator(new DenseDoubleMatrix2D(data.asArray()));
198         DoubleMatrix2D actuals = null;
199 
200         actuals = est.getMeanVariance();
201         assertEquals(15.33, actuals.get(0, 0), 0.01);
202         assertEquals(0.003220, actuals.get(0, 1), 0.001);
203         assertEquals(15.28, actuals.get(4, 0), 0.01);
204         assertEquals(0.001972, actuals.get(4, 1), 0.001);
205         assertEquals(15.97, actuals.get(18, 0), 0.01);
206         assertEquals(0.010072, actuals.get(18, 1), 0.001);
207 
208         actuals = est.getLoess();
209         assertEquals(15.24, actuals.get(0, 0), 0.01);
210         assertEquals(0.002620, actuals.get(0, 1), 0.001);
211         assertEquals(15.37, actuals.get(4, 0), 0.01);
212         assertEquals(0.006238, actuals.get(4, 1), 0.001);
213         assertEquals(16.72, actuals.get(18, 0), 0.01);
214         assertEquals(0.007320, actuals.get(18, 1), 0.001);
215 
216         actuals = est.getWeights();
217         assertNull(actuals);
218     }
219 
220 }