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

backup

Browse Source
This commit is contained in:
Haoyu Cheng
2019-11-19 01:32:39 -05:00
parent ada03ea388
commit e8ff0c1dbc
2 changed files with 465 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files

473
Correct.cpp
View File

@@ -8351,38 +8351,73 @@ void Preorder_Merge_Advance(uint32_t snpID, haplotype_evdience_alloc* hap, int p
int if_snp_vector_useful(haplotype_evdience_alloc* hap,
long long occ_0, long long occ_1, long long occ_1_low,
long long coverage, uint32_t* SNPs, long long SNPsLen, int roundID)
int if_snp_vector_useful_v2(haplotype_evdience_alloc* hap,
long long occ_0, long long occ_1, uint32_t* SNPs, long long SNPsLen)
{
double occ_1_coverage_low = coverage * 0.3;
if(occ_1 >= occ_1_low && occ_0 >= occ_1_low)
double occ_1_coverage_low = (occ_0 + occ_1) * 0.3;
if(occ_1 == 0 || occ_0 == 0)
{
if(occ_1 >= occ_1_coverage_low && occ_0 >= occ_1_coverage_low)
return 0;
}
if(occ_1 >= occ_1_coverage_low && occ_0 >= occ_1_coverage_low)
{
return 1;
}
else if(occ_1 >= 5 && occ_0 >= 5)
{
return 1;
}
else if(occ_1 >= 2 && occ_0 >= 2 && SNPsLen >= 2)
{
/**
int nearsnp;
int non_nearsnps;
count_nearby_snps(hap, SNPs, SNPsLen, &nearsnp, &non_nearsnps);
if(non_nearsnps > 0)
{
return 1;
}
else if(occ_1 >= 5 && occ_0 >= 5)
**/
return 1;
}
return 0;
}
int if_snp_vector_useful(haplotype_evdience_alloc* hap,
long long occ_0, long long occ_1, uint32_t* SNPs, long long SNPsLen)
{
double occ_1_coverage_low = (occ_0 + occ_1) * 0.3;
if(occ_1 == 0 || occ_0 == 0)
{
return 0;
}
if(occ_1 >= occ_1_coverage_low && occ_0 >= occ_1_coverage_low)
{
return 1;
}
else if(occ_1 >= 5 && occ_0 >= 5)
{
return 1;
}
else if(occ_1 >= 3 && occ_0 >= 3 && SNPsLen >= 2)
{
int nearsnp;
int non_nearsnps;
count_nearby_snps(hap, SNPs, SNPsLen, &nearsnp, &non_nearsnps);
if(non_nearsnps > 0)
{
return 1;
}
else if(occ_1 >= 2 && occ_0 >= 2 && SNPsLen >= 2)
{
int nearsnp;
int non_nearsnps;
count_nearby_snps(hap, SNPs, SNPsLen, &nearsnp, &non_nearsnps);
if(non_nearsnps > 0)
{
return 1;
}
}
else if(roundID == 1 && occ_0 >= occ_1_coverage_low)
{
return 1;
}
}
return 0;
@@ -8604,9 +8639,6 @@ void process_repeat_snps(haplotype_evdience_alloc* hap, int coverage, overlap_re
{
int i, snpID, vectorID, flag;
int8_t *vector;
long long occ_1_threshold_low;
long long occ_1_threshold_up;
occ_1_threshold_low = 0;
uint32_t* snp_ids;
@@ -8621,15 +8653,9 @@ void process_repeat_snps(haplotype_evdience_alloc* hap, int coverage, overlap_re
merge_SNP_Vectors(hap, snp_ids, length);
/**
if(if_snp_vector_useful(hap, hap->dp.SNP_IDs.IDs[i].occ_0, hap->dp.SNP_IDs.IDs[i].occ_1,
occ_1_threshold_low, coverage, snp_ids, length, 0))**/
if(if_snp_vector_useful(hap, hap->result_stat.occ_0, hap->result_stat.occ_1,
occ_1_threshold_low, coverage, snp_ids, length, 0))
if(if_snp_vector_useful(hap, hap->result_stat.occ_0, hap->result_stat.occ_1,
snp_ids, length))
{
//fprintf(stderr, "i: %d \n", i);
///remove_reads(hap, snp_ids, length, overlap_list);
try_to_remove_reads(Get_Result_SNP_Vector((*hap)), Get_SNP_Vector_Length((*hap)),
overlap_list, snp_ids, length, hap);
@@ -8639,7 +8665,6 @@ void process_repeat_snps(haplotype_evdience_alloc* hap, int coverage, overlap_re
{
hap->dp.SNP_IDs.IDs[i].is_remove = 0;
}
}
@@ -8694,8 +8719,8 @@ overlap_region_alloc* overlap_list, All_reads* R_INF)
merge_SNP_Vectors(hap, snp_ids, length);
if(if_snp_vector_useful(hap, hap->result_stat.occ_0, hap->result_stat.occ_1,
occ_1_threshold_low, coverage, snp_ids, length, 0))
if(if_snp_vector_useful(hap, hap->result_stat.occ_0, hap->result_stat.occ_1,
snp_ids, length))
{
@@ -8856,7 +8881,7 @@ void debug_repeat_vector(haplotype_evdience_alloc* hap)
}
int generate_haplotypes_DP(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF, long long rLen,
int generate_haplotypes_DP_back(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF, long long rLen,
int force_repeat)
{
int j, i;
@@ -9117,7 +9142,7 @@ int force_repeat)
}
/**
if(memcmp("m64011_190329_072846/59507330/ccs",
if(memcmp("m64016_190918_162737/49678749/ccs",
Get_NAME((*R_INF), overlap_list->list[0].x_id),
Get_NAME_LENGTH((*R_INF), overlap_list->list[0].x_id)) == 0)
{
@@ -9138,6 +9163,7 @@ int force_repeat)
if(overlap_list->mapped_overlaps_length > Coverage_Threshold(coverage, rLen) || force_repeat)
{
@@ -9407,6 +9433,258 @@ int force_repeat)
}
int generate_haplotypes_DP(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF, long long rLen,
int force_repeat)
{
int j, i;
int vectorID, vectorID2;
int diff_core_vector = 0;
int diff_vector_ID = -1;
int8_t *vector, *vector2;
if(hap->available_snp == 0)
{
return 0;
}
///if hap->available_snp == 1, the following codes would have bugs
///filter snps that are highly likly false
if(hap->available_snp > 1)
{
i = 0;
///if a snp is very near to others, it should not be a real snp
for (j = 0; j < hap->available_snp; j++)
{
if(j > 0 && j < hap->available_snp - 1)
{
if(hap->snp_stat[j].site != hap->snp_stat[j - 1].site + 1
&&
hap->snp_stat[j].site + 1 != hap->snp_stat[j + 1].site)
{
hap->snp_stat[i] = hap->snp_stat[j];
i++;
}
}
else if(j == 0)
{
if(hap->snp_stat[j].site + 1 != hap->snp_stat[j + 1].site)
{
hap->snp_stat[i] = hap->snp_stat[j];
i++;
}
}
else
{
if(hap->snp_stat[j].site != hap->snp_stat[j - 1].site + 1)
{
hap->snp_stat[i] = hap->snp_stat[j];
i++;
}
}
}
hap->available_snp = i;
}
int flag;
long long overlap_length, total_read, unuseful_read, last_j, last_j_ID, last_j_flag;
total_read = unuseful_read = 0;
///check if any read may be conflict with others
for (i = 0; i < overlap_list->length; i++)
{
overlap_length = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1;
if (overlap_list->list[i].is_match == 1)
{
total_read++;
flag = -1;
for (j = 0; j < hap->available_snp; j++)
{
vectorID = hap->snp_stat[j].id;
vector = Get_SNP_Vector((*hap), vectorID);
///flag == -1 means there are no useful signals yet
if (flag == -1)
{
if((vector[i] == 0 || vector[i] == 1 ))
{
flag = 0;
}
}///flag == 0 means there is at least one useful signal yet
else if (flag == 0)
{
if(vector[i] != 0 && vector[i] != 1)
{
flag = 2;
last_j = hap->snp_stat[j].site;
last_j_ID = j;
last_j_flag = vector[i];
}
}///flag == 0 means there is at least one useful signal first, and another unuseful signal after that
else if(flag == 2)
{
if((vector[i] == 0 || vector[i] == 1 ))
{
flag = 3;
break;
}
}
}
if(flag == 3)
{
unuseful_read++;
for (j = 0; j < hap->available_snp; j++)
{
vectorID = hap->snp_stat[j].id;
vector = Get_SNP_Vector((*hap), vectorID);
if(vector[i] == 0)
{
hap->snp_stat[j].occ_0--;
hap->snp_stat[j].occ_2++;
}
else if(vector[i] == 1)
{
hap->snp_stat[j].occ_1--;
hap->snp_stat[j].occ_2++;
}
else if(vector[i] != 2)
{
hap->snp_stat[j].occ_2++;
}
vector[i] = 2;
}
///this read may be unuseful
///overlap_list->list[i].is_match = 0;
///overlap_list->list[i].is_match = 2;
overlap_list->list[i].is_match = 4;
///overlap_list->mapped_overlaps--;
overlap_list->mapped_overlaps_length -= overlap_length;
}
}
}
/*******************************DP********************************/
init_DP_matrix(&(hap->dp), hap->available_snp);
long long equal_best = 0;
uint32_t* column;
long long column_length;
for (i = 0; i < hap->available_snp; i++)
{
///vector of snp i
vectorID = hap->snp_stat[i].id;
vector = Get_SNP_Vector((*hap), vectorID);
hap->dp.visit[i] = 0;
hap->dp.max[i] = 1;
hap->dp.backtrack_length[i] = 0;
equal_best = 0;
column = Get_DP_Backtrack_Column(hap->dp, i);
column_length = Get_DP_Backtrack_Column_Length(hap->dp, i);
for (j = 0; j < i; j++)
{
///vector of snp j
vectorID2 = hap->snp_stat[j].id;
vector2 = Get_SNP_Vector((*hap), vectorID2);
///vector is compatible with vector2
if(calculate_distance_snp_vector(vector, vector2, Get_SNP_Vector_Length((*hap))) == 0)
{
if(hap->dp.max[i] < hap->dp.max[j] + 1)
{
hap->dp.max[i] = hap->dp.max[j] + 1;
column[0] = j;
equal_best = 1;
}
else if(hap->dp.max[i] == hap->dp.max[j] + 1)
{
column[equal_best] = j;
equal_best++;
}
}
}
hap->dp.backtrack_length[i] = equal_best;
}
/*******************************DP********************************/
uint64_t tmp_mode = 0;
for (i = 0; i < hap->available_snp; i++)
{
tmp_mode = hap->dp.max[i];
tmp_mode = tmp_mode << 32;
tmp_mode = tmp_mode | (uint64_t)(i);
hap->dp.max_for_sort[i] = tmp_mode;
}
qsort(hap->dp.max_for_sort, hap->available_snp, sizeof(uint64_t), cmp_max_DP);
int snpID;
int group_num = 0;
///the minmum snp_num is 1
hap->dp.max_snp_num = 0;
hap->dp.max_score = -2;
for (i = 0; i < hap->available_snp; i++)
{
snpID = Get_Max_DP_ID(hap->dp.max_for_sort[i]);
if(hap->dp.visit[snpID] == 0)
{
hap->dp.current_snp_num = Get_Max_DP_Value(hap->dp.max_for_sort[i]);
Preorder_Merge_Advance_Repeat(snpID, hap, 0);
}
}
//if(hap->dp.max_snp_num > 0)
if(hap->available_snp > 0)
{
process_repeat_snps(hap, coverage, overlap_list);
return 1;
}
else
{
return 0;
}
}
inline int check_informative_site(haplotype_evdience_alloc* hap, SnpStats* snp)
{
long long vectorID = snp->id;
@@ -9536,38 +9814,6 @@ int force_repeat)
long long m, snp_occ;
/**
while (1)
{
m = 0;
for (i = 0; i < hap->available_snp; i++)
{
if(check_informative_site(hap, &(hap->snp_stat[i])))
{
hap->snp_stat[m] = hap->snp_stat[i];
m++;
}
}
if(m == hap->available_snp)
{
break;
}
hap->available_snp = m;
for (i = 0; i < overlap_list->length; i++)
{
if (overlap_list->list[i].is_match == 1)
{
snp_occ = snp_occ_in_one_read(hap, overlap_list, i);
if(snp_occ >= 1)
{
remove_read_from_snps(hap, overlap_list, i);
}
}
}
}
**/
if(hap->available_snp > 0)
{
///************************debug**************************///
@@ -10036,8 +10282,8 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF,
}
///debug_snp_matrix(hap);
///generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat);
generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat);
generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat);
///generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat);
///debug_snp_matrix(hap);
@@ -10145,7 +10391,7 @@ void correct_overlap_back(overlap_region_alloc* overlap_list, All_reads* R_INF,
void print_overlap(char* name, long long readID,
overlap_region_alloc* overlap_list, All_reads* R_INF)
overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads)
{
if(memcmp(name, Get_NAME((*R_INF), readID),
Get_NAME_LENGTH((*R_INF),readID)) == 0)
@@ -10212,7 +10458,80 @@ overlap_region_alloc* overlap_list, All_reads* R_INF)
overlap_list->list[i].strong);
}
}
if(output_reads)
{
UC_Read g_read;
init_UC_Read(&g_read);
recover_UC_Read(&g_read, R_INF, readID);
fprintf(stderr, "\n\nOutput all related reads\n");
fprintf(stderr, "ref_read:\n");
fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((*R_INF),readID), Get_NAME((*R_INF),readID));
fprintf(stderr, "%.*s\n", g_read.length, g_read.seq);
fprintf(stderr, "query_read:\n");
for (i = 0; i < overlap_list->length; i++)
{
fprintf(stderr, "i: %d\n", i);
recover_UC_Read(&g_read, R_INF, overlap_list->list[i].y_id);
fprintf(stderr, ">%.*s\n",
Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id),
Get_NAME((*R_INF),overlap_list->list[i].y_id));
fprintf(stderr, "%.*s\n", g_read.length, g_read.seq);
}
destory_UC_Read(&g_read);
for (i = 0; i < overlap_list->length; i++)
{
fprintf(stderr, "\ni: %d, %.*s, x_s: %d, x_e: %d, y_s: %d, y_end: %d, w_list_length: %d, dir: %d, strong: %d, is_match: %d\n",
i, Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id),
Get_NAME((*R_INF),overlap_list->list[i].y_id),
overlap_list->list[i].x_pos_s,
overlap_list->list[i].x_pos_e,
overlap_list->list[i].y_pos_s,
overlap_list->list[i].y_pos_e,
overlap_list->list[i].w_list_length,
overlap_list->list[i].y_pos_strand,
overlap_list->list[i].strong,
overlap_list->list[i].is_match);
for (j = 0; j < overlap_list->list[i].w_list_length; j++)
{
fprintf(stderr, "************************\ncigar_j: %d, x_s: %d, x_e: %d, y_s: %d, y_end: %d\n",
j, overlap_list->list[i].w_list[j].x_start,
overlap_list->list[i].w_list[j].x_end,
overlap_list->list[i].w_list[j].y_start,
overlap_list->list[i].w_list[j].y_end);
if(overlap_list->list[i].w_list[j].y_end == -1)
{
fprintf(stderr, "not match\n");
}
else
{
int cigar_i, operation, operationLen;
CIGAR* cigar = &(overlap_list->list[i].w_list[j].cigar);
fprintf(stderr, "length: %d\n", cigar->length);
for (cigar_i = 0; cigar_i < cigar->length; cigar_i++)
{
operation = cigar->C_C[cigar_i];
operationLen = cigar->C_L[cigar_i];
fprintf(stderr, "oper: %d, Len: %d\n", operation, operationLen);
}
}
}
}
}
}
}
@@ -10275,8 +10594,8 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat);
// print_overlap("m64011_190329_072846/59507330/ccs",
// overlap_list->list[0].x_id, overlap_list, R_INF);
print_overlap("m64016_190918_162737/53545052/ccs",
overlap_list->list[0].x_id, overlap_list, R_INF, 1);

93
Overlaps.cpp
View File

@@ -403,6 +403,8 @@ void normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources)
}
else
{
///must have this line
new_element.ml = 1;
set_reverse_overlap(&new_element, &(sources[i].buffer[j]));
add_ma_hit_t_alloc(&(sources[tn]), &new_element);
si_overlaps++;
@@ -5054,6 +5056,27 @@ long long weakID, uint32_t w_qs, uint32_t w_qe)
}
inline int check_weak_ma_hit_reverse(ma_hit_t_alloc* aim_paf, ma_hit_t_alloc* reverse_paf_list,
long long weakID)
{
long long i = 0;
long long strongID, index;
///all overlaps coming from another haplotye are strong
for (i = 0; i < aim_paf->length; i++)
{
strongID = Get_tn(aim_paf->buffer[i]);
index = get_specific_overlap
(&(reverse_paf_list[strongID]), strongID, weakID);
if(index != -1)
{
return 0;
}
}
return 1;
}
inline int check_weak_ma_hit_debug(ma_hit_t_alloc* aim_paf, ma_hit_t_alloc* reverse_paf_list,
long long weakID)
{
@@ -6897,6 +6920,17 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
uint32_t qn, tn;
ma_hit_t new_element;
long long qLen_0, qLen_1;
// if(memcmp("m64016_190918_162737/92668450/ccs", Get_NAME(R_INF, i),
// Get_NAME_LENGTH(R_INF, i)) == 0)
// {
// debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs",
// sources, reverse_sources, -1, "clean_weak_ma_hit_t");
// debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs",
// sources, reverse_sources, -1, "clean_weak_ma_hit_t");
// }
for (i = 0; i < num_sources; i++)
{
for (j = 0; j < sources[i].length; j++)
@@ -6907,20 +6941,13 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
//if this is a weak overlap
if(sources[i].buffer[j].ml == 0)
{
if(!check_weak_ma_hit(&(sources[qn]), reverse_sources, tn,
Get_qs(sources[i].buffer[j]), Get_qe(sources[i].buffer[j])))
if(
!check_weak_ma_hit(&(sources[qn]), reverse_sources, tn,
Get_qs(sources[i].buffer[j]), Get_qe(sources[i].buffer[j]))
/**
||
!check_weak_ma_hit_reverse(&(reverse_sources[qn]), sources, tn)**/)
{
/**
if(memcmp("m64016_190918_162737/76808505/ccs",
Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i)) == 0)
{
fprintf(stderr, "#### %.*s, %.*s\n",
Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn),
Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn));
}
**/
sources[i].buffer[j].bl = 0;
index = get_specific_overlap(&(sources[tn]), tn, qn);
sources[tn].buffer[index].bl = 0;
@@ -6929,6 +6956,8 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
}
}
long long m = 0;
long long pre_overlaps, current_overlaps, exact_overlaps;
exact_overlaps = pre_overlaps = current_overlaps = 0;
@@ -6954,10 +6983,17 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
current_overlaps += sources[i].length;
}
/**
fprintf(stdout, "pre_overlaps: %lld, current_overlaps: %lld, exact_overlaps: %lld\n",
pre_overlaps, current_overlaps, exact_overlaps);
**/
// if(memcmp("m64016_190918_162737/92668450/ccs", Get_NAME(R_INF, i),
// Get_NAME_LENGTH(R_INF, i)) == 0)
// {
// debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs",
// sources, reverse_sources, -1, "clean_weak_ma_hit_t");
// debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs",
// sources, reverse_sources, -1, "clean_weak_ma_hit_t");
// }
fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
}
@@ -7093,14 +7129,15 @@ ma_hit_t_alloc* reverse_sources, int id, char* command)
{
qn = Get_qn(sources[i].buffer[j]);
tn = Get_tn(sources[i].buffer[j]);
fprintf(stderr, "target: %.*s, qs: %d, qe: %d, ts: %d, te: %d, ml: %d, rev: %d\n",
fprintf(stderr, "target: %.*s, qs: %d, qe: %d, ts: %d, te: %d, ml: %d, rev: %d, el: %d\n",
Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn),
Get_qs(sources[i].buffer[j]),
Get_qe(sources[i].buffer[j]),
Get_ts(sources[i].buffer[j]),
Get_te(sources[i].buffer[j]),
sources[i].buffer[j].ml,
sources[i].buffer[j].rev);
sources[i].buffer[j].rev,
sources[i].buffer[j].el);
@@ -7419,16 +7456,24 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
}
// debug_info_of_specfic_read("m64016_190918_162737/72220752/ccs",
// debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs",
// sources, reverse_sources, -1, "init");
debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs",
sources, reverse_sources, -1, "init");
ma_sub_t* coverage_cut;
normalize_ma_hit_t(sources, n_read);
///normalize_ma_hit_t_single_side(sources, n_read);
///normalize_ma_hit_t(sources, n_read);
normalize_ma_hit_t_single_side(sources, n_read);
// debug_info_of_specfic_read("m64016_190918_162737/72220752/ccs",
// debug_info_of_specfic_read("m64016_190918_162737/92668450/ccs",
// sources, reverse_sources, -1, "normalize");
debug_info_of_specfic_read("m64016_190918_162737/53545052/ccs",
sources, reverse_sources, -1, "normalize");
@@ -7444,7 +7489,7 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
// debug_info_of_specfic_read("m64016_190918_162737/72220752/ccs",
// debug_info_of_specfic_read("m64016_190918_162737/49678749/ccs",
// sources, reverse_sources, -1, "clean");