debug

Assembly.cpp

@@ -16,6 +16,9 @@ All_reads R_INF;
 pthread_mutex_t statistics;
 long long total_matched_overlap_0 = 0;
 long long total_matched_overlap_1 = 0;
+long long total_potiental_matched_overlap_0 = 0;
+long long total_potiental_matched_overlap_1 = 0;
+
 long long complete_threads = 0;
 
 
@@ -689,6 +692,8 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     long long matched_overlap_0 = 0;
     long long matched_overlap_1 = 0;
+    long long potiental_matched_overlap_0 = 0;
+    long long potiental_matched_overlap_1 = 0;
     long long j;
 
     int thr_ID = *((int*)arg);
@@ -698,6 +703,10 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     UC_Read g_read;
     init_UC_Read(&g_read);
+
+    UC_Read overlap_read;
+    init_UC_Read(&overlap_read);
+
     HPC_seq HPC_read;
     Hash_code k_code;
     uint64_t code;
@@ -852,29 +861,30 @@ void* Overlap_calculate_heap_merge(void* arg)
         
         ///merge_k_mer_pos_list_alloc(&array_list, &l);
         merge_k_mer_pos_list_alloc_heap_sort(&array_list, &l, &heap);
-        /**
+        
-        if (array_list.length < 3)
+        // if (array_list.length < 3)
-        {
+        // {
-            merge_k_mer_pos_list_alloc(&array_list, &l);
+        //     merge_k_mer_pos_list_alloc(&array_list, &l);
-        }
+        // }
-        else
+        // else
-        {
+        // {
-            merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap);
+        //     merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap);
-        }
+        // }
-        **/
+        
 
        ///以x_pos_e，即结束位置为主元排序
        calculate_overlap_region(&l, &overlap_list, i, g_read.length, &R_INF);
 
 
-       correct_overlap(&overlap_list, &R_INF, &g_read, &correct);
+       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.6 <=  overlap_list.list[j].align_length)
+            if (Len_x * 0.9 <=  overlap_list.list[j].align_length)
             {
                 if (overlap_list.list[j].y_pos_strand == 0)
                 {
@@ -886,6 +896,7 @@ void* Overlap_calculate_heap_merge(void* arg)
                 }
             }
         }
+        **/
         
         /**
         POA_i = 0;
@@ -988,6 +999,8 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     destory_Graph(&POA_Graph);
     destory_UC_Read(&g_read);
+    destory_UC_Read(&overlap_read);
+    
 
     destory_Correct_dumy(&correct);
 
@@ -995,11 +1008,15 @@ void* Overlap_calculate_heap_merge(void* arg)
     pthread_mutex_lock(&statistics);
     total_matched_overlap_0 += matched_overlap_0;
     total_matched_overlap_1 += matched_overlap_1;
+    total_potiental_matched_overlap_0 += potiental_matched_overlap_0;
+    total_potiental_matched_overlap_1 += potiental_matched_overlap_1;
     complete_threads++;
     if(complete_threads == thread_num)
     {
         fprintf(stderr, "total_matched_overlap_0: %llu\n", total_matched_overlap_0);
         fprintf(stderr, "total_matched_overlap_1: %llu\n", total_matched_overlap_1);
+        fprintf(stderr, "total_potiental_matched_overlap_0: %llu\n", total_potiental_matched_overlap_0);
+        fprintf(stderr, "total_potiental_matched_overlap_1: %llu\n", total_potiental_matched_overlap_1);
     }
 	pthread_mutex_unlock(&statistics);
 }

Correct.cpp

@@ -674,7 +674,11 @@ char* r_string)
         
 }
 
-void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy)
+
+void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, 
+                        long long* matched_overlap_0, long long* matched_overlap_1, 
+                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1)
 {
     reverse_complement(g_read->seq, g_read->length);
 
@@ -733,6 +737,86 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re
     // T_total_match, T_total_unmatch, T_total_mis);
 	// pthread_mutex_unlock(&debug_statistics);
     /************需要注释掉********* */
+    /**
+    long long j;
+    long long Len_x;
+    int threshold;
+    long long y_start;
+    long long Len_y;
+    long long currentIDLen;
+    for (j = 0; j < overlap_list->length; j++)
+    {
+        Len_x = overlap_list->list[j].x_pos_e -  overlap_list->list[j].x_pos_s + 1;
+
+        if (Len_x * 0.6 <=  overlap_list->list[j].align_length)
+        {
+            threshold = Len_x * 0.04;
+
+            y_start = overlap_list->list[j].y_pos_s - threshold;
+            if (y_start < 0)
+            {
+                y_start = 0;
+            }
+
+            
+            ///当前y的长度
+            currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[j].y_id);
+            ///不能超过y的剩余长度
+            Len_y = MIN(Len_x + 2 * threshold, currentIDLen - y_start);
+            
+            
+            if (overlap_list->list[j].y_pos_strand == 0)
+            {
+                recover_UC_Read(overlap_read, R_INF, overlap_list->list[j].y_id);
+                (*matched_overlap_0)++;
+            }
+            else
+            {
+                recover_UC_Read_RC(overlap_read, R_INF, overlap_list->list[j].y_id);
+                (*matched_overlap_1)++;
+            }
+        
+
+            
+
+            EdlibAlignResult result = edlibAlign(g_read->seq + overlap_list->list[j].x_pos_s, Len_x, 
+            overlap_read->seq + y_start, Len_y, 
+            edlibNewAlignConfig(threshold, EDLIB_MODE_HW, EDLIB_TASK_PATH, NULL, 0));
+            
+
+            if (result.status == EDLIB_STATUS_OK && result.editDistance != -1) {
+                char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD);
+                int cigar_length = strlen(cigar);
+                free(cigar);
+
+                
+                if (overlap_list->list[j].y_pos_strand == 0)
+                {
+                    (*potiental_matched_overlap_0)++;
+                }
+                else
+                {
+                    (*potiental_matched_overlap_1)++;
+                }
+                
+                // if (overlap_list->list[j].shared_seed == 1 && Len_x >= 1000)
+                // {
+                //     (*matched_overlap_0)++;
+                // }
+
+                // if (overlap_list->list[j].shared_seed == 2 && Len_x >= 1000)
+                // {
+                //     (*matched_overlap_1)++;
+                // }
+            
+                
+            }
+
+            edlibFreeAlignResult(result);
+        }
+    }
+    **/
+    
 }
 
 

Correct.h

@@ -26,8 +26,10 @@ typedef struct
 } Correct_dumy;
 
 
-void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy);
+void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
-
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, 
+                        long long* matched_overlap_0, long long* matched_overlap_1, 
+                        long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1);
 void init_Correct_dumy(Correct_dumy* list);
 void destory_Correct_dumy(Correct_dumy* list);
 void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list);

Levenshtein_distance.h

@@ -228,6 +228,386 @@ inline int Reserve_Banded_BPM
 
 }
 
+
+
+
+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)
+{
+	(*return_err) = (unsigned int)-1;
+
+	Word Peq[256];
+
+	int band_length = (errthold << 1) + 1;
+	int i = 0;
+	Word tmp_Peq_1 = (Word)1;
+
+	Peq['A'] = (Word)0;
+	Peq['T'] = (Word)0;
+	Peq['G'] = (Word)0;
+	Peq['C'] = (Word)0;
+
+
+	Word Peq_A;
+	Word Peq_T;
+	Word Peq_C;
+	Word Peq_G;
+
+	///band_length = 2k + 1
+	for (i = 0; i<band_length; i++)
+	{
+		Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
+		tmp_Peq_1 = tmp_Peq_1 << 1;
+	}
+
+	///Peq['T'] = Peq['T'] | Peq['C'];
+
+	Peq_A = Peq['A'];
+	Peq_C = Peq['C'];
+	Peq_T = Peq['T'];
+	Peq_G = Peq['G'];
+
+
+	memset(Peq, 0, sizeof(Word)* 256);
+
+
+	Peq['A'] = Peq_A;
+	Peq['C'] = Peq_C;
+	Peq['T'] = Peq_T;
+	Peq['G'] = Peq_G;
+
+
+	
+
+	Word Mask = ((Word)1 << (errthold << 1));
+
+	Word VP = 0;
+	Word VN = 0;
+	Word X = 0;
+	Word D0 = 0;
+	Word HN = 0;
+	Word HP = 0;
+
+
+	i = 0;
+
+	Word column_start;
+
+	int err = 0;
+
+	Word err_mask = (Word)1;
+
+
+	///band_down = 2k
+	///i_bd = 2k
+	///int i_bd = i + band_down;
+	int i_bd = (errthold << 1);
+
+
+	int last_high = (errthold << 1);
+
+
+	/// t_length_1 = SEQ_LENGTH - 1
+	int t_length_1 = t_length - 1;
+	//while(i<t_length)
+
+	while (i<t_length_1)
+	{
+		///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
+		X = Peq[text[i]] | VN;
+
+		D0 = ((VP + (X&VP)) ^ VP) | X;
+
+		HN = VP&D0;
+		HP = VN | ~(VP | D0);
+
+		X = D0 >> 1;
+		VN = X&HP;
+		VP = HN | ~(X | HP);
+
+		if (!(D0&err_mask))
+		{
+			++err;
+
+			///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
+			if ((err - last_high)>errthold)
+			{
+				///fprintf(stderr, "0 ######, i: %u\n", i);
+				return -1;
+			}
+				
+		}
+
+
+		Peq['A'] = Peq['A'] >> 1;
+		Peq['C'] = Peq['C'] >> 1;
+		Peq['G'] = Peq['G'] >> 1;
+		Peq['T'] = Peq['T'] >> 1;
+
+
+		++i;
+		++i_bd;
+		Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
+
+
+		///Peq['T'] = Peq['T'] | Peq['C'];
+
+		column_start = i << 3;
+		matrix_bit[column_start] = D0;
+		matrix_bit[column_start + 1] = VP;
+		matrix_bit[column_start + 2] = VN;
+		matrix_bit[column_start + 3] = HP;
+		matrix_bit[column_start + 4] = HN;
+	}
+
+
+
+
+
+	X = Peq[text[i]] | VN;
+	D0 = ((VP + (X&VP)) ^ VP) | X;
+	HN = VP&D0;
+	HP = VN | ~(VP | D0);
+	X = D0 >> 1;
+	VN = X&HP;
+	VP = HN | ~(X | HP);
+	if (!(D0&err_mask))
+	{
+		++err;
+		if ((err - last_high)>errthold)
+			return -1;
+	}
+
+
+	column_start = (i + 1) << 3;
+	matrix_bit[column_start] = D0;
+	matrix_bit[column_start + 1] = VP;
+	matrix_bit[column_start + 2] = VN;
+	matrix_bit[column_start + 3] = HP;
+	matrix_bit[column_start + 4] = HN;
+
+
+	////fprintf(stderr, "sucess(2)\n");
+
+	/// last_high = 2k
+	/// site = (SEQ_LENGTH + 2k) - 2k -1
+	/// site = SEQ_LENGTH - 1
+	///int site = p_length - last_high - 1;
+	int site = t_length - 1;
+	int return_site = -1;
+	if ((err <= errthold) && (err<=*return_err))
+	{
+		*return_err = err;
+		return_site = site;
+	}
+	int i_last = i;
+	i = 0;
+
+
+
+
+	while (i<errthold)
+	{
+		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;
+		}
+	}
+
+
+	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;
+	}
+	
+
+
+	///end_site是正确的
+	int end_site = return_site;
+	int start_site = end_site;
+	///这个是各个bit-vector里面，end_site对应bit所在的位置
+	int back_track_site = band_length - (p_length - end_site);
+
+	Word v_value, h_value, delta_value, min_value, current_value;
+	Word direction, is_mismatch; ///0 is match, 1 is mismatch, 2 is up, 3 is left
+
+	///代表pattern到哪了，就是短的那个到哪了
+	i = t_length;
+	int path_length = 0;
+
+	///到0就结束了，后面的路径可以直接match
+	current_value = *return_err;
+
+
+	int low_bound = band_length - 1;
+
+
+	while (i>0)
+	{
+		if (current_value == 0)
+		{
+			break;
+		}
+
+		column_start = i << 3;
+
+		delta_value = current_value -
+			((~(matrix_bit[column_start] >> back_track_site))&err_mask);
+
+
+		if (back_track_site == 0)
+		{
+			///HP
+			h_value = current_value - ((matrix_bit[column_start + 3] >> back_track_site)&err_mask);
+			//HN
+			h_value = h_value + ((matrix_bit[column_start + 4] >> back_track_site)&err_mask);
+
+
+			min_value = delta_value;
+			direction = 0;
+
+			if (h_value < min_value)
+			{
+				min_value = h_value;
+				direction = 3;
+			}
+
+		}
+		else if (back_track_site == low_bound)
+		{
+			v_value = current_value - ((matrix_bit[column_start + 1] >> (back_track_site - 1))&err_mask);
+			v_value = v_value + ((matrix_bit[column_start + 2] >> (back_track_site - 1))&err_mask);
+
+			min_value = delta_value;
+			direction = 0;
+
+			if (v_value < min_value)
+			{
+				min_value = v_value;
+				direction = 2;
+			}
+
+		}
+		else
+		{
+
+			h_value = current_value - ((matrix_bit[column_start + 3] >> back_track_site)&(Word)1);
+
+
+			h_value = h_value + ((matrix_bit[column_start + 4] >> back_track_site)&(Word)1);
+
+
+			v_value = current_value - ((matrix_bit[column_start + 1] >> (back_track_site - 1))&err_mask);
+			v_value = v_value + ((matrix_bit[column_start + 2] >> (back_track_site - 1))&err_mask);
+
+
+			min_value = delta_value;
+			direction = 0;
+
+			if (v_value < min_value)
+			{
+				min_value = v_value;
+				direction = 2;
+			}
+
+
+			if (h_value < min_value)
+			{
+				min_value = h_value;
+				direction = 3;
+			}
+
+		}
+
+
+		if (direction == 0)
+		{
+
+			if (delta_value != current_value)
+			{
+				direction = 1;
+			}
+
+			i--;
+
+			start_site--;
+
+		}
+		if (direction == 2)///ru guo xiang shang yi dong, bing bu huan lie
+		{
+			back_track_site--;
+			start_site--;
+		}
+		else if (direction == 3)///ru guo xiang zuo yi dong
+		{
+			i--;
+			back_track_site++;
+		}
+
+
+		path[path_length++] = direction;
+
+
+		current_value = min_value;
+
+	}
+
+
+	if (i > 0)
+	{
+		memset(path + path_length, 0, i);
+		start_site = start_site - i;
+		direction = 0;
+		path_length = path_length + i;
+	}
+
+	if (direction != 3)
+	{
+		start_site++;
+	}
+	(*return_start_site) = start_site;
+	(*return_path_length) = path_length;
+
+
+
+
+	return return_site;
+
+}
+
 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)