X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=mothurmetastats.cpp;fp=mothurmetastats.cpp;h=31542e80fe9f2d2b01c9d07cd97501954e841bcf;hb=113672689e4c84383a0304bfd21fb9574a77866c;hp=9ba40502a2294d5ce42ba591c076da6fd3015979;hpb=30910bfac62a1299cd20d99683c4fae68bfc2c1b;p=mothur.git diff --git a/mothurmetastats.cpp b/mothurmetastats.cpp index 9ba4050..31542e8 100644 --- a/mothurmetastats.cpp +++ b/mothurmetastats.cpp @@ -204,6 +204,13 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector storage[i][7]=temp[i][1]; storage[i][8]=pvalues[i]; } + cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); + cout << "pvalues" << endl; + for (int i = 0; i < row; i++){ if (pvalues[i] < 0.0000001) { + pvalues[i] = 0.0; + }cout << pvalues[i] << '\t'; } + cout << endl; + calc_qvalues(pvalues); // BACKUP checks cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); @@ -246,6 +253,7 @@ int MothurMetastats::runMetastats(string outputFileName, vector< vector out << endl << endl; out.close(); + return 0; }catch(exception& e) { @@ -484,5 +492,325 @@ int MothurMetastats::calc_twosample_ts(vector& Pmatrix, int secondGroupi } } /***********************************************************/ +vector MothurMetastats::calc_qvalues(vector& pValues) { + try { + vector qvalues; + int numRows = pValues.size(); + + //fill lambdas with 0.00, 0.01, 0.02... 0.95 + vector lambdas(96, 0); + for (int i = 1; i < lambdas.size(); i++) { lambdas[i] = lambdas[i-1] + 0.01; } + + vector pi0_hat(lambdas.size(), 0); + + //calculate pi0_hat + for (int l = 0; l < lambdas.size(); l++){ // for each lambda value + int count = 0; + 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])); + } + } + + //vector f_spline = smoothSpline(lambdas, pi0_hat, 3); + //double pi0 = f_spline[(f_spline.size()-1)]; // this is the essential pi0_hat value + + //cubic Spline + double pi0, notsure; //this is the essential pi0_hat value + int notSure; + vector resultsSpline(lambdas.size(), 0.0); + spline(pi0_hat, lambdas, notSure, notSure, resultsSpline); + //some sort of loop to get best value?? + splint(pi0_hat, lambdas, notsure, pi0, resultsSpline); + + //order p-values + vector ordered_qs = qvalues; + vector ordered_ps(pValues.size(), 0); + for (int i = 1; i < ordered_ps.size(); i++) { ordered_ps[i] = ordered_ps[i-1] + 1; } + vector tempPvalues = pValues; + OrderPValues(0, numRows-1, tempPvalues, ordered_ps); + + ordered_qs[numRows-1] <- min((pValues[ordered_ps[numRows-1]]*pi0), 1.0); + for (int i = (numRows-2); i >= 0; i--){ + double p = pValues[ordered_ps[i]]; + double temp = p*numRows*pi0/(double)i; + + ordered_qs[i] = min(temp, ordered_qs[i+1]); + ordered_qs[i] = min(ordered_qs[i], 1.0); + } + + //re-distribute calculated qvalues to appropriate rows + for (int i = 0; i < numRows; i++){ + qvalues[ordered_ps[i]] = ordered_qs[i]; + } + + return qvalues; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "calc_qvalues"); + exit(1); + } +} +/***********************************************************/ +int MothurMetastats::OrderPValues(int low, int high, vector& p, vector& order) { + try { + + if (low < high) { + int i = low+1; + int j = high; + int pivot = (low+high) / 2; + + swapElements(low, pivot, p, order); //puts pivot in final spot + + /* compare value */ + double key = p[low]; + + /* partition */ + while(i <= j) { + /* find member above ... */ + while((i <= high) && (p[i] <= key)) { i++; } + + /* find element below ... */ + while((j >= low) && (p[j] > key)) { j--; } + + if(i < j) { + swapElements(i, j, p, order); + } + } + + swapElements(low, j, p, order); + + /* recurse */ + OrderPValues(low, j-1, p, order); + OrderPValues(j+1, high, p, order); + } + + return 0; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "OrderPValues"); + exit(1); + } +} +/***********************************************************/ +int MothurMetastats::swapElements(int i, int j, vector& p, vector& order) { + try { + + double z = p[i]; + p[i] = p[j]; + p[j] = z; + + int temp = order[i]; + order[i] = order[j]; + order[j] = temp; + + return 0; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "swapElements"); + exit(1); + } +} +/*********************************************************** +vector MothurMetastats::smoothSpline(vector x, vector y, int df) { + try { + + cout << "lambdas" << endl; + for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; } + cout << endl << "pi0_hat" << endl; + for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; } + cout << endl; + + //double low = -1.5; double high = 1.5; double tol = 1e-04; double eps = 2e-08; double maxit = 500; + int n = x.size(); + vector w(n, 1); + + //x <- signif(x, 6L) - I think this rounds to 6 decimals places + //ux <- unique(sort(x)) //x will be unique and sorted since we created it + //ox <- match(x, ux) //since its unique and sort it will match + + vector wbar(n, 0); + vector ybar(n, 0); + vector xbar(n, 0.0); + double yssw = 0.0; + for (int i = 0; i < n; i++) { + wbar[i] = w[i]; + ybar[i] = w[i]*y[i]; + yssw += (w[i] * y[i] * y[i]) - wbar[i] * (ybar[i] * ybar[i]); + xbar[i] = (x[i] - x[0]) / (x[n-1] - x[0]); + } + + vector knot = sknot1(xbar); + int nk = knot.size() - 4; + + //double ispar = 0.0; double spar = 0.0; double icrit = 3.0; double dofoff = 3.0; + + return y; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "smoothSpline"); + exit(1); + } +} +/***********************************************************/ +// This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479 +int MothurMetastats::spline(vector& x, vector& y, int yp1, int ypn, vector& y2) { + try { + + cout << "lambdas" << endl; + for (int l = 0; l < x.size(); l++){ cout << x[l] << '\t'; } + cout << endl << "pi0_hat" << endl; + for (int l = 0; l < y.size(); l++){ cout << y[l] << '\t'; } + cout << endl; + + double p, qn, sig, un; + + int n = y2.size(); + vector u(n-1, 0.0); + + if (yp1 > 0.99e30) { y2[0] = u[0] = 0.0; } + else { + y2[0] = -0.5; + u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1]-x[0]) - yp1); + } + + for (int i = 1; i < n-1; i++) { + sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]); + p = sig * y2[i-1] + 2.0; + y2[i] = (sig - 1.0)/p; + u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i] - y[i-1]) / (x[i]-x[i-1]); + u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p; + } + + if (ypn > 0.99e30) { qn=un=0.0; } + else { + qn=0.5; + un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); + } + + y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0); + + for (int k=n-2; k>=0; k--) { + y2[k]=y2[k]*y2[k+1]+u[k]; + } + + return 0; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "spline"); + exit(1); + } +} +/***********************************************************/ +// This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479 +int MothurMetastats::splint(vector& xa, vector& ya, double& x, double& y, vector& y2a) { + try { + int k; + double h,b,a; + + int n = xa.size(); + + int klo=0; + int khi=n-1; + while (khi-klo > 1) { + + if (m->control_pressed) { break; } + + k = (khi+klo) >> 1; + if (xa[k] > x) { khi=k; } + else { klo=k; } + } + + h=xa[khi]-xa[klo]; + if (h == 0.0) { m->mothurOut("[ERROR]: Bad xa input to splint routine."); m->mothurOutEndLine(); m->control_pressed = true; return 0; } + a=(xa[khi]-x)/h; + b=(x - xa[klo])/h; + y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; + + return 0; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "splint"); + exit(1); + } +} +/***********************************************************/ +vector MothurMetastats::sknot1(vector x) { + try { + vector knots; + int n = x.size(); + int nk = nkn(n); + + cout << nk << endl; + //R equivalent - rep(x[1L], 3L) + knots.push_back(x[0]); knots.push_back(x[0]); knots.push_back(x[0]); + + //generate a sequence of nk equally spaced values from 1 to n. R equivalent = seq.int(1, n, length.out = nk) + vector indexes = getSequence(0, n-1, nk); + for (int i = 0; i < indexes.size(); i++) { knots.push_back(x[indexes[i]]); } + + //R equivalent - rep(x[n], 3L) + knots.push_back(x[n-1]); knots.push_back(x[n-1]); knots.push_back(x[n-1]); + + return knots; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "sknot1"); + exit(1); + } +} +/***********************************************************/ +vector MothurMetastats::getSequence(int start, int end, int length) { + try { + vector sequence; + double increment = (end-start) / (double) (length-1); + + sequence.push_back(start); + for (int i = 1; i < length-1; i++) { + sequence.push_back(int(i*increment)); + } + sequence.push_back(end); + + return sequence; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "getSequence"); + exit(1); + } +} +/***********************************************************/ +//not right, havent fully figured out the variable types yet... +int MothurMetastats::nkn(int n) { + try { + + if (n < 50) { return n; } + else { + double a1 = log2(50); + double a2 = log2(100); + double a3 = log2(140); + double a4 = log2(200); + + if (n < 200) { + return (int)pow(2.0, (a1 + (a2 - a1) * (n - 50)/(float)150)); + }else if (n < 800) { + return (int)pow(2.0, (a2 + (a3 - a2) * (n - 200)/(float)600)); + }else if (n < 3200) { + return (int)pow(2.0, (a3 + (a4 - a3) * (n - 800)/(float)2400)); + }else { + return (int)pow((double)(200 + (n - 3200)), 0.2); + } + } + + return 0; + + }catch(exception& e) { + m->errorOut(e, "MothurMetastats", "nkn"); + exit(1); + } +} +/***********************************************************/ + + +