X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=mothurmetastats.cpp;h=f6a7800c3697833eeb9160330671804f13633bfb;hb=250e3b11b1c9c1e1ad458ab6c7e71ac2e67e11d9;hp=56789736e5b0d8d9bbf273baa4a07c36efb86ec8;hpb=91a27e0483827c06c21c4fe89558923bbfe86573;p=mothur.git diff --git a/mothurmetastats.cpp b/mothurmetastats.cpp index 5678973..f6a7800 100644 --- a/mothurmetastats.cpp +++ b/mothurmetastats.cpp @@ -31,7 +31,7 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector try { row = data.size(); //numBins column = data[0].size(); //numGroups in subset - secondGroupingStart = secGroupingStart; //g + secondGroupingStart = secGroupingStart; //g number of samples in group 1 vector< vector > Pmatrix; Pmatrix.resize(row); for (int i = 0; i < row; i++) { Pmatrix[i].resize(column, 0.0); } // the relative proportion matrix @@ -58,7 +58,6 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector for (int i = 0; i < column; i++) { for (int j = 0; j < row; j++) { Pmatrix[j][i] = data[j][i]/totals[i]; - } } @@ -70,7 +69,7 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector //#************************************************************ //# generate p values fisher's exact test //#************************************************************ - double total1, total2; + double total1, total2; total1 = 0; total2 = 0; //total for first grouping for (int i = 0; i < secondGroupingStart; i++) { total1 += totals[i]; } @@ -139,13 +138,20 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector T_statistics[i] = xbar_diff/denom; // calculate two sample t-statistic } - /*for (int i = 0; i < row; i++) { + /*for (int i = 0; i < row; i++) { for (int j = 0; j < 3; j++) { cout << "C1[" << i+1 << "," << j+1 << "]=" << C1[i][j] << ";" << endl; cout << "C2[" << i+1 << "," << j+1 << "]=" << C2[i][j] << ";" << endl; } cout << "T_statistics[" << i+1 << "]=" << T_statistics[i] << ";" << endl; + } + + for (int i = 0; i < row; i++) { + for (int j = 0; j < column; j++) { + cout << "Fmatrix[" << i+1 << "," << j+1 << "]=" << data[i][j] << ";" << endl; + } }*/ + //#************************************* //# generate initial permuted p-values //#************************************* @@ -155,12 +161,12 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector //# generate p values for sparse data //# using fisher's exact test //#************************************* - double total1, total2; + double total1, total2; total1 = 0; total2 = 0; //total for first grouping - for (int i = 0; i < secondGroupingStart; i++) { total1 += totals[i]; } + for (int i = 0; i < secondGroupingStart; i++) { total1 += totals[i]; } //total for second grouping - for (int i = secondGroupingStart; i < column; i++) { total2 += totals[i]; } + for (int i = secondGroupingStart; i < column; i++) { total2 += totals[i]; } vector fish; fish.resize(row, 0.0); vector fish2; fish2.resize(row, 0.0); @@ -169,8 +175,9 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector for(int j = 0; j < secondGroupingStart; j++) { fish[i] += data[i][j]; } for(int j = secondGroupingStart; j < column; j++) { fish2[i] += data[i][j]; } - - if ((fish[1] < secondGroupingStart) && (fish2[i] < (column-secondGroupingStart))) { + + if ((fish[i] < secondGroupingStart) && (fish2[i] < (column-secondGroupingStart))) { + double f11, f12, f21, f22; f11 = fish[i]; f12 = fish2[i]; @@ -222,7 +229,7 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector if (m->control_pressed) { out.close(); return 0; } //if there are binlabels use them otherwise count. - if (m->binLabelsInFile.size() == row) { out << m->binLabelsInFile[i] << '\t'; } + if (i < m->currentSharedBinLabels.size()) { out << m->currentSharedBinLabels[i] << '\t'; } else { out << (i+1) << '\t'; } out << C1[i][0] << '\t' << C1[i][1] << '\t' << C1[i][2] << '\t' << C2[i][0] << '\t' << C2[i][1] << '\t' << C2[i][2] << '\t' << pvalues[i] << '\t' << qvalues[i] << endl; @@ -269,7 +276,7 @@ vector MothurMetastats::permuted_pvalues(vector< vector >& Imatr for (int j = 0; j < numPermutations; j++) { cleanedpermuted_ttests[j].push_back(permuted_ttests[j][i]); } } } - + //#now for each taxa for (int i = 0; i < row; i++) { //number of cleanedpermuted_ttests greater than tstat[i] @@ -304,7 +311,11 @@ vector MothurMetastats::permute_and_calc_ts(vector< vector >& Im vector< vector > permutedMatrix = Imatrix; //randomize columns, ie group abundances. - for (int i = 0; i < permutedMatrix.size(); i++) { random_shuffle(permutedMatrix[i].begin(), permutedMatrix[i].end()); } + map randomMap; + vector randoms; + for (int i = 0; i < column; i++) { randoms.push_back(i); } + random_shuffle(randoms.begin(), randoms.end()); + for (int i = 0; i < randoms.size(); i++) { randomMap[i] = randoms[i]; } //calc ts vector< vector > C1; C1.resize(row); @@ -319,16 +330,16 @@ vector MothurMetastats::permute_and_calc_ts(vector< vector >& Im for(int i = 0; i < row; i++){ // for each taxa //# find the mean of each group double g1Total = 0.0; double g2Total = 0.0; - for (int j = 0; j < secondGroupingStart; j++) { g1Total += permutedMatrix[i][j]; } + for (int j = 0; j < secondGroupingStart; j++) { g1Total += permutedMatrix[i][randomMap[j]]; } C1[i][0] = g1Total/(double)(secondGroupingStart); - for (int j = secondGroupingStart; j < column; j++) { g2Total += permutedMatrix[i][j]; } + for (int j = secondGroupingStart; j < column; j++) { g2Total += permutedMatrix[i][randomMap[j]]; } C2[i][0] = g2Total/(double)(column-secondGroupingStart); //# find the variance of each group double g1Var = 0.0; double g2Var = 0.0; - for (int j = 0; j < secondGroupingStart; j++) { g1Var += pow((permutedMatrix[i][j]-C1[i][0]), 2); } + for (int j = 0; j < secondGroupingStart; j++) { g1Var += pow((permutedMatrix[i][randomMap[j]]-C1[i][0]), 2); } C1[i][1] = g1Var/(double)(secondGroupingStart-1); - for (int j = secondGroupingStart; j < column; j++) { g2Var += pow((permutedMatrix[i][j]-C2[i][0]), 2); } + for (int j = secondGroupingStart; j < column; j++) { g2Var += pow((permutedMatrix[i][randomMap[j]]-C2[i][0]), 2); } C2[i][1] = g2Var/(double)(column-secondGroupingStart-1); //# find the std error of each group -std err^2 (will change to std err at end) @@ -336,6 +347,8 @@ vector MothurMetastats::permute_and_calc_ts(vector< vector >& Im C2[i][2] = C2[i][1]/(double)(column-secondGroupingStart); } + + //#************************************* //# two sample t-statistics //#************************************* @@ -356,14 +369,7 @@ vector MothurMetastats::permute_and_calc_ts(vector< vector >& Im /***********************************************************/ vector MothurMetastats::calc_qvalues(vector& pValues) { try { - - /* cout << "x <- c(" << pValues[0]; - for (int l = 1; l < pValues.size(); l++){ - cout << ", " << pValues[l]; - } - cout << ")\n";*/ - - int numRows = pValues.size(); + int numRows = pValues.size(); vector qvalues(numRows, 0.0); //fill lambdas with 0.00, 0.01, 0.02... 0.95 @@ -378,7 +384,8 @@ vector MothurMetastats::calc_qvalues(vector& pValues) { for (int i = 0; i < numRows; i++){ // for each p-value in order if (pValues[i] > lambdas[l]){ count++; } } - pi0_hat[l] = count/(double)(numRows*(1-lambdas[l])); + pi0_hat[l] = count/(double)(numRows*(1.0-lambdas[l])); + //cout << lambdas[l] << '\t' << count << '\t' << numRows*(1.0-lambdas[l]) << '\t' << pi0_hat[l] << endl; } double pi0 = smoothSpline(lambdas, pi0_hat, 3);