]> git.donarmstrong.com Git - mothur.git/blobdiff - chimera.cpp
working on pam
[mothur.git] / chimera.cpp
index 83d40fe8342558eacb598525862364dba0f7e5ca..53e77322f52ba40eb478e0b385a84982de752d6b 100644 (file)
@@ -8,6 +8,7 @@
  */
 
 #include "chimera.h"
+#include "referencedb.h"
 
 //***************************************************************************************************************
 //this is a vertical soft filter
@@ -31,6 +32,8 @@ string Chimera::createFilter(vector<Sequence*> seqs, float t) {
                        if (m->control_pressed) { return filterString; }
                
                        string seqAligned = seqs[i]->getAligned();
+                       
+                       if (seqAligned.length() != filterString.length()) {  m->mothurOut(seqs[i]->getName() + " is not the same length as the template sequences. Aborting!\n");  exit(1); }
                
                        for (int j = 0; j < seqAligned.length(); j++) {
                                //if this spot is a gap
@@ -92,89 +95,163 @@ map<int, int> Chimera::runFilter(Sequence* seq) {
 //***************************************************************************************************************
 vector<Sequence*> Chimera::readSeqs(string file) {
        try {
-       
+               
                vector<Sequence*> container;
                int count = 0;
                length = 0;
                unaligned = false;
-
-               m->mothurOut("Reading sequences from " + file + "..."); cout.flush();
+               ReferenceDB* rdb = ReferenceDB::getInstance();
                
-               #ifdef USE_MPI  
-                       int pid;
-                       vector<long> positions;
-                       int numSeqs;
-               
-                       MPI_Status status; 
-                       MPI_File inMPI;
-                       MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
-
-                       char* inFileName = new char[file.length()];
-                       memcpy(inFileName, file.c_str(), file.length());
-       
-                       MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
-                       delete inFileName;
-
-                       if (pid == 0) {
-                               positions = setFilePosFasta(file, numSeqs); //fills MPIPos, returns numSeqs
-
-                               //send file positions to all processes
-                               MPI_Bcast(&numSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
-                               MPI_Bcast(&positions[0], (numSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos     
-                       }else{
-                               MPI_Bcast(&numSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
-                               positions.resize(numSeqs+1);
-                               MPI_Bcast(&positions[0], (numSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
-                       }
+               if (file == "saved") {
                        
-                       //read file 
-                       for(int i=0;i<numSeqs;i++){
                        
-                               if (m->control_pressed) { MPI_File_close(&inMPI); return container; }
-       
-                               //read next sequence
-                               int seqlength = positions[i+1] - positions[i];
-                               char* buf4 = new char[seqlength];
-
-                               MPI_File_read_at(inMPI, positions[i], buf4, seqlength, MPI_CHAR, &status);
+                       m->mothurOutEndLine();  m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory.");        m->mothurOutEndLine();
+                       
+                       for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
+                               Sequence* temp = new Sequence(rdb->referenceSeqs[i].getName(), rdb->referenceSeqs[i].getAligned());
+                               
+                               if (count == 0) {  length = temp->getAligned().length();  count++;  } //gets first seqs length
+                               else if (length != temp->getAligned().length()) {       unaligned = true;       }
                                
-                               string tempBuf = buf4;
-                               if (tempBuf.length() > seqlength) { tempBuf = tempBuf.substr(0, seqlength); }
-                               delete buf4;
+                               if (temp->getName() != "") {  container.push_back(temp);  }
+                       }
+                       
+                       templateFileName = rdb->getSavedReference();
+                       
+               }else {
+                       
+                       m->mothurOut("Reading sequences from " + file + "..."); cout.flush();
+                       
+                       #ifdef USE_MPI
+                               int pid, processors;
+                               vector<unsigned long long> positions;
+                               int numSeqs;
+                               int tag = 2001;
+                               MPI_File inMPI;
+                               MPI_Status status; 
+                       
+                               if (byGroup) {
+                                       char inFileName[1024];
+                                       strcpy(inFileName, file.c_str());
+                                       
+                                       MPI_File_open(MPI_COMM_SELF, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
+                                       
+                                       positions = m->setFilePosFasta(file, numSeqs); //fills MPIPos, returns numSeqs
+                                       
+                                       //read file 
+                                       for(int i=0;i<numSeqs;i++){
+                                               
+                                               if (m->control_pressed) { MPI_File_close(&inMPI); return container; }
+                                               
+                                               //read next sequence
+                                               int seqlength = positions[i+1] - positions[i];
+                                               char* buf4 = new char[seqlength];
+                                               
+                                               MPI_File_read_at(inMPI, positions[i], buf4, seqlength, MPI_CHAR, &status);
+                                               
+                                               string tempBuf = buf4;
+                                               if (tempBuf.length() > seqlength) { tempBuf = tempBuf.substr(0, seqlength); }
+                                               delete buf4;
+                                               
+                                               istringstream iss (tempBuf,istringstream::in);
+                                               
+                                               Sequence* current = new Sequence(iss);   
+                                               if (current->getName() != "") {
+                                                       if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
+                                                       else if (length != current->getAligned().length()) {    unaligned = true;       }
+                                                       
+                                                       container.push_back(current);  
+                                                       if (rdb->save) { rdb->referenceSeqs.push_back(*current); }
+                                               }
+                                       }
+                                       
+                                       MPI_File_close(&inMPI);
+                                       
+                               }else {                                 
+                                       
+                                       MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
+                                       MPI_Comm_size(MPI_COMM_WORLD, &processors);
+                                       
+                                       //char* inFileName = new char[file.length()];
+                                       //memcpy(inFileName, file.c_str(), file.length());
+                                       
+                                       char inFileName[1024];
+                                       strcpy(inFileName, file.c_str());
+                                       
+                                       MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
+                                       //delete inFileName;
+                                       
+                                       if (pid == 0) {
+                                               positions = m->setFilePosFasta(file, numSeqs); //fills MPIPos, returns numSeqs
+                                               
+                                               //send file positions to all processes
+                                               for(int i = 1; i < processors; i++) { 
+                                                       MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+                                                       MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
+                                               }
+                                       }else{
+                                               MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
+                                               positions.resize(numSeqs+1);
+                                               MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
+                                       }
+                                       
+                                       //read file 
+                                       for(int i=0;i<numSeqs;i++){
+                                               
+                                               if (m->control_pressed) { MPI_File_close(&inMPI); return container; }
+                                               
+                                               //read next sequence
+                                               int seqlength = positions[i+1] - positions[i];
+                                               char* buf4 = new char[seqlength];
+                                               
+                                               MPI_File_read_at(inMPI, positions[i], buf4, seqlength, MPI_CHAR, &status);
+                                               
+                                               string tempBuf = buf4;
+                                               if (tempBuf.length() > seqlength) { tempBuf = tempBuf.substr(0, seqlength); }
+                                               delete buf4;
+                                               
+                                               istringstream iss (tempBuf,istringstream::in);
+                                               
+                                               Sequence* current = new Sequence(iss);   
+                                               if (current->getName() != "") {
+                                                       if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
+                                                       else if (length != current->getAligned().length()) {    unaligned = true;       }
+                                                       
+                                                       container.push_back(current);  
+                                                       if (rdb->save) { rdb->referenceSeqs.push_back(*current); }
+                                               }
+                                       }
+                                       
+                                       MPI_File_close(&inMPI);
+                                       MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+                               }
+               #else
 
-                               istringstream iss (tempBuf,istringstream::in);
-               
-                               Sequence* current = new Sequence(iss);   
-                               if (current->getName() != "") {
-                                       if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
-                                       else if (length != current->getAligned().length()) {    unaligned = true;       }
+                       ifstream in;
+                       m->openInputFile(file, in);
                        
+                       //read in seqs and store in vector
+                       while(!in.eof()){
+                               
+                               if (m->control_pressed) { return container; }
+                               
+                               Sequence* current = new Sequence(in);  m->gobble(in);
+                               
+                               if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
+                               else if (length != current->getAligned().length()) {   unaligned = true;        }
+                                                       
+                               if (current->getName() != "") {  
                                        container.push_back(current);  
+                                       if (rdb->save) { rdb->referenceSeqs.push_back(*current); }
                                }
                        }
-                       
-                       MPI_File_close(&inMPI);
-       #else
-
-               ifstream in;
-               openInputFile(file, in);
+                       in.close();
+               #endif
                
-               //read in seqs and store in vector
-               while(!in.eof()){
+                       m->mothurOut("Done."); m->mothurOutEndLine();
                        
-                       if (m->control_pressed) { return container; }
-                       
-                       Sequence* current = new Sequence(in);  gobble(in);
-                       
-                       if (count == 0) {  length = current->getAligned().length();  count++;  } //gets first seqs length
-                       else if (length != current->getAligned().length()) {    unaligned = true;       }
-                                               
-                       if (current->getName() != "") {  container.push_back(current);  }
+                       filterString = (string(container[0]->getAligned().length(), '1'));
                }
-               in.close();
-       #endif
-       
-               m->mothurOut("Done."); m->mothurOutEndLine();
                
                return container;
        }
@@ -199,13 +276,16 @@ void Chimera::setMask(string filename) {
                        MPI_Offset size;
                        MPI_Status status;
                        
-                       char* inFileName = new char[filename.length()];
-                       memcpy(inFileName, filename.c_str(), filename.length());
+                       //char* inFileName = new char[filename.length()];
+                       //memcpy(inFileName, filename.c_str(), filename.length());
+                       
+                       char inFileName[1024];
+                       strcpy(inFileName, filename.c_str());
        
                        MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
                        MPI_File_get_size(inMPI, &size);
 
-                       delete inFileName;
+                       //delete inFileName;
                        
                        char* buffer = new char[size];
                        MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
@@ -228,7 +308,7 @@ void Chimera::setMask(string filename) {
        #else
        
                        ifstream infile;
-                       openInputFile(filename, infile);
+                       m->openInputFile(filename, infile);
                        
                        if (!infile.eof()) {
                                Sequence temp(infile);