build_ref_db: fixed bug erasing some of the higher LCAs (i.e. lowest

similarities)

python/obitools3/commands/annotate.pyx

@@ -190,58 +190,50 @@ def sequenceTaggerGenerator(config, taxo=None):
             seq['seq_rank']=counter[0]
     
         for i,v in toSet:
-            #try:
-            if taxo is not None:
-                environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
-            else:
-                environ = {'sequence':seq, 'counter':counter[0], 'math':math}
-            val = eval(v, environ, seq)
-            #except Exception,e:       # TODO discuss usefulness of this
-            #    if options.onlyValid:
-            #        raise e
-            #    val = v
+            try:
+                if taxo is not None:
+                    environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
+                else:
+                    environ = {'sequence':seq, 'counter':counter[0], 'math':math}
+                val = eval(v, environ, seq)
+            except Exception:       # set string if not a valid expression
+                val = v
             seq[i]=val
 
         if length:
             seq['seq_length']=len(seq)
 
         if newId is not None:
-#            try:
-            if taxo is not None:
-                environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
-            else:
-                environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
-            val = eval(newId, environ, seq)
-#            except Exception,e:
-#                if options.onlyValid:
-#                    raise e
-#                val = newId
+            try:
+                if taxo is not None:
+                    environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
+                else:
+                    environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
+                val = eval(newId, environ, seq)
+            except Exception:      # set string if not a valid expression
+                val = newId
             seq.id=val
 
         if newDef is not None:
-#            try:
-            if taxo is not None:
-                environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
-            else:
-                environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
-            val = eval(newDef, environ, seq)
-#            except Exception,e:
-#                if options.onlyValid:
-#                    raise e
-#                val = newDef
+            try:
+                if taxo is not None:
+                    environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
+                else:
+                    environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
+                val = eval(newDef, environ, seq)
+            except Exception:    # set string if not a valid expression
+                val = newDef
             seq.definition=val
-#             
+             
         if newSeq is not None:
-#            try:
-            if taxo is not None:
-                environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
-            else:
-                environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
-            val = eval(newSeq, environ, seq)
-#            except Exception,e:
-#                if options.onlyValid:
-#                    raise e
-#                val = newSeq
+            try:
+                if taxo is not None:
+                    environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
+                else:
+                    environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
+                val = eval(newSeq, environ, seq)
+            except Exception:    # set string if not a valid expression
+                val = newSeq
             seq.seq=val
             if 'seq_length' in seq:
                 seq['seq_length']=len(seq)
@@ -251,15 +243,14 @@ def sequenceTaggerGenerator(config, taxo=None):
                 seq.view.delete_column(QUALITY_COLUMN)
                     
         if run is not None:
-#            try:
-            if taxo is not None:
-                environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
-            else:
-                environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
-            eval(run, environ, seq)
-#            except Exception,e:
-#                if options.onlyValid:
-#                    raise e
+            try:
+                if taxo is not None:
+                    environ = {'taxonomy' : taxo, 'sequence':seq, 'counter':counter[0], 'math':math}
+                else:
+                    environ = {'sequence':seq, 'counter':counter[0], 'math':math}     
+                eval(run, environ, seq)
+            except Exception,e:
+                raise e
 
     return sequenceTagger
 

python/obitools3/commands/export.pyx

@@ -59,13 +59,23 @@ def run(config):
      # Check that the input view has the type NUC_SEQS if needed    # TODO discuss, maybe bool property
     if (output[2] == Nuc_Seq) and (iview.type != b"NUC_SEQS_VIEW") :  # Nuc_Seq_Stored? TODO
         raise Exception("Error: the view to export in fasta or fastq format is not a NUC_SEQS view")
-   
+    
+    if config['obi']['only'] is not None:
+        withoutskip = min(input[4], config['obi']['only'])
+    else:
+        withoutskip = input[4]
+     
+    if config['obi']['skip'] is not None:
+        skip = min(input[4], config['obi']['skip'])
+    else:
+        skip = 0
+    
     # Initialize the progress bar
     if config['obi']['noprogressbar']:
         pb = None
     else:
-        pb = ProgressBar(len(iview), config, seconde=5)
-    
+        pb = ProgressBar(withoutskip - skip, config, seconde=5)
+
     i=0
     for seq in iview :
         PyErr_CheckSignals()

python/obitools3/version.py

@@ -1,5 +1,5 @@
 major = 3
 minor = 0
-serial= '0-beta13'
+serial= '0-beta15'
 
 version ="%d.%02d.%s" % (major,minor,serial)

src/build_reference_db.c

@@ -157,7 +157,7 @@ int build_reference_db(const char* dms_name,
 	ecotx_t* lca_2 = NULL;
 	ecotx_t* lca = NULL;
 	index_t idx1, idx2;
-	index_t i, j, k;
+	index_t i, j, k, count;
 	int32_t taxid_array_length;
 	int32_t score_array_length;
 	int32_t taxid_array_writable_length;
@@ -185,6 +185,7 @@ int build_reference_db(const char* dms_name,
 	matrix_view_name = strcpy(matrix_view_name, o_view_name);
 	strcat(matrix_view_name, "_matrix");
 
+	fprintf(stderr, "Aligning queries with reference database...\n");
 	if (obi_lcs_align_one_column(dms_name,
 								 refs_view_name,
 								 "",
@@ -320,13 +321,19 @@ int build_reference_db(const char* dms_name,
 		return -1;
 	}
 
+	count = (matrix_with_lca_view->infos)->line_count;
+	fprintf(stderr, "Computing LCAs...\n");
+
 	// Compute all the LCAs
 		// For each pair
-	for (i=0; i<(matrix_with_lca_view->infos)->line_count; i++)
+	for (i=0; i<count; i++)
 	{
 		if (! keep_running)
 			return -1;
 
+		if (i%1000 == 0)
+			fprintf(stderr,"\rDone : %f %%       ", (i / (float) count)*100);
+
 		// Read all taxids associated with the first sequence and compute their LCA
 		// Read line index
 		idx1 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx1_column, i, 0);
@@ -363,6 +370,7 @@ int build_reference_db(const char* dms_name,
 			return -1;
 		}
 	}
+	fprintf(stderr,"\rDone : 100 %%           \n");
 
 	// Clone refs view, add 2 arrays columns for lca and score, compute and write them
 
@@ -442,13 +450,18 @@ int build_reference_db(const char* dms_name,
 		return -1;
 	}
 
+	fprintf(stderr, "Building LCA arrays...\n");
+
 	// For each sequence, look for all its alignments in the matrix, and for each different LCA taxid/score, order them and write them
 	// Going through matrix once, filling refs arrays on the go for efficiency
-	for (i=0; i<(matrix_with_lca_view->infos)->line_count; i++)
+	for (i=0; i<count; i++)
 	{
 		if (! keep_running)
 			return -1;
 
+		if (i%1000 == 0)
+			fprintf(stderr,"\rDone : %f %%       ", (i / (float) count)*100);
+
 		// Read ref line indexes
 		idx1 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx1_column, i, 0);
 		idx2 = obi_get_int_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_idx2_column, i, 0);
@@ -464,6 +477,8 @@ int build_reference_db(const char* dms_name,
 		// Read alignment score
 		score = obi_get_float_with_elt_idx_and_col_p_in_view(matrix_with_lca_view, matrix_score_column, i, 0);
 
+		//fprintf(stderr, "\n\ntaxid_lca=%d, score=%f, idx1=%d, idx2=%d", taxid_lca, score, idx1, idx2);
+
 		///////////////// Compute for first sequence \\\\\\\\\\\\\\\\\\\\\\\     (TODO function)
 
 		// Read arrays
@@ -480,9 +495,11 @@ int build_reference_db(const char* dms_name,
 //			return -1;
 //		}
 
+		//fprintf(stderr, "\n1st sequence");
 		// If empty, add values
 		if (taxid_array_length == 0)
 		{
+			//fprintf(stderr, "\nEmpty, add value");
 			if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx1, &taxid_lca, (uint8_t) (obi_sizeof(OBI_INT) * 8), 1) < 0)
 			{
 				obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@@ -496,6 +513,8 @@ int build_reference_db(const char* dms_name,
 		}
 		else
 		{
+			//fprintf(stderr, "\nNot empty");
+
 			j = 0;
 			modified = false;
 			while (j < taxid_array_length)
@@ -509,6 +528,9 @@ int build_reference_db(const char* dms_name,
 						memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 						modified = true;
 
+						//fprintf(stderr, "\nSame LCA, replace %d and %f with %d and %f", lca_taxid_array_writable[j],
+						//		score_array_writable[j], taxid_lca, score);
+
 						// Better score for the same LCA, replace this LCA/score pair
 						lca_taxid_array_writable[j] = taxid_lca;
 						score_array_writable[j] = score;
@@ -535,6 +557,8 @@ int build_reference_db(const char* dms_name,
 				{
 					if (score > score_array[j])
 					{
+						//fprintf(stderr, "\nInsert new");
+
 						memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
 						memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 						modified = true;
@@ -579,10 +603,15 @@ int build_reference_db(const char* dms_name,
 				memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 				modified = true;
 
+				//fprintf(stderr, "\nAppend at the end");
+
 				// Append LCA
 				lca_taxid_array_writable[taxid_array_writable_length] = taxid_lca;
 				score_array_writable[score_array_writable_length] = score;
 
+				taxid_array_writable_length++;
+				score_array_writable_length++;
+
 				// Remove the previous (children) LCAs from the array if their score is equal or lower
 				while ((j>0) && (score_array_writable[j-1] <= score))
 				{
@@ -603,6 +632,13 @@ int build_reference_db(const char* dms_name,
 			// Write new arrays
 			if (modified)
 			{
+//				fprintf(stderr, "\n\nnew array:");
+//				for (k=0;k<taxid_array_writable_length;k++)
+//				{
+//					lca = obi_taxo_get_taxon_with_taxid(tax, lca_taxid_array_writable[k]);
+//					fprintf(stderr, "\nLCA=%d, %s, score=%f", lca_taxid_array_writable[k], lca->name, score_array_writable[k]);
+//				}
+
 				if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx1, lca_taxid_array_writable, (uint8_t) (obi_sizeof(OBI_INT) * 8), taxid_array_writable_length) < 0)
 				{
 					obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@@ -632,9 +668,13 @@ int build_reference_db(const char* dms_name,
 //			return -1;
 //		}
 
+		//fprintf(stderr, "\n2nd sequence");
+
 		// If empty, add values
 		if (taxid_array_length == 0)
 		{
+			//fprintf(stderr, "\nEmpty, add value");
+
 			if (obi_set_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, idx2, &taxid_lca, (uint8_t) (obi_sizeof(OBI_INT) * 8), 1) < 0)
 			{
 				obidebug(1, "\nError setting a LCA taxid array in a column when building a reference database");
@@ -648,6 +688,8 @@ int build_reference_db(const char* dms_name,
 		}
 		else
 		{
+			//fprintf(stderr, "\nNot empty");
+
 			j = 0;
 			modified = false;
 			while (j < taxid_array_length)
@@ -661,6 +703,9 @@ int build_reference_db(const char* dms_name,
 						memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 						modified = true;
 
+						//fprintf(stderr, "\nSame LCA, replace %d and %f with %d and %f", lca_taxid_array_writable[j],
+						//		score_array_writable[j], taxid_lca, score);
+
 						// Better score for the same LCA, replace this LCA/score pair
 						lca_taxid_array_writable[j] = taxid_lca;
 						score_array_writable[j] = score;
@@ -687,6 +732,8 @@ int build_reference_db(const char* dms_name,
 				{
 					if (score > score_array[j])
 					{
+						//fprintf(stderr, "\nInsert new");
+
 						memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
 						memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 						modified = true;
@@ -727,6 +774,8 @@ int build_reference_db(const char* dms_name,
 
 			if (j == taxid_array_length) // same or parent LCA not found, need to be appended at the end
 			{
+				//fprintf(stderr, "\nAppend at the end");
+
 				memcpy(lca_taxid_array_writable, lca_taxid_array, taxid_array_length*sizeof(obiint_t));
 				memcpy(score_array_writable, score_array, score_array_length*sizeof(obifloat_t));
 				modified = true;
@@ -735,6 +784,9 @@ int build_reference_db(const char* dms_name,
 				lca_taxid_array_writable[taxid_array_writable_length] = taxid_lca;
 				score_array_writable[score_array_writable_length] = score;
 
+				taxid_array_writable_length++;
+				score_array_writable_length++;
+
 				// Remove the previous (children) LCAs from the array if their score is equal or lower
 				while ((j>0) && (score_array_writable[j-1] <= score))
 				{
@@ -769,11 +821,17 @@ int build_reference_db(const char* dms_name,
 			}
 		}
 	}
+	fprintf(stderr,"\rDone : 100 %%           \n");
 
+	fprintf(stderr, "Writing results...\n");
+	count = (o_view->infos)->line_count;
 	// Fill empty LCA informations (because filling from potentially sparse alignment matrix) with the sequence taxid
 	score=1.0;  // technically getting LCA of identical sequences
-	for (i=0; i<(o_view->infos)->line_count; i++)
+	for (i=0; i<count; i++)
 	{
+		if (i%1000 == 0)
+			fprintf(stderr,"\rDone : %f %%       ", (i / (float) count)*100);
+
 		obi_get_array_with_col_p_in_view(o_view, final_lca_taxid_a_column, i, &taxid_array_length);
 		if (taxid_array_length == 0)  // no LCA set
 		{
@@ -799,6 +857,7 @@ int build_reference_db(const char* dms_name,
 			}
 		}
 	}
+	fprintf(stderr,"\rDone : 100 %%           \n");
 
 	// Add information about the threshold used to build the DB
 	snprintf(threshold_str, 5, "%f", threshold);
@@ -858,7 +917,6 @@ int build_reference_db(const char* dms_name,
 	free(matrix_view_name);
 	free(matrix_with_lca_view_name);
 
-	fprintf(stderr,"\rDone : 100 %%           \n");
 	return 0;
 }
 

src/obi_ecopcr.c

@@ -1061,7 +1061,7 @@ int obi_ecopcr(const char* i_dms_name,
 											length = 0;
 											if (posj > posi)
 												length = posj - posi - o1->patlen - o2->patlen;
-											if (posj < posi)
+											else if (circular > 0)
 												length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
 											if ((length>0) &&	// For when primers touch or overlap
 												(!min_len || (length >= min_len)) &&
@@ -1151,7 +1151,7 @@ int obi_ecopcr(const char* i_dms_name,
 											length = 0;
 											if (posj > posi)
 												length = posj - posi + 1  - o2->patlen - o1->patlen; /* - o1->patlen : deleted by <EC> (prior to the OBITools3) */
-											if (posj < posi)
+											else if (circular > 0)
 												length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
 											if ((length>0) &&	// For when primers touch or overlap
 												(!min_len || (length >= min_len)) &&
@@ -1232,7 +1232,7 @@ int obi_ecopcr(const char* i_dms_name,
 		return -1;
 	}
 
-	fprintf(stderr,"\rDone : 100 %%           ");
+	fprintf(stderr,"\rDone : 100 %%           \n");
 	return 0;
 
 	return 0;

src/obi_ecotag.c

@@ -455,7 +455,7 @@ int obi_ecotag(const char* dms_name,
 
 	for (i=0; i < query_count; i++)
 	{
-		if (i%100 == 0)
+		if (i%1000 == 0)
 			fprintf(stderr,"\rDone : %f %%       ", (i / (float) query_count)*100);
 
 		best_match_count = 0;
@@ -562,7 +562,7 @@ int obi_ecotag(const char* dms_name,
 			score_array = obi_get_array_with_col_p_in_view(ref_view, score_a_column, best_match_idx, &lca_array_length);
 
 			k = 0;
-			while ((k < lca_array_length) && (score_array[k] >= ecotag_threshold))
+			while ((k < lca_array_length) && (score_array[k] >= best_score))
 				k++;
 
 			if (k>0)
@@ -570,12 +570,12 @@ int obi_ecotag(const char* dms_name,
 				lca_array = obi_get_array_with_col_p_in_view(ref_view, lca_taxid_a_column, best_match_idx, &lca_array_length);
 				if (j>0)
 				{
-					lca = obi_taxo_get_taxon_with_taxid(taxonomy, lca_taxid);
-					if (lca == NULL)
-					{
-						obidebug(1, "\nError getting a taxon from a taxid when doing taxonomic assignment");
-						return -1;
-					}
+//					lca = obi_taxo_get_taxon_with_taxid(taxonomy, lca_taxid);
+//					if (lca == NULL)
+//					{
+//						obidebug(1, "\nError getting a taxon from a taxid when doing taxonomic assignment");
+//						return -1;
+//					}
 					lca_in_array = obi_taxo_get_taxon_with_taxid(taxonomy, lca_array[k-1]);
 					if (lca_in_array == NULL)
 					{

src/obidms_taxonomy.c

@@ -2376,9 +2376,10 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
 				// and the deleted taxids with no current reference. An element of the list is composed of the taxid, and the index
 				// of the taxon in the taxa structure, or -1 for deleted taxids.
 				// Creating the merged list requires to merge the 3 ordered lists into one.
-				while (((nT < (tax->taxa)->count) && ((tax->taxa)->taxon[nT].taxid < old_taxid)) || ((nD >= 0) && (delnodes[nD] < old_taxid)))
+				while (((nT < (tax->taxa)->count) && ((tax->taxa)->taxon[nT].taxid < old_taxid)) ||
+						((nD >= 0) && (delnodes[nD] < old_taxid)))
 				{
-					if ((tax->taxa)->taxon[nT].taxid < delnodes[nD])
+					if ((nT < (tax->taxa)->count) && (tax->taxa)->taxon[nT].taxid < delnodes[nD])
 					{ // Add element from taxa list
 						// Enlarge structure if needed
 						if (n == buffer_size)
@@ -2401,7 +2402,7 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
 						nT++;
 						n++;
 					}
-					else if (delnodes[nD] < (tax->taxa)->taxon[nT].taxid)
+					else
 					{ // Add element from deleted taxids list
 						// Enlarge structure if needed
 						if (n == buffer_size)
@@ -3036,12 +3037,12 @@ OBIDMS_taxonomy_p obi_read_taxonomy(OBIDMS_p dms, const char* taxonomy_name, boo
 
 	strcpy(tax->tax_name, taxonomy_name);
 
-	buffer_size = 2048;
-
 	taxonomy_path = get_taxonomy_path(dms, taxonomy_name);
 	if (taxonomy_path == NULL)
 		return NULL;
 
+	buffer_size = strlen(taxonomy_path) + strlen(taxonomy_name) + 6;
+
 	// Read ranks
 	ranks_file_name = (char*) malloc(buffer_size*sizeof(char));
 	if (ranks_file_name == NULL)

src/obidmscolumn.c

@@ -1973,7 +1973,11 @@ int obi_enlarge_column(OBIDMS_column_p column)
 
 	// Calculate the new file size
 	old_line_count = (column->header)->line_count;
-	new_line_count = old_line_count * COLUMN_GROWTH_FACTOR;
+	new_line_count = ceil((double) old_line_count * (double) COLUMN_GROWTH_FACTOR);
+	if (new_line_count > old_line_count+100000)
+		new_line_count = old_line_count+100000;
+	else if (new_line_count < old_line_count+1000)
+		new_line_count = old_line_count+1000;
 
 	if (new_line_count > MAXIMUM_LINE_COUNT)
 	{
@@ -2381,6 +2385,54 @@ char* obi_get_elements_names(OBIDMS_column_p column)
 }
 
 
+char* obi_get_formatted_elements_names(OBIDMS_column_p column)
+{
+	char* elements_names;
+	int   i, j;
+	int   elt_idx;
+	int   len;
+
+	elements_names = (char*) malloc(((column->header)->elements_names_length + (column->header)->nb_elements_per_line) * sizeof(char));
+	if (elements_names == NULL)
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError allocating memory for elements names");
+		return NULL;
+	}
+
+	j = 0;
+	for (i=0; i < (column->header)->nb_elements_per_line; i++)
+	{
+		elt_idx = ((column->header)->elements_names_idx)[i];
+		len = strlen(((column->header)->elements_names)+elt_idx);
+		memcpy(elements_names+j, ((column->header)->elements_names)+elt_idx, len*sizeof(char));
+		j = j + len;
+		elements_names[j] = ';';
+		j++;
+		elements_names[j] = ' ';
+		j++;
+	}
+
+	elements_names[j - 1] = '\0';
+
+	return elements_names;
+}
+
+
+char* obi_column_formatted_infos(OBIDMS_column_p column)
+{
+	char* column_infos;
+	char* elt_names;
+
+	column_infos = malloc(1024 * sizeof(char));
+
+	elt_names = obi_get_formatted_elements_names(column);
+
+
+	free(elt_names);
+	return column_infos;
+}
+
 
 int obi_column_prepare_to_set_value(OBIDMS_column_p column, index_t line_nb, index_t elt_idx)
 {

src/obidmscolumn.h

@@ -505,6 +505,14 @@ index_t obi_column_get_element_index_from_name(OBIDMS_column_p column, const cha
 char* obi_get_elements_names(OBIDMS_column_p column);
 
 
+// TODO
+//char* obi_get_formatted_elements_names(OBIDMS_column_p column);
+
+
+// TODO
+//char* obi_column_formatted_infos(OBIDMS_column_p column);
+
+
 /**
  * @brief Prepares a column to set a value.
  *

src/sse_banded_LCS_alignment.c

@@ -686,6 +686,9 @@ int calculateSizeToAllocate(int maxLen, int LCSmin)
 	size *=  3;
 	size +=  16;
 
+	size += 10;  // band-aid for memory bug I don't understand (triggered on specific db on ubuntu)
+				 // bug might have to do with the way different systems behave when aligning the address in obi_get_memory_aligned_on_16
+
 	return(size*sizeof(int16_t));
 }