M50 improved

Assembly.cpp

@@ -3321,6 +3321,132 @@ void Output_PAF()
 }
 
 
+int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold)
+{
+    long long i, k;
+    uint32_t qn, tn;
+    long long T_edges, A_edges;
+    T_edges = A_edges = 0;
+    for (i = 0; i < listLen; i++)
+    {
+        qn = (uint32_t)list[i];
+        for (k = i + 1; k < listLen; k++)
+        {
+            tn = (uint32_t)list[k];
+            if(get_specific_overlap(&(paf[qn]), qn, tn) != -1)
+            {
+                A_edges++;
+            }
+            
+            if(get_specific_overlap(&(paf[tn]), tn, qn) != -1)
+            {
+                A_edges++;
+            }
+
+            T_edges = T_edges + 2;
+        }
+        
+    }
+
+    if(A_edges >= (T_edges*threshold))
+    {
+        return 1;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+
+
+void rescue_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, 
+long long readNum, long long rescue_threshold, float cluster_threshold)
+{
+    double startTime = Get_T();
+    long long i, j, revises = 0;
+    uint32_t qn, tn;
+
+    kvec_t(uint64_t) edge_vector;
+    kv_init(edge_vector);
+    kvec_t(uint64_t) edge_vector_index;
+    kv_init(edge_vector_index);
+    uint64_t flag;
+    int index;
+    
+    for (i = 0; i < readNum; i++)
+    {
+        edge_vector.n = 0;
+        edge_vector_index.n = 0;
+        for (j = 0; j < paf[i].length; j++)
+        {
+            qn = Get_qn(paf[i].buffer[j]);
+            tn = Get_tn(paf[i].buffer[j]);
+            index = get_specific_overlap(&(rev_paf[tn]), tn, qn);
+            if(index != -1)
+            {
+                flag = tn;
+                flag = flag << 32;
+                flag = flag | (uint64_t)(index);
+                kv_push(uint64_t, edge_vector, flag);
+                kv_push(uint64_t, edge_vector_index, j);
+            }
+        }
+
+        ///the read itself has these overlaps, but all related reads do not have 
+        ///we need to remove all overlaps from paf[i], and then add all overlaps to rev_paf[i]
+        if(edge_vector.n >= rescue_threshold && 
+        check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1)
+        {
+            fprintf(stderr,"\nremove following %d edges...\n", edge_vector.n);
+            print_revise_edges(&(paf[i]), edge_vector_index.a, edge_vector_index.n);
+
+            add_overlaps(&(paf[i]), &(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
+            remove_overlaps(&(paf[i]), edge_vector_index.a, edge_vector_index.n);
+            revises = revises + edge_vector.n;
+        }
+
+        edge_vector.n = 0;
+        edge_vector_index.n = 0;
+        for (j = 0; j < rev_paf[i].length; j++)
+        {
+            qn = Get_qn(rev_paf[i].buffer[j]);
+            tn = Get_tn(rev_paf[i].buffer[j]);
+            index = get_specific_overlap(&(paf[tn]), tn, qn);
+            if(index != -1)
+            {
+                flag = tn;
+                flag = flag << 32;
+                flag = flag | (uint64_t)(index);
+                kv_push(uint64_t, edge_vector, flag);
+                kv_push(uint64_t, edge_vector_index, j);
+            }
+        }
+
+        ///the read itself do not have these overlaps, but all related reads have
+        ///we need to remove all overlaps from rev_paf[i], and then add all overlaps to paf[i]
+        if(edge_vector.n >= rescue_threshold && 
+        check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1)
+        {
+            fprintf(stderr,"\nadd following %d edges...\n", edge_vector.n);
+            print_revise_edges(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
+
+
+            remove_overlaps(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n);
+            add_overlaps_from_different_sources(paf, &(paf[i]), edge_vector.a, edge_vector.n);
+            revises = revises + edge_vector.n;
+        }
+    }
+
+
+    kv_destroy(edge_vector);
+    kv_destroy(edge_vector_index);
+
+    fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises);
+}
+
+
 
 void generate_overlaps(int last_round)
 {
@@ -3350,6 +3476,8 @@ void generate_overlaps(int last_round)
         pthread_join(_r_threads[i], NULL);
     free(_r_threads);
 
+    ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985);
+
     fprintf(stdout, "Final overlaps have been calculated.\n");
     fprintf(stdout, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); 
 
@@ -3369,6 +3497,8 @@ void generate_overlaps(int last_round)
 
 
 
+
+
 void Correct_Reads(int last_round)
 {
     

Correct.h

@@ -1447,6 +1447,10 @@ void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos
 uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy);
 /**********************for prefilter************************ */
 
+int verify_single_window(long long x_start, long long x_end, 
+long long overlap_x_s, long long overlap_y_s, int x_id,
+int y_id, int y_strand, char* x_buffer, char* y_buffer,
+All_reads* R_INF);
 
 void init_Cigar_record_alloc(Cigar_record_alloc* x);
 void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size);

Hash_Table.h

@@ -28,7 +28,7 @@ typedef khash_t(POS64) Pos_Table;
 #define OVERLAP_THRESHOLD_FILTER 0.9
 #define WINDOW_MAX_SIZE WINDOW + TAIL_LENGTH + 3
 #define THRESHOLD_MAX_SIZE  31
-
+#define FINAL_OVERLAP_ERROR_RATE 0.03
 
 #define GROUP_SIZE 4
 ///最长是10M10D10M10D10M这种

Overlaps.h

@@ -330,4 +330,11 @@ int write);
 void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, int id, char* command);
 
+void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum);
+
+void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen);
+void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
+void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
+uint64_t* source_index, long long listLen);
+void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
 #endif

Process_Read.cpp

@@ -329,6 +329,27 @@ void destory_UC_Read(UC_Read* r)
 	free(r->seq);
 }
 
+void init_aux_table()
+{
+	if (bit_t_seq_table[0][0] == 0)
+	{
+		uint64_t i = 0;
+
+		for (i = 0; i < 256; i++)
+		{
+			bit_t_seq_table[i][0] = s_H[((i >> 6)&(uint64_t)3)];
+			bit_t_seq_table[i][1] = s_H[((i >> 4)&(uint64_t)3)];
+			bit_t_seq_table[i][2] = s_H[((i >> 2)&(uint64_t)3)];
+			bit_t_seq_table[i][3] = s_H[(i&(uint64_t)3)];
+
+			bit_t_seq_table_rc[i][0] = RC_CHAR(bit_t_seq_table[i][3]);
+			bit_t_seq_table_rc[i][1] = RC_CHAR(bit_t_seq_table[i][2]);
+			bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]);
+			bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]);
+		}
+		
+	}
+}
 
 void init_UC_Read(UC_Read* r)
 {

Process_Read.h

@@ -56,7 +56,7 @@ static char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'};
 
 #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]]
 
-
+void init_aux_table();
 void init_kseq(char* file);
 void destory_kseq();
 int get_read(kseq_t *s, int adapterLen);