working on pam

[mothur.git] / trimseqscommand.h

diff --git a/trimseqscommand.h b/trimseqscommand.h
index 4d0a6112cb6e288600ae944e9340563f5973da47..14eed17e236a333583d461d61482a1d541c6a9cf 100644 (file)
--- a/trimseqscommand.h
+++ b/trimseqscommand.h
@@ -55,7 +55,7 @@ private:
        bool abort, createGroup;
        string fastaFile, oligoFile, qFileName, groupfile, nameFile, countfile, outputDir;
        
        bool abort, createGroup;
        string fastaFile, oligoFile, qFileName, groupfile, nameFile, countfile, outputDir;
        

-       bool flip, allFiles, qtrim, keepforward, pairedOligos, reorient;
+       bool flip, allFiles, qtrim, keepforward, pairedOligos, reorient, logtransform;

        int numFPrimers, numRPrimers, numLinkers, numSpacers, maxAmbig, maxHomoP, minLength, maxLength, processors, tdiffs, bdiffs, pdiffs, ldiffs, sdiffs, comboStarts;
        int qWindowSize, qWindowStep, keepFirst, removeLast;
        double qRollAverage, qThreshold, qWindowAverage, qAverage;
        int numFPrimers, numRPrimers, numLinkers, numSpacers, maxAmbig, maxHomoP, minLength, maxLength, processors, tdiffs, bdiffs, pdiffs, ldiffs, sdiffs, comboStarts;
        int qWindowSize, qWindowStep, keepFirst, removeLast;
        double qRollAverage, qThreshold, qWindowAverage, qAverage;


@@ -98,7 +98,7 @@ struct trimData {
     vector<vector<string> > qualFileNames;
     vector<vector<string> > nameFileNames;
     unsigned long long lineStart, lineEnd, qlineStart, qlineEnd;
     vector<vector<string> > qualFileNames;
     vector<vector<string> > nameFileNames;
     unsigned long long lineStart, lineEnd, qlineStart, qlineEnd;

-    bool flip, allFiles, qtrim, keepforward, createGroup, pairedOligos, reorient;
+    bool flip, allFiles, qtrim, keepforward, createGroup, pairedOligos, reorient, logtransform;

        int numFPrimers, numRPrimers, numLinkers, numSpacers, maxAmbig, maxHomoP, minLength, maxLength, tdiffs, bdiffs, pdiffs, ldiffs, sdiffs;
        int qWindowSize, qWindowStep, keepFirst, removeLast, count;
        double qRollAverage, qThreshold, qWindowAverage, qAverage;
        int numFPrimers, numRPrimers, numLinkers, numSpacers, maxAmbig, maxHomoP, minLength, maxLength, tdiffs, bdiffs, pdiffs, ldiffs, sdiffs;
        int qWindowSize, qWindowStep, keepFirst, removeLast, count;
        double qRollAverage, qThreshold, qWindowAverage, qAverage;


@@ -121,7 +121,7 @@ struct trimData {
        trimData(string fn, string qn, string nf, string cf, string tn, string sn, string tqn, string sqn, string tnn, string snn, string tcn, string scn,string gn, vector<vector<string> > ffn, vector<vector<string> > qfn, vector<vector<string> > nfn, unsigned long long lstart, unsigned long long lend, unsigned long long qstart, unsigned long long qend,  MothurOut* mout,
                       int pd, int bd, int ld, int sd, int td, map<string, int> pri, map<string, int> bar, vector<string> revP, vector<string> li, vector<string> spa, map<int, oligosPair> pbr, map<int, oligosPair> ppr, bool po,
                       vector<string> priNameVector, vector<string> barNameVector, bool cGroup, bool aFiles, bool keepF, int keepfi, int removeL,
        trimData(string fn, string qn, string nf, string cf, string tn, string sn, string tqn, string sqn, string tnn, string snn, string tcn, string scn,string gn, vector<vector<string> > ffn, vector<vector<string> > qfn, vector<vector<string> > nfn, unsigned long long lstart, unsigned long long lend, unsigned long long qstart, unsigned long long qend,  MothurOut* mout,
                       int pd, int bd, int ld, int sd, int td, map<string, int> pri, map<string, int> bar, vector<string> revP, vector<string> li, vector<string> spa, map<int, oligosPair> pbr, map<int, oligosPair> ppr, bool po,
                       vector<string> priNameVector, vector<string> barNameVector, bool cGroup, bool aFiles, bool keepF, int keepfi, int removeL,

-                      int WindowStep, int WindowSize, int WindowAverage, bool trim, double Threshold, double Average, double RollAverage,
+                      int WindowStep, int WindowSize, int WindowAverage, bool trim, double Threshold, double Average, double RollAverage, bool lt,

                       int minL, int maxA, int maxH, int maxL, bool fli, bool reo, map<string, string> nm, map<string, int> ncount) {
         filename = fn;
         qFileName = qn;
                       int minL, int maxA, int maxH, int maxL, bool fli, bool reo, map<string, string> nm, map<string, int> ncount) {
         filename = fn;
         qFileName = qn;


@@ -174,6 +174,7 @@ struct trimData {
         qThreshold = Threshold;
         qAverage = Average;
         qRollAverage = RollAverage;
         qThreshold = Threshold;
         qAverage = Average;
         qRollAverage = RollAverage;

+        logtransform = lt;

         minLength = minL;
         maxAmbig = maxA;
         maxHomoP = maxH;
         minLength = minL;
         maxAmbig = maxA;
         maxHomoP = maxH;


@@ -261,7 +262,7 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                
         TrimOligos* trimOligos = NULL;
         int numBarcodes = pDataArray->barcodes.size();
                
         TrimOligos* trimOligos = NULL;
         int numBarcodes = pDataArray->barcodes.size();

-        if (pDataArray->pairedOligos)   {   trimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, 0, 0, pDataArray->pairedPrimers, pDataArray->pairedBarcodes);   numBarcodes = pDataArray->pairedBarcodes.size(); }
+        if (pDataArray->pairedOligos)   {   trimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, 0, 0, pDataArray->pairedPrimers, pDataArray->pairedBarcodes);   numBarcodes = pDataArray->pairedBarcodes.size(); pDataArray->numFPrimers = pDataArray->pairedPrimers.size(); }

         else                {   trimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, pDataArray->ldiffs, pDataArray->sdiffs, pDataArray->primers, pDataArray->barcodes, pDataArray->revPrimer, pDataArray->linker, pDataArray->spacer);  }
         
         TrimOligos* rtrimOligos = NULL;
         else                {   trimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, pDataArray->ldiffs, pDataArray->sdiffs, pDataArray->primers, pDataArray->barcodes, pDataArray->revPrimer, pDataArray->linker, pDataArray->spacer);  }
         
         TrimOligos* rtrimOligos = NULL;


@@ -276,6 +277,19 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                 oligosPair tempPair(trimOligos->reverseOligo((it->second).reverse), (trimOligos->reverseOligo((it->second).forward))); //reverseBarcode, rc ForwardBarcode
                 rpairedBarcodes[it->first] = tempPair;
             }
                 oligosPair tempPair(trimOligos->reverseOligo((it->second).reverse), (trimOligos->reverseOligo((it->second).forward))); //reverseBarcode, rc ForwardBarcode
                 rpairedBarcodes[it->first] = tempPair;
             }

+            
+            int index = rpairedBarcodes.size();
+            for (map<string, int>::iterator it = pDataArray->barcodes.begin(); it != pDataArray->barcodes.end(); it++) {
+                oligosPair tempPair("", trimOligos->reverseOligo((it->first))); //reverseBarcode, rc ForwardBarcode
+                rpairedBarcodes[index] = tempPair; index++;
+            }
+            
+            index = rpairedPrimers.size();
+            for (map<string, int>::iterator it = pDataArray->primers.begin(); it != pDataArray->primers.end(); it++) {
+                oligosPair tempPair("", trimOligos->reverseOligo((it->first))); //reverseBarcode, rc ForwardBarcode
+                rpairedPrimers[index] = tempPair; index++;
+            }
+

             rtrimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, 0, 0, rpairedPrimers, rpairedBarcodes); numBarcodes = rpairedBarcodes.size();
         }
         
             rtrimOligos = new TrimOligos(pDataArray->pdiffs, pDataArray->bdiffs, 0, 0, rpairedPrimers, rpairedBarcodes); numBarcodes = rpairedBarcodes.size();
         }
         


@@ -283,7 +297,7 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                for(int i = 0; i < pDataArray->lineEnd; i++){ //end is the number of sequences to process
                                   
                        if (pDataArray->m->control_pressed) {
                for(int i = 0; i < pDataArray->lineEnd; i++){ //end is the number of sequences to process
                                   
                        if (pDataArray->m->control_pressed) {

-                delete trimOligos;
+                delete trimOligos; if (pDataArray->reorient) { delete rtrimOligos; }

                                inFASTA.close(); trimFASTAFile.close(); scrapFASTAFile.close();
                                if ((pDataArray->createGroup) && (pDataArray->countfile == "")) {        outGroupsFile.close();   }
                 if(pDataArray->qFileName != "")        {       qFile.close();  scrapQualFile.close(); trimQualFile.close();    }
                                inFASTA.close(); trimFASTAFile.close(); scrapFASTAFile.close();
                                if ((pDataArray->createGroup) && (pDataArray->countfile == "")) {        outGroupsFile.close();   }
                 if(pDataArray->qFileName != "")        {       qFile.close();  scrapQualFile.close(); trimQualFile.close();    }


@@ -299,11 +313,14 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                        int currentSeqsDiffs = 0;
             
                        Sequence currSeq(inFASTA); pDataArray->m->gobble(inFASTA);
                        int currentSeqsDiffs = 0;
             
                        Sequence currSeq(inFASTA); pDataArray->m->gobble(inFASTA);

+            Sequence savedSeq(currSeq.getName(), currSeq.getAligned());

                        
                        

-                       QualityScores currQual;
+                       QualityScores currQual; QualityScores savedQual;

                        if(pDataArray->qFileName != ""){
                                currQual = QualityScores(qFile);  pDataArray->m->gobble(qFile);
                        if(pDataArray->qFileName != ""){
                                currQual = QualityScores(qFile);  pDataArray->m->gobble(qFile);

+                savedQual.setName(currQual.getName()); savedQual.setScores(currQual.getScores());

                        }
                        }

+              

                        
                        string origSeq = currSeq.getUnaligned();
                        if (origSeq != "") {
                        
                        string origSeq = currSeq.getUnaligned();
                        if (origSeq != "") {


@@ -353,13 +370,13 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                     int thisPrimerIndex = 0;
                     
                     if(numBarcodes != 0){
                     int thisPrimerIndex = 0;
                     
                     if(numBarcodes != 0){

-                        thisSuccess = rtrimOligos->stripBarcode(currSeq, currQual, thisBarcodeIndex);
+                        thisSuccess = rtrimOligos->stripBarcode(savedSeq, savedQual, thisBarcodeIndex);

                         if(thisSuccess > pDataArray->bdiffs)           {       thisTrashCode += 'b';   }
                         else{ thisCurrentSeqsDiffs += thisSuccess;  }
                     }
                     
                     if(pDataArray->numFPrimers != 0){
                         if(thisSuccess > pDataArray->bdiffs)           {       thisTrashCode += 'b';   }
                         else{ thisCurrentSeqsDiffs += thisSuccess;  }
                     }
                     
                     if(pDataArray->numFPrimers != 0){

-                        thisSuccess = rtrimOligos->stripForward(currSeq, currQual, thisPrimerIndex, pDataArray->keepforward);
+                        thisSuccess = rtrimOligos->stripForward(savedSeq, savedQual, thisPrimerIndex, pDataArray->keepforward);

                         if(thisSuccess > pDataArray->pdiffs)           {       thisTrashCode += 'f';   }
                         else{ thisCurrentSeqsDiffs += thisSuccess;  }
                     }
                         if(thisSuccess > pDataArray->pdiffs)           {       thisTrashCode += 'f';   }
                         else{ thisCurrentSeqsDiffs += thisSuccess;  }
                     }


@@ -372,7 +389,13 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                         currentSeqsDiffs = thisCurrentSeqsDiffs;
                         barcodeIndex = thisBarcodeIndex;
                         primerIndex = thisPrimerIndex;
                         currentSeqsDiffs = thisCurrentSeqsDiffs;
                         barcodeIndex = thisBarcodeIndex;
                         primerIndex = thisPrimerIndex;

-                    }
+                        savedSeq.reverseComplement();
+                        currSeq.setAligned(savedSeq.getAligned());
+                        if(pDataArray->qFileName != ""){
+                            savedQual.flipQScores();
+                            currQual.setScores(savedQual.getScores());
+                        }
+                    }else { trashCode += "(" + thisTrashCode + ")";  }

                 }
 
                 
                 }
 
                 


@@ -407,9 +430,9 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                                        int origLength = currSeq.getNumBases();
                                        
                                        if(pDataArray->qThreshold != 0)                 {       success = currQual.stripQualThreshold(currSeq, pDataArray->qThreshold);                 }
                                        int origLength = currSeq.getNumBases();
                                        
                                        if(pDataArray->qThreshold != 0)                 {       success = currQual.stripQualThreshold(currSeq, pDataArray->qThreshold);                 }

-                                       else if(pDataArray->qAverage != 0)              {       success = currQual.cullQualAverage(currSeq, pDataArray->qAverage);                              }
-                                       else if(pDataArray->qRollAverage != 0)  {       success = currQual.stripQualRollingAverage(currSeq, pDataArray->qRollAverage);  }
-                                       else if(pDataArray->qWindowAverage != 0){       success = currQual.stripQualWindowAverage(currSeq, pDataArray->qWindowStep, pDataArray->qWindowSize, pDataArray->qWindowAverage);       }
+                                       else if(pDataArray->qAverage != 0)              {       success = currQual.cullQualAverage(currSeq, pDataArray->qAverage, pDataArray->logtransform);                            }
+                                       else if(pDataArray->qRollAverage != 0)  {       success = currQual.stripQualRollingAverage(currSeq, pDataArray->qRollAverage, pDataArray->logtransform);        }
+                                       else if(pDataArray->qWindowAverage != 0){       success = currQual.stripQualWindowAverage(currSeq, pDataArray->qWindowStep, pDataArray->qWindowSize, pDataArray->qWindowAverage, pDataArray->logtransform);     }

                                        else                                            {       success = 1;                            }
                                        
                                        //you don't want to trim, if it fails above then scrap it
                                        else                                            {       success = 1;                            }
                                        
                                        //you don't want to trim, if it fails above then scrap it


@@ -567,12 +590,13 @@ static DWORD WINAPI MyTrimThreadFunction(LPVOID lpParam){
                        }
                        
                        //report progress
                        }
                        
                        //report progress

-                       if((i) % 1000 == 0){    pDataArray->m->mothurOut(toString(i)); pDataArray->m->mothurOutEndLine();               }
+                       if((pDataArray->count) % 1000 == 0){    pDataArray->m->mothurOut(toString(pDataArray->count)); pDataArray->m->mothurOutEndLine();               }

                        
                }
                //report progress
                if((pDataArray->count) % 1000 != 0){    pDataArray->m->mothurOut(toString(pDataArray->count)); pDataArray->m->mothurOutEndLine();               }
                
                        
                }
                //report progress
                if((pDataArray->count) % 1000 != 0){    pDataArray->m->mothurOut(toString(pDataArray->count)); pDataArray->m->mothurOutEndLine();               }
                

+        if (pDataArray->reorient) { delete rtrimOligos; }

                delete trimOligos;
                inFASTA.close();
                trimFASTAFile.close();
                delete trimOligos;
                inFASTA.close();
                trimFASTAFile.close();