Merge overlap

Assembly.cpp

@@ -33,6 +33,8 @@ void* Perform_Counting(void* arg)
 
     uint64_t code;
 
+    uint64_t end_pos;
+
     Hash_code k_code;
 
     int avalible_k = 0;
@@ -54,7 +56,7 @@ void* Perform_Counting(void* arg)
 
             avalible_k = 0;
 
-            while ((code = get_HPC_code(&HPC_read)) != 6)
+            while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
             {
                 if(code < 4)
                 {
@@ -90,7 +92,7 @@ void* Perform_Counting(void* arg)
 
             avalible_k = 0;
 
-            while ((code = get_HPC_code(&HPC_read)) != 6)
+            while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
             {
                 if(code < 4)
                 {
@@ -148,8 +150,9 @@ void* Build_hash_table(void* arg)
 	int file_flag = 1;
 
     uint64_t code;
+    uint64_t end_pos;
 
-    long long HPC_base;
+    ///long long HPC_base;
 
     Hash_code k_code;
 
@@ -171,9 +174,9 @@ void* Build_hash_table(void* arg)
 
             avalible_k = 0;
 
-            HPC_base = 0;
+            ///HPC_base = 0;
 
-            while ((code = get_HPC_code(&HPC_read)) != 6)
+            while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
             {
                 if(code < 4)
                 {
@@ -189,7 +192,7 @@ void* Build_hash_table(void* arg)
                        ///所以如果当前k-mer在PCB表中存在，则他的位置一定要加入到候选位置中去
                         ///insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1, FORWARD);
                         insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, 
-                            curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1);
+                            curr_sub_block.read[i].ID, end_pos);
                     }
                     
                 }
@@ -199,7 +202,7 @@ void* Build_hash_table(void* arg)
                     init_Hash_code(&k_code);
                 }
 
-                HPC_base++;
+                ///HPC_base++;
             }
 
             
@@ -222,7 +225,7 @@ void* Build_hash_table(void* arg)
 
             HPC_base = 0;
 
-            while ((code = get_HPC_code(&HPC_read)) != 6)
+            while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
             {
                 if(code < 4)
                 {
@@ -411,6 +414,145 @@ void Build_hash_table_multiple_thr()
     
 }
 
+
+void* Overlap_calculate(void* arg)
+{
+
+    int thr_ID = *((int*)arg);
+
+    long long i = 0;
+    int avalible_k = 0;
+
+    UC_Read g_read;
+    init_UC_Read(&g_read);
+    HPC_seq HPC_read;
+    Hash_code k_code;
+    uint64_t code;
+    uint64_t end_pos;
+    k_mer_pos* list;
+    uint64_t list_length;
+    uint64_t sub_ID;
+    Candidates_list l;
+
+    k_mer_pos_list* merge_list;
+
+    init_Candidates_list(&l);
+
+    for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
+    ///for (i = thr_ID; i < R_INF.total_reads/50; i = i + thread_num)
+    {
+        /**
+        if (i % 1000 == 0)
+        {
+            fprintf(stderr, "i: %llu\n", i);
+        }
+        **/
+        
+        
+
+        clear_Candidates_list(&l);
+        recover_UC_Read(&g_read, &R_INF, i);
+
+
+        ///forward strand
+        init_HPC_seq(&HPC_read, g_read.seq, g_read.length);
+        init_Hash_code(&k_code);
+        avalible_k = 0;
+
+        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
+        {
+            if(code < 4)
+            {
+                k_mer_append(&k_code,code,k_mer_length);
+                avalible_k++;
+                if (avalible_k>=k_mer_length)
+                {
+                    list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID);
+                    merge_Candidates_list(&l, list, list_length, end_pos, 0);
+                    //merge_Candidates_list_version(&l, list, list_length, end_pos, 0);
+                }
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+
+            ///HPC_base++;
+        }
+
+
+
+        ///reverse complement strand
+        reverse_complement(g_read.seq, g_read.length);
+        init_HPC_seq(&HPC_read, g_read.seq, g_read.length);
+        init_Hash_code(&k_code);
+        avalible_k = 0;
+
+        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
+        {
+            if(code < 4)
+            {
+                k_mer_append(&k_code,code,k_mer_length);
+                avalible_k++;
+                if (avalible_k>=k_mer_length)
+                {
+                    list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID);
+                    merge_Candidates_list(&l, list, list_length, end_pos, 1);
+                    //merge_Candidates_list_version(&l, list, list_length, end_pos, 1);
+                }
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+
+            ///HPC_base++;
+        }
+
+    }
+    
+    destory_Candidates_list(&l);
+}
+
+
+void Overlap_calculate_multipe_thr()
+{
+    double start_time = Get_T();
+
+
+    fprintf(stdout, "Begin Overlap Calculate ...... \n");
+
+    pthread_t *_r_threads;
+
+	_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*thread_num);
+
+    int i = 0;
+
+	for (i = 0; i < thread_num; i++)
+	{
+		int *arg = (int*)malloc(sizeof(*arg));
+		*arg = i;
+
+		pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg);
+
+	}
+    
+
+    for (i = 0; i<thread_num; i++)
+		pthread_join(_r_threads[i], NULL);
+
+    free(_r_threads);
+
+    
+
+    fprintf(stdout, "Finish Overlap Calculate.\n");
+
+    fprintf(stdout, "%-30s%18.2f\n\n", "Calculate Overlap time:", Get_T() - start_time);    
+}
+
+
 int load_pre_cauculated_index()
 {
     if(load_Total_Pos_Table(&PCB, read_file_name) && load_All_reads(&R_INF, read_file_name))
@@ -437,6 +579,7 @@ void Verify_Counting()
 
     HPC_seq HPC_read;
     uint64_t code;
+    uint64_t end_pos;
     Hash_code k_code;
     int avalible_k = 0;
 
@@ -450,7 +593,7 @@ void Verify_Counting()
 
         avalible_k = 0;
 
-        while ((code = get_HPC_code(&HPC_read)) != 6)
+        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
         {
             if(code < 4)
             {
@@ -502,10 +645,11 @@ void verify_Position_hash_table()
 
 
     long long read_number = 0;
-    long long HPC_base;
+    ///long long HPC_base;
 
     HPC_seq HPC_read;
     uint64_t code;
+    uint64_t end_pos;
     Hash_code k_code;
     int avalible_k = 0;
     k_mer_pos* list;
@@ -554,9 +698,9 @@ void verify_Position_hash_table()
         init_Hash_code(&k_code);
 
         avalible_k = 0;
-        HPC_base = 0;
+        ///HPC_base = 0;
 
-        while ((code = get_HPC_code(&HPC_read)) != 6)
+        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
         {
 
 
@@ -601,7 +745,7 @@ void verify_Position_hash_table()
                         {
                             for(j=0; j<count2; j++)
                             {   
-                                if(list[j].readID == read_number && list[j].offset == HPC_base - k_mer_length + 1)
+                                if(list[j].readID == read_number && list[j].offset == end_pos)
                                 {
                                     break;
                                 }
@@ -610,7 +754,7 @@ void verify_Position_hash_table()
                             if(j == count2)
                             {
                                 fprintf(stderr, "locate error, read_number: %llu, pos: %llu\n", 
-                                read_number, HPC_base - k_mer_length + 1);
+                                read_number, end_pos);
 
                                 for(j=0; j<count2; j++)
                                 {   
@@ -637,7 +781,7 @@ void verify_Position_hash_table()
                 init_Hash_code(&k_code);
             }
 
-            HPC_base++;
+            ///HPC_base++;
             
         }
 

Assembly.h

@@ -9,7 +9,7 @@
 void Counting_multiple_thr();
 void Build_hash_table_multiple_thr();
 int load_pre_cauculated_index();
-
+void Overlap_calculate_multipe_thr();
 /********************************for debug***************************************/
 void Verify_Counting();
 

Hash_Table.cpp

@@ -4,6 +4,39 @@
 #include "Hash_Table.h"
 
 
+
+
+void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list)
+{
+    list->size = 1000;
+    list->length = 0;
+    list->list = (k_mer_pos_list*)malloc(sizeof(k_mer_pos_list)*list->size);
+}
+
+void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list)
+{
+    free(list->list);
+}
+
+void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list)
+{
+    list->length = 0;
+}
+
+
+void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length)
+{
+    list->length++;
+    if (list->length > list->size)
+    {
+        list->size = list->size * 2;
+        list->list = (k_mer_pos_list*)realloc(list->list, sizeof(k_mer_pos_list)*list->size);
+    }
+    
+    list->list[list->length - 1].list = n_list;
+    list->list[list->length - 1].length = n_length;
+}
+
 void init_Count_Table(Count_Table** table)
 {
     *table = kh_init(COUNT64);
@@ -300,6 +333,386 @@ int cmp_k_mer_pos(const void * a, const void * b)
 }
 
 
+void init_Candidates_list(Candidates_list* l)
+{
+    l->length = 0;
+    l->size = 0;
+    l->list = NULL;
+    l->tmp = NULL;
+    l->foward_pos = 0;
+    l->rc_pos = 0;
+}
+
+void clear_Candidates_list(Candidates_list* l)
+{
+    l->length = 0;
+    l->foward_pos = 0;
+    l->rc_pos = 0;
+}
+
+void destory_Candidates_list(Candidates_list* l)
+{
+    free(l->list);
+    free(l->tmp);
+}
+
+///1是x小，2是y小，0是相等
+inline int cmp_pos(k_mer_hit* x, k_mer_pos* y, uint64_t y_pos, uint8_t strand)
+{
+    if (x->strand < strand)
+    {
+        return 1;
+    }
+    else if (x->strand > strand)
+    {
+        return 2;
+    }
+    else
+    {
+        if (x->readID < y->readID)
+        {
+            return 1;
+        }
+        else if (x->readID > y->readID)
+        {
+            return 2;
+        }
+        else
+        {
+            if (x->offset < y->offset)
+            {
+                return 1;
+            }
+            else if (x->offset > y->offset)
+            {
+                return 2;
+            }
+            else
+            {
+                if (x->self_offset < y_pos)
+                {
+                    return 1;
+                }
+                else  if (x->self_offset > y_pos)
+                {
+                    return 2;
+                }
+                else
+                {
+                    return 0;
+                }
+                
+            }   
+        }    
+    }
+    
+}
+
+void debug_merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t i,
+uint64_t i_1, uint64_t i_2, uint64_t len1, uint64_t len2, uint64_t strand)
+{
+    if (i != l->length)
+    {
+        fprintf(stderr, "ERROR\n");
+        fprintf(stderr, "i: %llu, l->length: %llu\n", i , l->length);
+        fprintf(stderr, "i_1: %llu, len1: %llu\n", i_1 , len1);
+        fprintf(stderr, "i_2: %llu, len2: %llu\n", i_2 , len2);
+        fprintf(stderr, "strand: %llu\n", strand);
+        
+    }
+
+    for (i = 1; i < l->length; i++)
+    {
+        if (l->list[i].strand < l->list[i-1].strand)
+        {
+            fprintf(stderr, "ERROR -1\n");
+        }
+        else if (l->list[i].strand == l->list[i-1].strand)
+        {
+            if (l->list[i].readID < l->list[i-1].readID)
+            {
+                fprintf(stderr, "ERROR 0\n");
+            }
+            else if (l->list[i].readID == l->list[i-1].readID)
+            {
+                if (l->list[i].offset < l->list[i-1].offset)
+                {
+                    fprintf(stderr, "ERROR 1\n");
+                }
+                else if (l->list[i].offset == l->list[i-1].offset)
+                {
+                    if (l->list[i].self_offset < l->list[i-1].self_offset)
+                    {
+                        fprintf(stderr, "ERROR 2\n");
+                    }
+                }
+            } 
+        }
+    }
+
+    for (i = 0; i < n_lengh; i++)
+    {
+        int j = 0;
+        for (j = 0; j < l->length; j++)
+        {
+            if (
+                n_list[i].offset == l->list[j].offset
+                &&
+                n_list[i].readID == l->list[j].readID
+                &&
+                end_pos == l->list[j].self_offset
+                &&
+                strand == l->list[j].strand
+            )
+            {
+                break;
+            }   
+        }
+
+        if (j == l->length)
+        {
+            fprintf(stderr, "ERROR 4\n");
+        }
+    }
+
+
+    for (i = 0; i < l->length - n_lengh; i++)
+    {
+        int j = 0;
+        for (j = 0; j < l->length; j++)
+        {
+            if (
+                l->tmp[i].offset == l->list[j].offset
+                &&
+                l->tmp[i].readID == l->list[j].readID
+                &&
+                l->tmp[i].self_offset == l->list[j].self_offset
+                &&
+                l->tmp[i].strand == l->list[j].strand
+            )
+            {
+                break;
+            }   
+        }
+
+        if (j == l->length)
+        {
+            fprintf(stderr, "ERROR 5\n");
+        }
+    }
+}
+
+void merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand)
+{
+    if (n_lengh == 0) ///不加这个就会巨慢
+    {
+        return;
+    }
+    
+
+    if(l->length + n_lengh > l->size)
+    {
+        l->size = l->length + n_lengh;
+        l->list = (k_mer_hit*)realloc(l->list, sizeof(k_mer_hit)*l->size);
+        l->tmp = (k_mer_hit*)realloc(l->tmp, sizeof(k_mer_hit)*l->size);
+    }
+
+    int flag;
+    uint64_t len1, len2, i_1, i_2, i;
+    k_mer_hit *l_1;
+    k_mer_pos *l_2;
+
+    l_1 = l->list;
+    if (strand == 0)
+    {
+        ///由于标号全是0，则标号为1的位置不受影响(数组后半部分)，可以直接复制过去
+        ///需要从l_1 + l->rc_pos开始，复制l->length-l->rc_pos个元素，到l->tmp + l->rc_pos + n_lengh处
+        memcpy(l->tmp + l->rc_pos + n_lengh, l_1 + l->rc_pos, sizeof(k_mer_hit)*(l->length-l->rc_pos));
+        ///l_1只需要从0扫描到l->rc_pos就好了
+        len1 = l->rc_pos; 
+        i_1 = 0;
+        ///此时目标数组没数据，
+        i = 0;
+        ///更新l->rc_pos
+        l->rc_pos = l->rc_pos + n_lengh;
+    }
+    else
+    {
+        ///如果strand==1，则标号为0的位置不受影响(数组前半部分)，可以直接复制过去
+        ///从l_1开始，复制l->rc_pos个元素，到l->tmp
+        memcpy(l->tmp, l_1, sizeof(k_mer_hit)*l->rc_pos);
+        ///l_1只需要从l->rc_pos扫描到末尾就好了
+        i_1 = l->rc_pos;
+        len1 = l->length; 
+        ///此时目标数组已经有了l->rc_pos
+        i = i_1;
+        ///l->rc_pos不用更新
+    }
+    
+    
+    
+
+    l_2 = n_list;
+    len2 = n_lengh;
+    i_2 = 0;
+
+
+    while (i_1<len1 && i_2<len2)
+    {
+        flag = cmp_pos(&l_1[i_1], &l_2[i_2], end_pos, strand); 
+        ///l_1小
+        if (flag == 1)
+        {
+            l->tmp[i].readID = l_1[i_1].readID;
+            l->tmp[i].offset = l_1[i_1].offset;
+            l->tmp[i].self_offset = l_1[i_1].self_offset;
+            l->tmp[i].strand = l_1[i_1].strand;
+            i_1++;
+            i++;
+        }
+        else  ///l_2小或者l_1 == l_2
+        {
+            l->tmp[i].readID = l_2[i_2].readID;
+            l->tmp[i].offset = l_2[i_2].offset;
+            l->tmp[i].self_offset = end_pos;
+            l->tmp[i].strand = strand;
+            i_2++;
+            i++;
+        }        
+    }
+
+
+    while (i_1<len1)
+    {
+        memcpy(l->tmp + i, l_1 + i_1, sizeof(k_mer_hit)*(len1-i_1));
+        i = i + len1-i_1;
+        i_1 = len1;
+    }
+
+    while (i_2<len2)
+    {
+        l->tmp[i].readID = l_2[i_2].readID;
+        l->tmp[i].offset = l_2[i_2].offset;
+        l->tmp[i].self_offset = end_pos;
+        l->tmp[i].strand = strand;
+        i_2++;
+        i++;
+    }
+    
+    
+
+    k_mer_hit* k;
+    k = l->list;
+    l->list = l->tmp;
+    l->tmp = k;
+    l->length = l->length + n_lengh;
+
+    
+    /**
+    if (strand == 0)
+    {
+        i = l->length;
+        len1 = l->length - n_lengh;
+    }
+    debug_merge_Candidates_list(l, n_list, n_lengh, end_pos, i, i_1, i_2, len1, len2, strand);
+    **/
+    
+}
+
+
+
+
+
+
+
+/********************************for debug***************************************/
+void merge_Candidates_list_version(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand)
+{
+    if (n_lengh == 0) ///不加这个就会巨慢
+    {
+        return;
+    }
+    
+
+    if(l->length + n_lengh > l->size)
+    {
+        l->size = l->length + n_lengh;
+        l->list = (k_mer_hit*)realloc(l->list, sizeof(k_mer_hit)*l->size);
+        l->tmp = (k_mer_hit*)realloc(l->tmp, sizeof(k_mer_hit)*l->size);
+    }
+
+    int flag;
+    uint64_t len1, len2, i_1, i_2, i;
+    k_mer_hit *l_1;
+    k_mer_pos *l_2;
+
+    l_1 = l->list;
+    len1 = l->length; 
+    i_1 = 0;
+
+    l_2 = n_list;
+    len2 = n_lengh;
+    i_2 = 0;
+
+    i = 0;
+
+    while (i_1<len1 && i_2<len2)
+    {
+        flag = cmp_pos(&l_1[i_1], &l_2[i_2], end_pos, strand); 
+        ///l_1小
+        if (flag == 1)
+        {
+            l->tmp[i].readID = l_1[i_1].readID;
+            l->tmp[i].offset = l_1[i_1].offset;
+            l->tmp[i].self_offset = l_1[i_1].self_offset;
+            l->tmp[i].strand = l_1[i_1].strand;
+            i_1++;
+            i++;
+        }
+        else  ///l_2小或者l_1 == l_2
+        {
+            l->tmp[i].readID = l_2[i_2].readID;
+            l->tmp[i].offset = l_2[i_2].offset;
+            l->tmp[i].self_offset = end_pos;
+            l->tmp[i].strand = strand;
+            i_2++;
+            i++;
+        }        
+    }
+
+
+    while (i_1<len1)
+    {
+        l->tmp[i].readID = l_1[i_1].readID;
+        l->tmp[i].offset = l_1[i_1].offset;
+        l->tmp[i].self_offset = l_1[i_1].self_offset;
+        l->tmp[i].strand = l_1[i_1].strand;
+        i_1++;
+        i++;
+    }
+
+    while (i_2<len2)
+    {
+        l->tmp[i].readID = l_2[i_2].readID;
+        l->tmp[i].offset = l_2[i_2].offset;
+        l->tmp[i].self_offset = end_pos;
+        l->tmp[i].strand = strand;
+        i_2++;
+        i++;
+    }
+    
+    k_mer_hit* k;
+    k = l->list;
+    l->list = l->tmp;
+    l->tmp = k;
+    l->length = l->length + n_lengh;
+
+    ///debug_merge_Candidates_list(l, n_list, n_lengh, end_pos, i, i_1, i_2, len1, len2, strand);
+}
+
+
+
+
 
 
 

Hash_Table.h

@@ -37,6 +37,39 @@ typedef struct
     uint64_t readID;
 } k_mer_pos;
 
+typedef struct
+{
+    k_mer_pos* list;
+    uint64_t length;
+} k_mer_pos_list;
+
+typedef struct
+{
+    k_mer_pos_list* list;
+    uint64_t size;
+    uint64_t length;
+} k_mer_pos_list_alloc;
+
+
+
+typedef struct
+{
+    uint64_t offset;
+    uint64_t readID;
+    uint64_t self_offset;
+    uint8_t strand;
+} k_mer_hit;
+
+typedef struct
+{
+    k_mer_hit* list;
+    k_mer_hit* tmp;
+    long long length;
+    long long size;
+    uint64_t foward_pos;
+    uint64_t rc_pos;
+} Candidates_list;
+
 typedef struct
 {
 	Pos_Table** sub_h;
@@ -273,8 +306,8 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in
         {
             qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos);
         }
-        
 
+        
         return 1;
     }
     else
@@ -300,6 +333,10 @@ int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name);
 
 void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq);
 
+void init_Candidates_list(Candidates_list* l);
+void clear_Candidates_list(Candidates_list* l);
+void destory_Candidates_list(Candidates_list* l);
+void merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand);
 
 
 
@@ -429,4 +466,8 @@ void test_COUNT64();
 void debug_mode(uint64_t d, uint64_t thread_ID, uint64_t thread_num);
 
 
+
+/********************************for debug***************************************/
+void merge_Candidates_list_version(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand);
+
 #endif

kmer.h

@@ -33,13 +33,15 @@ typedef struct
 } HPC_seq;
 
 
-inline uint64_t get_HPC_code(HPC_seq* seq)
+inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos)
 {
 
     if(seq->i < seq ->l)
     {
         uint8_t code = seq_nt6_table[(uint8_t)seq->str[seq->i]];
 
+        (*end_pos) = seq->i;
+
         for (; seq->i < seq->l; seq->i++)
         {
             ///统计N的个数

main.cpp

@@ -35,11 +35,9 @@ int main(int argc, char *argv[])
         Build_hash_table_multiple_thr();
     }
     
-
-    
-
     ///verify_Position_hash_table();
 
+    Overlap_calculate_multipe_thr();