5 * Created by westcott on 7/6/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "mothurmetastats.h"
11 #include "mothurfisher.h"
13 /***********************************************************/
14 MothurMetastats::MothurMetastats(double t, int n) {
16 m = MothurOut::getInstance();
20 }catch(exception& e) {
21 m->errorOut(e, "MothurMetastats", "MothurMetastats");
25 /***********************************************************/
26 MothurMetastats::~MothurMetastats() {}
27 /***********************************************************/
28 //main metastats function
29 int MothurMetastats::runMetastats(string outputFileName, vector< vector<double> >& data, int secondGroupingStart) {
32 row = data.size(); //numBins
33 column = data[0].size(); //numGroups in subset
34 int size = row*column;
36 //consistent with original, but this should never be true
37 if ((secondGroupingStart >= column) || (secondGroupingStart <= 0)) { m->mothurOut("[ERROR]: Check your g value."); m->mothurOutEndLine(); return 0; }
39 //Initialize the matrices
40 vector<double> pmatrix; pmatrix.resize(size, 0.0);
41 vector<double> permuted; permuted.resize(size, 0.0);
42 vector< vector<double> > storage; storage.resize(row);
43 for (int i = 0; i < storage.size(); i++) { storage[i].resize(9, 0.0); }
45 //Produces the sum of each column
46 vector<double> total; total.resize(column, 0.0);
47 vector<double> ratio; ratio.resize(column, 0.0);
48 double total1 = 0.0; double total2 = 0.0;
50 //total[i] = total abundance for group[i]
51 for (int i = 0; i < column; i++) {
52 for (int j = 0; j < row; j++) {
53 total[i] += data[j][i];
57 //total for first grouping
58 for (int i = 0; i < secondGroupingStart; i++) { total1 += total[i]; }
60 //total for second grouping
61 for (int i = secondGroupingStart; i < column; i++) { total2 += total[i]; }
63 //Creates the ratios by first finding the minimum of totals
64 double min = total[0];
65 for (int i = 0; i < total.size(); i++) {
66 if (total[i] < min) { min = total[i]; }
70 if (min <= 0.0) { m->mothurOut("[ERROR]: the sum of one of the columns <= 0."); m->mothurOutEndLine(); return 0; }
73 for(int i = 0; i < ratio.size(); i++){ ratio[i] = total[i] / min; }
75 //Change matrix into an array as received by R for compatibility - kept to be consistent with original
77 for(int i = 0; i < column; i++){
78 for(int j = 0; j < row; j++){
79 pmatrix[count]=data[j][i];
85 for (int i =0; i < column; i++){ pmatrix[i] /= ratio[i]; }
89 for (int i=0; i < size; i++) {
90 if (count % row == 0) { j++; }
91 pmatrix[i] /= ratio[j];
96 vector<double> permuted_ttests; permuted_ttests.resize(row, 0.0);
97 vector<double> pvalues; pvalues.resize(row, 0.0);
98 vector<double> tinitial; tinitial.resize(row, 0.0);
100 if (m->control_pressed) { return 1; }
102 //Find the initial values for the matrix.
103 start(pmatrix, secondGroupingStart, tinitial, storage);
105 if (m->control_pressed) { return 1; }
107 // Start the calculations.
108 if ( (column == 2) || (secondGroupingStart < 8) || ((column-secondGroupingStart) < 8) ){
110 vector<double> fish; fish.resize(row, 0.0);
111 vector<double> fish2; fish2.resize(row, 0.0);
113 for(int i = 0; i < row; i++){
115 for(int j = 0; j < secondGroupingStart; j++) { fish[i] += data[i][j]; }
116 for(int j = secondGroupingStart; j < column; j++) { fish2[i] += data[i][j]; }
118 //vector<double> tempData; tempData.resize(4, 0.0);
119 double f11, f12, f21, f22;
122 f21 = total1 - fish[i];
123 f22 = total2 - fish2[i];
128 pre = fisher.fexact(f11, f12, f21, f22);
130 if (m->control_pressed) { return 1; }
132 if (pre > 0.999999999) { pre = 1.0; }
139 testp(permuted_ttests, permuted, pmatrix, secondGroupingStart, tinitial, pvalues);
141 if (m->control_pressed) { return 1; }
143 // Checks to make sure the matrix isn't sparse.
144 vector<double> sparse; sparse.resize(row, 0.0);
145 vector<double> sparse2; sparse2.resize(row, 0.0);
149 for(int i = 0; i < row; i++){
151 for(int j = 0; j < secondGroupingStart; j++) { sparse[i] += data[i][j]; }
152 if(sparse[i] < (double)secondGroupingStart) { c++; }
155 for(int j = secondGroupingStart; j < column; j++) { sparse2[i] += data[i][j]; }
156 if( (sparse2[i] < (double)(column-secondGroupingStart))) { c++; }
160 double f11,f12,f21,f22;
162 f11=sparse[i]; sparse[i]=0;
163 f12=sparse2[i]; sparse2[i]=0;
170 pre = fisher.fexact(f11, f12, f21, f22);
172 if (m->control_pressed) { return 1; }
174 if (pre > 0.999999999){
186 // Calculates the mean of counts (not normalized)
187 vector< vector<double> > temp; temp.resize(row);
188 for (int i = 0; i < temp.size(); i++) { temp[i].resize(2, 0.0); }
190 for (int j = 0; j < row; j++){
191 if (m->control_pressed) { return 1; }
193 for (int i = 0; i < secondGroupingStart; i++){ temp[j][0] += data[j][i]; }
194 temp[j][0] /= (double)secondGroupingStart;
196 for(int i = secondGroupingStart; i < column; i++){ temp[j][1] += data[j][i]; }
197 temp[j][1] /= (double)(column-secondGroupingStart);
200 for(int i = 0; i < row; i++){
201 if (m->control_pressed) { return 1; }
203 storage[i][3]=temp[i][0];
204 storage[i][7]=temp[i][1];
205 storage[i][8]=pvalues[i];
209 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
210 for (int i = 0; i < row; i++){
212 if (m->control_pressed) { return 1; }
214 if(pvalues[i] < threshold){
215 m->mothurOut("Feature " + toString((i+1)) + " is significant, p = ");
217 m->mothurOutJustToLog(toString(pvalues[i])); m->mothurOutEndLine();
221 // And now we write the files to a text file.
223 time_t t; t = time(NULL);
224 local = localtime(&t);
227 m->openOutputFile(outputFileName, out);
228 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
230 out << "Local time and date of test: " << asctime(local) << endl;
231 out << "# rows = " << row << ", # col = " << column << ", g = " << secondGroupingStart << endl << endl;
232 if (bflag == 1){ out << numPermutations << " permutations" << endl << endl; }
234 //output column headings - not really sure... documentation labels 9 columns, there are 10 in the output file
235 //storage 0 = meanGroup1 - line 529, 1 = varGroup1 - line 532, 2 = err rate1 - line 534, 3 = mean of counts group1?? - line 291, 4 = meanGroup2 - line 536, 5 = varGroup2 - line 539, 6 = err rate2 - line 541, 7 = mean of counts group2?? - line 292, 8 = pvalues - line 293
236 out << "OTU\tmean(group1)\tvariance(group1)\tstderr(group1)\tmean_of_counts(group1)\tmean(group2)\tvariance(group2)\tstderr(group2)\tmean_of_counts(group1)\tp-value\n";
238 for(int i = 0; i < row; i++){
239 if (m->control_pressed) { out.close(); return 0; }
242 for(int j = 0; j < 9; j++){ out << '\t' << storage[i][j]; }
251 }catch(exception& e) {
252 m->errorOut(e, "MothurMetastats", "runMetastats");
256 /***********************************************************/
257 //Find the initial values for the matrix
258 int MothurMetastats::start(vector<double>& Imatrix, int secondGroupingStart, vector<double>& initial, vector< vector<double> >& storage) {
263 double xbardiff = 0.0; double denom = 0.0;
264 vector<double> store; store.resize(a, 0.0);
265 vector<double> tool; tool.resize(a, 0.0);
266 vector< vector<double> > C1; C1.resize(row);
267 for (int i = 0; i < C1.size(); i++) { C1[i].resize(3, 0.0); }
268 vector< vector<double> > C2; C2.resize(row);
269 for (int i = 0; i < C2.size(); i++) { C2[i].resize(3, 0.0); }
271 meanvar(Imatrix, secondGroupingStart, store);
273 if (m->control_pressed) { return 0; }
275 //copy store into tool
278 for (int i = 0; i < row; i++){
279 C1[i][0]=tool[i]; //mean group 1
280 storage[i][0]=C1[i][0];
281 C1[i][1]=tool[i+row+row]; // var group 1
282 storage[i][1]=C1[i][1];
283 C1[i][2]=C1[i][1]/(secondGroupingStart);
284 storage[i][2]=sqrt(C1[i][2]);
286 C2[i][0]=tool[i+row]; // mean group 2
287 storage[i][4]=C2[i][0];
288 C2[i][1]=tool[i+row+row+row]; // var group 2
289 storage[i][5]=C2[i][1];
290 C2[i][2]=C2[i][1]/(column-secondGroupingStart);
291 storage[i][6]=sqrt(C2[i][2]);
294 if (m->control_pressed) { return 0; }
296 for (int i = 0; i < row; i++){
297 xbardiff = C1[i][0]-C2[i][0];
298 denom = sqrt(C1[i][2]+C2[i][2]);
299 initial[i]=fabs(xbardiff/denom);
304 }catch(exception& e) {
305 m->errorOut(e, "MothurMetastats", "start");
309 /***********************************************************/
310 int MothurMetastats::meanvar(vector<double>& pmatrix, int secondGroupingStart, vector<double>& store) {
312 vector<double> temp; temp.resize(row, 0.0);
313 vector<double> temp2; temp2.resize(row, 0.0);
314 vector<double> var; var.resize(row, 0.0);
315 vector<double> var2; var2.resize(row, 0.0);
317 double a = secondGroupingStart;
318 double b = column - a;
320 int n = row * column;
322 for (int i = 0; i < m; i++) { temp[i%row] += pmatrix[i]; }
323 for (int i = 0; i < n; i++) { temp2[i%row]+= pmatrix[i]; }
324 for (int i = 0; i < row; i++) { temp2[i] -= temp[i]; }
325 for (int i = 0; i <= row-1;i++) {
326 store[i] = temp[i]/a;
327 store[i+row]=temp2[i]/b;
330 //That completes the mean calculations.
332 for (int i = 0; i < m; i++) { var[i%row] += pow((pmatrix[i]-store[i%row]),2); }
333 for (int i = m; i < n; i++) { var2[i%row]+= pow((pmatrix[i]-store[(i%row)+row]),2); }
334 for (int i = 0; i <= row-1; i++){
335 store[i+2*row]=var[i]/(a-1);
336 store[i+3*row]=var2[i]/(b-1);
339 // That completes var calculations.
343 }catch(exception& e) {
344 m->errorOut(e, "MothurMetastats", "meanvar");
348 /***********************************************************/
349 int MothurMetastats::testp(vector<double>& permuted_ttests, vector<double>& permuted, vector<double>& Imatrix, int secondGroupingStart, vector<double>& Tinitial, vector<double>& ps) {
352 vector<double> Tvalues; Tvalues.resize(row, 0.0);
353 vector<double> counter; counter.resize(row, 0.0);
360 for (int j = 1; j <= row; j++) {
361 if (m->control_pressed) { return 0; }
362 permute_matrix(Imatrix, permuted, secondGroupingStart, Tvalues, Tinitial, counter);
365 for(int j = 0; j < row; j++) {
366 if (m->control_pressed) { return 0; }
367 ps[j] = ((counter[j]+1)/(double)(a+1));
372 }catch(exception& e) {
373 m->errorOut(e, "MothurMetastats", "testp");
377 /***********************************************************/
378 int MothurMetastats::permute_matrix(vector<double>& Imatrix, vector<double>& permuted, int secondGroupingStart, vector<double>& trial_ts, vector<double>& Tinitial, vector<double>& counter1){
381 vector<int> y; y.resize(column, 0);
382 for (int i = 1; i <= column; i++){ y[i-1] = i; }
386 int f = 0; int c = 0; int k = 0;
387 for (int i = 0; i < column; i++){
389 if (m->control_pressed) { return 0; }
391 f = y[i]; //column number
395 if (f == 1){ c = 0; } // starting value position in the Imatrix
397 for(int j = 1; j <= row; j++){
398 permuted[k] = Imatrix[c];
403 calc_twosample_ts(permuted, secondGroupingStart, trial_ts, Tinitial, counter1);
407 }catch(exception& e) {
408 m->errorOut(e, "MothurMetastats", "permute_matrix");
412 /***********************************************************/
413 int MothurMetastats::permute_array(vector<int>& array) {
415 static int seeded = 0;
422 for (int i = 0; i < array.size(); i++) {
423 if (m->control_pressed) { return 0; }
425 int selection = rand() % (array.size() - i);
426 int tmp = array[i + selection];
427 array[i + selection] = array[i];
433 }catch(exception& e) {
434 m->errorOut(e, "MothurMetastats", "permute_array");
438 /***********************************************************/
439 int MothurMetastats::calc_twosample_ts(vector<double>& Pmatrix, int secondGroupingStart, vector<double>& Ts, vector<double>& Tinitial, vector<double>& counter) {
443 vector< vector<double> > C1; C1.resize(row);
444 for (int i = 0; i < C1.size(); i++) { C1[i].resize(3, 0.0); }
445 vector< vector<double> > C2; C2.resize(row);
446 for (int i = 0; i < C2.size(); i++) { C2[i].resize(3, 0.0); }
447 vector<double> storage; storage.resize(a, 0.0);
448 vector<double> tool; tool.resize(a, 0.0);
449 double xbardiff = 0.0; double denom = 0.0;
451 meanvar(Pmatrix, secondGroupingStart, storage);
453 for(int i = 0;i <= (a-1); i++) {
454 if (m->control_pressed) { return 0; }
455 tool[i] = storage[i];
458 for (int i = 0; i < row; i++){
459 if (m->control_pressed) { return 0; }
461 C1[i][1]=tool[i+row+row];
462 C1[i][2]=C1[i][1]/(secondGroupingStart);
464 C2[i][0]=tool[i+row];
465 C2[i][1]=tool[i+row+row+row]; // var group 2
466 C2[i][2]=C2[i][1]/(column-secondGroupingStart);
469 for (int i = 0; i < row; i++){
470 if (m->control_pressed) { return 0; }
471 xbardiff = C1[i][0]-C2[i][0];
472 denom = sqrt(C1[i][2]+C2[i][2]);
473 Ts[i]=fabs(xbardiff/denom);
474 if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
481 }catch(exception& e) {
482 m->errorOut(e, "MothurMetastats", "calc_twosample_ts");
486 /***********************************************************/