//**********************************************************************************************************************
PCACommand::PCACommand(){
try {
- abort = true;
- //initialize outputTypes
+ abort = true; calledHelp = true;
vector<string> tempOutNames;
outputTypes["pca"] = tempOutNames;
outputTypes["loadings"] = tempOutNames;
PCACommand::PCACommand(string option) {
try {
- abort = false;
+ abort = false; calledHelp = false;
globaldata = GlobalData::getInstance();
//allow user to run help
- if(option == "help") { help(); abort = true; }
+ if(option == "help") { help(); abort = true; calledHelp = true; }
else {
//valid paramters for this command
metric = m->isTrue(temp);
label = validParameter.validFile(parameters, "label", false);
- if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); }
+ if (label == "not found") { label = ""; labels = globaldata->labels; if(labels.size() == 0) { m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); } }
else { m->splitAtDash(label, labels); }
groups = validParameter.validFile(parameters, "groups", false);
//**********************************************************************************************************************
void PCACommand::help(){
try {
- m->mothurOut("The pca command can only be run after a successful read.otu command."); m->mothurOutEndLine();
+ m->mothurOut("The pca command can only be run after a successful read.otu command of a shared or relabund file."); m->mothurOutEndLine();
m->mothurOut("The pca command parameters are label, groups and metric. No parameters are required."); m->mothurOutEndLine();
- m->mothurOut("The label parameter is used to analyze specific labels in your input. Default is the first label in your shared or relabund file. Multpile labels may be separated by dashes.\n");
+ m->mothurOut("The label parameter is used to analyze specific labels in your input. Default is the first label in your shared or relabund file. Multiple labels may be separated by dashes.\n");
m->mothurOut("The groups parameter allows you to specify which groups you would like analyzed. Groupnames are separated by dashes.\n");
m->mothurOut("The metric parameter allows indicate you if would like the pearson correlation coefficient calculated. Default=True"); m->mothurOutEndLine();
m->mothurOut("Example pca(groups=yourGroups).\n");
int PCACommand::execute(){
try {
- if (abort == true) { return 0; }
+ if (abort == true) { if (calledHelp) { return 0; } return 2; }
cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::showpoint);
exit(1);
}
}
-//**********************************************************************************************************************
+/**********************************************************************************************************************
vector< vector<double> > PCACommand::createMatrix(vector<SharedRAbundFloatVector*> lookupFloat){
try {
vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
m->errorOut(e, "PCACommand", "createMatrix");
exit(1);
}
-}
+}*/
//**********************************************************************************************************************
int PCACommand::process(vector<SharedRAbundFloatVector*>& lookupFloat){
try {
vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
+ ofstream out;
+ string temp = outputDir + "matrix.transpose.out";
+ m->openOutputFile(temp, out);
+ out << "matrix" << endl;
+
//fill matrix with shared files relative abundances
for (int i = 0; i < lookupFloat.size(); i++) {
for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
matrix[i].push_back(lookupFloat[i]->getAbundance(j));
+ out << lookupFloat[i]->getAbundance(j) << '\t';
}
+ out << endl;
}
-
+ out << endl << endl << "transpose" << endl;
vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
for (int i = 0; i < transposeMatrix.size(); i++) {
for (int j = 0; j < matrix.size(); j++) {
transposeMatrix[i].push_back(matrix[j][i]);
+ out << matrix[j][i] << '\t';
}
+ out << endl;
}
- matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
+ //matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
+
+ out << endl << endl << "matrix mult" << endl;
+ for (int i = 0; i < matrix.size(); i++) {
+ for (int j = 0; j < matrix[i].size(); j++) {
+ out << matrix[i][j] << '\t';
+ }
+ out << endl;
+ }
+ out.close();
+
double offset = 0.0000;
vector<double> d;
vector<double> e;
vector<vector<double> > G = matrix;
- vector<vector<double> > copy_G;
+ //vector<vector<double> > copy_G;
for(int count=0;count<2;count++){
- linearCalc.tred2(G, d, e); if (m->control_pressed) { return 0; }
- linearCalc.qtli(d, e, G); if (m->control_pressed) { return 0; }
+ linearCalc.recenter(offset, matrix, G); if (m->control_pressed) { return 0; }
+ linearCalc.tred2(G, d, e); if (m->control_pressed) { return 0; }
+ linearCalc.qtli(d, e, G); if (m->control_pressed) { return 0; }
offset = d[d.size()-1];
if(offset > 0.0) break;
}
for (int i = 1; i < 4; i++) {
- vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pcoa file
+ vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pca file
if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
- double corr = linearCalc.calcPearson(EuclidDists, matrix); //G is the pcoa file, D is the users distance matrix
+ double corr = linearCalc.calcPearson(EuclidDists, matrix); //G is the pca file, D is the users distance matrix
m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
+ m->mothurOut("Rsq " + toString(i) + " axis: " + toString(corr * corr)); m->mothurOutEndLine();
+
if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
}
}
}
}
- ofstream pcaData((fnameRoot+".pca").c_str(), ios::trunc);
+ ofstream pcaData((fnameRoot+".pca.axes").c_str(), ios::trunc);
pcaData.setf(ios::fixed, ios::floatfield);
pcaData.setf(ios::showpoint);
- outputNames.push_back(fnameRoot+".pca");
- outputTypes["pca"].push_back(fnameRoot+".pca");
+ outputNames.push_back(fnameRoot+".pca.axes");
+ outputTypes["pca"].push_back(fnameRoot+".pca.axes");
ofstream pcaLoadings((fnameRoot+".pca.loadings").c_str(), ios::trunc);
pcaLoadings.setf(ios::fixed, ios::floatfield);
projects / mothur.git / blobdiff