resure alignment

Assembly.cpp

@@ -773,8 +773,8 @@ void* Overlap_calculate_heap_merge(void* arg)
                         /************需要注释掉**********/
                         // debug_overlap = debug_overlap + list_length;
                         /************需要注释掉**********/
-                        ///append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 0);
+                        append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 0);
-                        append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 0, &g_read, &R_INF, &correct);
+                        ///append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 0, &g_read, &R_INF, &correct);
                     }
                     ///merge_Candidates_list(&l, list, list_length, end_pos, 0);
                     ///merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 0);
@@ -833,8 +833,8 @@ void* Overlap_calculate_heap_merge(void* arg)
                         /************需要注释掉**********/
                         // debug_overlap = debug_overlap + list_length;
                         /************需要注释掉**********/
-                        ///append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 1);
+                        append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 1);
-                        append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 1, &g_read, &R_INF, &correct);
+                        ///append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 1, &g_read, &R_INF, &correct);
                     }
                     ///merge_Candidates_list(&l, list, list_length, end_pos, 1);
                     //merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 1);
@@ -879,24 +879,6 @@ void* Overlap_calculate_heap_merge(void* arg)
        correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read,
        &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1);
       
-        /**
-        for (j = 0; j < overlap_list.length; j++)
-        {
-            long long Len_x = overlap_list.list[j].x_pos_e -  overlap_list.list[j].x_pos_s + 1;
-
-            if (Len_x * 0.9 <=  overlap_list.list[j].align_length)
-            {
-                if (overlap_list.list[j].y_pos_strand == 0)
-                {
-                    matched_overlap_0++;
-                }
-                else
-                {
-                    matched_overlap_1++;
-                }
-            }
-        }
-        **/
         
         /**
         POA_i = 0;
@@ -994,8 +976,8 @@ void* Overlap_calculate_heap_merge(void* arg)
     //destory_Candidates_list(&debug_l);
 
     destory_Heap(&heap);
-    ///destory_k_mer_pos_list_alloc(&array_list);
+    destory_k_mer_pos_list_alloc(&array_list);
-    destory_k_mer_pos_list_alloc_prefilter(&array_list);
+    ///destory_k_mer_pos_list_alloc_prefilter(&array_list);
 
     destory_Graph(&POA_Graph);
     destory_UC_Read(&g_read);

Correct.h

@@ -4,9 +4,8 @@
 #include "Hash_Table.h"
 #include "Levenshtein_distance.h"
 
-#define WINDOW 350
+
-#define THRESHOLD  15
+
-#define GROUP_SIZE 4
 
 #define MAX(x, y) ((x >= y)?x:y)
 #define MIN(x, y) ((x <= y)?x:y)
@@ -22,6 +21,9 @@ typedef struct
     uint64_t start_i;
     char overlap_region[WINDOW + THRESHOLD*2 + 10];
     char overlap_region_group[GROUP_SIZE][WINDOW + THRESHOLD*2 + 10];
+    char path[WINDOW + THRESHOLD*2 + 10];
+    int path_length;
+    Word matrix_bit[((WINDOW + 10)<<3)];
     __m128i Peq_SSE[256];
 } Correct_dumy;
 

Hash_Table.cpp

@@ -212,15 +212,29 @@ void init_overlap_region_alloc(overlap_region_alloc* list)
 {
     list->size = 1000;
     list->length = 0;
-    list->list = (overlap_region*)malloc(sizeof(overlap_region)*list->size);
+    ///list->list = (overlap_region*)malloc(sizeof(overlap_region)*list->size);
+    list->list = (overlap_region*)calloc(list->size, sizeof(overlap_region));
 }
 void clear_overlap_region_alloc(overlap_region_alloc* list)
 {
     list->length = 0;
+    int i = 0;
+    for (i = 0; i < list->size; i++)
+    {   
+        list->list[i].w_list_length = 0;
+    }
 }
 
 void destory_overlap_region_alloc(overlap_region_alloc* list)
 {
+    int i = 0;
+    for (i = 0; i < list->size; i++)
+    {
+        if (list->list[i].w_list_size != 0)
+        {
+            free(list->list[i].w_list);
+        }
+    }
     free(list->list);
 }
 
@@ -235,6 +249,8 @@ void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp
     {
         list->size = list->size * 2;
         list->list = (overlap_region*)realloc(list->list, sizeof(overlap_region)*list->size);
+        ///新分配空间要初始化
+        memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2));
     }
 
     if (list->length!=0 && 
@@ -641,6 +657,28 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF)
 
 
 
+void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error)
+{
+    
+    long long length = region->x_pos_e - region->x_pos_s + 1;
+    long long num_windows = length / WINDOW + 2;
+
+    ///w_list_length会在clear_overlap_region_alloc中定时清空
+    if (num_windows > region->w_list_size)
+    {
+        region->w_list_size = num_windows;
+        region->w_list = (window_list*)realloc(region->w_list, region->w_list_size*sizeof(window_list));
+    }
+    
+
+    region->w_list[region->w_list_length].x_start = x_start;
+    region->w_list[region->w_list_length].x_end = x_end;
+    region->w_list[region->w_list_length].y_start = y_start;
+    region->w_list[region->w_list_length].y_end = y_end;
+    region->w_list[region->w_list_length].error = error;
+    region->w_list[region->w_list_length].cigar.length = -1;
+    region->w_list_length++;
+}
 
 
 

Hash_Table.h

@@ -13,6 +13,14 @@ typedef khash_t(POS64) Pos_Table;
 #define MAX_SUFFIX_BITS 64
 #define MODE_VALUE 101
 
+///#define WINDOW 350
+///#define THRESHOLD  14
+#define WINDOW 375
+#define THRESHOLD  15
+#define THRESHOLD_RATE 0.04
+#define GROUP_SIZE 4
+///最长是10M10D10M10D10M这种
+#define CIGAR_MAX_LENGTH THRESHOLD*2+2
 
 typedef struct
 {
@@ -53,6 +61,26 @@ typedef struct
     uint64_t length;
 } k_mer_pos_list_alloc;
 
+
+typedef struct
+{
+  int C_L[CIGAR_MAX_LENGTH];
+  char C_C[CIGAR_MAX_LENGTH];
+  int length;
+} CIGAR;
+
+typedef struct
+{
+  uint64_t x_start;
+  uint64_t x_end;
+  int y_end;
+  int y_start;
+  ///int y_pre_start;
+  ///error小于等于0都要重新算
+  int error;
+  CIGAR cigar;
+} window_list;
+
 typedef struct
 {
     uint64_t x_id;
@@ -67,6 +95,10 @@ typedef struct
 
     uint64_t shared_seed;
     uint64_t align_length;
+
+    window_list* w_list;
+    uint64_t w_list_size;
+    uint64_t w_list_length;
 } overlap_region;
 
 
@@ -406,8 +438,7 @@ void destory_overlap_region_alloc(overlap_region_alloc* list);
 void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF);
 void calculate_overlap_region(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
 uint64_t readID, uint64_t readLength, All_reads* R_INF);
-
+void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error);
-
 
 
 

Levenshtein_distance.h

@@ -172,6 +172,10 @@ inline int Reserve_Banded_BPM
 	///int site = p_length - last_high - 1;
 	int site = t_length - 1;
 	int return_site = -1;
+	///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
+	///也就是p_length - t_length
+	///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
+	int available_i = p_length - t_length;
 	if ((err <= errthold) && (err<=*return_err))
 	{
 		*return_err = err;
@@ -181,9 +185,7 @@ inline int Reserve_Banded_BPM
 	i = 0;
 
 
-
+	while (i < available_i)
-
-	while (i<errthold)
 	{
 		err = err + ((VP >> i)&(Word)1);
 		err = err - ((VN >> i)&(Word)1);
@@ -196,45 +198,85 @@ inline int Reserve_Banded_BPM
 		}
 	}
 
-
-	unsigned int ungap_err;
-	ungap_err = err;
-
-
-	while (i<last_high)
-	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
-		++i;
-
-		if ((err <= errthold) && (err<=*return_err))
-		{
-			*return_err = err;
-			return_site = site + i;
-		}
-
-
-
-
-	}
-
-
-	if ((ungap_err <= errthold) && (ungap_err == *return_err))
-	{
-		return_site = site + errthold;
-	}
-
 	return return_site;
 
 }
 
 
 
+inline int try_cigar(char *pattern, int p_length,
+	char *text, int t_length, int end_site, char* path,
+	int error, 
+	int* return_start_site, 
+	int* return_path_length)
+{
+	int i = 0;
+	int tmp_err = 0;
+	///y上的起始位置
+	int start_site = end_site - t_length + 1;
 
+	if (start_site >= 0)
+	{
+
+		for (i = 0; i < t_length; i++)
+		{
+			///path[i] = 0;
+			///path倒着存
+			path[t_length - i - 1] = 0;
+			if (text[i] != pattern[i + start_site])
+			{
+				path[t_length - i - 1] = 1;
+				tmp_err++;
+
+				if (tmp_err > error)
+				{
+					return 0;
+				}
+			}
+		}
+
+		if (tmp_err == error)
+		{
+			(*return_path_length) = t_length;
+
+			(*return_start_site) = start_site;
+
+			return 6;
+		}
+	}
+
+
+	return 0;
+}
+
+
+///p_length有可能不够，但是t_length总是够的
+///就是p_length有可能小于t_length + 2 * errthold
 inline int Reserve_Banded_BPM_PATH
 (char *pattern, int p_length, char *text, int t_length, unsigned short errthold, 
-		unsigned int* return_err, int* return_start_site, int* return_path_length, Word* matrix_bit, char* path)
+		unsigned int* return_err, int* return_start_site, int* return_path_length, Word* matrix_bit, char* path, 
+		int old_error, int old_end_site)
 {
+	if (old_error != -1 && old_end_site != -1)
+	{
+		if (old_error == 0)
+		{
+			///fprintf(stderr, "0 error\n");
+			(*return_err) = old_error;
+			(*return_start_site) = old_end_site - t_length + 1;
+			return old_end_site;
+		}
+
+		if (try_cigar(pattern, p_length, text, t_length, old_end_site, path,
+		old_error, return_start_site, return_path_length))
+		{
+			///fprintf(stderr, "no gap error\n");
+			(*return_err) = old_error;
+			return old_end_site;
+		}
+		
+	}
+	
 	(*return_err) = (unsigned int)-1;
 
 	Word Peq[256];
@@ -330,10 +372,8 @@ inline int Reserve_Banded_BPM_PATH
 		{
 			++err;
 
-			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
 			if ((err - last_high)>errthold)
 			{
-				///fprintf(stderr, "0 ######, i: %u\n", i);
 				return -1;
 			}
 				
@@ -396,6 +436,10 @@ inline int Reserve_Banded_BPM_PATH
 	///int site = p_length - last_high - 1;
 	int site = t_length - 1;
 	int return_site = -1;
+	///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
+	///也就是p_length - t_length
+	///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
+	int available_i = p_length - t_length;
 	if ((err <= errthold) && (err<=*return_err))
 	{
 		*return_err = err;
@@ -404,10 +448,7 @@ inline int Reserve_Banded_BPM_PATH
 	int i_last = i;
 	i = 0;
 
-
+	while (i < available_i)
-
-
-	while (i<errthold)
 	{
 		err = err + ((VP >> i)&(Word)1);
 		err = err - ((VN >> i)&(Word)1);
@@ -421,40 +462,15 @@ inline int Reserve_Banded_BPM_PATH
 	}
 
 
-	unsigned int ungap_err;
-	ungap_err = err;
-
-
-	while (i<last_high)
-	{
-		err = err + ((VP >> i)&(Word)1);
-		err = err - ((VN >> i)&(Word)1);
-		++i;
-
-		if ((err <= errthold) && (err<=*return_err))
-		{
-			*return_err = err;
-			return_site = site + i;
-		}
-
-
-
-
-	}
-
-
-	if ((ungap_err <= errthold) && (ungap_err == *return_err))
-	{
-		return_site = site + errthold;
-	}
-
 	if ((*return_err) == (unsigned int)-1)
 	{
 		return return_site;
 	}
 	
-
+    ////注意，这里p_length要矫正啊啊
-
+	///不矫正会出错
+	///因为p_length有可能不够
+	p_length = t_length + 2 * errthold;
 	///end_site是正确的
 	int end_site = return_site;
 	int start_site = end_site;
@@ -608,6 +624,8 @@ inline int Reserve_Banded_BPM_PATH
 
 }
 
+
+////这个p_length四个是一样的
 inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *pattern3, char *pattern4, int p_length, char *text, int t_length,
 	int* return_sites, unsigned int* return_sites_error, unsigned short errthold, __m128i* Peq_SSE)
 
@@ -829,8 +847,12 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
 	i = 0;
 
 
+	///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
+	///也就是p_length - t_length
+	///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
+	int available_i = p_length - t_length;
 
-	while (i<errthold)
+	while (i < available_i)
 	{
 		///err = err + ((VP >> i)&(Word_32)1);
 		tmp_process = _mm_srli_epi32(VP, i);
@@ -871,87 +893,7 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
 		}
 	}
 
-
-	unsigned int ungap_err1;
-	unsigned int ungap_err2;
-	unsigned int ungap_err3;
-	unsigned int ungap_err4;
-	ungap_err1 = err1;
-	ungap_err2 = err2;
-	ungap_err3 = err3;
-	ungap_err4 = err4;
-
-
-
-
-	while (i<last_high)
-	{
-		///err = err + ((VP >> i)&(Word_32)1);
-		tmp_process = _mm_srli_epi32(VP, i);
-		tmp_process = _mm_and_si128(tmp_process, err_mask);
-		Err_4 = _mm_add_epi32(Err_4, tmp_process);
-
-		///err = err - ((VN >> i)&(Word_32)1);
-		tmp_process1 = _mm_srli_epi32(VN, i);
-		tmp_process1 = _mm_and_si128(tmp_process1, err_mask);
-		Err_4 = _mm_sub_epi32(Err_4, tmp_process1);
-		++i;
-
-		err1 = _mm_extract_epi32(Err_4, 0);
-		err2 = _mm_extract_epi32(Err_4, 1);
-		err3 = _mm_extract_epi32(Err_4, 2);
-		err4 = _mm_extract_epi32(Err_4, 3);
-
-
-		if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
-		{
-			return_sites[0] = site + i;
-			return_sites_error[0] = err1;
-		}
-		if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
-		{
-			return_sites[1] = site + i;
-			return_sites_error[1] = err2;
-		}
-		if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
-		{
-			return_sites[2] = site + i;
-			return_sites_error[2] = err3;
-		}
-		if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
-		{
-			return_sites[3] = site + i;
-			return_sites_error[3] = err4;
-		}
-	}
-
-
-	if ((ungap_err1 <= errthold) && ungap_err1 == return_sites_error[0])
-	{
-		return_sites[0] = site + errthold;
-	}
-
-
-	if ((ungap_err2 <= errthold) && ungap_err2 == return_sites_error[1])
-	{
-		return_sites[1] = site + errthold;
-	}
-
-
-	if ((ungap_err3 <= errthold) && ungap_err3 == return_sites_error[2])
-	{
-		return_sites[2] = site + errthold;
-	}
-
-
-	if ((ungap_err4 <= errthold) && ungap_err4 == return_sites_error[3])
-	{
-		return_sites[3] = site + errthold;
-	}
-
 	return 1;
-	
-
 }
 
 

main.cpp

@@ -138,9 +138,11 @@ void debug_edit_distance()
     char* x = "TTCCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCACTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCAGTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCGATTCCATTCAAATCAATTCCGTTCTATTCAATTTGATTCAT";
     char* y = "CCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCATTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCACTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCCATTCCATTCAATTCAATTCCGTGCTATTCAATTTGATTCATTTCCATTTAATTCCATTCCATTAGATTCCATT";
     **/
-   unsigned short toold = 7;
+   unsigned short toold = 15;
-   char* x = "TTCCATACGATTCCATTCAATTCGAGACCATTCTATTCCT";
+   char* x 
-    char* y = "CCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGT";
+   = "GAAAGAGAATCAAATGGAATTGAATCGAATGGAATCGAATGGATTGGAAAGGAATAGAATGGAATGGAATGGAATTGACTCAAATGGAATGGACTAGAATGGAATGGATTCGAATGGAAGGCAAAGGAATGGAATCTATCGGAATGGACTGTAATGGAATGGAATGGAAGGGATTGGAATGGATTCGAATGTAATGGACTGCAATAGAAAGGATTCGAATGGAATGAAAAAGAATTGAATGGAATAGAACAGAATGGAATCAAATCGAAGGAAATGGAATGGAATAGAAAGGAATGGAATGAAATGGAATGGAAAGGATTCGAATGGAATGCAATCGAATGGAATGGAATCGAACGGAATGGAATAAAATGGAAG";
+    char* y = 
+     "GAAAGAGAATCAAATGGAATTGAATCGAATGGAATCGAATGGATTGGAAAGGAATAGAATGGAATGGAATGGAATGGACTCAAATGGAATGTACTAGAATGGAATGGATTCGAATGGAAGGCAAAGGAATGGAATCTATTGGAATGGACTGTAATGGAATGGAATGGAAGGGATTGGAATGGACTCGAATGGAATGGACTGCAATAGAAAGGATTCGAATGGAATGAAAAAGAATTGAATGGAATAGAACAGAATGGAATCAAATCGAATGAAATGGAATGGAATAGAAAGGAATGGAATGAAATGGAATGGAAAGGATTCGAATGGAATGCAATCGAATGGAATGGAATCGAACGGAATGGAATAAATTTTCTG";
     fprintf(stderr, "x_length: %u\n", strlen(x));
     fprintf(stderr, "y_length: %u\n", strlen(y));
 
@@ -204,6 +206,8 @@ int main(int argc, char *argv[])
 
     return 1;
     **/
+    
+
 
 
     fprintf(stdout, "k-mer length: %d\n",k_mer_length);