xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-15 12:47:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

r183: reorganized the major workflow

Browse Source
This commit is contained in:
Heng Li
2020-04-03 13:34:07 -04:00
parent ffba3ce0ef
commit d529dcea3f
5 changed files with 75 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

206
Assembly.cpp
View File

@@ -1,5 +1,6 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <assert.h>
#include <zlib.h> #include <zlib.h>
#include "Assembly.h" #include "Assembly.h"
#include "Process_Read.h" #include "Process_Read.h"
@@ -482,7 +483,6 @@ void* Overlap_calculate_heap_merge(void* arg)
fprintf(stderr, "total recorrected bases: %lld\n", asm_opt.num_recorrected_bases); fprintf(stderr, "total recorrected bases: %lld\n", asm_opt.num_recorrected_bases);
} }
pthread_mutex_unlock(&statistics); pthread_mutex_unlock(&statistics);
free(arg);
return NULL; return NULL;
} }
@@ -578,9 +578,6 @@ void* Output_related_reads(void* arg)
destoryHaplotypeEvdience(&hap); destoryHaplotypeEvdience(&hap);
destory_Round2_alignment(&second_round); destory_Round2_alignment(&second_round);
free(arg);
return NULL; return NULL;
} }
@@ -694,7 +691,6 @@ void* Save_corrected_reads(void* arg)
destory_UC_Read(&g_read); destory_UC_Read(&g_read);
free(first_round_read); free(first_round_read);
free(second_round_read); free(second_round_read);
free(arg);
return NULL; return NULL;
} }
@@ -722,84 +718,37 @@ void Output_corrected_reads()
fclose(output_file); fclose(output_file);
} }
void ha_overlap_and_correct(int round)
void Overlap_calculate_multipe_thr()
{ {
double start_time = Get_T(); int i, *args;
pthread_t *_r_threads;
MALLOC(_r_threads, asm_opt.thread_num);
args = (int*)alloca(sizeof(int) * asm_opt.thread_num);
fprintf(stderr, "Begin calculating overlaps... \n"); ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF); // build the index
for (i = 0; i < asm_opt.thread_num; i++) {
pthread_t *_r_threads; args[i] = i;
if (!asm_opt.required_read_name)
_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num); pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)&args[i]);
else
int i = 0; pthread_create(_r_threads + i, NULL, Output_related_reads, (void*)&args[i]);
for (i = 0; i < asm_opt.thread_num; i++)
{
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
if(!asm_opt.required_read_name)
{
pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)arg);
}
else
{
pthread_create(_r_threads + i, NULL, Output_related_reads, (void*)arg);
}
} }
for (i = 0; i < asm_opt.thread_num; i++)
for (i = 0; i < asm_opt.thread_num; i++)
pthread_join(_r_threads[i], NULL); pthread_join(_r_threads[i], NULL);
free(_r_threads);
if(asm_opt.required_read_name)
{
exit(1);
}
ha_pt_destroy(ha_idx); ha_pt_destroy(ha_idx);
ha_idx = 0; ha_idx = 0;
fprintf(stderr, "All overlaps have been calculated.\n"); if (asm_opt.required_read_name) exit(0); // for debugging only
fprintf(stderr, "%-30s%18.2f\n\n", "Overlap calculation time:", Get_T() - start_time); for (i = 0; i < asm_opt.thread_num; i++) {
args[i] = i;
start_time = Get_T(); pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)&args[i]);
_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; i++)
{
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)arg);
} }
for (i = 0; i < asm_opt.thread_num; i++)
for (i = 0; i < asm_opt.thread_num; i++)
pthread_join(_r_threads[i], NULL); pthread_join(_r_threads[i], NULL);
free(_r_threads); free(_r_threads);
fprintf(stderr, "%-30s%18.2f\n\n", "Corrected read saving time:", Get_T() - start_time);
///only the last round can output read to disk
if (asm_opt.roundID == asm_opt.number_of_round - 1)
{
start_time = Get_T();
Output_corrected_reads();
fprintf(stderr, "%-30s%18.2f\n\n", "Output time:", Get_T() - start_time);
}
} }
void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf,
UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact) UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
{ {
@@ -1159,7 +1108,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
} }
} }
pthread_mutex_unlock(&statistics); pthread_mutex_unlock(&statistics);
free(arg);
return NULL; return NULL;
} }
@@ -1331,80 +1279,60 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises); fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises);
} }
void ha_overlap_final(void)
void generate_overlaps(int last_round)
{ {
double start_time = Get_T(); int i, *args;
asm_opt.roundID = asm_opt.number_of_round - last_round; pthread_t *_r_threads;
fprintf(stderr, "Begin calculting final overlaps ...\n");
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF); MALLOC(_r_threads, asm_opt.thread_num);
args = (int*)alloca(sizeof(int) * asm_opt.thread_num);
pthread_t *_r_threads; ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF);
for (i = 0; i < asm_opt.thread_num; i++) {
_r_threads = (pthread_t *)malloc(sizeof(pthread_t) * asm_opt.thread_num); args[i] = i;
pthread_create(_r_threads + i, NULL, Final_overlap_calculate_heap_merge, (void*)&args[i]);
int i = 0; }
for (i = 0; i < asm_opt.thread_num; i++)
for (i = 0; i < asm_opt.thread_num; i++) pthread_join(_r_threads[i], NULL);
{ free(_r_threads);
int *arg = (int*)malloc(sizeof(*arg));
*arg = i;
pthread_create(_r_threads + i, NULL, Final_overlap_calculate_heap_merge, (void*)arg);
}
for (i = 0; i < asm_opt.thread_num; i++)
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);
ha_pt_destroy(ha_idx); ha_pt_destroy(ha_idx);
ha_idx = 0; ha_idx = 0;
fprintf(stderr, "Final overlaps have been calculated.\n"); ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985);
fprintf(stderr, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> final overlap\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
Output_PAF();
trio_partition();
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name,
asm_opt.large_pop_bubble_size, 0, 1);
} }
int ha_assemble(void)
void Correct_Reads(int last_round)
{ {
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) int r, ovlp_loaded = 0;
{ if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, ovlp_loaded = 1;
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded overlaps from disk\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, 0); }
exit(0); if (!ovlp_loaded) {
} // construct hash table for high occurrence k-mers
else if (!asm_opt.no_kmer_flt)
{ ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF);
///fprintf(stderr, "Cannot find overlap file. Please run the whole hifiasm.\n"); // error correction
} assert(asm_opt.number_of_round > 0);
for (r = 0; r < asm_opt.number_of_round; ++r) {
clear_opt(&asm_opt, last_round); clear_opt(&asm_opt, asm_opt.number_of_round - r); // this update asm_opt.roundID and a few other fields
ha_overlap_and_correct(r);
if(last_round == 0) fprintf(stderr, "[M::%s::%.3f*%.2f] ==> corrected reads for round %d\n", __func__, yak_realtime(), yak_cputime() / yak_realtime(), r + 1);
{ }
generate_overlaps(last_round); Output_corrected_reads();
return; fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
} // overlap between corrected reads
clear_opt(&asm_opt, 0);
fprintf(stderr, "Error correction: Start the %d-th round ...\n", asm_opt.roundID); ha_overlap_final();
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cputime() / yak_realtime());
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF); ha_ft_destroy(ha_flt_tab);
Overlap_calculate_multipe_thr(); Output_PAF();
ha_pt_destroy(ha_idx); trio_partition();
}
fprintf(stderr, "Error correction: The %d-th round has been completed.\n", asm_opt.roundID); build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
Correct_Reads(last_round - 1); asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
destory_All_reads(&R_INF);
return 0;
} }

4
Assembly.h
View File

@@ -8,6 +8,6 @@
#define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Type(RECORD) (RECORD&3)
#define Get_Cigar_Length(RECORD) (RECORD>>2) #define Get_Cigar_Length(RECORD) (RECORD>>2)
void Overlap_calculate_multipe_thr(); int ha_assemble(void);
void Correct_Reads(int last_round);
#endif #endif

2
CommandLines.h
View File

@@ -3,7 +3,7 @@
#include <pthread.h> #include <pthread.h>
#define HA_VERSION "0.3.0-dirty-r182" #define HA_VERSION "0.3.0-dirty-r183"
#define VERBOSE 0 #define VERBOSE 0
#define VERBOSE_GFA 0 #define VERBOSE_GFA 0

2
Process_Read.cpp
View File

@@ -114,7 +114,7 @@ void write_All_reads(All_reads* r, char* read_file_name)
int load_All_reads(All_reads* r, char* read_file_name) int load_All_reads(All_reads* r, char* read_file_name)
{ {
fprintf(stderr, "Loading reads from disk... \n"); //fprintf(stderr, "Loading reads from disk... \n");
char* index_name = (char*)malloc(strlen(read_file_name)+15); char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name); sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "r"); FILE* fp = fopen(index_name, "r");

12
main.cpp
View File

@@ -8,20 +8,16 @@
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
int i; int i, ret;
yak_reset_realtime();
init_opt(&asm_opt); init_opt(&asm_opt);
if (!CommandLine_process(argc, argv, &asm_opt)) return 1; if (!CommandLine_process(argc, argv, &asm_opt)) return 1;
yak_reset_realtime(); ret = ha_assemble();
if (!asm_opt.no_kmer_flt)
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF);
Correct_Reads(asm_opt.number_of_round);
ha_ft_destroy(ha_flt_tab);
destory_All_reads(&R_INF);
destory_opt(&asm_opt); destory_opt(&asm_opt);
fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION); fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);
fprintf(stderr, "[M::%s] CMD:", __func__); fprintf(stderr, "[M::%s] CMD:", __func__);
for (i = 0; i < argc; ++i) for (i = 0; i < argc; ++i)
fprintf(stderr, " %s", argv[i]); fprintf(stderr, " %s", argv[i]);
fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss() / 1024.0 / 1024.0 / 1024.0); fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss() / 1024.0 / 1024.0 / 1024.0);
return 0; return ret;
} }