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
Files
6de4e147825cb2537cf600e5c744eb3192d10b9e
hifiasm/Correct.h
chhylp123 6de4e14782 update before fly HPC
2024-10-29 14:59:52 -04:00

1415 lines
41 KiB
C
Raw Blame History

#ifndef __CORRECT__
#define __CORRECT__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Hash_Table.h"
#include "Levenshtein_distance.h"
#include "POA.h"
#include "Process_Read.h"
#include "Correct.h"
#include "kalloc.h"
//#define CORRECT_THRESHOLD 0.70
#define CORRECT_THRESHOLD 0.60
///#define CORRECT_THRESHOLD_SECOND 0.55
#define CORRECT_THRESHOLD_HOMOPOLYMER 0.515
#define MIN_COVERAGE_THRESHOLD 3
#define CORRECT_INDEL_LENGTH 2
#define MISMATCH 1
#define INSERTION 2
#define DELETION 3
#define ERROR_RATE 1.5
#define UL_TOPN 50
#define SGAP 16
#define MAX_LGAP(ql) ((((ql)*0.2)<256)?((ql)*0.2):256)
#define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff?
///#define FLAG_THRE 0
#define MAX(x, y) (((x) >= (y))?(x):(y))
#define MIN(x, y) (((x) <= (y))?(x):(y))
#define DIFF(x, y) ((MAX((x), (y))) - (MIN((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)])
#define Coverage_Threshold(coverage, r_len) (coverage*r_len*1.1)
#define Get_Max_DP_Value(RECORD) (RECORD>>32)
#define Get_Max_DP_ID(RECORD) (RECORD&(uint64_t)0xffffffff)
#define Adjust_Threshold(threshold, x_len) ((threshold == 0 && x_len >= 4)? 1: threshold)
typedef struct
{
long long read_length;
long long window_length;
long long window_num;
long long window_start;
long long window_end;
long long tail_length;
int terminal;
}Window_Pool;
inline void init_Window_Pool(Window_Pool* dumy, long long read_length, long long window_length, long long tail_length)
{
dumy->terminal = 0;
dumy->read_length = read_length;
dumy->window_length = window_length;
dumy->tail_length = tail_length;
dumy->window_start = 0;
dumy->window_end = dumy->window_length - 1;
if (dumy->window_end >= dumy->read_length)
{
dumy->window_end = dumy->read_length - 1;
}
dumy->window_num = (dumy->read_length + dumy->window_length - 1) / dumy->window_length;
}
inline int get_Window(Window_Pool* dumy, long long* w_beg, long long* w_end)
{
(*w_beg) = dumy->window_start;
(*w_end) = dumy->window_end;
if(dumy->window_end == dumy->read_length - 1)
{
if(dumy->terminal == 1 || dumy->read_length == 0)
{
return 0;
}
else if(dumy->terminal == 0)
{
dumy->terminal = 1;
}
}
dumy->window_start = dumy->window_start + dumy->window_length;
dumy->window_end = dumy->window_end + dumy->window_length;
if (dumy->window_end >= dumy->read_length)
{
dumy->window_end = dumy->read_length - 1;
}
return 1;
}
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
{
long long length;
long long size;
Cigar_record* buffer;
}Cigar_record_alloc;
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;
uint32_t cov;
///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;
uint32_t homopolymer_num;
uint32_t non_homopolymer_num;
int score;
////the position of snp in read itself
uint32_t site;
uint8_t is_homopolymer;
}
SnpStats;
typedef struct
{
uint32_t beg;
uint32_t end;
uint32_t occ_0;
uint32_t occ_1;
uint32_t homopolymer_num;
uint32_t non_homopolymer_num;
uint32_t is_remove;
}
Snp_ID_Vector;
typedef struct
{
long long IDs_size;
long long IDs_length;
long long max_snp_id;
Snp_ID_Vector* IDs;
long long buffer_size;
long long buffer_length;
uint32_t* buffer;
}
Snp_ID_Vector_Alloc;
#define Get_DP_Backtrack_Column(matrix, i) (matrix.backtrack + matrix.snp_num * i)
#define Get_DP_Backtrack_Column_Length(matrix, i) (matrix.snp_num)
typedef struct
{
// uint32_t snp_size;
// uint32_t snp_num;
// uint32_t* max;
// uint32_t* colum_len;
// uint32_t* colum;
// uint32_t matrix_size;
uint32_t snp_num;
uint8_t* visit;
uint32_t* max;
uint64_t* max_for_sort;
uint32_t snp_size;
uint32_t* backtrack_length;
uint32_t* backtrack;
uint32_t backtrack_size;
uint32_t* buffer;
uint32_t* max_buffer;
int max_snp_num;
///int max_snp_ID;
int max_score;
int current_snp_num;
Snp_ID_Vector_Alloc SNP_IDs;
}
DP_matrix;
#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)
#define Get_Result_SNP_Vector(matrix) (matrix.r_snp)
typedef struct
{
SnpStats* a;
size_t n,m;
}kv_SnpStats_t;
typedef struct
{
haplotype_evdience* list;
uint32_t sub_list_start;
uint32_t sub_list_length;
uint32_t length;
uint32_t size;
/****************************may have bugs********************************/
uint8_t flag[WINDOW_MAX_SIZE];
/****************************may have bugs********************************/
uint32_t core_snp;
uint32_t overlap;
int8_t *snp_matrix;
uint32_t snp_matrix_size;
int8_t *r_snp;
uint32_t r_snp_size;
SnpStats result_stat;
kv_SnpStats_t snp_stat;
uint32_t nn_snp;
kvec_t(uint64_t) snp_srt;
// SnpStats* snp_stat;
// uint32_t snp;
// uint32_t snp_stat_size;
// uint32_t available_snp;
DP_matrix dp;
}
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)
///if(available < threshold || total < 10)
{
return 0;
}
return 1;
}
inline void count_nearby_snps(haplotype_evdience_alloc* hap, uint32_t* SNPs, int SNPsLen, int* nearsnp, int* non_nearsnps)
{
long long i, current_id, large_id, small_id;
long long distance = 10;
if(SNPsLen == 1)
{
(*non_nearsnps) = 1;
(*nearsnp) = 0;
return;
}
if(SNPsLen == 0)
{
(*nearsnp) = 0;
(*non_nearsnps) = 0;
return;
}
(*nearsnp) = 0;
(*non_nearsnps) = 0;
for (i = 0; i < SNPsLen; i++)
{
if(i > 0 && i < SNPsLen - 1)
{
current_id= SNPs[i];
///since SNPs[i - 1].site is larger than SNPs[i]
large_id = SNPs[i - 1];
small_id = SNPs[i + 1];
if(hap->snp_stat.a[large_id].site - hap->snp_stat.a[current_id].site < distance
||
hap->snp_stat.a[current_id].site - hap->snp_stat.a[small_id].site < distance)
{
(*nearsnp)++;
}
else
{
(*non_nearsnps)++;
}
if(hap->snp_stat.a[current_id].site > hap->snp_stat.a[large_id].site ||
hap->snp_stat.a[current_id].site < hap->snp_stat.a[small_id].site)
{
fprintf(stderr, "error\n");
}
}
else if(i == 0)
{
current_id= SNPs[i];
small_id = SNPs[i + 1];
if(hap->snp_stat.a[current_id].site - hap->snp_stat.a[small_id].site < distance)
{
(*nearsnp)++;
}
else
{
(*non_nearsnps)++;
}
if(hap->snp_stat.a[current_id].site < hap->snp_stat.a[small_id].site)
{
fprintf(stderr, "error\n");
}
}
else
{
large_id = SNPs[i - 1];
current_id= SNPs[i];
if(hap->snp_stat.a[large_id].site - hap->snp_stat.a[current_id].site < distance)
{
(*nearsnp)++;
}
else
{
(*non_nearsnps)++;
}
if(hap->snp_stat.a[current_id].site > hap->snp_stat.a[large_id].site)
{
fprintf(stderr, "error\n");
}
}
}
if((*nearsnp) + (*non_nearsnps) != SNPsLen)
{
fprintf(stderr, "(*nearsnp): %d, (*non_nearsnps): %d, SNPsLen: %d\n", (*nearsnp), (*non_nearsnps), SNPsLen);
}
}
inline int filter_one_snp_advance_nearby(haplotype_evdience_alloc* hap, int occ_0, int occ_1, int total,
long long homopolymer_num, long long non_homopolymer_num,
uint32_t* SNPs, int SNPsLen)
{
double available;
if(occ_0 <= occ_1)
{
available = occ_0;
}
else
{
available = occ_1;
}
int min = available;
double threshold1 = 0.35;
double threshold2 = 0.24;
available = available/((double)(total));
///if((non_homopolymer_num > 0 && min >= 5) || (min >= 6))
if(min >= 5)
{
if(available < threshold2 || total < 10)
{
return 0;
}
}
else if(available < threshold1 || occ_0 < MIN_COVERAGE_THRESHOLD + 1 || total < 10)
{
return 0;
}
return 1;
}
inline int if_is_homopolymer_strict(long long site, char* read, long long read_length)
{
long long beg, end, i;
long long threshold = 3;
beg = site - threshold;
if(beg < 0)
{
beg = 0;
}
end = site + threshold;
if(end >= read_length)
{
end = read_length - 1;
}
char f_homopolymer_ch = 0;
long long f_homopolymer_len = 0;
for (i = site + 1; i <= end; i++)
{
if(f_homopolymer_ch == 0)
{
f_homopolymer_ch = read[i];
f_homopolymer_len = 1;
}
else
{
if(read[i] != f_homopolymer_ch)
{
break;
}
else
{
f_homopolymer_len++;
}
}
}
char b_homopolymer_ch = 0;
long long b_homopolymer_len = 0;
for (i = site - 1; i >= beg; i--)
{
if(b_homopolymer_ch == 0)
{
b_homopolymer_ch = read[i];
b_homopolymer_len = 1;
}
else
{
if(read[i] != b_homopolymer_ch)
{
break;
}
else
{
b_homopolymer_len++;
}
}
}
if(f_homopolymer_ch == read[site])
{
f_homopolymer_len++;
}
else if(b_homopolymer_ch == read[site])
{
b_homopolymer_len++;
}
if(f_homopolymer_len >= threshold || b_homopolymer_len >= threshold)
{
return 1;
}
if (read[site] == f_homopolymer_ch
&&
b_homopolymer_ch == f_homopolymer_ch
&&
(f_homopolymer_len + b_homopolymer_len >= threshold))
{
return 1;
}
return 0;
}
inline int if_is_homopolymer_repeat(long long site, char* read, long long read_length)
{
long long beg, end, i;
long long threshold = 3;
beg = site - threshold;
if(beg < 0)
{
beg = 0;
}
end = site + threshold;
if(end >= read_length)
{
end = read_length - 1;
}
char f_homopolymer_ch = 0;
long long f_homopolymer_len = 0;
for (i = site + 1; i <= end; i++)
{
if(f_homopolymer_ch == 0)
{
f_homopolymer_ch = read[i];
f_homopolymer_len = 1;
}
else
{
if(read[i] != f_homopolymer_ch)
{
break;
}
else
{
f_homopolymer_len++;
}
}
}
char b_homopolymer_ch = 0;
long long b_homopolymer_len = 0;
for (i = site - 1; i >= beg; i--)
{
if(b_homopolymer_ch == 0)
{
b_homopolymer_ch = read[i];
b_homopolymer_len = 1;
}
else
{
if(read[i] != b_homopolymer_ch)
{
break;
}
else
{
b_homopolymer_len++;
}
}
}
if(f_homopolymer_ch == read[site])
{
f_homopolymer_len++;
}
else if(b_homopolymer_ch == read[site])
{
b_homopolymer_len++;
}
if(f_homopolymer_len >= threshold || b_homopolymer_len >= threshold)
{
return 1;
}
if (read[site] == f_homopolymer_ch
&&
b_homopolymer_ch == f_homopolymer_ch
&&
(f_homopolymer_len + b_homopolymer_len >= threshold))
{
return 1;
}
return 0;
}
inline void InsertSNPVector(haplotype_evdience_alloc* h, haplotype_evdience* sub_list, long long sub_length, char misBase,
UC_Read* g_read)
{
if(sub_length <= 0)
return;
long long i = 0; SnpStats *p = NULL;
kv_pushp(SnpStats, h->snp_stat, &p);
// 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;
// h->snp_stat[h->available_snp].is_homopolymer =
// if_is_homopolymer_strict(h->snp_stat[h->available_snp].site, g_read->seq, g_read->length);
// int8_t* vector = Get_SNP_Vector((*h), h->available_snp);
p->id = h->snp_stat.n-1;
p->occ_0 = 0;
p->occ_1 = 0;
p->occ_2 = 0;
p->overlap_num = 0;
p->site = sub_list[0].site;
p->is_homopolymer = if_is_homopolymer_strict(p->site, g_read->seq, g_read->length);
int8_t* vector = Get_SNP_Vector((*h), p->id);
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++;
h->snp_stat.a[p->id].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++;
h->snp_stat.a[p->id].occ_1++;
}
else
{
vector[sub_list[i].overlapID] = 2;
// h->snp_stat[h->available_snp].occ_2++;
h->snp_stat.a[p->id].occ_2++;
}
// h->snp_stat[h->available_snp].overlap_num++;
h->snp_stat.a[p->id].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;
int new_occ_0 = h->snp_stat.a[p->id].occ_0 + 1;
int new_occ_1 = h->snp_stat.a[p->id].occ_1;
if(filter_snp(new_occ_0, new_occ_1, new_occ_0 + new_occ_1) == 0) ///Fix-attention:definitely wrong
{
// h->snp_stat[h->available_snp].score = -1;
h->snp_stat.a[p->id].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%///Fix-attention:definitely wrong
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->snp_stat.a[p->id].score = consensus;
}
// h->available_snp++;
}
inline int calculate_score(int new_occ_0, int new_occ_1)
{
if(new_occ_0 + new_occ_1 == 0)
{
return -1;
}
if(filter_snp(new_occ_0, new_occ_1, new_occ_0 + new_occ_1) == 0)
{
return -1;
}
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));
return consensus;
}
inline void SetSnpMatrix(haplotype_evdience_alloc* h, uint32_t *nn_snp, uint64_t *overlap_num, int32_t set_matrix, void *km)
{
if(nn_snp && overlap_num) {
if(!km) kv_resize(SnpStats, h->snp_stat, *nn_snp);
else kv_resize_km(km, SnpStats, h->snp_stat, *nn_snp);
h->snp_stat.n = 0; h->overlap = *overlap_num; h->core_snp = 0;
}
if(set_matrix) {
uint64_t n_snp = nn_snp? *nn_snp:h->snp_stat.n;
uint64_t n_ovlp = overlap_num? *overlap_num:h->overlap;
uint64_t new_size = n_snp* n_ovlp;
if(h->snp_matrix_size < new_size) {
h->snp_matrix_size = new_size;
if(!km) REALLOC(h->snp_matrix, h->snp_matrix_size);
else KREALLOC(km, h->snp_matrix, h->snp_matrix_size);
}
memset(h->snp_matrix, -1, n_snp * n_ovlp);
if(h->r_snp_size < n_ovlp) {
h->r_snp_size = n_ovlp;
if(!km) REALLOC(h->r_snp, h->r_snp_size);
else KREALLOC(km, h->r_snp, h->r_snp_size);
}
}
}
inline void init_SNP_IDs(Snp_ID_Vector_Alloc* SNP_IDs)
{
SNP_IDs->max_snp_id = -1;
SNP_IDs->buffer_length = 0;
SNP_IDs->buffer_size = 1000;
SNP_IDs->buffer = (uint32_t*)malloc(sizeof(uint32_t) * SNP_IDs->buffer_size);
SNP_IDs->IDs_length = 0;
SNP_IDs->IDs_size = 10;
SNP_IDs->IDs = (Snp_ID_Vector*)calloc(SNP_IDs->IDs_size, sizeof(Snp_ID_Vector));
}
inline void clear_SNP_IDs(Snp_ID_Vector_Alloc* SNP_IDs)
{
SNP_IDs->IDs_length = 0;
SNP_IDs->buffer_length = 0;
SNP_IDs->max_snp_id = -1;
}
inline void destory_SNP_IDs(Snp_ID_Vector_Alloc* SNP_IDs)
{
free(SNP_IDs->buffer);
free(SNP_IDs->IDs);
}
inline void insert_SNP_IDs_addition(Snp_ID_Vector_Alloc* SNP_IDs, uint32_t* IDs_vec, int IDs_vec_length,
uint32_t occ_0, uint32_t occ_1, uint32_t homopolymer_num, uint32_t non_homopolymer_num)
{
if(SNP_IDs->IDs_length + 1 > SNP_IDs->IDs_size)
{
SNP_IDs->IDs_size = SNP_IDs->IDs_size * 2;
SNP_IDs->IDs = (Snp_ID_Vector*)realloc(SNP_IDs->IDs, SNP_IDs->IDs_size * sizeof(Snp_ID_Vector));
}
SNP_IDs->IDs[SNP_IDs->IDs_length].beg = SNP_IDs->buffer_length;
SNP_IDs->IDs[SNP_IDs->IDs_length].end = SNP_IDs->buffer_length + IDs_vec_length - 1;
SNP_IDs->IDs[SNP_IDs->IDs_length].occ_0 = occ_0;
SNP_IDs->IDs[SNP_IDs->IDs_length].occ_1 = occ_1;
SNP_IDs->IDs[SNP_IDs->IDs_length].homopolymer_num = homopolymer_num;
SNP_IDs->IDs[SNP_IDs->IDs_length].non_homopolymer_num = non_homopolymer_num;
if(SNP_IDs->buffer_length + IDs_vec_length > SNP_IDs->buffer_size)
{
SNP_IDs->buffer_size = (SNP_IDs->buffer_length + IDs_vec_length) * 2;
SNP_IDs->buffer = (uint32_t*)realloc(SNP_IDs->buffer, SNP_IDs->buffer_size * sizeof(uint32_t));
}
memcpy(SNP_IDs->buffer + SNP_IDs->buffer_length, IDs_vec, IDs_vec_length*sizeof(uint32_t));
SNP_IDs->IDs_length++;
SNP_IDs->buffer_length += IDs_vec_length;
}
inline void insert_SNP_IDs_addition(Snp_ID_Vector_Alloc* SNP_IDs, uint32_t* IDs_vec, int IDs_vec_length)
{
if(SNP_IDs->IDs_length + 1 > SNP_IDs->IDs_size)
{
SNP_IDs->IDs_size = SNP_IDs->IDs_size * 2;
SNP_IDs->IDs = (Snp_ID_Vector*)realloc(SNP_IDs->IDs, SNP_IDs->IDs_size * sizeof(Snp_ID_Vector));
}
SNP_IDs->IDs[SNP_IDs->IDs_length].beg = SNP_IDs->buffer_length;
SNP_IDs->IDs[SNP_IDs->IDs_length].end = SNP_IDs->buffer_length + IDs_vec_length - 1;
if(SNP_IDs->buffer_length + IDs_vec_length > SNP_IDs->buffer_size)
{
SNP_IDs->buffer_size = (SNP_IDs->buffer_length + IDs_vec_length) * 2;
SNP_IDs->buffer = (uint32_t*)realloc(SNP_IDs->buffer, SNP_IDs->buffer_size * sizeof(uint32_t));
}
memcpy(SNP_IDs->buffer + SNP_IDs->buffer_length, IDs_vec, IDs_vec_length*sizeof(uint32_t));
SNP_IDs->IDs_length++;
SNP_IDs->buffer_length += IDs_vec_length;
}
inline void init_DP_matrix(DP_matrix* dp, uint32_t snp_num)
{
if(snp_num > dp->snp_size)
{
dp->snp_size = snp_num;
dp->max = (uint32_t*)realloc(dp->max, dp->snp_size * sizeof(uint32_t));
dp->max_for_sort = (uint64_t*)realloc(dp->max_for_sort, dp->snp_size * sizeof(uint64_t));
dp->visit = (uint8_t*)realloc(dp->visit, dp->snp_size * sizeof(uint8_t));
dp->buffer = (uint32_t*)realloc(dp->buffer, dp->snp_size * sizeof(uint32_t));
dp->max_buffer = (uint32_t*)realloc(dp->max_buffer, dp->snp_size * sizeof(uint32_t));
dp->backtrack_length = (uint32_t*)realloc(dp->backtrack_length, dp->snp_size * sizeof(uint32_t));
dp->backtrack_size = snp_num*snp_num;
dp->backtrack = (uint32_t*)realloc(dp->backtrack, dp->backtrack_size * sizeof(uint32_t));
}
dp->snp_num = snp_num;
clear_SNP_IDs(&(dp->SNP_IDs));
}
inline void InitHaplotypeEvdience_buf(haplotype_evdience_alloc* h, void *km)
{
memset(h, 0, sizeof(haplotype_evdience_alloc));
/****************************may have bugs********************************/
memset(h->flag, 0, WINDOW_MAX_SIZE * sizeof(uint8_t));
/****************************may have bugs********************************/
// init_SNP_IDs(&(h->dp.SNP_IDs));
memset(&(h->dp.SNP_IDs), 0, sizeof(h->dp.SNP_IDs));
h->dp.SNP_IDs.max_snp_id = -1;
}
inline void InitHaplotypeEvdience(haplotype_evdience_alloc* h)
{
h->overlap = 0;
h->snp_matrix_size = 0;
h->snp_matrix = NULL;
h->r_snp_size = 0;
h->r_snp = NULL;
kv_init(h->snp_stat);
kv_init(h->snp_srt);
h->nn_snp = 0;
// h->snp_stat = NULL;
// h->snp = 0;
// h->snp_stat_size = 0;
// h->available_snp = 0;
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));
/****************************may have bugs********************************/
memset(h->flag, 0, WINDOW_MAX_SIZE * sizeof(uint8_t));
/****************************may have bugs********************************/
// h->dp.max = NULL;
// h->dp.colum = NULL;
// h->dp.colum_len = NULL;
// h->dp.matrix_size = 0;
// h->dp.snp_num = 0;
// h->dp.snp_size = 0;
h->dp.snp_num = 0;
h->dp.max = NULL;
h->dp.max_for_sort = NULL;
h->dp.visit = NULL;
h->dp.snp_size = 0;
h->dp.backtrack = NULL;
h->dp.backtrack_size = 0;
h->dp.backtrack_length = NULL;
h->dp.buffer = NULL;
h->dp.max_buffer = NULL;
init_SNP_IDs(&(h->dp.SNP_IDs));
}
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);
kv_destroy(h->snp_stat);
kv_destroy(h->snp_srt);
free(h->snp_matrix);
free(h->r_snp);
free(h->dp.backtrack);
free(h->dp.max);
free(h->dp.max_for_sort);
free(h->dp.visit);
free(h->dp.backtrack_length);
free(h->dp.buffer);
free(h->dp.max_buffer);
destory_SNP_IDs(&(h->dp.SNP_IDs));
}
inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h)
{
// h->snp = 0;
h->nn_snp = 0;
h->snp_stat.n = 0;
h->length = 0;
h->sub_list_start = 0;
h->sub_list_length = 0;
/****************************may have bugs********************************/
memset(h->flag, 0, WINDOW_MAX_SIZE * sizeof(uint8_t));
/****************************may have bugs********************************/
}
inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h, long long beg, long long useful_length)
{
/****************************may have bugs********************************/
memset(h->flag + beg, 0, useful_length * sizeof(uint8_t));
/****************************may have bugs********************************/
}
inline void addHaplotypeEvdience(haplotype_evdience_alloc* h, haplotype_evdience* ev, void *km)
{
if(h->length + 1 > h->size){
h->size = h->length + 1;
kroundup32(h->size);
if(!km) REALLOC(h->list, h->size);
else KREALLOC(km, h->list, h->size);
// h->list = (haplotype_evdience*)realloc(h->list, sizeof(haplotype_evdience)*h->size);
}
h->list[h->length] = (*ev);
h->length++;
}
typedef struct
{
char* corrected_read;
long long corrected_read_length;
long long last_boundary_length;
long long corrected_read_size;
long long corrected_base;
uint64_t* overlapID;
uint64_t length;
uint64_t lengthNT;
uint64_t size;
uint64_t start_i;
/****************************may have bugs********************************/
// char overlap_region[WINDOW + THRESHOLD*2 + 10];
// char overlap_region_group[GROUP_SIZE][WINDOW + THRESHOLD*2 + 10];
// char path[WINDOW + THRESHOLD*2 + 10];
// Word matrix_bit[((WINDOW + 10)<<3)];
char overlap_region[WINDOW_MAX_SIZE + THRESHOLD_MAX_SIZE*2 + 10];
char overlap_region_group[GROUP_SIZE][WINDOW_MAX_SIZE + THRESHOLD_MAX_SIZE*2 + 10];
char path[WINDOW_MAX_SIZE + THRESHOLD_MAX_SIZE*2 + 10];
char path_fix[WINDOW_MAX_SIZE + THRESHOLD_MAX_SIZE*2 + 10];
char overlap_region_fix[WINDOW_MAX_SIZE + THRESHOLD_MAX_SIZE*2 + 10];
Word matrix_bit[((WINDOW_MAX_SIZE + 10)<<3)];
/****************************may have bugs********************************/
int path_length;
__m128i Peq_SSE[256];
} Correct_dumy;
typedef struct
{
Correct_dumy dumy;
Cigar_record cigar;
Cigar_record tmp_cigar;
long long obtained_cigar_length;
}
Round2_alignment;
void init_Round2_alignment_buf(Round2_alignment* h, void *km);
void init_Round2_alignment(Round2_alignment* h);
void destory_Round2_alignment(Round2_alignment* h);
void clear_Round2_alignment(Round2_alignment* h);
typedef struct {
int32_t c_qs, c_qe, c_ts, c_te; //[c_qs, c_qe) && [c_ts, c_te)
int32_t c_wsid, c_weid, c_wsii, c_weii;//[c_wsid, c_weid] && [c_wsii, c_weii)
int32_t rev, sfx_e, pfx_e, mid_e, oid;
} rtrace_t;
typedef struct {
rtrace_t *a;
size_t n, m;
} kv_rtrace_t;
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, Graph* DAGCon,
Cigar_record* current_cigar, haplotype_evdience_alloc* hap,
Round2_alignment* second_round, kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal);
void init_Correct_dumy_buf(Correct_dumy* list, void *km);
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, void *km);
void clear_Correct_dumy_pure(Correct_dumy* list);
void get_seq_from_Graph(Graph* backbone, Graph* DAGCon, Correct_dumy* dumy, Cigar_record* current_cigar, char* self_string,
char* r_string, long long r_string_length, long long r_string_site);
void init_Cigar_record(Cigar_record* dummy);
void init_Cigar_record_buf(Cigar_record* dummy, void *km);
void destory_Cigar_record(Cigar_record* dummy);
void clear_Cigar_record(Cigar_record* dummy);
void add_new_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type);
void add_existing_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type);
void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq);
void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq);
void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read,
Graph* g, Graph* DAGCon, Cigar_record* current_cigar,
haplotype_evdience_alloc* hap, Round2_alignment* second_round,
kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km);
void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt, char *qstr,
uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz,
haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, overlap_region *aux_o,
double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t hpc_k, st_mt_t *stb, idx_emask_t *mm, mask_ul_ov_t *mk, void *km);
void ul_lalign_old_ed(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr,
uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy,
haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx,
double e_rate, int64_t wl, uint64_t is_base, void *km);
void lchain_align(overlap_region_alloc* overlap_list, const ul_idx_t *uref,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read,
Graph* g, Graph* DAGCon, Cigar_record* current_cigar,
haplotype_evdience_alloc* hap, Round2_alignment* second_round,
kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal,
double max_ov_diff_ec, long long winLen, void *km);
void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt,
char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz,
overlap_region *aux_o, double e_rate, int64_t wl, int64_t sid, uint64_t khit, uint64_t chain_cut,
void *km);
/***
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)
{
if(type == 0)///match
{
///add to existing cell, just increase length
if((uint32_t)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((uint32_t)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((uint32_t)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((uint32_t)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;
}
int verify_cigar_2(char* x, int x_len, char* y, int y_len, Cigar_record* cigar, int error);
int verify_single_window(long long x_start, long long x_end,
long long overlap_x_s, long long overlap_y_s, int x_id,
int y_id, int y_strand, char* x_buffer, char* y_buffer,
All_reads* R_INF);
void init_Cigar_record_alloc(Cigar_record_alloc* x);
void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size);
void destory_Cigar_record_alloc(Cigar_record_alloc* x);
void afine_gap_alignment(const char *qseq, uint8_t* qnum, const int ql,
const char *tseq, uint8_t* tnum, const int tl, const uint8_t *c2n, const int strand,
int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus,
long long* max_q_pos, long long* max_t_pos, long long* global_score,
long long* extention_score, long long* q_boundary_score, long long* q_boundary_t_coordinate,
long long* t_boundary_score, long long* t_boundary_q_coordinate,
long long* droped, int mode);
// void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF,
// UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp);
inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
{
long long overlapLen = xEnd - xBeg + 1;
if(yEnd - yBeg + 1 == overlapLen)
{
if(memcmp(x+xBeg, y+yBeg, overlapLen)==0) return 1;
}
return 0;
}
inline void get_cigar_cell(window_list *idx, window_list_alloc *cc, uint32_t i, uint8_t *c, uint32_t *len)
{
uint16_t p = cc->c.a[idx->cidx+i];
(*c) = (uint8_t)(p>>14); (*len) = (p&((uint16_t)0x3fff));
}
inline void push_cigar_cell(window_list_alloc *res, uint8_t c, uint32_t len)
{
uint16_t p = c; p <<= 14; p += (uint16_t)len;
kv_push(uint16_t, res->c, p);
}
int64_t get_num_wins(int64_t s, int64_t e, int64_t block_s);
void append_unmatched_wins(overlap_region *z, int64_t block_s);
///[w_s, w_e]
inline int64_t get_win_id_by_s(overlap_region *z, int64_t w_s, int64_t block_s, int64_t *w_e)
{
int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_s-n_s)/block_s;
if(w_e) {
(*w_e) = n_s + (wid+1)*block_s - 1;
if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e;
}
return wid;
}
///[w_s, w_e]
inline int64_t get_win_id_by_e(overlap_region *z, int64_t w_e, int64_t block_s, int64_t *w_s)
{
int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_e-n_s)/block_s;
if(w_s) {
(*w_s) = n_s + wid*block_s;
if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s;
}
return wid;
}
///[w_s, w_e]
inline void get_win_se_by_normalize_xs(overlap_region *z, int64_t norm_w_s, int64_t block_s, int64_t *w_s, int64_t *w_e)
{
int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (norm_w_s-n_s)/block_s;
if(w_s) {
(*w_s) = n_s + wid*block_s;
if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s;
}
if(w_e) {
(*w_e) = n_s + (wid+1)*block_s - 1;
if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e;
}
}
void inline resize_UC_Read(UC_Read *z, int64_t s)
{
if(z->size < s) {
REALLOC(z->seq, s); z->size = s;
}
}
void update_sketch_trace(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt,
All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl,
kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate);
int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug);
void convert_ul_ov_t(ul_ov_t *des, overlap_region *src, ma_ug_t *ug);
uint64_t check_connect_ug(const ul_idx_t *uref, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq);
uint64_t check_connect_rg(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t uv, uint32_t uw, int64_t bw, double diff_ec_ul, int64_t dq);
uint32_t govlp_check(const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, ul_ov_t *li, ul_ov_t *lj);
void ul_rid_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt,
char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate,
int64_t wl, kv_ul_ov_t *aln, kv_ul_ov_t *cln, kv_rtrace_t *trace, int64_t sid, uint64_t khit, void *km);
#define copy_asg_arr(des, src) ((des).a = (src).a, (des).n = (src).n, (des).m = (src).m)
#define is_ualn_win(a) (((a).error==INT16_MAX)&&((a).clen==0)&&((a).extra_end<0))
#define is_exact_aln(a) (((a).error<INT16_MAX)&&((a).clen>0))
#define is_est_aln(a) (((a).error<INT16_MAX)&&((a).clen==0))
#define FORWARD_KSW 0
#define BACKWARD_KSW 1
#define MATCH_SCORE_KSW 2
#define MISMATCH_SCORE_KSW 4
#define GAP_OPEN_KSW 4
#define GAP_EXT_KSW 2
#define Z_DROP_KSW 400
#define BAND_KSW 500
#define set_bit_extz_t(x, z, id) do {\
(x).cigar.a = (z).w_list.c.a+(z).w_list.a[(id)].cidx;\
(x).cigar.n = (x).cigar.m = (z).w_list.a[(id)].clen;\
(x).ts = (z).w_list.a[(id)].x_start;\
(x).te = (z).w_list.a[(id)].x_end;\
(x).ps = (z).w_list.a[(id)].y_start;\
(x).pe = (z).w_list.a[(id)].y_end;\
(x).err = (z).w_list.a[(id)].error;\
(x).thre = (z).w_list.a[(id)].error;\
} while (0)
typedef struct {
int64_t k, q[2], t[2], cq[2], ct[2], ci[2], werr, werr0, cerr;
int64_t qoff, f, toff, coff, cur_qoff;
} rtrace_iter;
int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *trace, rtrace_t *tc,
const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o,
bit_extz_t *exz, double e_rate, int64_t qs);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs);
void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
void push_alnw(overlap_region *aux_o, bit_extz_t *exz);
void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez);
#define ovlp_id(x) ((x).tn)
#define ovlp_min_wid(x) ((x).ts)
#define ovlp_max_wid(x) ((x).te)
#define ovlp_cur_wid(x) ((x).qn)
#define ovlp_cur_xoff(x) ((x).qs)
#define ovlp_cur_yoff(x) ((x).ts)
#define ovlp_cur_ylen(x) ((x).te)
#define ovlp_cur_coff(x) ((x).qe)
#define ovlp_bd(x) ((x).sec)
#define HPC_PL 12
#define HPC_RR 4
#define HPC_CC 2
#endif
Reference in New Issue View Git Blame Copy Permalink