Fixed msa arr generation in lib code

TimoLassmann · Feb 4, 2023 · 90747e3 · 90747e3
1 parent 166fc7b
commit 90747e3
Show file tree

Hide file tree

Showing 3 changed files with 98 additions and 25 deletions.
diff --git a/lib/src/aln_wrap.c b/lib/src/aln_wrap.c
@@ -23,14 +23,17 @@
 
 int kalign(char **seq, int *len, int numseq,int n_threads, int type, float gpo, float gpe, float tgpe, char ***aligned, int *out_aln_len)
 {
-        struct msa *msa;
+        struct msa *msa = NULL;
         RUN(kalign_arr_to_msa(seq, len,numseq, &msa));
+
         msa->quiet = 1;
         if(n_threads < 1){
                 n_threads = 1;
         }
         RUN(kalign_run(msa,n_threads, type,  gpo, gpe, tgpe));
+
         RUN(kalign_msa_to_arr(msa, aligned, out_aln_len));
+
         kalign_free_msa(msa);
 
         return OK;

diff --git a/lib/src/msa_op.c b/lib/src/msa_op.c
@@ -307,36 +307,59 @@ int convert_msa_to_internal(struct msa* msa, int type)
 
 int kalign_msa_to_arr(struct msa* msa, char ***aligned, int *out_aln_len)
 {
-        int aln_len = 0;
-        char** out = NULL;
+        ASSERT(msa != NULL,"No MSA!");
+        ASSERT(msa->aligned == ALN_STATUS_FINAL,"Sequences are not finalized");
 
+        char** out = NULL;
         int numseq = msa->numseq;
-        for (int j = 0; j <= msa->sequences[0]->len;j++){
-                aln_len += msa->sequences[0]->gaps[j];
+        int aln_len = msa->alnlen;
+        MMALLOC(out, sizeof(char*) * numseq);
+        for(int i = 0 ; i < numseq;i++){
+                out[i] = NULL;
+                MMALLOC(out[i], sizeof(char) * (aln_len +1));
         }
-        aln_len += msa->sequences[0]->len;
-        aln_len += 1;
 
-        MMALLOC(out, sizeof(char*) * numseq);
-        for(int i = 0; i < numseq; i++){
-                out[i] = 0;
-                MMALLOC(out[i], sizeof(char) * (uint64_t)aln_len);
-                int pos = 0;
-                for(int j = 0;j < msa->sequences[i]->len;j++){
-                        for(int c = 0;c < msa->sequences[i]->gaps[j];c++){
-                                out[i][pos] = '-';
-                                pos++;
-                        }
-                        out[i][pos] = msa->sequences[i]->seq[j];
-                        pos++;
+        /* galloc(&out, numseq,aln_len+1); */
+        for(int i = 0 ; i < numseq;i++){
+                for(int j = 0; j < aln_len;j++){
+                        out[i][j] = msa->sequences[i]->seq[j];
                 }
-                for(int c = 0;c < msa->sequences[i]->gaps[ msa->sequences[i]->len];c++){
-                        out[i][pos] = '-';
-                        pos++;
-                }
-                out[i][pos] = 0;
+                out[i][aln_len] = 0;
+                /* fprintf(stdout,"IN msa to ALIGNED: %s    %d\n", out[i], msa->alnlen); */
         }
 
+        /* /\* msa->alnlen = aln_len; *\/ */
+        /* /\* msa->aligned = ALN_STATUS_FINAL; *\/ */
+        /* /\* int aln_len = 0; *\/ */
+
+
+
+        /* for (int j = 0; j <= msa->sequences[0]->len;j++){ */
+        /*         aln_len += msa->sequences[0]->gaps[j]; */
+        /* } */
+        /* aln_len += msa->sequences[0]->len; */
+        /* aln_len += 1; */
+
+        /* MMALLOC(out, sizeof(char*) * numseq); */
+        /* for(int i = 0; i < numseq; i++){ */
+        /*         out[i] = 0; */
+        /*         MMALLOC(out[i], sizeof(char) * (uint64_t)aln_len); */
+        /*         int pos = 0; */
+        /*         for(int j = 0;j < msa->sequences[i]->len;j++){ */
+        /*                 for(int c = 0;c < msa->sequences[i]->gaps[j];c++){ */
+        /*                         out[i][pos] = '-'; */
+        /*                         pos++; */
+        /*                 } */
+        /*                 out[i][pos] = msa->sequences[i]->seq[j]; */
+        /*                 pos++; */
+        /*         } */
+        /*         for(int c = 0;c < msa->sequences[i]->gaps[ msa->sequences[i]->len];c++){ */
+        /*                 out[i][pos] = '-'; */
+        /*                 pos++; */
+        /*         } */
+        /*         out[i][pos] = 0; */
+        /* } */
+
         *aligned = out;
         *out_aln_len = aln_len;
         return OK;

diff --git a/tests/kalign_lib_test.c b/tests/kalign_lib_test.c
@@ -2,7 +2,7 @@
 #include <stdlib.h>
 
 #include "kalign/kalign.h"
-
+int test_2(void);
 int main(int argc, char *argv[])
 {
         struct msa* msa = NULL;
@@ -23,5 +23,52 @@ int main(int argc, char *argv[])
         kalign_write_msa(msa, "test.fa", "fasta");
         /* cleaning up */
         kalign_free_msa(msa);
+        test_2();
         return EXIT_SUCCESS;
 }
+
+int test_2(void){
+
+        // Initialize array
+        char* inseq[] = {
+                "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTTGCATGCATGCATGCATGCATGCATGCA",
+                "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTTGCATGCATGCATGCATGCATGCATGCA",
+                "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTTGCATGCATGCATGCATGCATGCATGCA",
+                "ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT"
+        };
+
+        int numseq = 4;
+        int L[] = {76, 76, 76, 48};
+
+        char** aln = NULL;
+        int aln_len = 0;
+
+        kalign(inseq, L, numseq, 4 , KALIGN_TYPE_DNA_INTERNAL, -1, -1 , -1, &aln, &aln_len);
+
+        printf("Aligned:\n");
+        for(int i = 0; i < numseq; i++){
+                printf("%s\n", aln[i]);
+                /* printf("%s\n", inseq[i]); */
+        }
+
+        for (int i = 0; i < numseq; ++i) {
+                int num_aligned = 0;
+                for (int j = 0; j < aln_len; ++j) {
+                        if (aln[i][j] != '-' && aln[i][j] != '\0') {
+                                num_aligned++;
+                        }
+                }
+                if (num_aligned != L[i]) {
+                        fprintf(stderr, "error: sequence %d is not fully aligned\n", i);
+                        fprintf(stderr, "%s\n", inseq[i]);
+                        fprintf(stderr, "%s\n", aln[i]);
+                }
+        }
+
+        /* Free alignment  */
+        for(int i = 0; i < numseq;i++){
+                free(aln[i]);
+        }
+        free(aln);
+        return 0;
+}