added mantel command

[mothur.git] / chimeraslayer.cpp

diff --git a/chimeraslayer.cpp b/chimeraslayer.cpp
index ee449cc9fd30f7ec521a03311473a805739159d1..38f7abe49bfc7266e6eb4acd123ff057b8780c42 100644 (file)
--- a/chimeraslayer.cpp
+++ b/chimeraslayer.cpp
@@ -13,7 +13,7 @@
 #include "blastdb.hpp"
 
 //***************************************************************************************************************
-ChimeraSlayer::ChimeraSlayer(string file, string temp, string mode, int k, int ms, int mms, int win, float div, 
+ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, string mode, int k, int ms, int mms, int win, float div, 
 int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera()  {      
        try {
                fastafile = file;
@@ -33,6 +33,7 @@ int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int num
                increment = inc;
                numWanted = numw;
                realign = r; 
+               trimChimera = trim;
        
                decalc = new DeCalculator();    
                
@@ -44,7 +45,7 @@ int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int num
        }
 }
 //***************************************************************************************************************
-ChimeraSlayer::ChimeraSlayer(string file, string temp, string name, string mode, int k, int ms, int mms, int win, float div, 
+ChimeraSlayer::ChimeraSlayer(string file, string temp, bool trim, string name, string mode, string abunds, int k, int ms, int mms, int win, float div, 
                                                         int minsim, int mincov, int minbs, int minsnp, int par, int it, int inc, int numw, bool r) : Chimera()  {      
        try {
                fastafile = file; templateSeqs = readSeqs(fastafile);
@@ -64,6 +65,8 @@ ChimeraSlayer::ChimeraSlayer(string file, string temp, string name, string mode,
                increment = inc;
                numWanted = numw;
                realign = r; 
+               includeAbunds = abunds;
+               trimChimera = trim;
                
                //read name file and create nameMapRank
                readNameFile(name);
@@ -283,10 +286,26 @@ vector<Sequence*> ChimeraSlayer::getTemplate(Sequence* q) {
                //create list of names we want to put into the template
                set<string> namesToAdd;
                for (itRank = nameMapRank.begin(); itRank != nameMapRank.end(); itRank++) {
-                       if ((itRank->second).size() > thisRank) {
-                               //you are more abundant than me
-                               for (int i = 0; i < (itRank->second).size(); i++) {
-                                       namesToAdd.insert((itRank->second)[i]);
+                       if (itRank->first != thisName) {
+                               if (includeAbunds == "greaterequal") {
+                                       if ((itRank->second).size() >= thisRank) {
+                                               //you are more abundant than me or equal to my abundance
+                                               for (int i = 0; i < (itRank->second).size(); i++) {
+                                                       namesToAdd.insert((itRank->second)[i]);
+                                               }
+                                       }
+                               }else if (includeAbunds == "greater") {
+                                       if ((itRank->second).size() > thisRank) {
+                                               //you are more abundant than me
+                                               for (int i = 0; i < (itRank->second).size(); i++) {
+                                                       namesToAdd.insert((itRank->second)[i]);
+                                               }
+                                       }
+                               }else if (includeAbunds == "all") {
+                                       //add everyone
+                                       for (int i = 0; i < (itRank->second).size(); i++) {
+                                               namesToAdd.insert((itRank->second)[i]);
+                                       }
                                }
                        }
                }
@@ -397,8 +416,11 @@ void ChimeraSlayer::printHeader(ostream& out) {
        out << "Name\tLeftParent\tRightParent\tDivQLAQRB\tPerIDQLAQRB\tBootStrapA\tDivQLBQRA\tPerIDQLBQRA\tBootStrapB\tFlag\tLeftWindow\tRightWindow\n";
 }
 //***************************************************************************************************************
-int ChimeraSlayer::print(ostream& out, ostream& outAcc) {
+Sequence* ChimeraSlayer::print(ostream& out, ostream& outAcc) {
        try {
+               Sequence* trim = NULL;
+               if (trimChimera) { trim = trimQuery; }
+               
                if (chimeraFlags == "yes") {
                        string chimeraFlag = "no";
                        if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
@@ -410,6 +432,21 @@ int ChimeraSlayer::print(ostream& out, ostream& outAcc) {
                                if ((chimeraResults[0].bsa >= minBS) || (chimeraResults[0].bsb >= minBS)) {
                                        m->mothurOut(querySeq->getName() + "\tyes"); m->mothurOutEndLine();
                                        outAcc << querySeq->getName() << endl;
+                                       
+                                       if (trimChimera) {  
+                                               int lengthLeft = spotMap[chimeraResults[0].winLEnd] - spotMap[chimeraResults[0].winLStart];
+                                               int lengthRight = spotMap[chimeraResults[0].winREnd] - spotMap[chimeraResults[0].winRStart];
+                                               
+                                               string newAligned = trim->getAligned();
+
+                                               if (lengthLeft > lengthRight) { //trim right
+                                                       for (int i = (spotMap[chimeraResults[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
+                                               }else { //trim left
+                                                       for (int i = 0; i < spotMap[chimeraResults[0].winLEnd]; i++) { newAligned[i] = '.'; }
+                                               }
+                                               trim->setAligned(newAligned);
+                                       }
+                                               
                                }
                        }
                        
@@ -417,7 +454,7 @@ int ChimeraSlayer::print(ostream& out, ostream& outAcc) {
                        out << endl;
                }else {  out << querySeq->getName() << "\tno" << endl;  }
                
-               return 0;
+               return trim;
                
        }
        catch(exception& e) {
@@ -427,13 +464,16 @@ int ChimeraSlayer::print(ostream& out, ostream& outAcc) {
 }
 #ifdef USE_MPI
 //***************************************************************************************************************
-int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
+Sequence* ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
        try {
                MPI_Status status;
                bool results = false;
                string outAccString = "";
                string outputString = "";
                
+               Sequence* trim = NULL;
+               if (trimChimera) { trim = trimQuery; }
+               
                if (chimeraFlags == "yes") {
                        string chimeraFlag = "no";
                        if(  (chimeraResults[0].bsa >= minBS && chimeraResults[0].divr_qla_qrb >= divR)
@@ -454,6 +494,19 @@ int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
                                
                                        MPI_File_write_shared(outAcc, buf2, length, MPI_CHAR, &status);
                                        delete buf2;
+                                       
+                                       if (trimChimera) {  
+                                               int lengthLeft = spotMap[chimeraResults[0].winLEnd] - spotMap[chimeraResults[0].winLStart];
+                                               int lengthRight = spotMap[chimeraResults[0].winREnd] - spotMap[chimeraResults[0].winRStart];
+                                               
+                                               string newAligned = trim->getAligned();
+                                               if (lengthLeft > lengthRight) { //trim right
+                                                       for (int i = (spotMap[chimeraResults[0].winRStart]-1); i < newAligned.length(); i++) { newAligned[i] = '.'; }
+                                               }else { //trim left
+                                                       for (int i = 0; i < (spotMap[chimeraResults[0].winLEnd]-1); i++) { newAligned[i] = '.'; }
+                                               }
+                                               trim->setAligned(newAligned);   
+                                       }
                                }
                        }
                        
@@ -481,7 +534,7 @@ int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
                }
                
                
-               return results;
+               return trim;
        }
        catch(exception& e) {
                m->errorOut(e, "ChimeraSlayer", "print");
@@ -493,6 +546,8 @@ int ChimeraSlayer::print(MPI_File& out, MPI_File& outAcc) {
 //***************************************************************************************************************
 int ChimeraSlayer::getChimeras(Sequence* query) {
        try {
+               if (trimChimera) { trimQuery = new Sequence(query->getName(), query->getAligned());  }
+               
                chimeraFlags = "no";
 
                //filter query