#include "parsimonycommand.h"
#include "unifracunweightedcommand.h"
#include "unifracweightedcommand.h"
+#include "libshuffcommand.h"
#include "mothur.h"
else if(commandName == "summary.shared") { command = new SummarySharedCommand(); }
else if(commandName == "unifrac.weighted") { command = new UnifracWeightedCommand(); }
else if(commandName == "unifrac.unweighted") { command = new UnifracUnweightedCommand(); }
- else if(commandName == "get.group") { command = new GetgroupCommand(); }
- else if(commandName == "get.label") { command = new GetlabelCommand(); }
- else if(commandName == "get.line") { command = new GetlineCommand(); }
+ else if(commandName == "get.group") { command = new GetgroupCommand(); }
+ else if(commandName == "get.label") { command = new GetlabelCommand(); }
+ else if(commandName == "get.line") { command = new GetlineCommand(); }
+ else if(commandName == "libshuff") { command = new LibShuffCommand(); }
else { command = new NoCommand(); }
return command;
#include "coverage.h"
//**********************************************************************************************************************
-vector< vector<float> > Coverage::getValues(FullMatrix*, float) {
- try {
+Coverage::Coverage() {
globaldata = GlobalData::getInstance();
- numGroups = globaldata->Groups.size();
+ numUserGroups = globaldata->Groups.size();
+ numGroups = globaldata->gGroupmap->getNumGroups();
+}
+//**********************************************************************************************************************
+void Coverage::getValues(FullMatrix* matrix, float d, vector< vector<float> >& data) {
+ try {
+ vector<float> min;
+ vector<string> groups;
//initialize data
- data.resize(numGroups);
+ data.resize(numUserGroups);
for (int l = 0; l < data.size(); l++) {
data[l].push_back(0.0);
}
/**************************************/
- //get the minumums for each comparision
+ //get the minimums for each comparision
/**************************************/
- return data;
+ int count = 0;
+ int a = 0;
+ int b = 0;
+ for (int i = 0; i < numGroups; i++) {
+ for (int j = 0; j < numGroups; j++) {
+
+ //is this "box" one hte user wants analyzed?
+ if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
+
+ min = matrix->getMins(count); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
+
+ //find the coverage at this distance
+ sort(min.begin(), min.end());
+//cout << "minvector for : " << globaldata->gGroupmap->namesOfGroups[i] + globaldata->gGroupmap->namesOfGroups[j] << endl;
+//for(int h = 0; h<min.size(); h++) {
+// cout << min[h] << " ";
+//}
+//cout << endl;
+ int index = -1;
+ //find index in min where value is higher than d
+ for (int m = 0; m < min.size(); m++) {
+ if (min[m] > d) { index = m; break; }
+ }
+
+ // if you don't find one than all the mins are less than d
+ if (index == -1) { index = min.size(); }
+
+ //save value in data
+ data[a][b] = 1.0 - ((min.size()-index)/(float)min.size());
+//cout << "D = " << d << "data " << a << b << " = " << data[a][b] << endl;
+ if (b < numUserGroups-1) { b++; }
+ else{ //you are moving to a new row of "boxes"
+ b = 0;
+ a++;
+ }
+ }
+ count++;
+ }
+ }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the Coverage class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+//For the random matrices
+void Coverage::getValues(FullMatrix* matrix, float d, vector< vector<float> >& data, string r) {
+ try {
+ vector<float> min;
+ vector<string> groups;
+
+ //initialize data
+ data.resize(numUserGroups);
+ for (int l = 0; l < data.size(); l++) {
+ data[l].push_back(0.0);
+ }
+
+ /**************************************/
+ //get the minimums for each comparision
+ /**************************************/
+ int count = 0;
+ int a = 0;
+ int b = 0;
+ for (int i = 0; i < numGroups; i++) {
+ for (int j = 0; j < numGroups; j++) {
+
+ //is this "box" one hte user wants analyzed?
+ if ((inUsersGroups(globaldata->gGroupmap->namesOfGroups[i], globaldata->Groups) == true) && (inUsersGroups(globaldata->gGroupmap->namesOfGroups[j], globaldata->Groups) == true)) {
+cout << "combo " << a << b << endl;
+cout << "original matrix mins4rows " << endl;
+matrix->printMinsForRows(cout);
+ //create random matrix for this comparison
+ matrix->shuffle(count);
+
+ min = matrix->getMins(count); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
+cout << "shuffled matrix mins4rows " << endl;
+matrix->printMinsForRows(cout);
+
+ //find the coverage at this distance
+ sort(min.begin(), min.end());
+
+ int index = -1;
+ //find index in min where value is higher than d
+ for (int m = 0; m < min.size(); m++) {
+ if (min[m] > d) { index = m; break; }
+ }
+
+ // if you don't find one than all the mins are less than d
+ if (index == -1) { index = min.size(); }
+
+ //save value in data
+ data[a][b] = 1.0 - ((min.size()-index)/(float)min.size());
+cout << "D = " << d << "data " << a << b << " = " << data[a][b] << endl;
+ if (b < numUserGroups-1) { b++; }
+ else{ //you are moving to a new row of "boxes"
+ b = 0;
+ a++;
+ }
+ }
+ count++;
+
+ //restore matrix to original form for next shuffle
+ matrix->restore();
+min = matrix->getMins(count-1);
+cout << "restored matrix mins4rows " << endl;
+matrix->printMinsForRows(cout);
+ }
+ }
}
catch(exception& e) {
cout << "Standard Error: " << e.what() << " has occurred in the Coverage class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
class Coverage {
public:
- Coverage(){};
+ Coverage();
~Coverage(){};
- vector< vector<float> > getValues(FullMatrix*, float);
+ void getValues(FullMatrix*, float, vector< vector<float> >&);
+ void getValues(FullMatrix*, float, vector< vector<float> >&, string); //for random matrices
private:
GlobalData* globaldata;
- vector< vector<float> > data;
- int numGroups;
+ int numGroups, numUserGroups;
};
}else if (listfile != "") { //you want to do single commands
validateReadFiles();
validateReadPhil();
- }else {//you are reading a shared file
+ }else if ((listfile == "") && (sharedfile == "")) {
+ cout << "You must enter either a listfile or a sharedfile with the read.otu command. " << endl; return false;
+ }else{//you are reading a shared file
validateReadFiles();
}
}else if (commandName == "read.tree") {
errorFree = false;
}
+ if ((commandName == "libshuff") && (globaldata->gMatrix == NULL)) {
+ cout << "You must read in a matrix before you use the libshuff command. " << endl; return false;
+ }
+
if (commandName == "parsimony") {
//are you trying to use parsimony without reading a tree or saying you want random distribution
if (randomtree == "") {
//unable to open
if (ableToOpen == 1) { errorFree = false; }
else { globaldata->inputFileName = phylipfile; }
- //are we reading a phylipfile
+ //are we reading a columnfile
}else if (columnfile != "") {
ableToOpen = openInputFile(columnfile, filehandle);
filehandle.close();
ifstream filehandle;
int ableToOpen;
+ if (groupfile != "") {
+ ableToOpen = openInputFile(groupfile, filehandle);
+ filehandle.close();
+ //unable to open
+ if (ableToOpen == 1) { errorFree = false; }
+ }
+
if ((phylipfile == "") && (columnfile == "")) { cout << "When executing a read.dist you must enter a phylip or a column." << endl; errorFree = false; }
else if ((phylipfile != "") && (columnfile != "")) { cout << "When executing a read.dist you must enter ONLY ONE of the following: phylip or column." << endl; errorFree = false; }
reading = new Progress("Reading matrix: ", numSeqs * numSeqs);
int count = 0;
- float distance;
string group, name;
if(group == "not found") { cout << "Error: Sequence '" << name << "' was not found in the group file, please correct." << endl; exit(1); }
for(int j=0;j<numSeqs;j++){
- filehandle >> distance;
-
- matrix[i][j] = distance;
+ filehandle >> matrix[i][j];
+
count++;
reading->update(count);
}
int i = low;
int j = high;
- int y = 0;
+ float y = 0;
string name;
/* compare value */
for (int i = 0; i < numSeqs; i++) {
out << "row " << i << " group = " << index[i].groupname << " name = " << index[i].seqName << endl;
for (int j = 0; j < numSeqs; j++) {
- out << matrix[i][j] << " ";
+ //out << matrix[i][j] << " ";
}
out << endl;
}
}
/**************************************************************************/
-void FullMatrix::getMinsForRowsVectors(){
+void FullMatrix::setBounds(){
try{
numGroups = globaldata->gGroupmap->namesOfGroups.size();
//sort globaldata->gGroupmap.namesOfGroups so that it will match the matrix
sort(globaldata->gGroupmap->namesOfGroups.begin(), globaldata->gGroupmap->namesOfGroups.end());
+ //one for each comparision
+ //minsForRows.resize(numGroups*numGroups);
+
/*************************************************/
//find where in matrix each group starts and stops
/*************************************************/
- vector<int> bounds; //bounds[1] = starting row in matrix from group B, bounds[2] = starting row in matrix from group C, bounds[3] = no need to find upper bound of C because its numSeqs.
bounds.resize(numGroups);
bounds[0] = 0;
- bounds[numGroups] = numSeqs-1;
+ bounds[numGroups] = numSeqs;
+
//for each group find bounds of subgroup/comparison
for (int i = 1; i < numGroups; i++) {
- getBounds(bounds[i], globaldata->gGroupmap->namesOfGroups[i]);
+ getBounds(bounds[i], globaldata->gGroupmap->namesOfGroups[i-1]);
}
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the FullMatrix class function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+
+}
+/**************************************************************************/
+vector<float> FullMatrix::getMins(int x) {
+ try{
+ //clear out old data
+ minsForRows.clear();
+
/************************************************************/
- //fill the minsForRows vectors for each group the user wants
+ //fill the minsForRows vector for the box the user wants
/************************************************************/
- int countx = bounds[1]; //where second group starts
- int county = bounds[1];
+ int count = 0;
+ int lowBoundx = bounds[0]; //where first group starts
+ int lowBoundy = bounds[0];
+ int highBoundx = bounds[1]; //where second group starts
+ int highBoundy = bounds[1];
- //go through the entire matrix
- for (int x = 0; x < numSeqs; x++) {
- for (int y = 0; y < numSeqs; y++) {
- //if have not changed groups
- if ((x < countx) && (y < county)) {
-
- }
+ int countx = 1; //index in bound
+ int county = 1; //index in bound
+
+ //find the bounds for the box the user wants
+ for (int i = 0; i < (numGroups * numGroups); i++) {
+
+ //are you at the box?
+ if (count == x) { break; }
+ else { count++; }
+
+ //move to next box
+ if (county < numGroups) {
+ county++;
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
+ }else{ //you are moving to a new row of "boxes"
+ county = 1;
+ countx++;
+ highBoundx = bounds[countx];
+ lowBoundx = bounds[countx-1];
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
}
}
-
+ //each row in the box
+ for (int x = lowBoundx; x < highBoundx; x++) {
+ float min4Row = 100000.0;
+ //each entry in that row
+ for (int y = lowBoundy; y < highBoundy; y++) {
+ //if you are not on the diagonal and you are less than previous minimum
+ if ((x != y) && (matrix[x][y] < min4Row)) {
+ min4Row = matrix[x][y];
+ }
+ }
+ //save minimum value for that row in minsForRows vector of vectors
+ minsForRows.push_back(min4Row);
+ }
-
+ return minsForRows;
}
catch(exception& e) {
- cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function getMins. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
catch(...) {
- cout << "An unknown error has occurred in the FullMatrix class function getMinsForRowsVectors. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ cout << "An unknown error has occurred in the FullMatrix class function getMins. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
exit(1);
}
-
}
-
/**************************************************************************/
void FullMatrix::getBounds(int& higher, string group) {
try{
//for each group find bounds of subgroup/comparison
for (it = index.begin(); it != index.end(); it++) {
if (it->second.groupname == group) {
- if (gotLower != true) { gotLower = true; }
+ gotLower = true;
}else if ((gotLower == true) && (it->second.groupname != group)) { higher = it->first; break; }
}
}
+/**************************************************************************/
+//print out matrix
+void FullMatrix::printMinsForRows(ostream& out) {
+ try{
+ for (int j = 0; j < minsForRows.size(); j++) {
+ out << minsForRows[j] << " ";
+ }
+ out << endl;
+
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the FullMatrix class function printMatrix. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+
+}
+/**************************************************************************/
+//shuffles the sequences in the 2 groups passed in.
+void FullMatrix::shuffle(int box){
+ try{
+ vector<int> rows2Swap;
+ vector<int> shuffled;
+ float y = 0;
+ string name = "";
+
+ /****************************/
+ //find the box the user wants
+ /****************************/
+ int count = 0;
+ int lowBoundy = bounds[0]; //where first group starts
+ int highBoundy = bounds[1]; //where second group starts
+ int county = 1; //index in bound
+
+ //find the bounds for the box the user wants
+ for (int i = 0; i < (numGroups * numGroups); i++) {
+
+ //are you at the box?
+ if (count == box) { break; }
+ else { count++; }
+
+ //move to next box
+ if (county < numGroups) {
+ county++;
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
+ }else{ //you are moving to a new row of "boxes"
+ county = 1;
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
+ }
+ }
+
+ /************************/
+ //save its rows locations
+ /************************/
+ //go through the matrix map to find the rows from groups you want to randomize
+ for (int y = lowBoundy; y < highBoundy; y++) {
+ rows2Swap.push_back(y);
+ shuffled.push_back(y);
+ }
+
+ //randomize rows to shuffle in shuffled
+ random_shuffle(shuffled.begin(), shuffled.end());
+
+ /***************************************/
+ //swap rows and columns to randomize box
+ /***************************************/
+ for (int i = 0; i < shuffled.size(); i++) {
+ //record the swaps you are making so you can undo them in restore function
+ restoreIndex[i].a = shuffled[i];
+ restoreIndex[i].b = rows2Swap[i];
+
+ /* swap rows*/
+ for (int h = 0; h < numSeqs; h++) {
+ y = matrix[shuffled[i]][h];
+ matrix[shuffled[i]][h] = matrix[rows2Swap[i]][h];
+ matrix[rows2Swap[i]][h] = y;
+ }
+
+ /* swap columns */
+ for (int b = 0; b < numSeqs; b++) {
+ y = matrix[b][shuffled[i]];
+ matrix[b][shuffled[i]] = matrix[b][rows2Swap[i]];
+ matrix[b][rows2Swap[i]] = y;
+ }
+
+ //swap map elements
+ name = index[shuffled[i]].seqName;
+ index[shuffled[i]].seqName = index[rows2Swap[i]].seqName;
+ index[rows2Swap[i]].seqName = name;
+ }
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function shuffle. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the FullMatrix class function shuffle. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/**************************************************************************/
+//unshuffles the matrix.
+void FullMatrix::restore(){
+ try{
+ float y = 0;
+ string name = "";
+
+ //reverse iterate through swaps and undo them to restore original matrix and index map.
+ for(it2 = restoreIndex.rbegin(); it2 != restoreIndex.rend(); it2++) {
+ /* swap rows */
+ for (int h = 0; h < numSeqs; h++) {
+ y = matrix[it2->second.a][h];
+ matrix[it2->second.a][h] = matrix[it2->second.b][h];
+ matrix[it2->second.b][h] = y;
+ }
+
+ /* swap columns */
+ for (int b = 0; b < numSeqs; b++) {
+ y = matrix[b][it2->second.a];
+ matrix[b][it2->second.a] = matrix[b][it2->second.b];
+ matrix[b][it2->second.b] = y;
+ }
+
+
+ //swap map elements
+ name = index[it2->second.a].seqName;
+ index[it2->second.a].seqName = index[it2->second.b].seqName;
+ index[it2->second.b].seqName = name;
+ }
+
+ //clear restore for next shuffle
+ restoreIndex.clear();
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the FullMatrix class function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+/**************************************************************************/
+void FullMatrix::getDist(vector<float>& distances) {
+ try{
+ map<float, float> dist; //holds the distances for the integral form
+ map<float, float>::iterator it;
+
+ /************************************************************/
+ //fill the minsForRows vectors for each group the user wants
+ /************************************************************/
+ int lowBoundx = bounds[0]; //where first group starts
+ int lowBoundy = bounds[0];
+ int highBoundx = bounds[1]; //where second group starts
+ int highBoundy = bounds[1];
+
+ int countx = 1; //index in bound
+ int county = 1; //index in bound
+
+ //go through each "box" in the matrix
+ for (int i = 0; i < (numGroups * numGroups); i++) {
+ //each row in the box
+ for (int x = lowBoundx; x < highBoundx; x++) {
+ float min4Row = 100000.0;
+ //each entry in that row
+ for (int y = lowBoundy; y < highBoundy; y++) {
+ //if you are not on the diagonal and you are less than previous minimum
+ if ((x != y) && (matrix[x][y] < min4Row)){
+ min4Row = matrix[x][y];
+ }
+ }
+ //save minimum value
+ dist[min4Row] = min4Row;
+ }
+
+ //****** reset bounds to process next "box" ********
+ //if you still have more "boxes" in that row
+ if (county < numGroups) {
+ county++;
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
+ }else{ //you are moving to a new row of "boxes"
+ county = 1;
+ countx++;
+ highBoundx = bounds[countx];
+ lowBoundx = bounds[countx-1];
+ highBoundy = bounds[county];
+ lowBoundy = bounds[county-1];
+ }
+ }
+
+ //store distances in users vector
+ for (it = dist.begin(); it != dist.end(); it++) {
+ distances.push_back(it->first);
+ }
+
+ }
+ catch(exception& e) {
+ cout << "Standard Error: " << e.what() << " has occurred in the FullMatrix class Function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+ catch(...) {
+ cout << "An unknown error has occurred in the FullMatrix class function restore. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+ exit(1);
+ }
+}
+
string seqName;
};
+struct Swap {
+ int a;
+ int b;
+};
+
class FullMatrix {
int getNumSeqs();
void printMatrix(ostream&);
- void getMinsForRowsVectors(); //requires globaldata->Groups to be filled
-
+ void setBounds(); //requires globaldata->Groups to be filled
+ vector<float> getMins(int); //returns vector of mins for "box" requested ie. groups A, B, 0 = AA, 1 = AB, 2 = BA, 3 = BB;
+ void getDist(vector<float>&); //fills a vector with the valid distances for the integral form.
+ void shuffle(int); //shuffles the sequences in the box passed in.
+ void restore(); //unshuffles the matrix.
+
+
+ void printMinsForRows(ostream&);
private:
void sortGroups(int, int); //this function sorts the sequences within the matrix.
void getBounds(int&, string);
void readSquareMatrix(ifstream&);
void readLTMatrix(ifstream&);
- vector< vector<float> > matrix; //a 2D distance matrix of all the sequences and their distances to eachother.
- vector< vector<float> > minsForRows; //vector< minimum distance for that subrow> -one for each comparison.
+
map<int, Names> index; // row in vector, sequence group. need to know this so when we sort it can be updated.
+ map<int, Swap> restoreIndex; //a map of the swaps made so you can undo them in restore.
map<int, Names>::iterator it;
+ map<int, Swap>::reverse_iterator it2;
+
+ vector< vector<float> > matrix; //a 2D distance matrix of all the sequences and their distances to eachother.
+ vector<float> minsForRows; //vector< minimum distance for that subrow> - one for each comparison.
+ vector<int> bounds; //bounds[1] = starting row in matrix from group B, bounds[2] = starting row in matrix from group C, bounds[3] = no need to find upper bound of C because its numSeqs.
+
GroupMap* groupmap; //maps sequences to groups they belong to.
GlobalData* globaldata;
int numSeqs, numGroups, numUserGroups;
helpRequest = optionText;
}
+ if (commandName == "libshuff") {
+ iters = "10000";
+ cutoff = "1.0";
+ }
+
string key, value;
//reads in parameters and values
if((optionText != "") && (commandName != "help")){
if (key == "fileroot" ) { fileroot = value; }
if (key == "abund" ) { abund = value; }
if (key == "random" ) { randomtree = value; }
- if (key == "calc") { calc = value; }
+ if (key == "calc") { calc = value; }
+ if (key == "step") { step = value; }
+ if (key == "form") { form = value; }
+
if (key == "line") {//stores lines to be used in a set
if (key == "abund" ) { abund = value; }
if (key == "random" ) { randomtree = value; }
if (key == "calc") { calc = value; }
-
+ if (key == "step") { step = value; }
+ if (key == "form") { form = value; }
if (key == "line") {//stores lines to be used in a vector
lines.clear();
string GlobalData::getAbund() { return abund; }
string GlobalData::getRandomTree() { return randomtree; }
string GlobalData::getGroups() { return groups; }
+string GlobalData::getStep() { return step; }
+string GlobalData::getForm() { return form; }
void GlobalData::setListFile(string file) { listfile = file; inputFileName = file;}
void GlobalData::setRabundFile(string file) { rabundfile = file; inputFileName = file;}
void GlobalData::setSabundFile(string file) { sabundfile = file; inputFileName = file;}
method = "furthest";
fileroot = "";
abund = "10";
+ step = "0.01";
+ form = "integral";
}
//*******************************************************/
freq = "100";
method = "furthest";
calc = "";
- abund = "10";
+ abund = "10";
+ step = "0.01";
+ form = "integral";
}
/*******************************************************/
string getAbund();
string getRandomTree();
string getGroups();
+ string getStep();
+ string getForm();
void setListFile(string);
void setPhylipFile(string);
private:
string phylipfile, columnfile, listfile, rabundfile, sabundfile, namefile, groupfile, orderfile, fastafile, treefile, sharedfile, line, label, randomtree, groups;
- string cutoff, format, precision, method, fileroot, iters, jumble, freq, calc, abund;
+ string cutoff, format, precision, method, fileroot, iters, jumble, freq, calc, abund, step, form;
static GlobalData* _uniqueInstance;
GlobalData( const GlobalData& ); // Disable copy constructor
globaldata = GlobalData::getInstance();
convert(globaldata->getCutOff(), cutOff);
convert(globaldata->getIters(), iters);
+ convert(globaldata->getStep(), step);
+ form = globaldata->getForm();
matrix = globaldata->gMatrix;
coverageFile = getRootName(globaldata->getPhylipFile()) + "coverage";
summaryFile = getRootName(globaldata->getPhylipFile()) + "slsummary";
try {
//deltaValues[0] = scores for the difference between AA and AB.
//cValues[0][0] = AA, cValues[0][1] = AB, cValues[0][2] = AC, cValues[1][0] = BA, cValues[1][1] = BB...
+ vector<float> dist;
+ int next;
sumDelta.resize(numComp-numGroups, 0.0);
-
- float D = 0.0;
+ float D = 0.0;
+
/*****************************/
//get values for users matrix
/*****************************/
- matrix->getMinsForRowsVectors();
+ matrix->setBounds();
+ if (form != "discrete") { matrix->getDist(dist); next = 1; }
+//cout << "Distances" << endl;
+//for (int i = 0; i < dist.size(); i++) { cout << dist[i] << " "; }
+//cout << endl;
+
//get values for users matrix
while (D <= cutOff) {
- cValues = coverage->getValues(matrix, D);
+ //clear out old Values
+ deltaValues.clear();
+ coverage->getValues(matrix, D, cValues);
//find delta values
int count = 0;
}
}
- D += 0.01;
-
printCoverageFile(D);
- //clear out old Values
- cValues.clear();
- deltaValues.clear();
+ //check form
+ if (form != "discrete") {
+ if (next == dist.size()) { break; }
+ else { D = dist[next]; next++; }
+ }else { D += step; }
+
+
}
+ //output sum Deltas
+ for (int i = 0; i < numGroups; i++) {
+ for (int j = 0; j < numGroups; j++) {
+ //don't output AA to AA
+ if (i != j) {
+ cout << "Delta " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t';
+ }
+ }
+ }
+ cout << endl;
+
+ for (int i = 0; i < sumDelta.size(); i++) {
+ cout << setprecision(6) << sumDelta[i] << '\t';
+ }
+ cout << endl;
+
/*******************************************************************************/
//create and score random matrixes finding the sumDelta values for summary file
/******************************************************************************/
}
}
+
for (int m = 0; m < iters; m++) {
//generate random matrix in getValues
//values for random matrix
+ cout << "Iteration " << m+1 << endl;
+ D = 0.0;
+ next = 1;
+
while (D <= cutOff) {
- cValues = coverage->getValues(matrix, D);
+ coverage->getValues(matrix, D, cValues, "random");
//find delta values
int count = 0;
//don't save AA to AA
if (i != j) {
//(Caa - Cab)^2
- rsumDelta[count][m] += (abs(cValues[i][i]-cValues[i][j]) * abs(cValues[i][i]-cValues[i][j]));
+ rsumDelta[count][m] += ((abs(cValues[i][i]-cValues[i][j]) * abs(cValues[i][i]-cValues[i][j])));
+//cout << "iter " << m << " box " << i << j << " delta = " << ((abs(cValues[i][i]-cValues[i][j]) * abs(cValues[i][i]-cValues[i][j]))) << endl;
count++;
}
}
}
-
- D += 0.01;
+
+ //check form
+ if (form != "discrete") {
+ if (next == dist.size()) { break; }
+ else { D = dist[next]; next++; }
+ }else { D += step; }
+
//clear out old Values
cValues.clear();
}
+cout << "random sum delta for iter " << m << endl;
+for (int i = 0; i < rsumDelta.size(); i++) {
+ cout << setprecision(6) << rsumDelta[i][m] << '\t';
+}
+cout << endl;
+
}
/**********************************************************/
//format output
out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
- out << setprecision(globaldata->getIters().length()) << d << '\t';
+ out << setprecision(6) << d << '\t';
//print out coverage values
for (int i = 0; i < numGroups; i++) {
for (int j = 0; j < numGroups; j++) {
//don't output AA to AA
if (i != j) {
- outSum << "Delta" + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'<< "DeltaSig" + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t';
+ outSum << "Delta " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'<< "DeltaSig " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t';
+ cout << "Delta " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t'<< "DeltaSig " + globaldata->Groups[i] + "-" + globaldata->Groups[j] << '\t';
}
}
}
outSum << endl;
+ cout << endl;
//print out delta values
for (int i = 0; i < sumDelta.size(); i++) {
outSum << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << sumDeltaSig[i] << '\t';
+ cout << setprecision(6) << sumDelta[i] << '\t' << setprecision(globaldata->getIters().length()) << sumDeltaSig[i] << '\t';
}
outSum << endl;
+ cout << endl;
}
catch(exception& e) {
}
}
+ //sort so labels match
+ sort(globaldata->Groups.begin(), globaldata->Groups.end());
+
// number of comparisons i.e. with groups A,B,C = AA, AB, AC, BA, BB, BC...;
for (int i=0; i<numGroups; i++) {
for (int l = 0; l < numGroups; l++) {
GlobalData* globaldata;
Coverage* coverage;
FullMatrix* matrix;
- float cutOff;
+ float cutOff, step;
int numGroups, numComp, iters;
- string coverageFile, summaryFile;
+ string coverageFile, summaryFile, form;
ofstream out, outSum;
commands["summary.shared"] = "summary.shared";
commands["unifrac.weighted"] = "unifrac.weighted";
commands["unifrac.unweighted"] = "unifrac.unweighted";
+ commands["libshuff"] = "libshuff";
}
parameters["random"] = "random";
parameters["groups"] = "groups";
parameters["calc"] = "calc";
+ parameters["step"] = "step";
+ parameters["form"] = "form";
}
catch(exception& e) {
projects / mothur.git / commitdiff