- vector<string> subclasses;
- map<string, string> subclass2Class;
- map<string, int> subclassCounts;
- map<string, vector<int> > subClass2GroupIndex; //maps subclass name to vector of indexes in lookup from that subclass. old -> 1,2,3 means groups in location 1,2,3 of lookup are from old. Saves time below.
- bool error = false;
- for (int j = 0; j < lookup.size(); j++) {
- string group = lookup[j]->getGroup();
- string treatment = designMap.get(group, mclass); //get value for this group in this category
- string thisSub = designMap.get(group, subclass);
- map<string, string>::iterator it = subclass2Class.find(thisSub);
- if (it == subclass2Class.end()) {
- subclass2Class[thisSub] = treatment;
- subclassCounts[thisSub] = 1;
- vector<int> temp; temp.push_back(j);
- subClass2GroupIndex[thisSub] = temp;
+
+ int numBins = lookup[0]->getNumBins();
+ for (int i = 0; i < numBins; i++) {
+ if (m->control_pressed) { break; }
+
+ it = bins.find(i);
+ if (it != bins.end()) { //flagged in Kruskal Wallis
+
+ vector<float> abunds; for (int j = 0; j < lookup.size(); j++) { abunds.push_back(lookup[j]->getAbundance(i)); }
+
+ bool sig = testOTUWilcoxon(class2SubClasses, abunds, subClass2GroupIndex, subclass2Class);
+ if (sig) { significantOtuLabels[i] = it->second; }
+
+ }//bins flagged from kw
+ }//for bins
+
+ return significantOtuLabels;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LefseCommand", "runWilcoxon");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+//lefse.py - test_rep_wilcoxon_r function
+bool LefseCommand::testOTUWilcoxon(map<string, set<string> >& class2SubClasses, vector<float> abunds, map<string, vector<int> >& subClass2GroupIndex, map<string, string> subclass2Class) {
+ try {
+ int totalOk = 0;
+ double alphaMtc = wilcoxonAlpha;
+ vector< set<string> > allDiffs;
+ LinearAlgebra linear;
+
+ //for each subclass comparision
+ map<string, set<string> >::iterator itB;
+ for(map<string, set<string> >::iterator it=class2SubClasses.begin();it!=class2SubClasses.end();it++){
+ itB = it;itB++;
+ for(itB;itB!=class2SubClasses.end();itB++){
+ if (m->control_pressed) { return false; }
+ bool first = true;
+ int dirCmp = 0; // not set?? dir_cmp = "not_set" # 0=notset or none, 1=true, 2=false.
+ int curv_sign = 0;
+ int ok = 0;
+ int count = 0;
+ for (set<string>::iterator itClass1 = (it->second).begin(); itClass1 != (it->second).end(); itClass1++) {
+ bool br = false;
+ for (set<string>::iterator itClass2 = (itB->second).begin(); itClass2 != (itB->second).end(); itClass2++) {
+ string subclass1 = *itClass1;
+ string subclass2 = *itClass2;
+ count++;
+
+ if (m->debug) { m->mothurOut( "[DEBUG comparing " + it->first + "-" + *itClass1 + " to " + itB->first + "-" + *itClass2 + "\n"); }
+
+ string treatment1 = subclass2Class[subclass1];
+ string treatment2 = subclass2Class[subclass2];
+ int numSubs1 = class2SubClasses[treatment1].size();
+ int numSubs2 = class2SubClasses[treatment2].size();
+
+ //if mul_cor != 0: alpha_mtc = th*l_subcl1*l_subcl2 if mul_cor == 2 else 1.0-math.pow(1.0-th,l_subcl1*l_subcl2)
+ if (strict != 0) { alphaMtc = wilcoxonAlpha * numSubs1 * numSubs2 ; }
+ if (strict == 2) {}else{ alphaMtc = 1.0-pow((1.0-wilcoxonAlpha),(double)(numSubs1 * numSubs2)); }
+
+ //fill x and y with this comparisons data
+ vector<double> x; vector<double> y;
+
+ //fill x and y
+ vector<int> xIndexes = subClass2GroupIndex[subclass1]; //indexes in lookup for this subclass
+ vector<int> yIndexes = subClass2GroupIndex[subclass2]; //indexes in lookup for this subclass
+ for (int k = 0; k < yIndexes.size(); k++) { y.push_back(abunds[yIndexes[k]]); }
+ for (int k = 0; k < xIndexes.size(); k++) { x.push_back(abunds[xIndexes[k]]); }
+
+ // med_comp = False
+ //if len(cl1) < min_c or len(cl2) < min_c:
+ //med_comp = True
+ bool medComp = false; // are there enough samples per subclass
+ if ((xIndexes.size() < minC) || (yIndexes.size() < minC)) { medComp = true; }
+
+ double sx = m->median(x);
+ double sy = m->median(y);
+
+ //if cl1[0] == cl2[0] and len(set(cl1)) == 1 and len(set(cl2)) == 1:
+ //tres, first = False, False
+ double pValue = 0.0;
+ double H = 0.0;
+ bool tres = true; //don't think this is set in the python source. Not sure how that is handled, but setting it here.
+ if ((x[0] == y[0]) && (x.size() == 1) && (y.size() == 1)) { tres = false; first = false; }
+ else if (!medComp) {
+ H = linear.calcWilcoxon(x, y, pValue);
+ if (pValue < (alphaMtc*2.0)) { tres = true; }
+ else { tres = false; }
+ }
+ /*if first:
+ first = False
+ if not curv and ( med_comp or tres ):
+ dir_cmp = sx < sy
+ if sx == sy: br = True
+ elif curv:
+ dir_cmp = None
+ if med_comp or tres:
+ curv_sign += 1
+ dir_cmp = sx < sy
+ else: br = True
+ elif not curv and med_comp:
+ if ((sx < sy) != dir_cmp or sx == sy): br = True
+ elif curv:
+ if tres and dir_cmp == None:
+ curv_sign += 1
+ dir_cmp = sx < sy
+ if tres and dir_cmp != (sx < sy):
+ br = True
+ curv_sign = -1
+ elif not tres or (sx < sy) != dir_cmp or sx == sy: br = True
+ */
+ int sxSy = 2; //false
+ if (sx<sy) { sxSy = 1; } //true
+
+ if (first) {
+ first = false;
+ if ((!curv) && (medComp || tres)) {
+ dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
+ if (sx == sy) { br = true; }
+ }else if (curv) {
+ dirCmp = 0;
+ if (medComp || tres) {
+ curv_sign++;
+ dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
+ }
+ }else { br = true; }
+ }else if (!curv && medComp) {
+ if (sxSy != dirCmp || sx == sy) { br = true; }
+ }else if (curv) {
+ if (tres && dirCmp == 0) { curv_sign++; }
+ dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
+ if (tres && dirCmp != sxSy) { //if tres and dir_cmp != (sx < sy):
+ br = true;
+ curv_sign = -1;
+ }
+ }else if (!tres || sxSy != dirCmp || sx == sy) { br = true; } //elif not tres or (sx < sy) != dir_cmp or sx == sy: br = True
+ if (br) { break; }
+ ok++;
+ }//for class2 subclasses
+ if (br) { break; }
+ }//for class1 subclasses
+ bool diff = false;
+ if (curv) { diff = false; if (curv_sign > 0) { diff = true; } } //if curv: diff = curv_sign > 0
+ else { //else: diff = (ok == len(cl_hie[pair[1]])*len(cl_hie[pair[0]]))
+ diff = false;
+ if (ok == count) { diff = true; }
+ }
+ if (diff) { totalOk++; }
+ if (!diff && (multiClassStrat == "onevone")) { return false; }
+ if (diff && (multiClassStrat == "onevall")) { //all_diff.append(pair)
+ set<string> pair; pair.insert(it->first); pair.insert(itB->first);
+ allDiffs.push_back(pair);
+ }
+ }//classes
+ }//classes
+
+ if (multiClassStrat == "onevall") {
+ int tot_k = class2SubClasses.size();
+ for(map<string, set<string> >::iterator it=class2SubClasses.begin();it!=class2SubClasses.end();it++){
+ if (m->control_pressed) { return false; }
+ int nk = 0;
+ //is this class okay in all comparisons
+ for (int h = 0; h < allDiffs.size(); h++) {
+ if (allDiffs[h].count(it->first) != 0) { nk++; }
+ }
+ if (nk == (tot_k-1)) { return true; }//if nk == tot_k-1: return True