merge branch

Author: y3tseng

README.md

@@ -20,6 +20,8 @@
 
 ## <a name="intro"></a> Introduction
 
+**\**This experimental version is still in development and has not been fully debugged. It may contain bugs. Please refer to the main branch for a stable version.\****
+
 TWILIGHT (**T**all and **Wi**de A**lig**nments at **H**igh **T**hroughput) is a tool designed for ultrafast and ultralarge multiple sequence alignment. It is able to scale to millions of long nucleotide sequences (>10000 bases). TWILIGHT can run on CPU-only platforms (Linux/Mac) or take advantage of CUDA-capable GPUs for further acceleration.
 
 By default, TWILIGHT requires an unaligned sequence file in FASTA format and an input guide tree in Newick format to generate the output alignment in FASTA format (Fig. 1a). When a guide tree is unavailable, users can utilize the iterative mode, which provides a Snakemake workflow to estimate guide trees using external tools (Fig 1b.).

src/TALCO-XDrop.cpp

@@ -171,6 +171,17 @@ void Talco_xdrop::Traceback(
             ref_idx--;
         }
         aln.push_back(dir);   
+        if (firstTile &&  (ref_idx < 0 || query_idx < 0)) break;
+    }
+    if (firstTile) {
+        while (ref_idx > -1) {
+            aln.push_back(2);
+            ref_idx--;
+        }
+        while (query_idx > -1) {
+            aln.push_back(1);
+            query_idx--;
+        }
     }
 }
 
@@ -215,8 +226,8 @@ void Talco_xdrop::Tile (
         // int8_t tb_state_marker = 0;
 
         float denominator = refNum * qryNum;
-                    
 
+        
         int32_t L[3], U[3];
         int32_t *S[3], *I[2], *D[2];
         int32_t *CS[3], *CI[2], *CD[2];
@@ -330,10 +341,9 @@ void Talco_xdrop::Tile (
                 int32_t offsetUp = L[k%3]-L[(k+2)%3]+offset;
                 int32_t offsetLeft = L[k%3]-L[(k+2)%3]+offset-1;
 
-
 
                 if ((k==0) || 
-                    ((offsetDiag >= 0) && (offsetDiag <= U[(k+1)%3]-L[(k+1)%3])) || 
+                    ((offsetDiag >= 0) && (offsetDiag <= U[(k+1)%3]-L[(k+1)%3])) ||
                     (tile == 0 && (i == 0 || j == 0 ))) {
                     int32_t similarScore = 0;
                     float numerator = 0;
@@ -345,10 +355,10 @@ void Talco_xdrop::Tile (
                         }
                     }  
 
-                    similarScore = static_cast<int32_t>(std::nearbyint(numerator/denominator));
-                    if  (tile == 0 && (i == 0 || j == 0 )) match = similarScore + param.gapExtend * std::max(reference_idx + j, query_idx + i);
-                    else if (offsetDiag < 0) match = similarScore;
-                    else                match = S[(k+1)%3][offsetDiag] + similarScore;
+                    similarScore = static_cast<int32_t>(roundeven(numerator/denominator));
+                    if  (tile == 0 && (i == 0 || j == 0 )) match = similarScore + param.gapBoundary * std::max(reference_idx + j, query_idx + i);
+                    else if (offsetDiag < 0)               match = similarScore;
+                    else                                   match = S[(k+1)%3][offsetDiag] + similarScore;
                 }
 
                 int32_t pos_gapOpen_ref =   static_cast<int32_t>(std::nearbyint(gapOp[0][reference_idx+j]));
@@ -357,29 +367,24 @@ void Talco_xdrop::Tile (
                 int32_t pos_gapExtend_qry = static_cast<int32_t>(std::nearbyint(gapEx[1][query_idx+i]));
 
                 if (query_idx + i == query.size() - 1) {
-                    pos_gapOpen_ref = param.gapExtend;
-                    pos_gapExtend_ref = param.gapExtend;
+                    pos_gapOpen_ref = param.gapBoundary;
+                    pos_gapExtend_ref = param.gapBoundary;
                 }
                 if (reference_idx + j == reference.size() - 1) {
-                    pos_gapOpen_qry = param.gapExtend;
-                    pos_gapExtend_qry = param.gapExtend;
+                    pos_gapOpen_qry = param.gapBoundary;
+                    pos_gapExtend_qry = param.gapBoundary;
                 }
-                
+
                 if ((offsetUp >= 0) && (offsetUp <= U[(k+2)%3]-L[(k+2)%3])) {
-                    // delOp = S[(k+2)%3][offsetUp] + gapOpen;
                     delOp = S[(k+2)%3][offsetUp] + pos_gapOpen_ref;
-                    // delExt = D[(k+1)%2][offsetUp] + gapExtend;
                     delExt = D[(k+1)%2][offsetUp] + pos_gapExtend_ref;
                 }
 
                 if ((offsetLeft >= 0) && (offsetLeft <= U[(k+2)%3]-L[(k+2)%3])) {
-                    // insOp = S[(k+2)%3][offsetLeft] + gapOpen;
                     insOp = S[(k+2)%3][offsetLeft] + pos_gapOpen_qry;
-                    // insExt = I[(k+1)%2][offsetLeft] + gapExtend;
                     insExt = I[(k+1)%2][offsetLeft] + pos_gapExtend_qry;
                 }
 
-                
                 I[k%2][offset] = insOp;
                 D[k%2][offset] = delOp;
 
@@ -395,8 +400,8 @@ void Talco_xdrop::Tile (
                     Dptr = true;
                 }
 
-                if (match > I[k%2][offset]) {
-                    if (match > D[k%2][offset]) {
+                if (match >= I[k%2][offset]) {
+                    if (match >= D[k%2][offset]) {
                         S[k%3][offset] = match;
                         ptr = 0;
                     }
@@ -548,7 +553,7 @@ void Talco_xdrop::Tile (
                 tb_start_ftr = (tb_state == 3) ? ftr_length.size() - 2: ftr_length.size() - 1;
             }
         }
-
+            
         reference_idx += conv_ref_idx;
         query_idx += conv_query_idx;
 

src/TALCO-XDrop.hpp

@@ -40,11 +40,12 @@ namespace Talco_xdrop {
     struct Params {
         float scoreMatrix [5][5];
         float gapOpen;
-        float gapClose;
+        float gapBoundary;
         float gapExtend;
         int32_t xdrop;
         int32_t fLen;
         int32_t marker;
+        int type;
 
         void updateXDrop(int32_t new_xdrop) { 
             this->xdrop = new_xdrop;
@@ -59,7 +60,7 @@ namespace Talco_xdrop {
             }
             this->gapOpen = t_param[25];
             this->gapExtend = t_param[26];
-            this->gapClose = t_param[27];
+            this->gapBoundary = t_param[27];
             this->xdrop = static_cast<int32_t> (1000 * -1 * t_param[26]);
             this->fLen = (1 << 12);
             this->marker = (1 << 10); //reduce this value to save memory