xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-04 13:38:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

phasing

Browse Source
This commit is contained in:
Haoyu Cheng
2019-09-02 10:30:14 -04:00
parent 3c0c201cdf
commit d5327fe9ee
15 changed files with 3396 additions and 656 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.vscode/settings.json
Vendored
+6 -1
View File
@@ -10,6 +10,11 @@
"*.tcc": "cpp",
"bitset": "cpp",
"algorithm": "cpp",
"hashtable": "cpp"
"hashtable": "cpp",
"string_view": "cpp",
"list": "cpp",
"string": "cpp",
"array": "cpp",
"utility": "cpp"
}
}
Assembly.cpp
+724 -62
View File
@@ -21,6 +21,7 @@ long long total_potiental_matched_overlap_0 = 0;
long long total_potiental_matched_overlap_1 = 0;
long long total_num_read_base = 0;
long long total_num_correct_base = 0;
int roundID = 0;
long long complete_threads = 0;
@@ -139,12 +140,146 @@ void* Perform_Counting(void* arg)
}
destory_R_buffer_block(&curr_sub_block);
free(arg);
///free(arg);
}
void* Perform_Counting_non_first(void* arg)
{
int thr_ID = *((int*)arg);
int i = 0;
HPC_seq HPC_read;
long long read_number = 0;
long long select_k_mer_number = 0 ;
long long k_mer_number = 0 ;
int file_flag = 1;
uint64_t code;
uint64_t end_pos;
Hash_code k_code;
int avalible_k = 0;
UC_Read g_read;
init_UC_Read(&g_read);
for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
{
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)
{
///插入
if(insert_Total_Count_Table(&TCB, &k_code, k_mer_length))
{
select_k_mer_number++;
}
k_mer_number++;
}
}
else
{
avalible_k = 0;
init_Hash_code(&k_code);
}
}
}
destory_UC_Read(&g_read);
}
void* Build_hash_table_non_first(void* arg)
{
int thr_ID = *((int*)arg);
int i = 0;
HPC_seq HPC_read;
int file_flag = 1;
uint64_t code;
uint64_t end_pos;
///long long HPC_base;
Hash_code k_code;
int avalible_k = 0;
UC_Read g_read;
init_UC_Read(&g_read);
for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
{
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;
///HPC_base = 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)
{
///选取的k-mer满足两个要求
///1. hash(k-mer) % 101 <= 3
///2. occ(k-mer)要满足范围
///TCB表中的元素仅满足第一个要求，而PCB表中的元素满足两个要求
///所以如果当前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,
i, end_pos);
}
}
else
{
avalible_k = 0;
init_Hash_code(&k_code);
}
///HPC_base++;
}
}
destory_UC_Read(&g_read);
}
@@ -275,7 +410,7 @@ void* Build_hash_table(void* arg)
}
destory_R_buffer_block(&curr_sub_block);
free(arg);
///free(arg);
@@ -296,20 +431,21 @@ void Counting_multiple_thr()
fprintf(stdout, "Begin Counting ...... \n");
init_kseq(read_file_name);
init_All_reads(&R_INF);
init_Total_Count_Table(k_mer_length, &TCB);
pthread_t inputReadsHandle;
init_R_buffer(thread_num);
int *is_insert = (int*)malloc(sizeof(*is_insert));
*is_insert = 1;
pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
if (roundID == 0)
{
init_kseq(read_file_name);
init_All_reads(&R_INF);
init_R_buffer(thread_num);
pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
}
pthread_t *_r_threads;
@@ -321,14 +457,19 @@ void Counting_multiple_thr()
{
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
pthread_create(_r_threads + i, NULL, Perform_Counting, (void*)arg);
if (roundID == 0)
{
pthread_create(_r_threads + i, NULL, Perform_Counting, (void*)arg);
}
else
{
pthread_create(_r_threads + i, NULL, Perform_Counting_non_first, (void*)arg);
}
}
pthread_join(inputReadsHandle, NULL);
for (i = 0; i<thread_num; i++)
pthread_join(_r_threads[i], NULL);
@@ -338,8 +479,11 @@ void Counting_multiple_thr()
///destory_R_buffer();
///destory_Total_Count_Table(&TCB);
destory_kseq();
if (roundID == 0)
{
pthread_join(inputReadsHandle, NULL);
destory_kseq();
}
fprintf(stdout, "Finish Counting ...... \n");
@@ -368,21 +512,24 @@ void Build_hash_table_multiple_thr()
fprintf(stdout, "%-30s%18.2f\n\n", "Traverse time:", Get_T() - T_start_time);
fflush(stdout);
init_kseq(read_file_name);
clear_R_buffer();
malloc_All_reads(&R_INF);
///at this moment, TCB can be free
destory_Total_Count_Table(&TCB);
pthread_t inputReadsHandle;
int *is_insert = (int*)malloc(sizeof(*is_insert));
*is_insert = 0;
pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
if (roundID == 0)
{
init_kseq(read_file_name);
clear_R_buffer();
malloc_All_reads(&R_INF);
pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
}
@@ -397,36 +544,49 @@ void Build_hash_table_multiple_thr()
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
pthread_create(_r_threads + i, NULL, Build_hash_table, (void*)arg);
if (roundID == 0)
{
pthread_create(_r_threads + i, NULL, Build_hash_table, (void*)arg);
}
else
{
pthread_create(_r_threads + i, NULL, Build_hash_table_non_first, (void*)arg);
}
}
pthread_join(inputReadsHandle, NULL);
if (roundID == 0)
{
pthread_join(inputReadsHandle, NULL);
}
for (i = 0; i<thread_num; i++)
pthread_join(_r_threads[i], NULL);
free(_r_threads);
destory_kseq();
fprintf(stdout, "Finish Building hash table ...... \n");
fprintf(stdout, "%-30s%18.2f\n\n", "Build hash table time:", Get_T() - start_time);
///destory_Total_Count_Table(&TCB);
if (write_index_to_disk)
if (roundID == 0)
{
write_Total_Pos_Table(&PCB, read_file_name);
///destory_Total_Pos_Table(&PCB);
///load_Total_Pos_Table(&PCB, read_file_name);
write_All_reads(&R_INF, read_file_name);
///destory_All_reads(&R_INF);
///load_All_reads(&R_INF, read_file_name);
destory_kseq();
destory_R_buffer();
///destory_Total_Count_Table(&TCB);
if (write_index_to_disk)
{
write_Total_Pos_Table(&PCB, read_file_name);
///destory_Total_Pos_Table(&PCB);
///load_Total_Pos_Table(&PCB, read_file_name);
write_All_reads(&R_INF, read_file_name);
///destory_All_reads(&R_INF);
///load_All_reads(&R_INF, read_file_name);
}
}
destory_Total_Count_Table(&TCB);
///destory_Total_Count_Table(&TCB);
free(is_insert);
@@ -787,6 +947,294 @@ int get_required_read(const char *required_name, long long RID, All_reads* R_INF
}
void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
char* new_read, int* new_length)
{
int i, j;
int pre_i, new_i;
int operation, operation_length;
pre_i = new_i = 0;
int diff_char_i = 0;
for (i = 0; i < cigar->length; i++)
{
operation = Get_Cigar_Type(cigar->record[i]);
operation_length = Get_Cigar_Length(cigar->record[i]);
if (operation == 0)
{
memcpy(new_read + new_i, pre_read + pre_i, operation_length);
pre_i = pre_i + operation_length;
new_i = new_i + operation_length;
}
else if (operation == 1)
{
for (j = 0; j < operation_length; j++)
{
new_read[new_i] = Get_MisMatch_Base(cigar->lost_base[diff_char_i]);
new_i++;
diff_char_i++;
}
pre_i = pre_i + operation_length;
/**
pre_i = pre_i + operation_length;
memcpy(new_read + new_i, cigar->lost_base + diff_char_i, operation_length);
new_i = new_i + operation_length;
diff_char_i = diff_char_i + operation_length;
**/
}
else if (operation == 3)
{
pre_i = pre_i + operation_length;
diff_char_i = diff_char_i + operation_length;
}
else if (operation == 2)
{
memcpy(new_read + new_i, cigar->lost_base + diff_char_i, operation_length);
new_i = new_i + operation_length;
diff_char_i = diff_char_i + operation_length;
}
}
*new_length = new_i;
}
void get_uncorrected_read_from_cigar(Cigar_record* cigar, char* new_read, int new_length, char* pre_read, int* pre_length)
{
int i, j;
int pre_i, new_i;
int operation, operation_length;
pre_i = new_i = 0;
int diff_char_i = 0;
for (i = 0; i < cigar->length; i++)
{
operation = Get_Cigar_Type(cigar->record[i]);
operation_length = Get_Cigar_Length(cigar->record[i]);
if (operation == 0)
{
memcpy(pre_read + pre_i, new_read + new_i, operation_length);
pre_i = pre_i + operation_length;
new_i = new_i + operation_length;
}
else if (operation == 1)
{
for (j = 0; j < operation_length; j++)
{
pre_read[pre_i] = Get_Match_Base(cigar->lost_base[diff_char_i]);
pre_i++;
diff_char_i++;
}
new_i = new_i + operation_length;
}
else if (operation == 3)
{
memcpy(pre_read + pre_i, cigar->lost_base + diff_char_i, operation_length);
pre_i = pre_i + operation_length;
diff_char_i = diff_char_i + operation_length;
}
else if (operation == 2)
{
new_i = new_i + operation_length;
diff_char_i = diff_char_i + operation_length;
}
}
*pre_length = pre_i;
}
int debug_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
char* new_read, int new_length, int correct_base)
{
int i;
int total_errors = 0;
for (i = 0; i < cigar->length; i++)
{
if (Get_Cigar_Type(cigar->record[i]) > 0)
{
total_errors = total_errors + Get_Cigar_Length(cigar->record[i]);
}
}
if(total_errors!=correct_base)
{
fprintf(stderr, "total_errors: %d, correct_base: %d\n", total_errors, correct_base);
}
int pre_i, new_i;
int operation, operation_length;
pre_i = new_i = 0;
for (i = 0; i < cigar->length; i++)
{
operation = Get_Cigar_Type(cigar->record[i]);
operation_length = Get_Cigar_Length(cigar->record[i]);
if (operation == 0)
{
pre_i = pre_i + operation_length;
new_i = new_i + operation_length;
}
if (operation == 1)
{
pre_i = pre_i + operation_length;
new_i = new_i + operation_length;
}
if (operation == 3)
{
pre_i = pre_i + operation_length;
}
if (operation == 2)
{
new_i = new_i + operation_length;
}
}
/**
fprintf(stderr, "total_errors: %d, correct_base: %d, length: %d, lost_base_length: %d\n",
total_errors, correct_base, cigar->length, cigar->lost_base_length);
**/
if (pre_i != pre_length || new_i != new_length)
{
fprintf(stderr, "pre_i: %d, pre_length: %d\n", pre_i, pre_length);
fprintf(stderr, "new_i: %d, new_length: %d\n", new_i, new_length);
}
char* tmp_seq = (char*)malloc(new_length + pre_length);
int tmp_length;
get_corrected_read_from_cigar(cigar, pre_read, pre_length, tmp_seq, &tmp_length);
if(tmp_length != new_length)
{
fprintf(stderr, "tmp_length: %d, new_length: %d\n", tmp_length, new_length);
}
if(memcmp(new_read, tmp_seq, new_length)!=0)
{
fprintf(stderr, "error new string\n");
}
get_uncorrected_read_from_cigar(cigar, new_read, new_length, tmp_seq, &tmp_length);
if(tmp_length != pre_length)
{
fprintf(stderr, "tmp_length: %d, pre_length: %d\n", tmp_length, pre_length);
}
if(memcmp(pre_read, tmp_seq, pre_length)!=0)
{
fprintf(stderr, "error pre string\n");
}
free(tmp_seq);
if(cigar->new_read_length != new_length)
{
fprintf(stderr, "cigar->new_read_length: %d, new_length: %d\n", cigar->new_read_length, new_length);
}
}
inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_record* input)
{
if (input->length > records[ID].size)
{
records[ID].size = input->length;
records[ID].record = (uint32_t*)realloc(records[ID].record, records[ID].size*sizeof(uint32_t));
}
records[ID].length = input->length;
memcpy(records[ID].record, input->record, input->length*sizeof(uint32_t));
if (input->lost_base_length > records[ID].lost_base_size)
{
records[ID].lost_base_size = input->lost_base_length;
records[ID].lost_base = (char*)realloc(records[ID].lost_base, records[ID].lost_base_size);
}
records[ID].lost_base_length = input->lost_base_length;
memcpy(records[ID].lost_base, input->lost_base, input->lost_base_length);
records[ID].new_length = input->new_read_length;
}
void just_debug(overlap_region_alloc* overlap_list, All_reads* R_INF)
{
int j;
int i;
int y_id;
int y_strand;
int y_readLen;
int overlap_length;
int high_quality_overlaps = 0;
UC_Read g_read;
init_UC_Read(&g_read);
for (j = 0; j < overlap_list->length; j++)
{
y_id = overlap_list->list[j].y_id;
y_strand = overlap_list->list[j].y_pos_strand;
y_readLen = Get_READ_LENGTH((*R_INF), y_id);
overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
if (overlap_length * OVERLAP_THRESHOLD <= overlap_list->list[j].align_length)
{
high_quality_overlaps++;
fprintf(stderr, "i: %d, overlap_length: %d, align_length: %d, y_strand: %d\n",
high_quality_overlaps, overlap_length, overlap_list->list[j].align_length, y_strand);
fprintf(stderr, "x_pos_s: %d, x_pos_e: %d\n",
overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e);
fprintf(stderr, "%.*s\n\n", Get_NAME_LENGTH((*R_INF), y_id), Get_NAME((*R_INF), y_id));
}
}
recover_UC_Read_RC(&g_read, R_INF, overlap_list->list[0].x_id);
for (i = 0; i < g_read.length && i < 81; i++)
{
fprintf(stderr, "%c", g_read.seq[i]);
}
fprintf(stderr, "\n");
fprintf(stderr, "x_length: %d\n", g_read.length);
fprintf(stderr, "high_quality_overlaps: %d, overlap_list->length: %d\n", high_quality_overlaps,overlap_list->length);
destory_UC_Read(&g_read);
}
void* Overlap_calculate_heap_merge(void* arg)
{
/************需要注释掉**********/
@@ -850,17 +1298,19 @@ void* Overlap_calculate_heap_merge(void* arg)
init_buffer_sub_block(&current_sub_buffer);
Cigar_record current_cigar;
init_Cigar_record(&current_cigar);
haplotype_evdience_alloc hap;
InitHaplotypeEvdience(&hap);
for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
{
/**
if(!get_required_read("m54238_180916_191625\/43253934\/ccs", i, &R_INF))
{
continue;
}
**/
clear_Cigar_record(&current_cigar);
clear_Heap(&heap);
clear_Candidates_list(&l);
///clear_Candidates_list(&debug_l);
@@ -996,14 +1446,38 @@ void* Overlap_calculate_heap_merge(void* arg)
///clear_Graph(&POA_Graph);
correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph,
&matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1);
&matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
&current_cigar, &hap);
num_read_base = num_read_base + g_read.length;
num_correct_base = num_correct_base + correct.corrected_base;
output_read_to_buffer(i, &R_INF, correct.corrected_read, correct.corrected_read_length, &current_sub_buffer);
///output_read_to_buffer(i, &R_INF, g_read.seq, g_read.length, &current_sub_buffer);
push_cigar(R_INF.cigars, i, &current_cigar);
/**
if(memcmp("m54334_180926_225337/39780640/ccs", Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i)) == 0)
{
fprintf(stderr, "i: %d\n", i);
///fprintf(stderr, "overlap_list.length: %d\n", overlap_list.length);
just_debug(&overlap_list, &R_INF);
}
**/
///output_read_to_buffer(i, &R_INF, correct.corrected_read, correct.corrected_read_length, &current_sub_buffer);
/**
debug_cigar(&current_cigar, g_read.seq, g_read.length, correct.corrected_read, correct.corrected_read_length,
correct.corrected_base);
**/
/**
POA_i = 0;
@@ -1109,10 +1583,12 @@ void* Overlap_calculate_heap_merge(void* arg)
destory_Graph(&POA_Graph);
destory_UC_Read(&g_read);
destory_UC_Read(&overlap_read);
destory_Cigar_record(&current_cigar);
destory_Correct_dumy(&correct);
destoryHaplotypeEvdience(&hap);
pthread_mutex_lock(&statistics);
total_matched_overlap_0 += matched_overlap_0;
@@ -1133,9 +1609,115 @@ void* Overlap_calculate_heap_merge(void* arg)
fprintf(stderr, "total_num_correct_base: %llu\n", total_num_correct_base);
}
pthread_mutex_unlock(&statistics);
free(arg);
}
void* Save_corrected_reads(void* arg)
{
int thr_ID = *((int*)arg);
long long i, j;
UC_Read g_read;
init_UC_Read(&g_read);
int new_read_size = 10000;
char* new_read = (char*)malloc(new_read_size);
Cigar_record cigar;
int new_read_length;
uint64_t N_occ;
for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
{
recover_UC_Read(&g_read, &R_INF, i);
if(R_INF.cigars[i].new_length>new_read_size)
{
new_read_size = R_INF.cigars[i].new_length;
new_read = (char*)realloc(new_read, new_read_size);
}
cigar.length = R_INF.cigars[i].length;
cigar.lost_base_length = R_INF.cigars[i].lost_base_length;
cigar.record = R_INF.cigars[i].record;
cigar.lost_base = R_INF.cigars[i].lost_base;
get_corrected_read_from_cigar(&cigar, g_read.seq, g_read.length, new_read, &new_read_length);
///need modification
reverse_complement(new_read, new_read_length);
N_occ = 0;
for (j = 0; j < new_read_length; j++)
{
if(new_read[j] == 'N')
{
N_occ++;
}
}
if(R_INF.read_size[i] < new_read_length)
{
R_INF.read_size[i] = new_read_length;
R_INF.read_sperate[i] = (uint8_t*)realloc(R_INF.read_sperate[i], R_INF.read_size[i]/4+1);
}
R_INF.read_length[i] = new_read_length;
compress_base(Get_READ(R_INF, i),
new_read, new_read_length,
&R_INF.N_site[i], N_occ);
}
destory_UC_Read(&g_read);
free(new_read);
free(arg);
}
void Output_corrected_reads()
{
long long i, j;
UC_Read g_read;
init_UC_Read(&g_read);
FILE* output_file = fopen(output_file_name, "w");
if(number_of_round % 2 == 0)
{
for (i = 0; i < R_INF.total_reads; i++)
{
recover_UC_Read(&g_read, &R_INF, i);
fwrite(">", 1, 1, output_file);
fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
fwrite("\n", 1, 1, output_file);
fwrite(g_read.seq, 1, g_read.length, output_file);
fwrite("\n", 1, 1, output_file);
}
}
else
{
for (i = 0; i < R_INF.total_reads; i++)
{
recover_UC_Read_RC(&g_read, &R_INF, i);
fwrite(">", 1, 1, output_file);
fwrite(Get_NAME(R_INF, i), 1, Get_NAME_LENGTH(R_INF, i), output_file);
fwrite("\n", 1, 1, output_file);
fwrite(g_read.seq, 1, g_read.length, output_file);
fwrite("\n", 1, 1, output_file);
}
}
destory_UC_Read(&g_read);
fclose(output_file);
}
@@ -1145,20 +1727,19 @@ void* Overlap_calculate_heap_merge(void* arg)
void Overlap_calculate_multipe_thr()
{
fprintf(stderr, "k_mer_min_freq: %lld, CORRECT_THRESHOLD: %f\n", k_mer_min_freq, CORRECT_THRESHOLD);
double start_time = Get_T();
init_output_buffer(thread_num);
pthread_t outputResultSinkHandle;
/**
if (roundID == number_of_round - 1)
{
init_output_buffer(thread_num);
pthread_create(&outputResultSinkHandle, NULL, pop_buffer, NULL);
}
**/
pthread_create(&outputResultSinkHandle, NULL, pop_buffer, NULL);
fprintf(stdout, "R_INF.total_reads: %llu\n", R_INF.total_reads);
fprintf(stdout, "Begin Overlap Calculate ...... \n");
@@ -1174,7 +1755,6 @@ void Overlap_calculate_multipe_thr()
*arg = i;
pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)arg);
//pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg);
}
@@ -1182,17 +1762,55 @@ void Overlap_calculate_multipe_thr()
for (i = 0; i<thread_num; i++)
pthread_join(_r_threads[i], NULL);
pthread_join(outputResultSinkHandle, NULL);
/**
if (roundID == number_of_round - 1)
{
pthread_join(outputResultSinkHandle, NULL);
destory_output_buffer();
///destory_All_reads(&R_INF);
}
**/
free(_r_threads);
destory_output_buffer();
destory_Total_Pos_Table(&PCB);
fprintf(stdout, "Finish Overlap Calculate.\n");
fprintf(stdout, "%-30s%18.2f\n\n", "Calculate Overlap time:", Get_T() - start_time);
start_time = Get_T();
_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*thread_num);
for (i = 0; i < thread_num; i++)
{
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)arg);
//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, "%-30s%18.2f\n\n", "Save corrected read time:", Get_T() - start_time);
///fflush(stdout);
///only the last round can output read to disk
if (roundID == number_of_round - 1)
{
start_time = Get_T();
Output_corrected_reads();
fprintf(stdout, "%-30s%18.2f\n\n", "Output time:", Get_T() - start_time);
}
}
@@ -1438,6 +2056,50 @@ void verify_Position_hash_table()
}
void Correct_Reads(int last_round)
{
complete_threads = 0;
total_matched_overlap_0 = 0;
total_matched_overlap_1 = 0;
total_potiental_matched_overlap_0 = 0;
total_potiental_matched_overlap_1 = 0;
total_num_read_base = 0;
total_num_correct_base = 0;
if(last_round == 0)
return;
roundID = number_of_round - last_round;
fprintf(stdout, "Error correction: start the %d-th round ...\n", roundID);
///only the first round correction can load index from disk
if (roundID == 0 & load_index_from_disk && load_pre_cauculated_index())
{
;
}
else
{
Counting_multiple_thr();
Build_hash_table_multiple_thr();
}
fprintf(stdout, "Total pos in hash tabe: %d\n", PCB.total_occ);
fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq);
fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq);
///verify_Position_hash_table();
Overlap_calculate_multipe_thr();
fprintf(stdout, "Error correction: the %d-th round has been completed.\n", roundID);
Correct_Reads(last_round - 1);
}
Assembly.h
+4 -1
View File
@@ -4,12 +4,15 @@
#define FORWARD 0
#define REVERSE_COMPLEMENT (0x8000000000000000)
#define Get_Cigar_Type(RECORD) (RECORD&3)
#define Get_Cigar_Length(RECORD) (RECORD>>2)
void Counting_multiple_thr();
void Build_hash_table_multiple_thr();
int load_pre_cauculated_index();
void Overlap_calculate_multipe_thr();
void Correct_Reads(int last_round);
/********************************for debug***************************************/
void Verify_Counting();
CommandLines.cpp
+4 -1
View File
@@ -15,6 +15,7 @@ int k_mer_min_freq = 3;
int k_mer_max_freq = 66;
int load_index_from_disk = 0;
int write_index_to_disk = 0;
int number_of_round = 1;
double Get_T(void)
@@ -42,19 +43,21 @@ int CommandLine_process (int argc, char *argv[])
{ "thread", ko_required_argument, 103},
{ "k_mer_min_freq", ko_required_argument, 104},
{ "k_mer_max_freq", ko_required_argument, 105},
{ "round", ko_required_argument, 106},
{ NULL, 0, 0 }
};
ketopt_t opt = KETOPT_INIT;
int i, c;
while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:", longopts)) >= 0) {
while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:r:", longopts)) >= 0) {
if (c == 100 || c == 'h') Print_H();
else if (c == 103 || c == 't') thread_num = atoi(opt.arg);
else if (c == 102 || c == 'o') output_file_name = opt.arg;
else if (c == 101 || c == 'q') read_file_name = opt.arg;
else if (c == 104 || c == 'n') k_mer_min_freq = atoi(opt.arg);
else if (c == 105 || c == 'm') k_mer_max_freq = atoi(opt.arg);
else if (c == 106 || c == 'r') number_of_round = atoi(opt.arg);
else if (c == 'k') k_mer_length = atoi(opt.arg);
else if (c == 'l') load_index_from_disk = 1;
else if (c == 'w') write_index_to_disk = 1;
CommandLines.h
+1
View File
@@ -12,6 +12,7 @@ extern int k_mer_min_freq;
extern int k_mer_max_freq;
extern int load_index_from_disk;
extern int write_index_to_disk;
extern int number_of_round;
int CommandLine_process (int argc, char *argv[]);
Correct.cpp
+1956 -565
View File
File diff suppressed because it is too large Load Diff
Correct.h
+502 -1
View File
@@ -4,6 +4,7 @@
#include "Hash_Table.h"
#include "Levenshtein_distance.h"
#include "POA.h"
#include "Process_Read.h"
#define CORRECT_THRESHOLD 0.7
#define MIN_COVERAGE_THRESHOLD 4
@@ -12,11 +13,324 @@
#define INSERTION 2
#define DELETION 3
#define FLAG_THRE 1
#define MAX(x, y) ((x >= y)?x:y)
#define MIN(x, y) ((x <= y)?x:y)
#define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1)
///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1
#define Get_MisMatch_Base(RECORD) (s_H[(RECORD>>3)])
#define Get_Match_Base(RECORD) (s_H[(RECORD&7)])
typedef struct
{
/**[0-1] bits are type:**/
/**[2-31] bits are length**/
char current_operation;
int current_operation_length;
uint32_t* record;
uint64_t size;
uint64_t length;
uint32_t new_read_length;
char* lost_base;
uint64_t lost_base_size;
uint64_t lost_base_length;
}Cigar_record;
typedef struct
{
////the position of snp in read itself
uint32_t site;
////the overlapID
uint32_t overlapID;
////the position of snp in that overlap
uint32_t overlapSite;
///there are several types: 0: equal to read 1: not equal to read, but it is a mismatch 2: is a gap
uint8_t type;
///misbase
char misBase;
}haplotype_evdience;
typedef struct
{
///the id of this snp
uint32_t id;
uint32_t overlap_num;
uint32_t occ_0;
uint32_t occ_1;
uint32_t occ_2;
int score;
////the position of snp in read itself
uint32_t site;
}
SnpStats;
#define Get_SNP_Martix_Size(matrix) (matrix.snp * matrix.overlap)
#define Get_SNP_Vector(matrix, i) (matrix.snp_matrix + matrix.overlap * i)
#define Get_SNP_Vector_Length(matrix) (matrix.overlap)
#define Get_Result_SNP_Vector(matrix) (matrix.snp_matrix + matrix.overlap*matrix.snp)
typedef struct
{
haplotype_evdience* list;
uint32_t sub_list_start;
uint32_t sub_list_length;
uint32_t length;
uint32_t size;
uint32_t flag[WINDOW];
uint32_t available_snp;
uint32_t core_snp;
uint32_t snp;
uint32_t overlap;
int8_t* snp_matrix;
uint32_t snp_matrix_size;
SnpStats* snp_stat;
SnpStats result_stat;
uint32_t snp_stat_size;
}
haplotype_evdience_alloc;
inline int filter_snp(int x, int y, int total)
{
double available;
if(x <= y)
{
available = x;
}
else
{
available = y;
}
double threshold = 0.30;
available = available/((double)(total));
if(available <= threshold && available < 6)
{
return 0;
}
return 1;
}
inline int filter_one_snp(int occ_0, int occ_1, int total)
{
double available;
if(occ_0 <= occ_1)
{
available = occ_0;
}
else
{
available = occ_1;
}
double threshold = 0.35;
available = available/((double)(total));
if(available < threshold || occ_0 < MIN_COVERAGE_THRESHOLD + 1 || total < 10)
{
return 0;
}
return 1;
}
inline void InsertSNPVector(haplotype_evdience_alloc* h, haplotype_evdience* sub_list, long long sub_length, char misBase)
{
if(sub_length <= 0)
return;
long long i = 0;
h->snp_stat[h->available_snp].id = h->available_snp;
h->snp_stat[h->available_snp].occ_0 = 0;
h->snp_stat[h->available_snp].occ_1 = 0;
h->snp_stat[h->available_snp].occ_2 = 0;
h->snp_stat[h->available_snp].overlap_num = 0;
h->snp_stat[h->available_snp].site = sub_list[0].site;
int8_t* vector = Get_SNP_Vector((*h), h->available_snp);
for (i = 0; i < sub_length; i++)
{
if(sub_list[i].type == 0)
{
vector[sub_list[i].overlapID] = 0;
h->snp_stat[h->available_snp].occ_0++;
}
else if(sub_list[i].type == 1 && sub_list[i].misBase == misBase)
{
vector[sub_list[i].overlapID] = 1;
h->snp_stat[h->available_snp].occ_1++;
}
else
{
vector[sub_list[i].overlapID] = 2;
h->snp_stat[h->available_snp].occ_2++;
}
h->snp_stat[h->available_snp].overlap_num++;
}
int new_occ_0 = h->snp_stat[h->available_snp].occ_0 + 1;
int new_occ_1 = h->snp_stat[h->available_snp].occ_1;
if(filter_snp(new_occ_0, new_occ_1, new_occ_0 + new_occ_1) == 0)
{
h->snp_stat[h->available_snp].score = -1;
}
else
{
h->core_snp++;
double consensus = new_occ_0 + new_occ_1 - abs(new_occ_0 - new_occ_1);
consensus = consensus /((double)(new_occ_0 + new_occ_1));
///50% vs 50%
if(new_occ_0 == new_occ_1)
{
consensus = consensus + 0.25;
}
else if(consensus >= 0.8)
{
consensus = consensus + 0.2;
}
else if(consensus >= 0.6)
{
consensus = consensus + 0.15;
}
else if(consensus >= 0.4)
{
consensus = consensus + 0.1;
}
else if(consensus >= 0.2)
{
consensus = consensus + 0.05;
}
consensus= consensus*((double)(new_occ_0 + new_occ_1));
h->snp_stat[h->available_snp].score = consensus;
}
h->available_snp++;
}
inline void SetSnpMatrix(haplotype_evdience_alloc* h, long long snp_num, long long overlap_num)
{
long long new_size = (snp_num + 1)* overlap_num;
if(h->snp_matrix_size < new_size)
{
h->snp_matrix_size = new_size;
h->snp_matrix = (int8_t*)realloc(h->snp_matrix, h->snp_matrix_size);
}
if(h->snp_stat_size < snp_num)
{
h->snp_stat_size = snp_num;
h->snp_stat = (SnpStats*)realloc(h->snp_stat, h->snp_stat_size * sizeof(SnpStats));
}
///h->snp may be different with the number of snp vector
///since some snps have been filtered
h->snp = snp_num;
h->overlap = overlap_num;
h->available_snp = 0;
h->core_snp = 0;
memset(h->snp_matrix, -1, h->snp * h->overlap);
}
inline void InitHaplotypeEvdience(haplotype_evdience_alloc* h)
{
h->snp = 0;
h->available_snp = 0;
h->overlap = 0;
h->snp_matrix_size = 0;
h->snp_stat_size = 0;
h->snp_matrix = NULL;
h->snp_stat = NULL;
h->sub_list_start = 0;
h->sub_list_length = 0;
h->length = 0;
h->size = 100;
h->list = (haplotype_evdience*)calloc(h->size, sizeof(haplotype_evdience));
memset(h->flag, 0, WINDOW * sizeof(uint32_t));
}
inline void StarSubListHaplotypeEvdience(haplotype_evdience_alloc* h)
{
h->sub_list_start = h->length;
}
inline void EndSubListHaplotypeEvdience(haplotype_evdience_alloc* h)
{
h->sub_list_length = h->length - h->sub_list_start;
}
inline void destoryHaplotypeEvdience(haplotype_evdience_alloc* h)
{
free(h->list);
free(h->snp_stat);
free(h->snp_matrix);
}
inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h)
{
h->snp = 0;
h->length = 0;
h->sub_list_start = 0;
h->sub_list_length = 0;
memset(h->flag, 0, WINDOW * sizeof(uint32_t));
}
inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h)
{
memset(h->flag, 0, WINDOW * sizeof(uint32_t));
}
inline void addHaplotypeEvdience(haplotype_evdience_alloc* h, haplotype_evdience* ev)
{
uint32_t new_length = h->length + 1;
if(new_length > h->size)
{
h->size = h->size * 2;
if(h->size < new_length)
{
h->size = new_length;
}
h->list = (haplotype_evdience*)realloc(h->list, sizeof(haplotype_evdience)*h->size);
}
h->list[h->length] = (*ev);
h->length++;
}
typedef struct
{
char* corrected_read;
@@ -41,7 +355,8 @@ typedef struct
void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g,
long long* matched_overlap_0, long long* matched_overlap_1,
long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1);
long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
Cigar_record* current_cigar, haplotype_evdience_alloc* hap);
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);
@@ -51,6 +366,192 @@ void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, u
All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length);
void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy);
void init_Cigar_record(Cigar_record* dummy);
void destory_Cigar_record(Cigar_record* dummy);
void clear_Cigar_record(Cigar_record* dummy);
inline void add_new_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type)
{
uint32_t tmp;
tmp = len;
tmp = tmp << 2;
tmp = tmp | type;
dummy->length++;
if(dummy->length > dummy->size)
{
dummy->size = dummy->size * 2;
dummy->record = (uint32_t*)realloc(dummy->record, dummy->size*sizeof(uint32_t));
}
dummy->record[dummy->length - 1] = tmp;
}
inline void add_existing_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type)
{
uint32_t tmp;
tmp = dummy->record[dummy->length - 1] >> 2;
tmp = tmp + len;
tmp = tmp << 2;
tmp = tmp | type;
dummy->record[dummy->length - 1] = tmp;
}
inline void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq)
{
uint32_t tmp;
tmp = len;
tmp = tmp << 2;
tmp = tmp | type;
dummy->length++;
if(dummy->length > dummy->size)
{
dummy->size = dummy->size * 2;
dummy->record = (uint32_t*)realloc(dummy->record, dummy->size*sizeof(uint32_t));
}
dummy->record[dummy->length - 1] = tmp;
if (dummy->lost_base_length + len> dummy->lost_base_size)
{
dummy->lost_base_size = (dummy->lost_base_length + len) * 2;
dummy->lost_base = (char*)realloc(dummy->lost_base, dummy->lost_base_size*sizeof(char));
}
int i = 0;
for (i = 0; i < len; i++, dummy->lost_base_length++)
{
dummy->lost_base[dummy->lost_base_length] = seq[i];
}
}
inline void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq)
{
uint32_t tmp;
tmp = dummy->record[dummy->length - 1] >> 2;
tmp = tmp + len;
tmp = tmp << 2;
tmp = tmp | type;
dummy->record[dummy->length - 1] = tmp;
if (dummy->lost_base_length + len> dummy->lost_base_size)
{
dummy->lost_base_size = (dummy->lost_base_length + len) * 2;
dummy->lost_base = (char*)realloc(dummy->lost_base, dummy->lost_base_size*sizeof(char));
}
int i = 0;
for (i = 0; i < len; i++, dummy->lost_base_length++)
{
dummy->lost_base[dummy->lost_base_length] = seq[i];
}
}
/***
type:
0. match
1. mismatch
2. insertion
3. deletion
***/
inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint32_t type)
{
uint32_t tmp;
if(type == 0)///match
{
///add to existing cell, just increase length
if(dummy->current_operation == type)
{
add_existing_cell_to_cigar_record(dummy, len, type);
}
else ///add to new cell
{
add_new_cell_to_cigar_record(dummy, len, type);
}
dummy->new_read_length += len;
}
else if(type == 1)///mismatch
{
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
else
{
add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
dummy->new_read_length += len;
}
else if(type == 3)///deletion, the bases in previous read will be removed
{
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
else
{
add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
}
else if(type == 2)///insertion
{
/**
///add to existing cell, just increase length
if(dummy->current_operation == type)
{
add_existing_cell_to_cigar_record(dummy, len, type);
}
else ///add to new cell
{
add_new_cell_to_cigar_record(dummy, len, type);
}
**/
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
else
{
add_new_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
dummy->new_read_length += len;
}
dummy->current_operation = type;
}
/**********************for prefilter************************ */
void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list);
void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length,
Hash_Table.cpp
+4 -1
View File
@@ -658,7 +658,8 @@ 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)
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end)
{
long long length = region->x_pos_e - region->x_pos_s + 1;
@@ -678,6 +679,8 @@ void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end
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[region->w_list_length].extra_begin = extra_begin;
region->w_list[region->w_list_length].extra_end = extra_end;
region->w_list_length++;
}
Hash_Table.h
+6 -1
View File
@@ -82,10 +82,13 @@ typedef struct
typedef struct
{
///the begining and end of a window, instead of the whole overlap
uint64_t x_start;
uint64_t x_end;
int y_end;
int y_start;
int extra_begin;
int extra_end;
///int y_pre_start;
///error小于等于0都要重新算
int error;
@@ -95,6 +98,7 @@ typedef struct
typedef struct
{
uint64_t x_id;
///the begining and end of the whole overlap
uint64_t x_pos_s;
uint64_t x_pos_e;
uint64_t x_pos_strand;
@@ -450,7 +454,8 @@ 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);
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end);
POA.cpp
-2
View File
@@ -256,8 +256,6 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
int last_operation = -1;
///note that node 0 is the start node
///0 is match, 1 is mismatch, 2 is up, 3 is left
///2是x缺字符（y多字符），而3是y缺字符（x多字符）
Process_Read.cpp
+142
View File
@@ -186,6 +186,9 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
fread(r->read_length, sizeof(uint64_t), r->total_reads, fp);
r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
/**********should remove**********/
///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5));
///fread(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp);
@@ -205,6 +208,19 @@ int load_All_reads(All_reads* r, char* read_file_name)
fread(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
{
r->cigars[i].size = 0;
r->cigars[i].length = 0;
r->cigars[i].record = NULL;
r->cigars[i].lost_base_size = 0;
r->cigars[i].lost_base_length = 0;
r->cigars[i].lost_base = NULL;
}
free(index_name);
fclose(fp);
fprintf(stdout, "Reads has been loaded.\n");
@@ -246,6 +262,9 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name)
void malloc_All_reads(All_reads* r)
{
r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
///必须加r->total_reads
/**********should remove**********/
///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5));
@@ -257,7 +276,17 @@ void malloc_All_reads(All_reads* r)
r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
}
r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
{
r->cigars[i].size = 0;
r->cigars[i].length = 0;
r->cigars[i].record = NULL;
r->cigars[i].lost_base_size = 0;
r->cigars[i].lost_base_length = 0;
r->cigars[i].lost_base = NULL;
}
r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -298,6 +327,119 @@ void init_UC_Read(UC_Read* r)
}
void recover_UC_Read_sub_region_begin_end
(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, int extra_begin, int extra_end)
{
long long readLen = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
long long i;
long long copyLen;
long long end_pos = start_pos + length - 1;
if (strand == 0)
{
i = start_pos;
copyLen = 0;
long long initLen = start_pos % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= start_pos && R_INF->N_site[ID][i] <= end_pos)
{
r[R_INF->N_site[ID][i] - start_pos] = 'N';
}
else if(R_INF->N_site[ID][i] > end_pos)
{
break;
}
}
}
}
else
{
start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1;
///start_pos > end_pos
i = start_pos;
copyLen = 0;
long long initLen = (start_pos + 1) % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
copyLen = copyLen + initLen;
i = i - initLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
copyLen = copyLen + 4;
i = i - 4;
}
if (R_INF->N_site[ID])
{
long long offset = readLen - start_pos - 1;
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= end_pos && R_INF->N_site[ID][i] <= start_pos)
{
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
}
else if(R_INF->N_site[ID][i] > start_pos)
{
break;
}
}
}
}
}
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
{
Process_Read.h
+18 -2
View File
@@ -48,8 +48,8 @@ static uint8_t seq_nt6_table[256] = {
static char bit_t_seq_table[256][4] = {0};
static char bit_t_seq_table_rc[256][4] = {0};
static char s_H[4] = {'A', 'C', 'G', 'T'};
static char rc_Table[4] = {'T', 'G', 'C', 'A'};
static char s_H[5] = {'A', 'C', 'G', 'T', 'N'};
static char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'};
#define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]]
@@ -58,6 +58,19 @@ void init_kseq(char* file);
void destory_kseq();
int get_read(kseq_t *s);
typedef struct
{
/**[0-1] bits are type:**/
/**[2-31] bits are length**/
uint32_t* record;
uint32_t length;
uint32_t size;
char* lost_base;
uint32_t lost_base_length;
uint32_t lost_base_size;
uint32_t new_length;
}Compressed_Cigar_record;
typedef struct
@@ -69,6 +82,7 @@ typedef struct
uint8_t** read_sperate;
uint64_t* read_length;
uint64_t* read_size;
///seq start pos in uint8_t* read
///do not need it
@@ -83,6 +97,8 @@ typedef struct
uint64_t total_reads_bases;
uint64_t total_name_length;
Compressed_Cigar_record* cigars;
} All_reads;
extern All_reads R_INF;
debug_assembly.sh
+17
View File
@@ -0,0 +1,17 @@
#!/bin/bash
# My first script
if [ $# -eq 1 ]
then
echo "../minimap2/minimap2 -ax asm20 -t 32 ../minimap2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz "$1" >"$1".sam"
../minimap2/minimap2 -ax asm20 -t 32 ../minimap2/Homo_sapiens.GRCh38.dna.primary_assembly.fa $1 >$1.sam
echo "samtools view -Sb "$1".sam >"$1".bam"
samtools view -Sb $1.sam >$1.bam
echo "samtools sort "$1".bam sort_"$1
samtools sort $1.bam sort_$1
echo "rm sort_"$1".bam.bai"
rm sort_$1.bam.bai
echo "samtools index sort_"$1".bam"
samtools index sort_$1.bam
else
echo "debug_assembly.sh intput.fa"
fi
main.cpp
+2 -19
View File
@@ -206,6 +206,7 @@ int main(int argc, char *argv[])
return 1;
**/
fprintf(stdout, "Will perform %d round of error correction...\n", number_of_round);
fprintf(stdout, "defined k_mer_min_freq by user: %d\n", k_mer_min_freq);
fprintf(stdout, "defined k_mer_max_freq by user: %d\n", k_mer_max_freq);
@@ -213,25 +214,7 @@ int main(int argc, char *argv[])
fprintf(stdout, "k-mer length: %d\n",k_mer_length);
if (load_index_from_disk && load_pre_cauculated_index())
{
;
}
else
{
Counting_multiple_thr();
Build_hash_table_multiple_thr();
}
fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq);
fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq);
///verify_Position_hash_table();
Overlap_calculate_multipe_thr();
Correct_Reads(number_of_round);
sort_fasta.sh
+10
View File
@@ -0,0 +1,10 @@
#!/bin/bash
# My first script
if [ $# -eq 1 ]
then
echo "input file is: "$1
echo $1" | paste - - | sort -k1,1 -t \" \" | tr \"\t\" \"\n\" > sort_"$1
cat $1 | paste - - | sort -k1,1 -t " " | tr "\t" "\n" > sort_$1
else
echo "./sort_fasta.sh file_name"
fi