//**********************************************************************************************************************
string CooccurrenceCommand::getHelpString(){
try {
- string helpString = "help!";
-
+ string helpString = "The cooccurrence command calculates four metrics and tests their significance to assess whether presence-absence patterns are different than what one would expect by chance.";
+ helpString += "The cooccurrence command parameters are shared, metric, matrixmodel, iters, label and groups.";
+ helpString += "The matrixmodel parameter options are sim1, sim2, sim3, sim4, sim5, sim6, sim7, sim8 and sim9. Default=sim2";
+ helpString += "The metric parameter options are cscore, checker, combo and vratio. Default=cscore";
+ helpString += "The label parameter is used to analyze specific labels in your input.\n";
+ helpString += "The groups parameter allows you to specify which of the groups you would like analyzed.\n";
+ helpString += "The cooccurrence command should be in the following format: \n";
+ helpString += "cooccurrence(shared=yourSharedFile) \n";
+ helpString += "Example cooccurrence(shared=final.an.shared).\n";
+ helpString += "Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourShared).\n";
return helpString;
}
catch(exception& e) {
m->mothurOut("[ERROR]: " + metric + " is not a valid metric option for the cooccurrence command. Choices are cscore, checker, combo, vratio."); m->mothurOutEndLine(); abort = true;
}
- matrix = validParameter.validFile(parameters, "matrix", false); if (matrix == "not found") { matrix = "sim2"; }
+ matrix = validParameter.validFile(parameters, "matrixmodel", false); if (matrix == "not found") { matrix = "sim2"; }
if ((matrix != "sim1") && (matrix != "sim2") && (matrix != "sim3") && (matrix != "sim4") && (matrix != "sim5" ) && (matrix != "sim6" ) && (matrix != "sim7" ) && (matrix != "sim8" ) && (matrix != "sim9" )) {
m->mothurOut("[ERROR]: " + matrix + " is not a valid matrix option for the cooccurrence command. Choices are sim1, sim2, sim3, sim4, sim5, sim6, sim7, sim8, sim9."); m->mothurOutEndLine(); abort = true;
//reset groups parameter
delete input;
+ m->clearGroups();
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
int CooccurrenceCommand::getCooccurrence(vector<SharedRAbundVector*>& thisLookUp, ofstream& out){
try {
int numOTUS = thisLookUp[0]->getNumBins();
- vector< vector<int> > initmatrix (thisLookUp.size());
- vector< vector<int> > co_matrix (thisLookUp[0]->getNumBins());
+ vector< vector<int> > initmatrix; initmatrix.resize(thisLookUp.size());
+ vector< vector<int> > co_matrix; co_matrix.resize(thisLookUp[0]->getNumBins());
for (int i = 0; i < thisLookUp[0]->getNumBins(); i++) { co_matrix[i].resize((thisLookUp.size()), 0); }
for (int i = 0; i < thisLookUp.size(); i++) { initmatrix[i].resize((thisLookUp[i]->getNumBins()), 0); }
- vector<int> columntotal(thisLookUp.size(), 0);
- vector<int> rowtotal(numOTUS, 0);
+ vector<int> columntotal; columntotal.resize(thisLookUp.size(), 0);
+ vector<int> rowtotal; rowtotal.resize(numOTUS, 0);
- int rowcount = 0;
for (int i = 0; i < thisLookUp.size(); i++) {
for (int j = 0; j < thisLookUp[i]->getNumBins(); j++) {
if (m->control_pressed) { return 0; }
if(abund > 0) {
initmatrix[i][j] = 1;
co_matrix[j][i] = 1;
- rowcount++;
- columntotal[j]++;
+ rowtotal[j]++;
+ columntotal[i]++;
}
}
- rowtotal[i] = rowcount;
- rowcount = 0;
}
//nrows is ncols of inital matrix. All the functions need this value. They assume the transposition has already taken place and nrows and ncols refer to that matrix.
//comatrix and initmatrix are still vectors of vectors of ints as in the original script. The abundancevector is only what was read in ie not a co-occurrence matrix!
- int ncols = numOTUS;//rows of inital matrix
- int nrows = thisLookUp.size();//OTUs
+ int nrows = numOTUS;//rows of inital matrix
+ int ncols = thisLookUp.size();//groups
double initscore = 0.0;
- //transpose matrix
- int newmatrows = ncols;
- int newmatcols = nrows;
-
- //swap for transposed matrix
- nrows = newmatrows;//ncols;
- ncols = newmatcols;//nrows;
-
- vector<int> initcolumntotal(ncols, 0);
- vector<int> initrowtotal(nrows, 0);
+
+ vector<int> columntotal; columntotal.resize(ncols, 0);
+ vector<int> rowtotal; rowtotal.resize(nrows, 0);
vector<double> stats;
+ double probabilityMatrix[ncols * nrows];
+ vector<vector<int> > nullmatrix(nrows, vector<int>(ncols, 0));
TrialSwap2 trial;
- initcolumntotal = rowtotal;
- initrowtotal = columntotal;
- trial.update_row_col_totals(co_matrix, rowtotal, columntotal);
+ n = accumulate( columntotal.begin(), columntotal.end(), 0 );
+
+ //============================================================
- if (metric == "cscore") { initscore = trial.calc_c_score(co_matrix, rowtotal); }
- else if (metric == "checker") { initscore = trial.calc_checker(co_matrix, rowtotal); }
- else if (metric == "vratio") { initscore = trial.calc_vratio(rowtotal, columntotal); }
- else if (metric == "combo") { initscore = trial.calc_combo(co_matrix); }
+ //generate a probability matrix. Only do this once.
+ float start = 0.0;
+
+ if (matrix == "sim1") {
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + 1/double(nrows*ncols);
+ start = start + 1/double(nrows*ncols);
+ }
+ }
+ }
+ else if (matrix == "sim2") {
+ for(int i=0;i<nrows;i++) {
+ start = 0.0;
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + 1/double(ncols);
+ start = start + 1/double(ncols);
+ }
+ }
+ }
+
+ else if (matrix == "sim3") {
+ for(int j=0;j<ncols;j++) {
+ start = 0.0;
+ for(int i=0;i<nrows;i++) {
+ probabilityMatrix[ncols * i + j] = start + 1/double(nrows);
+ start = start + 1/double(nrows);
+ }
+ }
+ }
+
+ else if (matrix == "sim4") {
+ for(int i=0;i<nrows;i++) {
+ start = 0.0;
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n);
+ start = start + columntotal[j]/double(n);
+ }
+ }
+ }
+
+ else if (matrix == "sim5") {
+ for(int j=0;j<ncols;j++) {
+ start = 0.0;
+ for(int i=0;i<nrows;i++) {
+ probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n);
+ start = start + rowtotal[i]/double(n);
+ }
+ }
+ }
+
+ else if (matrix == "sim6") {
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + columntotal[j]/double(n*nrows);
+ start = start + columntotal[j]/double(n*nrows);
+ }
+ }
+ }
+
+
+ else if (matrix == "sim7") {
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + rowtotal[i]/double(n*ncols);
+ start = start + rowtotal[i]/double(n*ncols);
+ }
+ }
+ }
+
+ else if (matrix == "sim8") {
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ probabilityMatrix[ncols * i + j] = start + (rowtotal[i]*columntotal[j])/double(n*n);
+ start = start + (rowtotal[i]*columntotal[j])/double(n*n);
+ }
+ }
+ }
+
+ else {
+ if(sim != 9) {
+ m->mothurOut("[ERROR]: No model selected! \n");
+ m->control_pressed = true;
+ }
+ }
+
+
+
+ if (metric == "cscore") { initscore = trial.calc_c_score(initmatrix, rowtotal, ncols, nrows); }
+ else if (metric == "checker") { initscore = trial.calc_checker(initmatrix, rowtotal, ncols, nrows); }
+ else if (metric == "vratio") { initscore = trial.calc_vratio(nrows, ncols, rowtotal, columntotal); }
+ else if (metric == "combo") { initscore = trial.calc_combo(nrows, ncols, initmatrix); }
else { m->mothurOut("[ERROR]: No metric selected!\n"); m->control_pressed = true; return 1; }
m->mothurOut("Initial c score: " + toString(initscore)); m->mothurOutEndLine();
- //nullmatrix burn in
- for(int i=0;i<10000;i++) {
- if (m->control_pressed) { return 0; }
- if (matrix == "sim1") {
- trial.sim1(co_matrix);
- }else if (matrix == "sim2") {
- trial.sim2(co_matrix);
- }else if (matrix == "sim3") {
- trial.sim3(initmatrix);
- co_matrix = initmatrix;
- }else if (matrix == "sim4") {
- trial.sim4(columntotal, rowtotal, co_matrix);
- }else if (matrix == "sim5") {
- trial.sim5(initcolumntotal, initrowtotal, initmatrix);
- trial.transpose_matrix(initmatrix,co_matrix);
- }else if (matrix == "sim6") {
- trial.sim6(columntotal, co_matrix);
- }else if (matrix == "sim7") {
- trial.sim7(initcolumntotal, initmatrix);
- co_matrix = initmatrix;
- }else if (matrix == "sim8") {
- trial.sim8(columntotal, rowtotal, co_matrix);
- }else if (matrix == "sim9") {
- trial.swap_checkerboards (co_matrix);
- }else{
- m->mothurOut("[ERROR]: No model selected! \n");
- m->control_pressed = true;
+ double previous;
+ double current;
+ double randnum;
+ int count;
+
+ //burn-in
+ for(int i=0;i<10000;i++){
+ nullmatrix.clear();
+ //zero-fill the null matrix
+ nullmatrix.assign(nrows, vector<int>(ncols, 0));
+
+ if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
+ count = 0;
+ while(count < n) {
+ nextnum:
+ previous = 0.0;
+ randnum = rand() / double(RAND_MAX);
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous) {
+ nullmatrix[i][j] = 1;
+ count++;
+ if (count > n) break;
+ else
+ goto nextnum;
+ }
+ previous = current;
+ }
+ }
+ }
+ }
+
+ else if(matrix == "sim2" || matrix == "sim4") {
+ for(int i=0;i<nrows;i++) {
+ previous = 0.0;
+ count = 0;
+ while(count < rowtotal[i]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
+ }
+ previous = current;
+ }
+ }
+ }
}
+
+ else if(matrix == "sim3" || matrix == "sim5") {
+ //columns
+ for(int j=0;j<ncols;j++) {
+ count = 0;
+ while(count < columntotal[j]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int i=0;i<nrows;i++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
+ }
+ previous = current;
+ }
+ }
+ }
+ }
+
}
-
- //run
- for(int i=0;i<runs;i++) {
- if (m->control_pressed) { return 0; }
- //calc metric of nullmatrix
- if (matrix == "sim1") {
- trial.sim1(co_matrix);
- }else if (matrix == "sim2") {
- trial.sim2(co_matrix);
- }else if (matrix == "sim3") {
- trial.sim3(initmatrix);
- co_matrix = initmatrix;
- }else if (matrix == "sim4") {
- trial.sim4(columntotal, rowtotal, co_matrix);
- }else if (matrix == "sim5") {
- trial.sim5(initcolumntotal, initrowtotal, initmatrix);
- trial.transpose_matrix(initmatrix,co_matrix);
- }else if (matrix == "sim6") {
- trial.sim6(columntotal, co_matrix);
- }else if (matrix == "sim7") {
- trial.sim7(initcolumntotal, initmatrix);
- co_matrix = initmatrix;
- }else if (matrix == "sim8") {
- trial.sim8(columntotal, rowtotal, co_matrix);
- }else if (matrix == "sim9") {
- trial.swap_checkerboards (co_matrix);
- }else{
- m->mothurOut("[ERROR]: No model selected! \n");
- m->control_pressed = true;
+
+ //populate null matrix from probability matrix, do this a lot.
+ for(int i=0;i<runs;i++){
+ nullmatrix.clear();
+ //zero-fill the null matrix
+ nullmatrix.assign(nrows, vector<int>(ncols, 0));
+
+ if(matrix == "sim1" || matrix == "sim6" || matrix == "sim8" || matrix == "sim7") {
+ count = 0;
+ while(count < n) {
+ nextnum:
+ previous = 0.0;
+ randnum = rand() / double(RAND_MAX);
+ for(int i=0;i<nrows;i++) {
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous) {
+ nullmatrix[i][j] = 1;
+ count++;
+ if (count > n) break;
+ else
+ goto nextnum;
+ }
+ previous = current;
+ }
+ }
+ }
+ }
+
+ else if(matrix == "sim2" || matrix == "sim4") {
+ for(int i=0;i<nrows;i++) {
+ previous = 0.0;
+ count = 0;
+ while(count < rowtotal[i]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int j=0;j<ncols;j++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
+ }
+ previous = current;
+ }
+ }
+ }
+ }
+
+ else if(matrix == "sim3" || matrix == "sim5") {
+ //columns
+ for(int j=0;j<ncols;j++) {
+ count = 0;
+ while(count < columntotal[j]) {
+ randnum = rand() / double(RAND_MAX);
+ for(int i=0;i<nrows;i++) {
+ current = probabilityMatrix[ncols * i + j];
+ if(randnum <= current && randnum > previous && nullmatrix[i][j] != 1) {
+ nullmatrix[i][j] = 1;
+ count++;
+ previous = 0.0;
+ break;
+ }
+ previous = current;
+ }
+ }
+ }
}
- //
- //
- trial.update_row_col_totals(co_matrix, rowtotal, columntotal);
- if (metric == "cscore") {
- stats.push_back(trial.calc_c_score(co_matrix, rowtotal));
- }else if (metric == "checker") {
- stats.push_back(trial.calc_checker(co_matrix, rowtotal));
- }else if (metric == "vratio") {
- stats.push_back(trial.calc_vratio(rowtotal, columntotal));
- }else if (metric == "combo") {
- stats.push_back(trial.calc_combo(co_matrix));
- }else {
- m->mothurOut("[ERROR]: No metric selected!\n");
- m->control_pressed = true;
+ //swap_checkerboards takes the original matrix and swaps checkerboards
+ else if(matrix == "sim9") {
+ trial.swap_checkerboards (initmatrix, rowtotal, columntotal, ncols, nrows);
+ }
+ else {
+ cout << "[ERROR]: No null model selected!\n" << endl;
return 1;
}
-
+
+ //run metric on null matrix and add score to the stats vector
+ if (metric == "cscore"){
+ stats.push_back(trial.calc_c_score(nullmatrix, rowtotal, ncols, nrows));
+ }
+ else if (metric == "checker") {
+ stats.push_back(trial.calc_checker(nullmatrix, rowtotal, ncols, nrows));
+ }
+ else if (metric == "vratio") {
+ stats.push_back(trial.calc_vratio(nrows, ncols, rowtotal, columntotal));
+ }
+ else if (metric == "combo") {
+ stats.push_back(trial.calc_combo(nrows, ncols, nullmatrix));
+ }
+ else {
+ cout << "[ERROR]: No metric selected!\n" << endl;
+ return 1;
+ }
+
}
+
+
double total = 0.0;
for (int i=0; i<stats.size();i++) { total+=stats[i]; }
projects / mothur.git / blobdiff