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

clean code

Browse Source
This commit is contained in:
Haoyu Cheng
2019-12-26 02:26:24 -05:00
parent 7f9c73507d
commit 4ac9907e49
26 changed files with 1585 additions and 19325 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+205 -1793
View File
File diff suppressed because it is too large Load Diff
Assembly.h
+1 -6
View File
@@ -1,5 +1,6 @@
#ifndef __ASSEMBLY__
#define __ASSEMBLY__
#include "CommandLines.h"
#define FORWARD 0
#define REVERSE_COMPLEMENT (0x8000000000000000)
@@ -12,10 +13,4 @@ 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();
/********************************for debug*****************************************/
void verify_Position_hash_table();
#endif
CommandLines.cpp
+178 -61
View File
@@ -4,23 +4,9 @@
#include "ketopt.h"
#include <sys/time.h>
#define VERSION "0.0.0.1"
char* read_file_name = NULL;
char* output_file_name = NULL;
int thread_num = 1;
int k_mer_length = 40;
int coverage = -1;
//int k_mer_min_freq = 9;
int k_mer_min_freq = 3;
//int k_mer_min_freq = 2;
int k_mer_max_freq = 66;
int load_index_from_disk = 0;
int write_index_to_disk = 0;
int number_of_round = 1;
int read_graph = 0;
int c_round = 4;
int adapterLen = 0;
char* required_read_name = NULL;
hifiasm_opt_t asm_opt;
double Get_T(void)
{
@@ -29,66 +15,197 @@ double Get_T(void)
return t.tv_sec+t.tv_usec/1000000.0;
}
void Print_H()
void Print_H(hifiasm_opt_t* asm_opt)
{
fprintf(stderr, "Incorrect options.\n");
fprintf(stderr, "./ccs_assembly -w -l -q NA12878_chr1_10M.fq -o output_NA12878_v11_2.fa -c 31 -k 40 -t 32 -r 2 -a 4 -z 0\n");
fprintf(stderr, "Usage: hifiasm [options] -q <input.fa> -o <output_asm>\n");
fprintf(stderr, "Options:\n");
fprintf(stderr, " -q FILE input in the fastq(.gz)/fasta(.gz) formats\n");
fprintf(stderr, " -k FILE output assembly (in gfa format) and corrected reads (in fasta format)\n");
fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num);
fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round);
fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round);
fprintf(stderr, " -k INT k-mer length [%d] (must be < 64)\n", asm_opt->k_mer_length);
fprintf(stderr, " -w write all overlaps to disk, can accelerate assembly next time\n");
fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation\n");
fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
fprintf(stderr, " -p INT size of popped bubbles [%lld]\n", asm_opt->pop_bubble_size);
fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
fprintf(stderr, " -v show version number\n");
fprintf(stderr, " -h show help information\n");
fprintf(stderr, "Example: ./hifiasm -w -l -q NA12878.fq.gz -o NA12878.asm -k 40 -t 32 -r 2 -a 4 -z 0\n");
}
int CommandLine_process (int argc, char *argv[])
void init_opt(hifiasm_opt_t* asm_opt)
{
asm_opt->read_file_name = NULL;
asm_opt->output_file_name = NULL;
asm_opt->required_read_name = NULL;
asm_opt->thread_num = 1;
asm_opt->k_mer_length = 40;
asm_opt->k_mer_min_freq = 3;
asm_opt->k_mer_max_freq = 66;
asm_opt->load_index_from_disk = 0;
asm_opt->write_index_to_disk = 0;
asm_opt->number_of_round = 2;
asm_opt->adapterLen = 0;
asm_opt->clean_round = 4;
asm_opt->complete_threads = 0;
asm_opt->pop_bubble_size = 100000;
asm_opt->min_drop_rate = 0.2;
asm_opt->max_drop_rate = 0.8;
}
static ko_longopt_t longopts[] = {
{ "help", ko_no_argument, 100},
{ "seq", ko_required_argument, 101},
{ "output", ko_required_argument, 102},
{ "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 }
};
void clear_opt(hifiasm_opt_t* asm_opt, int last_round)
{
asm_opt->complete_threads = 0;
asm_opt->num_bases = 0;
asm_opt->num_corrected_bases = 0;
asm_opt->num_recorrected_bases = 0;
asm_opt->roundID = asm_opt->number_of_round - last_round;
}
ketopt_t opt = KETOPT_INIT;
int i, c;
while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:r:c:a:b:z:", 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;
else if (c == 'c') coverage = atoi(opt.arg);
else if (c == 'a') c_round = atoi(opt.arg);
else if (c == 'z') adapterLen = atoi(opt.arg);
else if (c == 'b') required_read_name = opt.arg;
else if (c == '?') printf("unknown opt: -%c\n", opt.opt? opt.opt : ':');
else if (c == ':') printf("missing arg: -%c\n", opt.opt? opt.opt : ':');
int check_option(hifiasm_opt_t* asm_opt)
{
if(asm_opt->read_file_name == NULL)
{
fprintf(stderr, "[ERROR] missing input: please specify a read file\n");
return 0;
}
if(asm_opt->output_file_name == NULL)
{
fprintf(stderr, "[ERROR] missing output: please specify the output name\n");
return 0;
}
if(asm_opt->thread_num < 1)
{
fprintf(stderr, "[ERROR] the number of threads must be > 0\n");
return 0;
}
if(asm_opt->number_of_round < 1)
{
fprintf(stderr, "[ERROR] the number of rounds for correction must be > 0\n");
return 0;
}
if(asm_opt->clean_round < 1)
{
fprintf(stderr, "[ERROR] the number of rounds for assembly cleaning must be > 0\n");
return 0;
}
if(asm_opt->adapterLen < 0)
{
fprintf(stderr, "[ERROR] the length of removed adapters must be >= 0\n");
return 0;
}
if(asm_opt->k_mer_length >= 64)
{
fprintf(stderr, "[ERROR] the length of k_mer must be < 64\n");
return 0;
}
if(asm_opt->max_drop_rate < 0 || asm_opt->max_drop_rate >= 1 )
{
fprintf(stderr, "[ERROR] max overlap drop ratio must be [0.0, 1.0)\n");
return 0;
}
if(asm_opt->min_drop_rate < 0 || asm_opt->min_drop_rate >= 1)
{
fprintf(stderr, "[ERROR] min overlap drop ratio must be [0.0, 1.0)\n");
return 0;
}
if(asm_opt->max_drop_rate <= asm_opt->min_drop_rate)
{
fprintf(stderr, "[ERROR] min overlap drop ratio must be less than max overlap drop ratio\n");
return 0;
}
if(asm_opt->pop_bubble_size < 0)
{
fprintf(stderr, "[ERROR] the size of popped bubbles must be >= 0\n");
return 0;
}
// fprintf(stderr, "input file: %s\n", asm_opt->read_file_name);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
// fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
// fprintf(stderr, "number of rounds for correction: %d\n", asm_opt->number_of_round);
// fprintf(stderr, "number of rounds for assembly cleaning: %d\n", asm_opt->clean_round);
// fprintf(stderr, "length of removed adapters: %d\n", asm_opt->adapterLen);
// fprintf(stderr, "length of k_mer: %d\n", asm_opt->k_mer_length);
// fprintf(stderr, "min overlap drop ratio: %.2g\n", asm_opt->min_drop_rate);
// fprintf(stderr, "max overlap drop ratio: %.2g\n", asm_opt->max_drop_rate);
// fprintf(stderr, "size of popped bubbles: %lld\n", asm_opt->pop_bubble_size);
return 1;
}
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
{
ketopt_t opt = KETOPT_INIT;
int c;
while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:q:k:lwm:n:r:a:b:z:x:y:p:", 0)) >= 0) {
if (c == 'h')
{
Print_H(asm_opt);
return 0;
}
else if (c == 'v')
{
fprintf(stderr, "[Version] %s\n", VERSION);
return 0;
}
else if (c == 't') asm_opt->thread_num = atoi(opt.arg);
else if (c == 'o') asm_opt->output_file_name = opt.arg;
else if (c == 'q') asm_opt->read_file_name = opt.arg;
else if (c == 'n') asm_opt->k_mer_min_freq = atoi(opt.arg);
else if (c == 'm') asm_opt->k_mer_max_freq = atoi(opt.arg);
else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg);
else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg);
else if (c == 'l') asm_opt->load_index_from_disk = 1;
else if (c == 'w') asm_opt->write_index_to_disk = 1;
else if (c == 'a') asm_opt->clean_round = atoi(opt.arg);
else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg);
else if (c == 'b') asm_opt->required_read_name = opt.arg;
else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg);
else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg);
else if (c == 'p') asm_opt->pop_bubble_size = atoll(opt.arg);
else if (c == ':')
{
fprintf(stderr, "[ERROR] missing option argument in \"%s\"\n", argv[opt.i - 1]);
return 0;
}
else if (c == '?')
{
fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[opt.i - 1]);
return 0;
}
}
if (argc == 1)
{
Print_H();
return 0;
}
if(coverage == -1)
{
fprintf(stdout, "Please set -c!\n");
Print_H();
Print_H(asm_opt);
return 0;
}
return 1;
return check_option(asm_opt);
}
CommandLines.h
+28 -16
View File
@@ -3,25 +3,37 @@
#include <pthread.h>
#define VERBOSE 0
extern char* read_file_name;
extern char* output_file_name;
extern int thread_num;
extern int k_mer_length;
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;
extern int roundID;
extern int coverage;
extern int read_graph;
extern int c_round;
extern char* required_read_name;
extern int adapterLen;
typedef struct {
char* read_file_name;
char* output_file_name;
char* required_read_name;
int thread_num;
int k_mer_length;
int k_mer_min_freq;
int k_mer_max_freq;
int load_index_from_disk;
int write_index_to_disk;
int number_of_round;
int adapterLen;
int clean_round;
int complete_threads;
int roundID;
float min_drop_rate;
float max_drop_rate;
long long pop_bubble_size;
long long num_bases;
long long num_corrected_bases;
long long num_recorrected_bases;
} hifiasm_opt_t;
int CommandLine_process (int argc, char *argv[]);
extern hifiasm_opt_t asm_opt;
void init_opt(hifiasm_opt_t* asm_opt);
void clear_opt(hifiasm_opt_t* asm_opt, int last_round);
int CommandLine_process (int argc, char *argv[], hifiasm_opt_t* asm_opt);
double Get_T(void);
#endif
Correct.cpp
+431 -6841
View File
File diff suppressed because it is too large Load Diff
Correct.h
+13 -309
View File
@@ -16,12 +16,11 @@
#define INSERTION 2
#define DELETION 3
#define MIN(x,y) ((x)<=(y)?(x):(y))
///#define FLAG_THRE 0
#define MAX(x, y) ((x >= y)?x:y)
#define MIN(x, y) ((x <= y)?x:y)
#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
@@ -37,8 +36,6 @@
#define Adjust_Threshold(threshold, x_len) ((threshold == 0 && x_len >= 4)? 1: threshold)
typedef struct
{
long long read_length;
@@ -67,7 +64,6 @@ inline void init_Window_Pool(Window_Pool* dumy, long long read_length, long long
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;
@@ -92,28 +88,10 @@ inline int get_Window(Window_Pool* dumy, long long* w_beg, long long* w_end)
{
dumy->window_end = dumy->read_length - 1;
}
// else if (dumy->read_length - dumy->window_end - 1 <= dumy->tail_length)
// {
// dumy->window_end = dumy->read_length - 1;
// }
/**
if((*w_end) - (*w_beg) + 1 < 375 && (*w_end) + 1 != dumy->read_length)
{
fprintf(stderr, "(*w_beg):%d, (*w_end): %d, dumy->read_length: %d\n",
(*w_beg), (*w_end), dumy->read_length);
}
**/
return 1;
}
typedef struct
{
/**[0-1] bits are type:**/
@@ -326,45 +304,6 @@ inline int filter_one_snp(int occ_0, int occ_1, int total)
return 1;
}
inline int filter_one_snp_advance_back(int occ_0, int occ_1, int total, int group_size,
long long homopolymer_num, long long non_homopolymer_num)
{
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 && group_size > 1)
if(non_homopolymer_num > 0 && 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 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;
@@ -502,96 +441,6 @@ uint32_t* SNPs, int SNPsLen)
return 1;
}
/**
inline int if_is_homopolymer(long long site, char* read, long long read_length)
{
long long beg, end, i;
beg = site - 10;
if(beg < 0)
{
beg = 0;
}
end = site + 10;
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 >= 5 || b_homopolymer_len >= 5)
{
return 1;
}
if (b_homopolymer_ch == f_homopolymer_ch
&&
(f_homopolymer_len + b_homopolymer_len >= 5))
{
return 1;
}
return 0;
}
**/
inline int if_is_homopolymer_strict(long long site, char* read, long long read_length)
{
@@ -666,22 +515,6 @@ inline int if_is_homopolymer_strict(long long site, char* read, long long read_l
b_homopolymer_len++;
}
/**
fprintf(stderr, "site: %d, beg: %d, end: %d\n", site, beg, end);
for (i = beg; i <= end; i++)
{
if (i == site)
{
fprintf(stderr, "|%c|", read[i]);
}
else
{
fprintf(stderr, "%c", read[i]);
}
}
fprintf(stderr, "\n");
**/
if(f_homopolymer_len >= threshold || b_homopolymer_len >= threshold)
{
return 1;
@@ -774,22 +607,6 @@ inline int if_is_homopolymer_repeat(long long site, char* read, long long read_l
b_homopolymer_len++;
}
/**
fprintf(stderr, "site: %d, beg: %d, end: %d\n", site, beg, end);
for (i = beg; i <= end; i++)
{
if (i == site)
{
fprintf(stderr, "|%c|", read[i]);
}
else
{
fprintf(stderr, "%c", read[i]);
}
}
fprintf(stderr, "\n");
**/
if(f_homopolymer_len >= threshold || b_homopolymer_len >= threshold)
{
return 1;
@@ -825,8 +642,6 @@ UC_Read* g_read)
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);
///fprintf(stderr, "is_homopolymer: %d\n", h->snp_stat[h->available_snp].is_homopolymer);
int8_t* vector = Get_SNP_Vector((*h), h->available_snp);
for (i = 0; i < sub_length; i++)
{
@@ -1053,7 +868,7 @@ inline void insert_SNP_IDs_addition(Snp_ID_Vector_Alloc* SNP_IDs, uint32_t* IDs_
}
inline void init_DP_matrix(DP_matrix* dp, int32_t snp_num)
inline void init_DP_matrix(DP_matrix* dp, uint32_t snp_num)
{
if(snp_num > dp->snp_size)
@@ -1147,9 +962,6 @@ inline void destoryHaplotypeEvdience(haplotype_evdience_alloc* h)
free(h->dp.buffer);
free(h->dp.max_buffer);
destory_SNP_IDs(&(h->dp.SNP_IDs));
}
inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h)
@@ -1237,115 +1049,22 @@ void clear_Round2_alignment(Round2_alignment* h);
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,
long long* matched_overlap_0, long long* matched_overlap_1,
long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1,
Cigar_record* current_cigar, haplotype_evdience_alloc* hap,
Round2_alignment* second_round, int force_repeat, int is_consensus,
int* fully_cov, int* abnormal, uint8_t* c2n);
int* fully_cov, int* abnormal);
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 clear_Correct_dumy_pure(Correct_dumy* list);
void pre_filter_by_nearby(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read,
All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length);
void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read,
All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length);
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 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];
}
}
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);
/***
@@ -1356,15 +1075,11 @@ inline void add_existing_cell_to_cigar_record_with_different_base(Cigar_record*
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)
if((uint32_t)dummy->current_operation == type)
{
add_existing_cell_to_cigar_record(dummy, len, type);
}
@@ -1379,7 +1094,7 @@ inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint3
{
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
if((uint32_t)dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
@@ -1394,7 +1109,7 @@ inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint3
{
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
if((uint32_t)dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
@@ -1418,7 +1133,7 @@ inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint3
**/
///add to existing cell, just increase length
///and add different bases
if(dummy->current_operation == type)
if((uint32_t)dummy->current_operation == type)
{
add_existing_cell_to_cigar_record_with_different_base(dummy, len, type, seq);
}
@@ -1431,21 +1146,10 @@ inline void add_cigar_record(char* seq, uint32_t len, Cigar_record* dummy, uint3
}
dummy->current_operation = type;
}
int verify_cigar_2(char* x, int x_len, char* y, int y_len, Cigar_record* cigar, int error);
/**********************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,
uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy);
/**********************for prefilter************************ */
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,
Hash_Table.cpp
+104 -2547
View File
File diff suppressed because it is too large Load Diff
Hash_Table.h
+17 -211
View File
@@ -31,7 +31,7 @@ typedef khash_t(POS64) Pos_Table;
#define FINAL_OVERLAP_ERROR_RATE 0.03
#define GROUP_SIZE 4
///最长是10M10D10M10D10M这种
///the max cigar likes 10M10D10M10D10M
///#define CIGAR_MAX_LENGTH THRESHOLD*2+2
#define CIGAR_MAX_LENGTH 31*2+4
@@ -93,8 +93,6 @@ typedef struct
int extra_begin;
int extra_end;
int error_threshold;
///int y_pre_start;
///error小于等于0都要重新算
int error;
CIGAR cigar;
} window_list;
@@ -218,20 +216,7 @@ typedef struct
} Total_Pos_Table;
/********************************for debug***************************************/
inline void print_64bit(uint64_t x)
{
int i;
for(i = 63; i >= 0; i--)
{
if(x & ((1ULL<<i)))
fprintf(stderr, "1");
else
fprintf(stderr, "0");
}
fprintf(stderr, "\n");
}
inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d)
{
@@ -243,22 +228,7 @@ inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d)
return result;
}
inline int if_k_mer_available(Hash_code* code, int k)
{
uint64_t h_key, low_key;
///k有可能是64，所以可能会有问题
///low_key = code->x[0] | (code->x[1] << k);
low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k));
//k不可能为0, 所以这个右移不会有问题
h_key = code->x[1] >> (64 - k);
if(mod_d(h_key, low_key, MODE_VALUE) > 3)
{
return 0;
}
return 1;
}
////suffix_bits = 64 in default
inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code,
@@ -278,6 +248,8 @@ uint64_t suffix_mode, int suffix_bits, int k)
code->x[1] = h_key << (64 - k);
code->x[1] = code->x[1] | (low_key >> SAFE_SHIFT(k));
return 1;
}
///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k)
@@ -285,10 +257,10 @@ inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t s
Hash_code* code, int k)
{
uint64_t h_key, low_key;
///k有可能是64，所以可能会有问题
///k might be 64，so it is unsafe
///low_key = code->x[0] | (code->x[1] << k);
low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k));
//k不可能为0, 所以这个右移不会有问题
//k cannot be 0, so this shift is safe
h_key = code->x[1] >> (64 - k);
if(mod_d(h_key, low_key, MODE_VALUE) > 3)
@@ -296,25 +268,13 @@ Hash_code* code, int k)
return 0;
}
///注意suffix_bits最大就是64
///前一个右移不安全，因为TCB->suffix_bits有可能为64
///后一个左移安全，因为TCB->suffix_bits不可能为0
//uint64_t sub_ID = (low_key >> TCB->suffix_bits) | (h_key << (64 - TCB->suffix_bits));
uint64_t sub_ID = (low_key >> SAFE_SHIFT(suffix_bits)) | (h_key << (64 - suffix_bits));
uint64_t sub_key = (low_key & suffix_mode);
*get_sub_ID = sub_ID;
*get_sub_key = sub_key;
// Hash_code de_code;
// recover_hash_code(sub_ID, sub_key, &de_code, suffix_mode, suffix_bits, k);
///if(de_code.x[0] != (*code).x[0] || de_code.x[1] != (*code).x[1]) fprintf(stderr, "hehe\n");
///if(de_code.x[0] == (*code).x[0] || de_code.x[1] == (*code).x[1]) fprintf(stderr, "hehe\n");
return 1;
}
@@ -326,7 +286,7 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int
return 0;
}
khint_t t; ///这就是个迭代器
khint_t t;
int absent;
@@ -340,7 +300,7 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int
{
kh_value(TCB->sub_h[sub_ID], t) = 1;
}
else ///哈希表中已有的元素
else
{
//kh_value(TCB->sub_h[sub_ID], t) = kh_value(TCB->sub_h[sub_ID], t) + 1;
kh_value(TCB->sub_h[sub_ID], t)++;
@@ -360,10 +320,9 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
return 0;
}
khint_t t; ///这就是个迭代器
int absent;
khint_t t;
///查询哈希表，key为k
///query hash table，key is k
t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key);
if (t != kh_end(TCB->sub_h[sub_ID]))
@@ -389,10 +348,9 @@ inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k
return (uint64_t)-1;
}
khint_t t; ///这就是个迭代器
int absent;
khint_t t;
///查询哈希表，key为k
///query hash table，key is k
t = kh_get(POS64, PCB->sub_h[sub_ID], sub_key);
if (t != kh_end(PCB->sub_h[sub_ID]))
@@ -441,7 +399,6 @@ inline uint64_t locate_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, k_
int cmp_k_mer_pos(const void * a, const void * b);
//inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos, uint64_t direction)
inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos)
{
k_mer_pos* list;
@@ -474,8 +431,7 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in
__sync_lock_release(&PCB->sub_h_lock[sub_ID].lock);
///当所有位置都存好后，不会再有其他线程修改该list
///所以可以在临界区外排序
//if all pos has been saved, it is safe to sort
if (flag && occ>1)
{
qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos);
@@ -510,7 +466,6 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k
void init_Candidates_list(Candidates_list* l);
void clear_Candidates_list(Candidates_list* l);
void destory_Candidates_list(Candidates_list* l);
void merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand);
void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
@@ -521,9 +476,7 @@ uint64_t n_end_pos, uint8_t n_direction);
void merge_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, Candidates_list* candidates);
void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT);
void merge_k_mer_pos_list_alloc_heap_sort_advance(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT);
void Init_Heap(HeapSq* HBT);
void destory_Heap(HeapSq* HBT);
@@ -532,48 +485,24 @@ void clear_Heap(HeapSq* HBT);
void init_overlap_region_alloc(overlap_region_alloc* list);
void clear_overlap_region_alloc(overlap_region_alloc* list);
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,
int extra_begin, int extra_end, int error_threshold);
void insert_kv_list_to_candidates(k_v* list, long long occ, long long y_id, long long y_offset, long long y_strand,
Candidates_list* candidates);
void overlap_region_sort_y_id(overlap_region *a, long long n);
void calculate_inexact_overlap_region(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF);
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
static const char LogTable256[256] = {
#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
-1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
};
static inline int ilog2_32(uint32_t v)
{
uint32_t t, tt;
if ((tt = v>>16)) return (t = tt>>8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v];
}
void init_fake_cigar(Fake_Cigar* x);
void destory_fake_cigar(Fake_Cigar* x);
void clear_fake_cigar(Fake_Cigar* x);
void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift);
void resize_fake_cigar(Fake_Cigar* x, long long size);
void resize_fake_cigar(Fake_Cigar* x, uint64_t size);
int get_fake_gap_pos(Fake_Cigar* x, int index);
int get_fake_gap_shift(Fake_Cigar* x, int index);
inline long long y_start_offset(long long x_start, Fake_Cigar* o)
@@ -585,7 +514,7 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
long long i;
for (i = 0; i < o->length; i++)
for (i = 0; i < (long long)o->length; i++)
{
if(x_start < get_fake_gap_pos(o, i))
{
@@ -593,7 +522,7 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
}
}
if(i == 0 || i == o->length)
if(i == 0 || i == (long long)o->length)
{
fprintf(stderr, "ERROR\n");
exit(0);
@@ -606,139 +535,16 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
inline void print_fake_gap(Fake_Cigar* o)
{
long long i;
for (i = 0; i < o->length; i++)
for (i = 0; i < (long long)o->length; i++)
{
fprintf(stderr, "**i: %d, gap_pos_in_x: %d, gap_shift: %d\n",
fprintf(stderr, "**i: %lld, gap_pos_in_x: %d, gap_shift: %d\n",
i, get_fake_gap_pos(o, i),
get_fake_gap_shift(o, i));
}
}
/********************************for debug***************************************/
inline int verify_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{
uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k))
{
return 0;
}
khint_t t; ///这就是个迭代器
int absent;
///查询哈希表，key为k
t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key);
if (t != kh_end(TCB->sub_h[sub_ID]))
{
kh_value(TCB->sub_h[sub_ID], t)--;
if (kh_value(TCB->sub_h[sub_ID], t)<0)
{
return -1;
}
else
{
return 1;
}
}
else
{
return -1;
}
}
/********************************for debug***************************************/
inline int Traverse_Total_Count_Table(Total_Count_Table* TCB)
{
int i;
Count_Table* h;
khint_t k;
long long non_empty_k_mer = 0;
for (i = 0; i < TCB->size; i++)
{
h = TCB->sub_h[i];
for (k = kh_begin(h); k != kh_end(h); ++k)
{
if (kh_exist(h, k)) // test if a bucket contains data
{
non_empty_k_mer++;
if (kh_value(h, k)!= 0)
{
fprintf(stderr, "ERROR when Traversing!\n");
}
}
}
}
fprintf(stdout, "non_empty_k_mer: %lld\n", non_empty_k_mer);
}
/********************************for debug***************************************/
void test_COUNT64();
/********************************for debug***************************************/
void debug_mode(uint64_t d, uint64_t thread_ID, uint64_t thread_num);
/********************************for debug***************************************/
void merge_Candidates_list_version(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand);
void sort_candidates(Candidates_list* candidates, long long readID,
overlap_region_alloc* overlap_list, All_reads* R_INF);
void append_overlap_region_alloc_from_existing(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF);
int cmp_by_x_pos_s(const void * a, const void * b);
void resize_Chain_Data(Chain_Data* x, long long size);
void init_window_list_alloc(window_list_alloc* x);
void clear_window_list_alloc(window_list_alloc* x);
void destory_window_list_alloc(window_list_alloc* x);
Levenshtein_distance.cpp
-805
View File
@@ -1,806 +1 @@
#include "Levenshtein_distance.h"
void output_bit_myers(Word x, int length)
{
int i = 0;
while (i < length)
{
fprintf(stderr, "%u", (x >> i) & ((Word)1));
i++;
}
fprintf(stderr, "\n");
}
void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i)
{
i++;
int k = 0;
int j = i;
int diff;
while (k < length)
{
int x_p = (VP >> k) & ((Word)1);
int x_n = (VN >> k) & ((Word)1);
if (x_p == 1 && x_n == 1)
{
fprintf(stderr, "error\n");
}
if (x_p == 1)
{
diff = 1;
///fprintf(stderr, "[+1]");
}
if (x_n == 1)
{
diff = -1;
///fprintf(stderr, "[-1]");
}
if (x_p == 0 && x_n == 0)
{
diff = 0;
///fprintf(stderr, "[+0]");
}
j++;
if (matrix[i][j] - matrix[i][j - 1] != diff)
{
fprintf(stderr, "*************V(k): %u\n", k);
///return;
}
k++;
}
///fprintf(stderr, "\n");
}
void prase_D0(Word D0, int length, int matrix[1000][1000], int i)
{
i++;
int k = 0;
int j = i;
int diff;
while (k < length)
{
int diff = (D0 >> k) & ((Word)1);
if(diff==matrix[i][j] - matrix[i-1][j-1])
{
fprintf(stderr, "*************D(k): %u\n", k);
fprintf(stderr, "diff: %u, matrix[i][j]: %u, matrix[i-1][j-1]: %u\n", diff, matrix[i][j], matrix[i-1][j-1]);
///return;
}
j++;
k++;
}
///fprintf(stderr, "\n");
}
void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i)
{
i++;
int k = 0;
int j = i;
int diff;
while (k < length)
{
int x_p = (HP >> k) & ((Word)1);
int x_n = (HN >> k) & ((Word)1);
if (x_p == 1 && x_n == 1)
{
fprintf(stderr, "error\n");
}
if (x_p == 1)
{
diff = 1;
///fprintf(stderr, "[+1]");
}
if (x_n == 1)
{
diff = -1;
///fprintf(stderr, "[-1]");
}
if (x_p == 0 && x_n == 0)
{
diff = 0;
///fprintf(stderr, "[+0]");
}
if(diff!=matrix[i][j] - matrix[i-1][j])
{
fprintf(stderr, "*************H(k): %u\n", k);
///return;
}
j++;
k++;
}
///fprintf(stderr, "\n");
}
/**
pattern是长的那个，是y
p_length是长的那个的长度, p_length实际没用
text是短的那个，是x
t_length是短的那个的长度
errthold是阈值
return_err是编辑距离
返回值是结束位置
**/
int Reserve_Banded_BPM_debug
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err, int matrix[1000][1000])
{
(*return_err) = (unsigned int)-1;
Word Peq[256];
int band_length = (errthold << 1) + 1;
int i = 0;
Word tmp_Peq_1 = (Word)1;
Peq['A'] = (Word)0;
Peq['T'] = (Word)0;
Peq['G'] = (Word)0;
Peq['C'] = (Word)0;
Word Peq_A;
Word Peq_T;
Word Peq_C;
Word Peq_G;
///band_length = 2k + 1
for (i = 0; i<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
///Peq['T'] = Peq['T'] | Peq['C'];
Peq_A = Peq['A'];
Peq_C = Peq['C'];
Peq_T = Peq['T'];
Peq_G = Peq['G'];
memset(Peq, 0, sizeof(Word)* 256);
Peq['A'] = Peq_A;
Peq['C'] = Peq_C;
Peq['T'] = Peq_T;
Peq['G'] = Peq_G;
Word Mask = ((Word)1 << (errthold << 1));
Word VP = 0;
Word VN = 0;
Word X = 0;
Word D0 = 0;
Word HN = 0;
Word HP = 0;
i = 0;
int err = 0;
Word err_mask = (Word)1;
///band_down = 2k
///i_bd = 2k
///int i_bd = i + band_down;
int i_bd = (errthold << 1);
int last_high = (errthold << 1);
/// t_length_1 = SEQ_LENGTH - 1
int t_length_1 = t_length - 1;
//while(i<t_length)
while (i<t_length_1)
{
fprintf(stderr, "i: %u, j_s: %u, j_e: %u\n", i + 1, i + 1, i + band_length + 1);
///if (i >= 6)
if (i >= 0)
{
/**
fprintf(stderr, "VP:\n");
output_bit_myers(VP, band_length);
fprintf(stderr, "VN:\n");
output_bit_myers(VN, band_length);
fprintf(stderr, "HP:\n");
output_bit_myers(HP, band_length);
fprintf(stderr, "HN:\n");
output_bit_myers(HN, band_length);
fprintf(stderr, "D0:\n");
output_bit_myers(D0, band_length);
fprintf(stderr, "text[i]: %c\n", text[i]);
fprintf(stderr, "Peq[text[i]]:\n");
output_bit_myers(Peq[text[i]], band_length);
for (size_t j = i; j < i + band_length; j++)
{
fprintf(stderr, "%c", pattern[j]);
}
fprintf(stderr, "\n");
fprintf(stderr, "Previous begin.\n");
prase_vertical(VP, VN, band_length, matrix, i-1);
prase_D0(D0, band_length, matrix, i-1);
prase_H(HP, HN, band_length, matrix, i-1);
fprintf(stderr, "Previous test done.\n");
**/
X = Peq[text[i]] | VN;
/**
fprintf(stderr, "#X:\n");
output_bit_myers(X, band_length);
**/
D0 = ((VP + (X&VP)) ^ VP) | X;
/**
fprintf(stderr, "#(X&VP):\n");
output_bit_myers((X&VP), band_length);
fprintf(stderr, "#(VP + (X&VP)):\n");
output_bit_myers((VP + (X&VP)), band_length);
fprintf(stderr, "#((VP + (X&VP)) ^ VP):\n");
output_bit_myers(((VP + (X&VP)) ^ VP), band_length);
fprintf(stderr, "#D0:\n");
output_bit_myers(D0, band_length);
**/
HN = VP&D0;
HP = VN | ~(VP | D0);
X = D0 >> 1;
VN = X&HP;
VP = HN | ~(X | HP);
}
else
{
///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
X = Peq[text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
X = D0 >> 1;
VN = X&HP;
VP = HN | ~(X | HP);
}
/**
for (size_t j = i + 1; j <= i + band_length + 1; j++)
{
fprintf(stderr, "[%u]", matrix[i + 1][j]);
}
fprintf(stderr, "\n");
**/
prase_vertical(VP, VN, band_length, matrix, i);
prase_D0(D0, band_length, matrix, i);
prase_H(HP, HN, band_length, matrix, i);
/**
fprintf(stderr, "VP:\n");
output_bit_myers(VP, band_length);
fprintf(stderr, "VN:\n");
output_bit_myers(VN, band_length);
**/
if (!(D0&err_mask))
{
++err;
///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
if ((err - last_high)>errthold)
{
///fprintf(stderr, "0 ######, i: %u\n", i);
return -1;
}
}
Peq['A'] = Peq['A'] >> 1;
Peq['C'] = Peq['C'] >> 1;
Peq['G'] = Peq['G'] >> 1;
Peq['T'] = Peq['T'] >> 1;
++i;
++i_bd;
Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
///Peq['T'] = Peq['T'] | Peq['C'];
}
X = Peq[text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
X = D0 >> 1;
VN = X&HP;
VP = HN | ~(X | HP);
if (!(D0&err_mask))
{
++err;
if ((err - last_high)>errthold)
return -1;
}
prase_vertical(VP, VN, band_length, matrix, i);
prase_D0(D0, band_length, matrix, i);
prase_H(HP, HN, band_length, matrix, i);
fprintf(stderr, "err: %d, matrix[][]: %d\n", err, matrix[i+1][i+1]);
fprintf(stderr, "VP:\n");
output_bit_myers(VP, band_length);
fprintf(stderr, "VN:\n");
output_bit_myers(VN, band_length);
////fprintf(stderr, "sucess(2)\n");
/// last_high = 2k
/// site = (SEQ_LENGTH + 2k) - 2k -1
/// site = SEQ_LENGTH - 1
///int site = p_length - last_high - 1;
int site = t_length - 1;
int return_site = -1;
if ((err <= errthold) && (err<=*return_err))
{
*return_err = err;
return_site = site;
}
int i_last = i;
i = 0;
while (i<errthold)
{
err = err + ((VP >> i)&(Word)1);
err = err - ((VN >> i)&(Word)1);
++i;
fprintf(stderr, "*i: %u, err: %d\n", i, err);
if ((err <= errthold) && (err <= *return_err))
{
*return_err = err;
return_site = site + i;
}
}
unsigned int ungap_err;
ungap_err = err;
while (i<last_high)
{
err = err + ((VP >> i)&(Word)1);
err = err - ((VN >> i)&(Word)1);
++i;
fprintf(stderr, "*i: %u, err: %d\n", i, err);
if ((err <= errthold) && (err<=*return_err))
{
*return_err = err;
return_site = site + i;
}
}
if ((ungap_err <= errthold) && (ungap_err == *return_err))
{
return_site = site + errthold;
}
return return_site;
}
int BS_Reserve_Banded_BPM
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err)
{
(*return_err) = (unsigned int)-1;
///Õâ¸öÊÇÄǸöÐèÒªÔ¤´¦ÀíµÄÏòÁ¿
Word Peq[256];
int band_length = (errthold << 1) + 1;
int i = 0;
Word tmp_Peq_1 = (Word)1;
Peq['A'] = (Word)0;
Peq['T'] = (Word)0;
Peq['G'] = (Word)0;
Peq['C'] = (Word)0;
Word Peq_A;
Word Peq_T;
Word Peq_C;
Word Peq_G;
///band_length = 2k + 1
///ÕâÊÇ°ÑpatternµÄÇ°2k + 1¸ö×Ö·ûÔ¤´¦Àí
///pattern[0]¶ÔÓ¦Peq[0]
///pattern[2k]¶ÔÓ¦Peq[2k]
for (i = 0; i<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
Peq['T'] = Peq['T'] | Peq['C'];
Peq_A = Peq['A'];
Peq_C = Peq['C'];
Peq_T = Peq['T'];
Peq_G = Peq['G'];
memset(Peq, 0, sizeof(Word)* 256);
Peq['A'] = Peq_A;
Peq['C'] = Peq_C;
Peq['T'] = Peq_T;
Peq['G'] = Peq_G;
Word Mask = ((Word)1 << (errthold << 1));
Word VP = 0;
Word VN = 0;
Word X = 0;
Word D0 = 0;
Word HN = 0;
Word HP = 0;
i = 0;
int err = 0;
Word err_mask = (Word)1;
///band_down = 2k
///i_bd = 2k
///int i_bd = i + band_down;
int i_bd = (errthold << 1);
int last_high = (errthold << 1);
/// t_length_1 = SEQ_LENGTH - 1
int t_length_1 = t_length - 1;
//while(i<t_length)
while (i<t_length_1)
{
///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
X = Peq[text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
X = D0 >> 1;
VN = X&HP;
VP = HN | ~(X | HP);
///Èç¹ûб¶Ô½ÇÏß·½ÏòÆ¥ÅäÔòD0ÊÇ1
///Èç¹û²»Æ¥ÅäÔòD0ÊÇ0
///Õâ¸öÒâ˼ÊÇÈç¹û×îÉÏÃæÄÇÌõ¶Ô½ÇÏßÉϵÄб¶Ô½ÇÏß·½Ïò·¢ÉúÎóÅä,ÔòÖ´ÐÐÄÚ²¿³ÌÐò
///
if (!(D0&err_mask))
{
++err;
///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
if ((err - last_high)>errthold)
return -1;
}
///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
//ÓÒÒÆʵ¼ÊÉÏÊÇ°Ñpattern[0]ÒƵôÁË
Peq['A'] = Peq['A'] >> 1;
Peq['C'] = Peq['C'] >> 1;
Peq['G'] = Peq['G'] >> 1;
Peq['T'] = Peq['T'] >> 1;
++i;
++i_bd;
///ÕâÊÇ°ÑеÄpattern[2k]¼Ó½øÀ´, ÕâòËÆÊǼӵ½Peq[2k]ÉÏÁË
Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
Peq['T'] = Peq['T'] | Peq['C'];
}
///fprintf(stderr, "sucess(1)\n");
///Õâ¸öÑ­»·ÄóöÀ´ÊÇΪÁË·ÀÖ¹ÄÚ´æй¶
///ÆäʵҲ¾ÍÊÇÑ­»·ÀïµÄ×îºóÒ»ÐÐÓï¾ä°É
///ÍêÈ«¿ÉÒÔ°ÑpatternÔö´óһλ
///²»¹ýÕâÑùÒ²ºÃ£¬¿ÉÒÔ¼õÉÙ¼ÆË㿪Ïú
X = Peq[text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
X = D0 >> 1;
VN = X&HP;
VP = HN | ~(X | HP);
if (!(D0&err_mask))
{
++err;
if ((err - last_high)>errthold)
return -1;
}
////fprintf(stderr, "sucess(2)\n");
/// last_high = 2k
/// site = (SEQ_LENGTH + 2k) - 2k -1
/// site = SEQ_LENGTH - 1
///´ËʱÕâ¸ösiteòËÆÊÇ×îÉÏÃæÄÇÌõ¶Ô½ÇÏßµÄλÖÃ
///int site = p_length - last_high - 1;
int site = t_length - 1;
int return_site = -1;
if ((err <= errthold) && (err<=*return_err))
{
*return_err = err;
return_site = site;
}
int i_last = i;
i = 0;
while (i<errthold)
{
err = err + ((VP >> i)&(Word)1);
err = err - ((VN >> i)&(Word)1);
++i;
if ((err <= errthold) && (err <= *return_err))
{
*return_err = err;
return_site = site + i;
}
}
unsigned int ungap_err;
ungap_err = err;
while (i<last_high)
{
err = err + ((VP >> i)&(Word)1);
err = err - ((VN >> i)&(Word)1);
++i;
if ((err <= errthold) && (err<=*return_err))
{
*return_err = err;
return_site = site + i;
}
}
if ((ungap_err <= errthold) && (ungap_err == *return_err))
{
return_site = site + errthold;
}
return return_site;
}
int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length,
unsigned short errthold,unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id)
{
Word Peq[128];
char Peq_index[4]= {'A','C','G','T'};
int symbol = 0;
int r;
Word tmp_Peq_1=(Word)1;
Peq['A']=(Word)0;
Peq['T']=(Word)0;
Peq['G']=(Word)0;
Peq['C']=(Word)0;
Word Peq_A;
Word Peq_T;
Word Peq_C;
Word Peq_G;
for (r =0; r<band_length; r++)
{
Peq[pattern[r]]=Peq[pattern[r]]|tmp_Peq_1;
tmp_Peq_1=tmp_Peq_1<<1;
}
Peq_A=Peq['A'];
Peq_C=Peq['C'];
Peq_T=Peq['T'];
Peq_G=Peq['G'];
for(symbol = 0; symbol < 128; symbol++)
{
Peq[symbol]=(Word)0;
//jump_c[symbol]=128;
}
Peq['A']=Peq_A;
Peq['C']=Peq_C;
Peq['T']=Peq_T;
Peq['G']=Peq_G;
Word Mask_Pre=(Word)1<<(band_length-2);
Word Mask=(Word)1<<(band_length-1);
Word VP=0;
Word VN=0;
Word X=0;
Word D0=0;
Word HN=0;
Word HP=0;
int s=0;
int i = 0;
int j=0;
int bound=band_length-2-band_down;
int err=0;
Word err_mask=(Word)1;
int s1=band_length-2;
int i_bd=i+band_down;
int last_high=band_length-t_length+p_length-band_down-1;
int t_length_1=t_length-1;
//while(i<t_length)
while(i<t_length_1)
{
X=Peq[text[i]]|VN;
D0=((VP+(X&VP))^VP)|X;
HN=VP&D0;
HP=VN|~(VP|D0);
X=D0>>1;
VN=X&HP;
VP=HN|~(X|HP);
if(!(D0&err_mask))
{
++err;
if((err-last_high)>errthold)
return -1;
}
Peq['A']=Peq['A']>>1;
Peq['C']=Peq['C']>>1;
Peq['G']=Peq['G']>>1;
Peq['T']=Peq['T']>>1;
++i;
++i_bd;
Peq[pattern[i_bd]]=Peq[pattern[i_bd]]|Mask;
}
///这个循环拿出来是为了防止内存泄露
X=Peq[text[i]]|VN;
D0=((VP+(X&VP))^VP)|X;
HN=VP&D0;
HP=VN|~(VP|D0);
X=D0>>1;
VN=X&HP;
VP=HN|~(X|HP);
if(!(D0&err_mask))
{
++err;
if((err-last_high)>errthold)
return -1;
}
int site=p_length-last_high-1;
int return_site=-1;
if((err<=errthold)&&(err<*return_err))
{
*return_err=err;
return_site=site;
}
int i_last=i;
i=0;
while(i<last_high)
{
err=err+((VP>>i)&(Word)1);
err=err-((VN>>i)&(Word)1);
++i;
if((err<=errthold)&&(err<*return_err))
{
*return_err=err;
return_site=site+i;
}
}
return return_site;
}
Levenshtein_distance.h
+58 -234
View File
@@ -12,152 +12,18 @@
typedef uint64_t Word;
typedef uint32_t Word_32;
typedef struct
{
/**
* vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0]
* num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0
**/
Word* vec;
Word num_words;
Word num_bits;
///for vec[num_words-1]
Word last_bit_mode;
Word last_bit_shift;
} bit_vectors;
inline void init_bit_vector(bit_vectors* b, Word num_bits)
{
b->num_bits = num_bits;
b->num_words = ((b->num_bits)>>6);
b->last_bit_mode = (Word)-1;
b->last_bit_shift = b->num_bits&63;
if(b->last_bit_shift != 0)
{
b->num_words++;
b->last_bit_mode = b->last_bit_mode >> (64 - b->last_bit_shift);
}
b->vec = (Word*)calloc(b->num_words, sizeof(Word));
}
inline void destory_bit_vector(bit_vectors* b, Word num_bits)
{
free(b->vec);
}
///&
inline int AND_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] & y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///|
inline int OR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] | y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///^
inline int XOR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] ^ y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///~
inline int NOT_bit_vector(bit_vectors* source, bit_vectors* dest)
{
if(source->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < source->num_words; i++)
{
dest->vec[i] = ~(source->vec[i]);
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///<< 1
inline int L_shift_1_bit_vector(bit_vectors* source, bit_vectors* dest)
{
/**
* vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0]
* num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0
**/
if(source->num_bits != dest->num_bits || source->num_words < 1) return 0;
Word i;
for (i = source->num_words - 1; i >= 1; i--)
{
dest->vec[i] = (source->vec[i])<<1;
dest->vec[i] = dest->vec[i] | ((source->vec[i-1])>>63);
}
dest->vec[0] = (source->vec[0])<<1;
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
return 1;
}
inline void get_error(int t_length, int errthold, int init_err, Word VP, Word VN,
unsigned int* return_err, int* back_site)
{
(*return_err) = (unsigned int)-1;
int site = t_length - 1;
int return_site = -1;
///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
///也就是p_length - t_length
///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
///in most cases, p_length should be t_length + 2 * errthold
///int available_i = p_length - t_length;
int available_i = 2 * errthold;
if ((init_err <= errthold) && (init_err <= (*return_err)))
if ((init_err <= errthold) && ((unsigned int)init_err <= (*return_err)))
{
(*return_err) = init_err;
return_site = site;
@@ -173,7 +39,7 @@ unsigned int* return_err, int* back_site)
init_err = init_err - ((VN >> i)&(Word)1);
++i;
if ((init_err <= errthold) && (init_err <= *return_err))
if ((init_err <= errthold) && ((unsigned int)init_err <= *return_err))
{
*return_err = init_err;
return_site = site + i;
@@ -188,7 +54,7 @@ unsigned int* return_err, int* back_site)
}
/****************************may have bugs********************************/
if((ungap_error<=errthold) && (ungap_error == (*return_err)))
if((ungap_error<=(unsigned int)errthold) && (ungap_error == (*return_err)))
{
return_site = site + errthold;
}
@@ -227,7 +93,7 @@ unsigned int* return_err, int* return_p_end, int* return_t_end)
///band_length = 2k + 1
for (i = 0; i<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
Peq[(uint8_t)pattern[i]] = Peq[(uint8_t)pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
@@ -273,7 +139,7 @@ unsigned int* return_err, int* return_p_end, int* return_t_end)
while (i<t_length_1)
{
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
@@ -308,14 +174,14 @@ unsigned int* return_err, int* return_p_end, int* return_t_end)
++i;
++i_bd;
Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
Peq[(uint8_t)pattern[i_bd]] = Peq[(uint8_t)pattern[i_bd]] | Mask;
}
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
@@ -405,13 +271,7 @@ int* return_t_end, int* return_aligned_t_len)
/**
pattern是长的那个，是y
p_length是长的那个的长度, p_length实际没用
text是短的那个，是x
t_length是短的那个的长度
errthold是阈值
return_err是编辑距离
返回值是结束位置
pattern is the longer one, while text is the shorter one
**/
inline int Reserve_Banded_BPM
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err)
@@ -438,7 +298,7 @@ inline int Reserve_Banded_BPM
///band_length = 2k + 1
for (i = 0; i<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
Peq[(uint8_t)pattern[i]] = Peq[(uint8_t)pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
@@ -493,8 +353,7 @@ inline int Reserve_Banded_BPM
while (i<t_length_1)
{
///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
@@ -509,10 +368,8 @@ inline int Reserve_Banded_BPM
{
++err;
///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k
if ((err - last_high)>errthold)
if ((err - last_high)>(int)errthold)
{
///fprintf(stderr, "0 ######, i: %u\n", i);
return -1;
}
@@ -527,17 +384,14 @@ inline int Reserve_Banded_BPM
++i;
++i_bd;
Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
///Peq['T'] = Peq['T'] | Peq['C'];
Peq[(uint8_t)pattern[i_bd]] = Peq[(uint8_t)pattern[i_bd]] | Mask;
}
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
@@ -560,16 +414,13 @@ inline int Reserve_Banded_BPM
///int site = p_length - last_high - 1;
int site = t_length - 1;
int return_site = -1;
///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
///也就是p_length - t_length
///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
///in most cases, p_lengthshould be t_length + 2 * errthold
int available_i = p_length - t_length;
if ((err <= errthold) && (err<=*return_err))
if ((err <= errthold) && ((unsigned int)err<=*return_err))
{
*return_err = err;
return_site = site;
}
int i_last = i;
i = 0;
/****************************may have bugs********************************/
@@ -582,14 +433,14 @@ inline int Reserve_Banded_BPM
err = err - ((VN >> i)&(Word)1);
++i;
if ((err <= errthold) && (err <= *return_err))
if ((err <= (int)errthold) && ((unsigned int)err <= *return_err))
{
*return_err = err;
return_site = site + i;
}
/****************************may have bugs********************************/
if(i == errthold)
if(i == (int)errthold)
{
ungap_error = err;
}
@@ -617,7 +468,7 @@ inline int try_cigar(char *pattern, int p_length,
{
int i = 0;
int tmp_err = 0;
///y上的起始位置
///start pos of y
int start_site = end_site - t_length + 1;
if (start_site >= 0)
@@ -626,7 +477,7 @@ inline int try_cigar(char *pattern, int p_length,
for (i = 0; i < t_length; i++)
{
///path[i] = 0;
///path倒着存
///path is saved backwards
path[t_length - i - 1] = 0;
if (text[i] != pattern[i + start_site])
{
@@ -655,8 +506,8 @@ inline int try_cigar(char *pattern, int p_length,
}
///p_length有可能不够，但是t_length总是够的
///就是p_length有可能小于t_length + 2 * errthold
///p_length might be samller than t_length + 2 * errthold
inline int Reserve_Banded_BPM_PATH
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold,
unsigned int* return_err, int* return_start_site, int* return_path_length, Word* matrix_bit, char* path,
@@ -666,7 +517,6 @@ inline int Reserve_Banded_BPM_PATH
{
if (old_error == 0)
{
///fprintf(stderr, "0 error\n");
(*return_err) = old_error;
(*return_start_site) = old_end_site - t_length + 1;
return old_end_site;
@@ -675,7 +525,6 @@ inline int Reserve_Banded_BPM_PATH
if (try_cigar(pattern, p_length, text, t_length, old_end_site, path,
old_error, return_start_site, return_path_length))
{
///fprintf(stderr, "no gap error\n");
(*return_err) = old_error;
return old_end_site;
}
@@ -704,7 +553,7 @@ inline int Reserve_Banded_BPM_PATH
///band_length = 2k + 1
for (i = 0; i<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
Peq[(uint8_t)pattern[i]] = Peq[(uint8_t)pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
@@ -762,7 +611,7 @@ inline int Reserve_Banded_BPM_PATH
while (i<t_length_1)
{
///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0]
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
@@ -777,7 +626,7 @@ inline int Reserve_Banded_BPM_PATH
{
++err;
if ((err - last_high)>errthold)
if ((err - last_high)>(int)errthold)
{
return -1;
}
@@ -793,7 +642,7 @@ inline int Reserve_Banded_BPM_PATH
++i;
++i_bd;
Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask;
Peq[(uint8_t)pattern[i_bd]] = Peq[(uint8_t)pattern[i_bd]] | Mask;
///Peq['T'] = Peq['T'] | Peq['C'];
@@ -810,7 +659,7 @@ inline int Reserve_Banded_BPM_PATH
X = Peq[text[i]] | VN;
X = Peq[(uint8_t)text[i]] | VN;
D0 = ((VP + (X&VP)) ^ VP) | X;
HN = VP&D0;
HP = VN | ~(VP | D0);
@@ -820,7 +669,7 @@ inline int Reserve_Banded_BPM_PATH
if (!(D0&err_mask))
{
++err;
if ((err - last_high)>errthold)
if ((err - last_high)>(int)errthold)
return -1;
}
@@ -846,16 +695,13 @@ inline int Reserve_Banded_BPM_PATH
unsigned int ungap_error = (unsigned int)-1;
/****************************may have bugs********************************/
///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
///也就是p_length - t_length
///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
///in most cases, p_length should be t_length + 2 * errthold
int available_i = p_length - t_length;
if ((err <= errthold) && (err<=*return_err))
if ((err <= (int)errthold) && ((unsigned int)err<=*return_err))
{
*return_err = err;
return_site = site;
}
int i_last = i;
i = 0;
while (i < available_i)
@@ -864,14 +710,14 @@ inline int Reserve_Banded_BPM_PATH
err = err - ((VN >> i)&(Word)1);
++i;
if ((err <= errthold) && (err <= *return_err))
if ((err <= (int)errthold) && ((unsigned int)err <= *return_err))
{
*return_err = err;
return_site = site + i;
}
/****************************may have bugs********************************/
if(i == errthold)
if(i == (int)errthold)
{
ungap_error = err;
}
@@ -894,24 +740,19 @@ inline int Reserve_Banded_BPM_PATH
}
/****************************may have bugs********************************/
////注意，这里p_length要矫正啊啊
///不矫正会出错
///因为p_length有可能不够
///need to correct p_length here, since p_length might be smaller than t_length + 2* err_threashlod
p_length = t_length + 2 * errthold;
///end_site是正确的
///end_site is always correct
int end_site = return_site;
int start_site = end_site;
///这个是各个bit-vector里面，end_site对应bit所在的位置
int back_track_site = band_length - (p_length - end_site);
Word v_value, h_value, delta_value, min_value, current_value;
Word direction, is_mismatch; ///0 is match, 1 is mismatch, 2 is up, 3 is left
///代表pattern到哪了，就是短的那个到哪了
///Word direction; ///0 is match, 1 is mismatch, 2 is up, 3 is left
Word direction = 0; ///0 is match, 1 is mismatch, 2 is up, 3 is left
i = t_length;
int path_length = 0;
///到0就结束了，后面的路径可以直接match
current_value = *return_err;
@@ -1044,15 +885,11 @@ inline int Reserve_Banded_BPM_PATH
(*return_start_site) = start_site;
(*return_path_length) = path_length;
return return_site;
}
////这个p_length四个是一样的
////four patterns have the same p_length
inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *pattern3, char *pattern4, int p_length, char *text, int t_length,
int* return_sites, unsigned int* return_sites_error, unsigned short errthold, __m128i* Peq_SSE)
@@ -1076,10 +913,10 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
for (i = 0; i<band_length; i++)
{
Peq[pattern1[i]][0] = Peq[pattern1[i]][0] | tmp_Peq_1;
Peq[pattern2[i]][1] = Peq[pattern2[i]][1] | tmp_Peq_1;
Peq[pattern3[i]][2] = Peq[pattern3[i]][2] | tmp_Peq_1;
Peq[pattern4[i]][3] = Peq[pattern4[i]][3] | tmp_Peq_1;
Peq[(uint8_t)pattern1[i]][0] = Peq[(uint8_t)pattern1[i]][0] | tmp_Peq_1;
Peq[(uint8_t)pattern2[i]][1] = Peq[(uint8_t)pattern2[i]][1] | tmp_Peq_1;
Peq[(uint8_t)pattern3[i]][2] = Peq[(uint8_t)pattern3[i]][2] | tmp_Peq_1;
Peq[(uint8_t)pattern4[i]][3] = Peq[(uint8_t)pattern4[i]][3] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
@@ -1132,7 +969,7 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
while (i<t_length_1)
{
///X = Peq[text[i]] | VN;
X = _mm_or_si128(Peq_SSE[text[i]], VN);
X = _mm_or_si128(Peq_SSE[(uint8_t)text[i]], VN);
@@ -1188,17 +1025,17 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
++i;
++i_bd;
Peq_SSE[pattern1[i_bd]] = _mm_or_si128(Mask1, Peq_SSE[pattern1[i_bd]]);
Peq_SSE[pattern2[i_bd]] = _mm_or_si128(Mask2, Peq_SSE[pattern2[i_bd]]);
Peq_SSE[pattern3[i_bd]] = _mm_or_si128(Mask3, Peq_SSE[pattern3[i_bd]]);
Peq_SSE[pattern4[i_bd]] = _mm_or_si128(Mask4, Peq_SSE[pattern4[i_bd]]);
Peq_SSE[(uint8_t)pattern1[i_bd]] = _mm_or_si128(Mask1, Peq_SSE[(uint8_t)pattern1[i_bd]]);
Peq_SSE[(uint8_t)pattern2[i_bd]] = _mm_or_si128(Mask2, Peq_SSE[(uint8_t)pattern2[i_bd]]);
Peq_SSE[(uint8_t)pattern3[i_bd]] = _mm_or_si128(Mask3, Peq_SSE[(uint8_t)pattern3[i_bd]]);
Peq_SSE[(uint8_t)pattern4[i_bd]] = _mm_or_si128(Mask4, Peq_SSE[(uint8_t)pattern4[i_bd]]);
///Peq_SSE['T'] = _mm_or_si128(Peq_SSE['T'], Peq_SSE['C']);
}
///X = Peq[text[i]] | VN;
X = _mm_or_si128(Peq_SSE[text[i]], VN);
X = _mm_or_si128(Peq_SSE[(uint8_t)text[i]], VN);
/*************D0 = ((VP + (X&VP)) ^ VP) | X*********************/
///X&VP
@@ -1249,22 +1086,22 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
err4 = _mm_extract_epi32(Err_4, 3);
if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
if ((err1 <= (int)errthold) && ((unsigned int)err1 <= return_sites_error[0]))
{
return_sites[0] = site;
return_sites_error[0] = err1;
}
if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
if ((err2 <= (int)errthold) && ((unsigned int)err2 <= return_sites_error[1]))
{
return_sites[1] = site;
return_sites_error[1] = err2;
}
if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
if ((err3 <= (int)errthold) && ((unsigned int)err3 <= return_sites_error[2]))
{
return_sites[2] = site;
return_sites_error[2] = err3;
}
if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
if ((err4 <= (int)errthold) && ((unsigned int)err4 <= return_sites_error[3]))
{
return_sites[3] = site;
return_sites_error[3] = err4;
@@ -1281,9 +1118,7 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
/****************************may have bugs********************************/
///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
///也就是p_length - t_length
///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
///in most cases, p_length should be t_length + 2 * errthold
int available_i = p_length - t_length;
while (i < available_i)
@@ -1305,29 +1140,29 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
err4 = _mm_extract_epi32(Err_4, 3);
if ((err1 <= errthold) && (err1 <= return_sites_error[0]))
if ((err1 <= (int)errthold) && ((unsigned int)err1 <= return_sites_error[0]))
{
return_sites[0] = site + i;
return_sites_error[0] = err1;
}
if ((err2 <= errthold) && (err2 <= return_sites_error[1]))
if ((err2 <= (int)errthold) && ((unsigned int)err2 <= return_sites_error[1]))
{
return_sites[1] = site + i;
return_sites_error[1] = err2;
}
if ((err3 <= errthold) && (err3 <= return_sites_error[2]))
if ((err3 <= (int)errthold) && ((unsigned int)err3 <= return_sites_error[2]))
{
return_sites[2] = site + i;
return_sites_error[2] = err3;
}
if ((err4 <= errthold) && (err4 <= return_sites_error[3]))
if ((err4 <= (int)errthold) && ((unsigned int)err4 <= return_sites_error[3]))
{
return_sites[3] = site + i;
return_sites_error[3] = err4;
}
/****************************may have bugs********************************/
if(i == errthold)
if(i == (int)errthold)
{
ungap_error1 = err1;
ungap_error2 = err2;
@@ -1363,15 +1198,4 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p
}
void output_bit_myers(Word x, int length);
void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i);
void prase_D0(Word D0, int length, int matrix[1000][1000], int i);
void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i);
int Reserve_Banded_BPM_debug(char *pattern, int p_length, char *text, int t_length, unsigned short errthold,
unsigned int* return_err, int matrix[1000][1000]);
int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length,unsigned short errthold,
unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id);
int BS_Reserve_Banded_BPM
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err);
#endif
Makefile
+6 -7
View File
@@ -1,10 +1,10 @@
CXX= g++
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -Wall
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-Winline
CPPFLAGS=
INCLUDES=
OBJS= Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \
POA.o Correct.o Levenshtein_distance.o edlib.o Overlaps.o ksw2_extz2_sse.o
EXE= ccs_assembly
POA.o Correct.o Levenshtein_distance.o Overlaps.o #ksw2_extz2_sse.o
EXE= hifiasm
LIBS= -lz -lpthread -lm
ifneq ($(asan),)
@@ -37,7 +37,7 @@ Assembly.o: Levenshtein_distance.h Output.h
CommandLines.o: CommandLines.h ketopt.h
Correct.o: Correct.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h
Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h Levenshtein_distance.h
Correct.o: POA.h edlib.h Assembly.h ksw2.h
Correct.o: POA.h Assembly.h #ksw2.h
Hash_Table.o: Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h
Hash_Table.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
Hash_Table.o: POA.h ksort.h
@@ -48,8 +48,7 @@ Overlaps.o: CommandLines.h
POA.o: POA.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h
POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
edlib.o: edlib.h
kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
main.o: Assembly.h Levenshtein_distance.h edlib.h
ksw2_extz2_sse.o: ksw2.h
main.o: Assembly.h Levenshtein_distance.h
#ksw2_extz2_sse.o: ksw2.h
Output.cpp
+5 -4
View File
@@ -150,18 +150,18 @@ inline void push_single_buffer(Output_buffer_sub_block* curr_sub_block)
void* pop_buffer(void*)
{
FILE* output_file = fopen(output_file_name, "w");
FILE* output_file = fopen(asm_opt.output_file_name, "w");
init_buffer_sub_block(&tmp_buffer_sub_block);
while (buffer_out.all_buffer_end < thread_num)
while (buffer_out.all_buffer_end < asm_opt.thread_num)
{
pthread_mutex_lock(&o_queueMutex);
while (if_empty_buffer() && (buffer_out.all_buffer_end < thread_num))
while (if_empty_buffer() && (buffer_out.all_buffer_end < asm_opt.thread_num))
{
pthread_cond_signal(&o_stallCond);
pthread_cond_wait(&o_flushCond, &o_queueMutex);
@@ -196,6 +196,7 @@ void* pop_buffer(void*)
fclose(output_file);
return NULL;
}
@@ -225,7 +226,7 @@ void finish_output_buffer()
buffer_out.all_buffer_end++;
if (buffer_out.all_buffer_end == thread_num)
if (buffer_out.all_buffer_end == asm_opt.thread_num)
{
pthread_cond_signal(&o_flushCond);
}
Overlaps.cpp
+374 -2873
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+11 -34
View File
@@ -3,7 +3,6 @@
#include <stdint.h>
#include "kvec.h"
#include "kdq.h"
///#include "Hash_Table.h"
///#define MIN_OVERLAP_LEN 2000
///#define MIN_OVERLAP_LEN 500
@@ -68,16 +67,7 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n);
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources,
char* output_file_name);
typedef struct {
ma_hit_t_alloc overlaps;
} Assembly_Graph;
void init_Assembly_Graph(Assembly_Graph* x);
void destory_Assembly_Graph(Assembly_Graph* x);
void collect_ma_hit_t(ma_hit_t_alloc* dest, ma_hit_t_alloc* sources, long long num_sources);
void normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources);
void debug_normalize_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources);
typedef struct {
@@ -93,7 +83,6 @@ void ma_hit_flt(ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *cover
int max_hang, int min_ovlp);
long long get_specific_overlap(ma_hit_t_alloc* x, uint32_t qn, uint32_t tn);
void debug_cut_ma_hit_t(ma_hit_t_alloc* sources, long long num_sources, ma_sub_t *coverage_cut);
typedef struct {
uint64_t ul;
@@ -139,7 +128,7 @@ static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, fl
///ext5 and ext3 is the hang on left side and right side, respectively
ext5 = qs < tl5? qs : tl5;
ext3 = ql - h->qe < tl3? ql - h->qe : tl3;
ext3 = ql - (int)h->qe < tl3? ql - (int)h->qe : tl3;
/**
@@ -179,11 +168,11 @@ static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, fl
********************************target-to-query overlap****************************
**/
if (qs <= tl5 && ql - h->qe <= tl3) return MA_HT_QCONT; // query contained in target
else if (qs >= tl5 && ql - h->qe >= tl3) return MA_HT_TCONT; // target contained in query
if (qs <= tl5 && ql - (int)h->qe <= tl3) return MA_HT_QCONT; // query contained in target
else if (qs >= tl5 && ql - (int)h->qe >= tl3) return MA_HT_TCONT; // target contained in query
else if (qs > tl5) u = 0, v = !!h->rev, l = qs - tl5; ///u = 0 means query-to-target overlap, l is the length of node in string graph (not the overlap length)
else u = 1, v = !h->rev, l = (ql - h->qe) - tl3; ///u = 1 means target-to-query overlaps, l is the length of node in string graph (not the overlap length)
if (h->qe - qs + ext5 + ext3 < min_ovlp || h->te - h->ts + ext5 + ext3 < min_ovlp) return MA_HT_SHORT_OVLP; // short overlap
if ((int)h->qe - qs + ext5 + ext3 < min_ovlp || (int)h->te - (int)h->ts + ext5 + ext3 < min_ovlp) return MA_HT_SHORT_OVLP; // short overlap
///u = 0 / 1 means query-to-target / target-to-query overlaps,
///l is the length of node in string graph (not the overlap length between two reads)
u |= h->qns>>32<<1, v |= h->tn<<1;
@@ -206,11 +195,6 @@ static inline int ma_hit2arc(const ma_hit_t *h, int ql, int tl, int max_hang, fl
}
void build_string_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen,
long long mini_overlap_length, long long max_hang_length,
long long clean_round, float min_ovlp_drop_ratio, float max_ovlp_drop_ratio,
float final_ovlp_drop_ratio, char* output_file_name, long long bubble_dist);
#define asg_arc_len(arc) ((uint32_t)(arc).ul)
#define asg_arc_n(g, v) ((uint32_t)(g)->idx[(v)])
@@ -296,13 +280,7 @@ typedef struct {
// count the number of outgoing arcs, including reduced arcs
static inline int count_out_with_del(const asg_t *g, uint32_t v)
{
uint32_t i, n, nv = asg_arc_n(g, v);
const asg_arc_t *av = asg_arc_a(g, v);
/**
for (i = n = 0; i < nv; ++i)
if (!av[i].del) ++n;
return n;
**/
uint32_t nv = asg_arc_n(g, v);
return nv;
}
@@ -318,14 +296,13 @@ static inline int count_out_without_del(const asg_t *g, uint32_t v)
return n;
}
void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, int id, char* fun);
void build_string_graph_without_clean(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen,
long long mini_overlap_length, long long max_hang_length,
long long clean_round, float min_ovlp_drop_ratio, float max_ovlp_drop_ratio,
float corase_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph,
int write);
void build_string_graph_without_clean(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long n_read, uint64_t* readLen, long long mini_overlap_length,
long long max_hang_length, long long clean_round, long long pop_bubble_size,
float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, int write);
void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, int id, char* command);
POA.cpp
+9 -583
View File
@@ -191,7 +191,7 @@ void init_Node_alloc(Node_alloc* list)
list->sort.iterative_buffer_visit = NULL;
long long i;
uint64_t i;
for (i = 0; i < list->size; i++)
{
list->list[i].insertion_edges.list=NULL;
@@ -238,10 +238,8 @@ uint64_t append_Node_alloc(Node_alloc* list, char base)
if (list->length + 1 > list->size)
{
long long i = list->size;
uint64_t i = list->size;
///list->topo_order这里用不到，所以不用先分配空间
///但是还是一起分配了吧，免得麻烦
list->size = list->size * 2;
list->list = (Node*)realloc(list->list, sizeof(Node)*list->size);
///list->topo_order = (uint64_t*)realloc(list->topo_order, sizeof(uint64_t)*list->size);
@@ -336,8 +334,6 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
for (i = 0; i < g_read_length; i++)
{
nodeID = add_Node_Graph(g, g_read_seq[i]);
////fprintf(stderr, "nodeID: %llu\n", nodeID);
if (firstID == -1)
{
@@ -345,12 +341,7 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
}
if (lastID != -1)
{
/**
///0是match边
add_Edge_Graph(g, lastID, nodeID, 0);
**/
///只有match边长度是0
///mismatch边长度都是1
///the legnth of match edge is 0, while the length of musmatch is 1
append_Edge_alloc(&(g->g_nodes.list[lastID].mismatch_edges), lastID, nodeID, 1, 0);
}
@@ -367,287 +358,6 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
inline void add_insertionEdge_weight_print(Graph* g, long long alignNodeID, char* insert, long long insert_length)
{
long long nodeID;
long long edgeID;
Edge_alloc* edge = &(g->g_nodes.list[alignNodeID].insertion_edges);
/******************************for homopolymer*************************/
long long i = 0;
char hom;
if (insert_length > 0)
{
hom = insert[0];
}
for (i = 0; i < insert_length; i++)
{
if(insert[i] != hom)
{
break;
}
}
fprintf(stderr, "###insert_length: %d\n", insert_length);
///if it is a homopolymer
if(i == insert_length)
{
///single base
edgeID = getEdge(g, edge, 1, insert[0]);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else ///不存在这么一条边
{
nodeID = add_Node_Graph(g, insert[0]);
append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
///将新加入的节点通过insertion_edges接回backbone上
///应该连回到原节点，而不是原节点的下一个节点
///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
}
///multiple bases
for (i = 1; i < insert_length; i++)
{
edgeID = get_insertion_Edges(g, edge, i + 1, insert);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else
{
create_insertion_Edges(g, alignNodeID, i + 1, insert);
}
}
return;
}
/******************************for homopolymer*************************/
fprintf(stderr, "###not homopolymer: %d\n", insert_length);
if (insert_length == 1)
{
edgeID = getEdge(g, edge, 1, insert[0]);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else ///不存在这么一条边
{
nodeID = add_Node_Graph(g, insert[0]);
append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
///将新加入的节点通过insertion_edges接回backbone上
///应该连回到原节点，而不是原节点的下一个节点
///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
}
}
else if (insert_length == 2)
{
/*******************第0个字符********************* */
edgeID = getEdge(g, edge, 1, insert[0]);
fprintf(stderr, "edgeID[0]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else ///不存在这么一条边
{
nodeID = add_Node_Graph(g, insert[0]);
append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
///将新加入的节点通过insertion_edges接回backbone上
///应该连回到原节点，而不是原节点的下一个节点
///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
}
fprintf(stderr, "edge->length: %d\n", edge->length);
/*******************第0个字符********************* */
/*******************第1个字符********************* */
if (insert[1] != insert[0])
{
edgeID = getEdge(g, edge, 1, insert[1]);
fprintf(stderr, "edgeID[1]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else ///不存在这么一条边
{
nodeID = add_Node_Graph(g, insert[1]);
append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
///将新加入的节点通过insertion_edges接回backbone上
///应该连回到原节点，而不是原节点的下一个节点
///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
}
fprintf(stderr, "edge->length: %d\n", edge->length);
}
/*******************第1个字符********************* */
/**********************两个字符******************* */
edgeID = get_insertion_Edges(g, edge, 2, insert);
fprintf(stderr, "edgeID[len2]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else
{
create_insertion_Edges(g, alignNodeID, insert_length, insert);
}
fprintf(stderr, "edge->length: %d\n", edge->length);
for (i = 0; i < edge->length; i++)
{
fprintf(stderr, "edgeID[%d].length: %d\n", i, edge->list[i].length);
}
/**********************两个字符******************* */
}
else if (insert_length > 2)
{
////fprintf(stderr, "too long insertion\n");
/*************************大于2个字符************************** */
edgeID = get_insertion_Edges(g, edge, insert_length, insert);
if (edgeID != -1)
{
///这条路均只有一个出度
edge->list[edgeID].weight++;
}
else
{
create_insertion_Edges(g, alignNodeID, insert_length, insert);
}
}
}
void addmatchedSeqToGraph_print(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end)
{
int x_i, y_i, cigar_i;
x_i = 0;
y_i = 0;
cigar_i = 0;
int operation;
int operationLen;
int i;
int last_operation = -1;
fprintf(stderr,"*******\n");
///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多字符）
while (cigar_i < cigar->length)
{
operation = cigar->C_C[cigar_i];
operationLen = cigar->C_L[cigar_i];
// fprintf(stderr, "operation: %d, operationLen: %d\n",
// operation, operationLen);
///这种情况代表匹配和mismatch
if (operation == 0 || operation == 1)
{
for (i = 0; i < operationLen; i++)
{
//backbone->g_nodes.list[currentNodeID].weight++;
///前面是插入，后面有可能是误配，也有可能是匹配
add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_operation);
x_i++;
y_i++;
currentNodeID++;
}
}///insertion
else if (operation == 2)
{
///cigar的起始和结尾不可能是2，所以这里-1没问题
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
add_insertionEdge_weight_print(backbone, currentNodeID, y_string + y_i, operationLen);
backbone->g_nodes.list[currentNodeID].num_insertions++;
}
///fprintf(stderr, "y_string: %.*s\n", operationLen, y_string+y_i);
y_i += operationLen;
}
else if (operation == 3)
{
///3是y缺字符（x多字符），也就是backbone多字符
///这个相当于在backbone对应字符处变成了‘——’
///因此可以用mismatch类似的方法处理
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
///add_deletion_to_backbone(backbone, &currentNodeID, operationLen);
///在编辑距离中，前面是个插入，后面是个删除，这种情况是不存在的
///为了保险要不还给他加上吧
///先不加
add_deletionEdge_weight(backbone, currentNodeID, operationLen);
}
currentNodeID += operationLen;
x_i += operationLen;
}
last_operation = operation;
cigar_i++;
}
/**
///cigar的起始和结尾不可能是2
if (cigar->C_C[0] == 2 || cigar->C_C[cigar->length - 1] == 2)
{
fprintf(stderr, "error\n");
}
if (x_i != x_length)
{
fprintf(stderr, "x_i: %d, x_length: %d\n", x_i, x_length);
}
if (y_i != y_length)
{
fprintf(stderr, "y_i: %d, y_length: %d\n", y_i, y_length);
}
**/
}
void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end)
{
@@ -664,20 +374,19 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
///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多字符）
///2 mean y has more bases, while 3 means x has more bases
while (cigar_i < cigar->length)
{
operation = cigar->C_C[cigar_i];
operationLen = cigar->C_L[cigar_i];
///这种情况代表匹配和mismatch
///match/mismatch
if (operation == 0 || operation == 1)
{
for (i = 0; i < operationLen; i++)
{
//backbone->g_nodes.list[currentNodeID].weight++;
///前面是插入，后面有可能是误配，也有可能是匹配
///if the previous node is insertion, this node might be mismatch/match
add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_operation);
x_i++;
y_i++;
@@ -686,7 +395,7 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
}///insertion
else if (operation == 2)
{
///cigar的起始和结尾不可能是2，所以这里-1没问题
///the begin and end of cigar cannot be 2, so -1 is right here
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
add_insertionEdge_weight(backbone, currentNodeID, y_string + y_i, operationLen);
@@ -696,15 +405,10 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
}
else if (operation == 3)
{
///3是y缺字符（x多字符），也就是backbone多字符
///这个相当于在backbone对应字符处变成了‘——’
///因此可以用mismatch类似的方法处理
///3 means x has more bases, that means backbone has more bases
///like a mismatch (-)
///if (operationLen <= CORRECT_INDEL_LENGTH)
{
///add_deletion_to_backbone(backbone, &currentNodeID, operationLen);
///在编辑距离中，前面是个插入，后面是个删除，这种情况是不存在的
///为了保险要不还给他加上吧
///先不加
add_deletionEdge_weight(backbone, currentNodeID, operationLen);
}
@@ -720,284 +424,6 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
}
void debug_graph(Graph* g, long long backbone_length)
{
long long i = 0;
if (g->s_start_nodeID != 0 || g->s_end_nodeID != backbone_length)
{
fprintf(stderr, "error\n");
}
for (i = g->s_start_nodeID; i <= g->s_end_nodeID; i++)
{
if(g->g_nodes.list[i].weight != 1)
{
fprintf(stderr, "error node weight\n");
}
if(g->g_nodes.list[i].mismatch_edges.length > 4)
{
fprintf(stderr, "error mismatch_edges\n");
}
if(g->g_nodes.list[i].mismatch_edges.length < 1 && i != g->s_end_nodeID)
{
fprintf(stderr, "i: %d, error mismatch_edges: %d\n", i, g->g_nodes.list[i].mismatch_edges.length);
}
}
for (i = 0; i < g->g_nodes.length; i++)
{
if(g->g_nodes.list[i].ID < g->s_start_nodeID || g->g_nodes.list[i].ID > g->s_end_nodeID)
{
if (g->g_nodes.list[i].deletion_edges.length +
g->g_nodes.list[i].insertion_edges.length +
g->g_nodes.list[i].mismatch_edges.length
!= 1)
{
fprintf(stderr, "g->s_start_nodeID: %lld\n",
g->s_start_nodeID);
fprintf(stderr, "g->s_end_nodeID: %lld\n",
g->s_end_nodeID);
fprintf(stderr, "deletion_edges_length: %lld\n",
g->g_nodes.list[i].deletion_edges.length);
fprintf(stderr, "insertion_edges_length: %lld, \n",
g->g_nodes.list[i].insertion_edges.length);
fprintf(stderr, "g->g_nodes.list[i].insertion_edges.list[0].length: %lld, \n",
g->g_nodes.list[i].insertion_edges.list[0].length);
fprintf(stderr, "g->g_nodes.list[i].insertion_edges.list[1].length: %lld, \n",
g->g_nodes.list[i].insertion_edges.list[1].length);
fprintf(stderr, "mismatch_edges_length: %lld\n",
g->g_nodes.list[i].mismatch_edges.length);
}
else
{
///不是0肯定是1
if (g->g_nodes.list[i].deletion_edges.length != 0)
{
long long step = g->g_nodes.list[i].deletion_edges.list[0].length;
long long nodeID = i;
for (int j = 0; j < step; j++)
{
nodeID = g->g_nodes.list[nodeID].deletion_edges.list[0].out_node;
}
nodeID = g->g_nodes.list[nodeID].deletion_edges.list[0].out_node;
if (nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID)
{
fprintf(stderr, "error\n");
}
}
if (g->g_nodes.list[i].insertion_edges.length != 0)
{
long long step = g->g_nodes.list[i].insertion_edges.list[0].length;
long long nodeID = i;
for (int j = 0; j < step; j++)
{
nodeID = g->g_nodes.list[nodeID].insertion_edges.list[0].out_node;
}
nodeID = g->g_nodes.list[nodeID].insertion_edges.list[0].out_node;
if ((nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID))
{
fprintf(stderr, "error: step: %d\n", step);
}
}
if (g->g_nodes.list[i].mismatch_edges.length != 0)
{
long long step = g->g_nodes.list[i].mismatch_edges.list[0].length;
long long nodeID = i;
for (int j = 0; j < step; j++)
{
nodeID = g->g_nodes.list[nodeID].mismatch_edges.list[0].out_node;
}
nodeID = g->g_nodes.list[nodeID].mismatch_edges.list[0].out_node;
if (nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID)
{
fprintf(stderr, "error\n");
}
}
}
}
}
}
void Graph_debug(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end)
{
/**
int x_i, y_i, cigar_i;
x_i = 0;
y_i = 0;
cigar_i = 0;
int operation;
int operationLen;
int i;
///0 is match, 1 is mismatch, 2 is up, 3 is left
///2是x缺字符（y多字符），而3是y缺字符（x多字符）
///while (x_i < x_len && y_i < y_len && cigar_i < cigar->length)
while (cigar_i < cigar->length)
{
operation = cigar->C_C[cigar_i];
operationLen = cigar->C_L[cigar_i];
///这种情况代表匹配
if (operation == 0)
{
for (i = 0; i < operationLen; i++)
{
if (backbone->g_nodes.list[currentNodeID].base != y_string[y_i])
{
fprintf(stderr, "error match\n");
}
backbone->g_nodes.list[currentNodeID].weight--;
x_i++;
y_i++;
currentNodeID++;
}
}
else if (operation == 1)
{
for (i = 0; i < operationLen; i++)
{
if (backbone->g_nodes.list[currentNodeID].base == y_string[y_i])
{
fprintf(stderr, "error mismatch 1\n");
}
long long mismatchID = get_alignToNode(backbone, currentNodeID, y_string[y_i]);
if(mismatchID == -1)
{
fprintf(stderr, "error mismatch 2\n");
}
else
{
backbone->g_nodes.list[mismatchID].weight--;
}
x_i++;
y_i++;
currentNodeID++;
}
}
else if (operation == 2)
{
long long nodeID = currentNodeID - 1;
long long mismatchID;
for (i = 0; i < operationLen; i++)
{
mismatchID = get_insertion_Node(backbone, nodeID, y_string[y_i]);
if (mismatchID == -1)
{
fprintf(stderr, "error insertion 1, i: %d\n", i);
}
else
{
backbone->g_nodes.list[mismatchID].weight--;
}
nodeID = mismatchID;
y_i++;
}
///注意这里是x_string[x_i]而不是x_string[currentNodeID]
mismatchID = get_insertion_Node(backbone, nodeID, x_string[x_i]);
if (mismatchID == -1)
{
fprintf(stderr, "error insertion 2, i: %d, x_i: %d\n", i, x_i);
}
if (mismatchID != currentNodeID)
{
fprintf(stderr, "error insertion 3, i: mismatchID: %d, currentNodeID: %d\n", mismatchID, currentNodeID);
}
}
else if (operation == 3)
{
for (i = 0; i < operationLen; i++)
{
long long mismatchID = get_alignToNode(backbone, currentNodeID, 'D');
if(mismatchID == -1)
{
fprintf(stderr, "error deletion 2\n");
}
else
{
backbone->g_nodes.list[mismatchID].weight--;
}
x_i++;
currentNodeID++;
}
}
cigar_i++;
}
if (cigar->C_C[0] == 2 || cigar->C_C[cigar->length - 1] == 2)
{
fprintf(stderr, "error\n");
}
if (x_i != x_length)
{
fprintf(stderr, "x_i: %d, x_length: %d\n", x_i, x_length);
}
if (y_i != y_length)
{
fprintf(stderr, "y_i: %d, y_length: %d\n", y_i, y_length);
}
**/
}
POA.h
+45 -341
View File
@@ -4,27 +4,6 @@
#include "Hash_Table.h"
#include "Process_Read.h"
/**
1. 单个节点信息
(1) ID
(2) base
(3) 入边信息
(4) 出边信息
(5) 比对到什么节点
2. 各个节点信息，用数组下标组织，数组下标就是节点ID; 还要存拓扑排序后的下标和节点ID的对应关系
3. 边
(1) 边的起始
(2) 边的结束节点
(3) 过这条边的序列的label，也就是名称
4. 各个序列信息
(1) 这个序列本身
(2) 这个序列的name或者ID
(3) 这个序列的在图中对应的起始和结束节点ID
5. 有两个回溯矩阵，一个是graph的，一个是seq的
**/
typedef struct
{
long long beg;
@@ -141,10 +120,9 @@ typedef struct
{
uint64_t in_node;
uint64_t out_node;
///0是match，1是mismatch，2是x缺字符（y多字符），而3是y缺字符（x多字符）
///0 is match，1 is mismatch，2 means y has more bases, 3 means x has more bases
uint64_t weight;
uint64_t num_insertions;
///这条路径上到backbone节点之前总共有多少节点
uint64_t length;
uint64_t self_edge_ID;
uint64_t reverse_edge_ID;
@@ -180,7 +158,7 @@ typedef struct
{
uint64_t ID;
uint64_t weight;
///记录的是以当前节点为尾的deletion个数
///number of deletion end with current node
uint64_t num_insertions;
char base;
Edge_alloc mismatch_edges;
@@ -245,18 +223,20 @@ inline int Pop_Node(Graph* DAGCon, Node** node)
inline int Push_Node(Graph* DAGCon, Node** node)
{
push_to_Queue(&(DAGCon->node_q), (**node).ID);
return 1;
}
inline int getInputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
{
if(set->index >= Input_Edges(*node).length)
if(set->index >= (long long)Input_Edges(*node).length)
{
return 0;
}
///skip all deleted edges
while (
set->index < Input_Edges(*node).length
set->index < (long long)Input_Edges(*node).length
&&
!(If_Edge_Exist(Input_Edges(*node).list[set->index]))
)
@@ -266,7 +246,7 @@ inline int getInputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
if(
set->index < Input_Edges(*node).length
set->index < (long long)Input_Edges(*node).length
&&
If_Edge_Exist(Input_Edges(*node).list[set->index])
)
@@ -285,14 +265,14 @@ inline int getInputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
inline int getInputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge)
{
if(set->index >= Input_Edges(*node).length)
if(set->index >= (long long)Input_Edges(*node).length)
{
return 0;
}
///skip all deleted edges
while (
set->index < Input_Edges(*node).length
set->index < (long long)Input_Edges(*node).length
&&
!(If_Edge_Exist(Input_Edges(*node).list[set->index]))
)
@@ -302,7 +282,7 @@ inline int getInputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge)
if(
set->index < Input_Edges(*node).length
set->index < (long long)Input_Edges(*node).length
&&
If_Edge_Exist(Input_Edges(*node).list[set->index])
)
@@ -320,14 +300,14 @@ inline int getInputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge)
inline int getOutputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
{
if(set->index >= Output_Edges(*node).length)
if(set->index >= (long long)Output_Edges(*node).length)
{
return 0;
}
///skip all deleted edges
while (
set->index < Output_Edges(*node).length
set->index < (long long)Output_Edges(*node).length
&&
!(If_Edge_Exist(Output_Edges(*node).list[set->index]))
)
@@ -335,7 +315,7 @@ inline int getOutputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
set->index++;
}
if(set->index < Output_Edges(*node).length &&
if(set->index < (long long)Output_Edges(*node).length &&
If_Edge_Exist(Output_Edges(*node).list[set->index]))
{
(*get_Node) = &(G_Node((*graph), Output_Edges(*node).list[set->index].out_node));
@@ -352,14 +332,14 @@ inline int getOutputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node)
inline int getOutputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge)
{
if(set->index >= Output_Edges(*node).length)
if(set->index >= (long long)Output_Edges(*node).length)
{
return 0;
}
///skip all deleted edges
while (
set->index < Output_Edges(*node).length
set->index < (long long)Output_Edges(*node).length
&&
!(If_Edge_Exist(Output_Edges(*node).list[set->index]))
)
@@ -367,7 +347,7 @@ inline int getOutputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge)
set->index++;
}
if(set->index < Output_Edges(*node).length &&
if(set->index < (long long)Output_Edges(*node).length &&
If_Edge_Exist(Output_Edges(*node).list[set->index]))
{
(*get_Edge) = &(Output_Edges(*node).list[set->index]);
@@ -389,9 +369,9 @@ inline void get_bi_direction_edges(Graph* DAGCon, Edge* edge, Edge** e_forward,
if(
edge->self_edge_ID < Output_Edges(G_Node(*DAGCon, in_node)).length
&&
Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].in_node == in_node
(long long)Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].in_node == in_node
&&
Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].out_node == out_node
(long long)Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].out_node == out_node
)
{
(*e_forward) = &(Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID]);
@@ -456,10 +436,7 @@ void destory_Graph(Graph* g);
void clear_Graph(Graph* g);
void Perform_POA(Graph* g, overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read);
void Graph_debug(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end);
void debug_graph(Graph* g, long long backbone_length);
uint64_t inline add_Node_Graph(Graph* g, char base)
{
@@ -503,26 +480,28 @@ uint64_t inline delete_Node_DAGCon(Graph* g, Node* node)
clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].insertion_edges));
clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].mismatch_edges));
clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].deletion_edges));
return 1;
}
///仅仅用于误配边
///just for mimatch edges
inline void add_mismatchEdge_weight(Graph* g, uint64_t in_node, char base, int last_operation)
{
long long i = 0;
long long nodeID;
Edge_alloc* edge = &(g->g_nodes.list[in_node].mismatch_edges);
for (i = 0; i < edge->length; i++)
for (i = 0; i < (long long)edge->length; i++)
{
nodeID = edge->list[i].out_node;
if(g->g_nodes.list[nodeID].base == base)
{
edge->list[i].weight++;
///如果上一个操作是insertion
///if last operation is insertion
if (last_operation == 2)
{
edge->list[i].num_insertions++;
@@ -532,25 +511,22 @@ inline void add_mismatchEdge_weight(Graph* g, uint64_t in_node, char base, int l
}
}
///说明不存在这么一条边
if (i == edge->length)
///there are no such edge
if (i == (long long)edge->length)
{
nodeID = add_Node_Graph(g, base);
///只有match边长度是0
///mismatch边长度都是1
///the length of match edge is 0, while the length of mismatch edge is 1
append_Edge_alloc(edge, in_node, nodeID, 1, 1);
///如果上一个操作是insertion
///if last operation is insertion
if (last_operation == 2)
{
edge->list[edge->length - 1].num_insertions++;
}
///将新节点的mismatch_edges连到backbone上
///add the mismatch_edges of new node to the backbone
append_Edge_alloc(&(g->g_nodes.list[nodeID].mismatch_edges), nodeID, in_node + 1, 1, 0);
}
///获得节点的mismatch_edges长度为1，其他均为0
}
@@ -561,7 +537,7 @@ inline void add_single_deletionEdge_weight(Graph* g, long long alignNodeID, long
long long nodeID;
Edge_alloc* edge = &(g->g_nodes.list[alignNodeID].deletion_edges);
for (i = 0; i < edge->length; i++)
for (i = 0; i < (long long)edge->length; i++)
{
nodeID = edge->list[i].out_node;
if(nodeID == nextNodeID)
@@ -571,8 +547,8 @@ inline void add_single_deletionEdge_weight(Graph* g, long long alignNodeID, long
}
}
///说明不存在这么一条边
if (i == edge->length)
///there are no such edge
if (i == (long long)edge->length)
{
append_Edge_alloc(edge, alignNodeID, nextNodeID, 1, edge_length);
}
@@ -580,23 +556,6 @@ inline void add_single_deletionEdge_weight(Graph* g, long long alignNodeID, long
inline void add_deletionEdge_weight(Graph* g, long long alignNodeID, long long deletion_length)
{
/**
if (deletion_length == 1)
{
add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 1, 0);
}
else if (deletion_length == 2)
{
add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 1, 0);
add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 2, 0);
add_single_deletionEdge_weight(g, alignNodeID + 1, alignNodeID + 2, 0);
}
else if (deletion_length > 2)
{
///fprintf(stderr, "too long deletion!\n");
add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + deletion_length, 0);
}
**/
long long i;
for (i = 0; i < deletion_length; i++)
{
@@ -611,7 +570,7 @@ inline int getEdge(Graph* g, Edge_alloc* edge, uint64_t edge_length, char base)
long long i = 0;
long long nodeID;
for (i = 0; i < edge->length; i++)
for (i = 0; i < (long long)edge->length; i++)
{
if (edge->list[i].length == edge_length)
{
@@ -651,7 +610,7 @@ inline int get_insertion_Edges(Graph* g, Edge_alloc* edge, uint64_t edge_length,
Edge_alloc* new_edge = edge;
for (i = 1; i < edge_length; i++)
for (i = 1; i < (long long)edge_length; i++)
{
nodeID = new_edge->list[edgeID].out_node;
new_edge = &(g->g_nodes.list[nodeID].insertion_edges);
@@ -662,7 +621,6 @@ inline int get_insertion_Edges(Graph* g, Edge_alloc* edge, uint64_t edge_length,
}
}
/****************************may have bugs********************************/
///return edgeID;
return return_edgeID;
/****************************may have bugs********************************/
}
@@ -673,7 +631,7 @@ inline int create_insertion_Edges(Graph* g, long long alignNodeID, uint64_t edge
{
long long i = 0;
long long nodeID;
///最后应该连回原节点
///should link back to the intial node
///long long backboneID = alignNodeID + 1;
long long backboneID = alignNodeID;
@@ -685,28 +643,29 @@ inline int create_insertion_Edges(Graph* g, long long alignNodeID, uint64_t edge
nodeID = add_Node_Graph(g, bases[0]);
///将新加入的节点通过insertion_edges接到alignNodeID上
///add the new node to alignNodeID by insertion_edges
append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, nodeID, 1, edge_length);
alignNodeID = nodeID;
for (i = 1; i < edge_length; i++)
for (i = 1; i < (long long)edge_length; i++)
{
nodeID = add_Node_Graph(g, bases[i]);
///将新加入的节点通过insertion_edges接到alignNodeID上
///add the new node to alignNodeID by insertion_edges
append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, nodeID, 1, edge_length - i);
alignNodeID = nodeID;
}
append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, backboneID, 1, 0);
return 1;
}
inline void extract_path(Graph* backbone, int debug_node_in_backbone, int path_i, char* pre)
{
int step = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length;
int string_i, preNode, j;
int string_i = 0, preNode = 0, j = 0;
if(step != 0)
{
string_i = 0;
@@ -723,65 +682,17 @@ inline void extract_path(Graph* backbone, int debug_node_in_backbone, int path_i
}
inline void extract_path_debug(Graph* backbone, int debug_node_in_backbone, int path_i, char* pre)
{
int step = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length;
int string_i, preNode, preEdge, j;
if(step != 0)
{
string_i = 0;
preNode = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].out_node;
preEdge = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length;
for (j = 0; j < step; j++)
{
///pre[string_i++] = G_Node(*backbone, preNode).base;
fprintf(stderr, "j: %d (%c%d), ", j, G_Node(*backbone, preNode).base, preEdge);
preEdge = G_Node(*backbone, preNode).insertion_edges.list[0].length;
preNode = G_Node(*backbone, preNode).insertion_edges.list[0].out_node;
}
}
fprintf(stderr, "\n");
///pre[string_i] = '\0';
}
inline int getEdge_DEBUG(Graph* g, Edge_alloc* edge, uint64_t edge_length, char base)
{
long long i = 0;
long long nodeID;
for (i = 0; i < edge->length; i++)
{
///fprintf(stderr, "************i:%d, edge->list[i].length: %d, edge_length: %d\n",i, edge->list[i].length, edge_length);
if (edge->list[i].length == edge_length)
{
nodeID = edge->list[i].out_node;
fprintf(stderr, "########i:%d, edge->list[i].length: %d, edge_length: %d, nodeID: %d, list[nodeID].base: %c, base: %c\n",
i, edge->list[i].length, edge_length, nodeID, g->g_nodes.list[nodeID].base, base);
if(g->g_nodes.list[nodeID].base == base)
{
return i;
}
}
}
return -1;
}
inline int get_insertion_Edges_new(Graph* backbone, int debug_node_in_backbone, uint64_t edge_length, char* bases)
{
int path_i, j, step, preNode;
for (path_i = 0; path_i < G_Node(*backbone, debug_node_in_backbone).insertion_edges.length; path_i++)
for (path_i = 0; path_i < (long long)G_Node(*backbone, debug_node_in_backbone).insertion_edges.length; path_i++)
{
step = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length;
if(step != edge_length)
if(step != (long long)edge_length)
{
continue;
}
@@ -793,8 +704,6 @@ inline int get_insertion_Edges_new(Graph* backbone, int debug_node_in_backbone,
for (j = 0; j < step; j++)
{
///pre[string_i++] = G_Node(*backbone, preNode).base;
///fprintf(stderr, "path_i: %d, ID: %d\n", path_i, G_Node(*backbone, preNode).ID);
if(G_Node(*backbone, preNode).base != bases[j])
{
break;
@@ -815,53 +724,7 @@ inline int get_insertion_Edges_new(Graph* backbone, int debug_node_in_backbone,
}
inline int get_insertion_Edges_debug(Graph* g, Edge_alloc* edge, uint64_t edge_length, char* bases)
{
long long i = 0;
long long nodeID;
long long edgeID;
if (edge_length < 1)
{
return -1;
}
///fprintf(stderr, "edge_length: %d, edge: %.*s\n", edge_length, edge_length, bases);
edgeID = getEdge_DEBUG(g, edge, edge_length, bases[0]);
fprintf(stderr, "i: %d, edgeID: %d, edge_length - i: %d\n", i, edgeID, edge_length);
long long return_edgeID = edgeID;
if(edgeID == -1)
{
return -1;
}
Edge_alloc* new_edge = edge;
for (i = 1; i < edge_length; i++)
{
nodeID = new_edge->list[edgeID].out_node;
new_edge = &(g->g_nodes.list[nodeID].insertion_edges);
edgeID = getEdge_DEBUG(g, new_edge, edge_length - i, bases[i]);
fprintf(stderr, "i: %d, edgeID: %d, edge_length - i: %d\n", i, edgeID, edge_length - i);
if(edgeID == -1)
{
return -1;
}
}
/****************************may have bugs********************************/
///return edgeID;
return return_edgeID;
/****************************may have bugs********************************/
}
inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* insert, long long insert_length)
{
@@ -873,20 +736,16 @@ inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* inse
if (insert_length == 1)
{
edgeID = getEdge(g, edge, 1, insert[0]);
// if(edgeID != get_insertion_Edges_new(g, alignNodeID, insert_length, insert))
// {
// fprintf(stderr, "error\n");
// }
if (edgeID != -1)
{
edge->list[edgeID].weight++;
}
else ///不存在这么一条边
else ///there is no such edge
{
nodeID = add_Node_Graph(g, insert[0]);
append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
///将新加入的节点通过insertion_edges接回backbone上
///应该连回到原节点，而不是原节点的下一个节点
///add the new node to alignNodeID by insertion_edges
//should link to the initial node, instead of the next node of the initial node
///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
}
@@ -897,7 +756,7 @@ inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* inse
edgeID = get_insertion_Edges_new(g, alignNodeID, insert_length, insert);
if (edgeID != -1)
{
///这条路均只有一个出度
///just one outdegree
edge->list[edgeID].weight++;
}
else
@@ -905,163 +764,8 @@ inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* inse
create_insertion_Edges(g, alignNodeID, insert_length, insert);
}
}
// /******************************for homopolymer*************************/
// long long i = 0;
// char hom;
// if (insert_length > 0)
// {
// hom = insert[0];
// }
// for (i = 0; i < insert_length; i++)
// {
// if(insert[i] != hom)
// {
// break;
// }
// }
// ///if it is a homopolymer
// if(i == insert_length)
// {
// ///single base
// edgeID = getEdge(g, edge, 1, insert[0]);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else ///不存在这么一条边
// {
// nodeID = add_Node_Graph(g, insert[0]);
// append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
// ///将新加入的节点通过insertion_edges接回backbone上
// ///应该连回到原节点，而不是原节点的下一个节点
// ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
// append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
// }
// ///multiple bases
// for (i = 1; i < insert_length; i++)
// {
// edgeID = get_insertion_Edges(g, edge, i + 1, insert);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else
// {
// create_insertion_Edges(g, alignNodeID, i + 1, insert);
// }
// }
// return;
// }
// /******************************for homopolymer*************************/
// if (insert_length == 1)
// {
// edgeID = getEdge(g, edge, 1, insert[0]);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else ///不存在这么一条边
// {
// nodeID = add_Node_Graph(g, insert[0]);
// append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
// ///将新加入的节点通过insertion_edges接回backbone上
// ///应该连回到原节点，而不是原节点的下一个节点
// ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
// append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
// }
// }
// else if (insert_length == 2)
// {
// /*******************第0个字符********************* */
// edgeID = getEdge(g, edge, 1, insert[0]);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else ///不存在这么一条边
// {
// nodeID = add_Node_Graph(g, insert[0]);
// append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
// ///将新加入的节点通过insertion_edges接回backbone上
// ///应该连回到原节点，而不是原节点的下一个节点
// ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
// append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
// }
// /*******************第0个字符********************* */
// /*******************第1个字符********************* */
// if (insert[1] != insert[0])
// {
// edgeID = getEdge(g, edge, 1, insert[1]);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else ///不存在这么一条边
// {
// nodeID = add_Node_Graph(g, insert[1]);
// append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1);
// ///将新加入的节点通过insertion_edges接回backbone上
// ///应该连回到原节点，而不是原节点的下一个节点
// ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0);
// append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0);
// }
// }
// /*******************第1个字符********************* */
// /**********************两个字符******************* */
// edgeID = get_insertion_Edges(g, edge, 2, insert);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else
// {
// create_insertion_Edges(g, alignNodeID, insert_length, insert);
// }
// /**********************两个字符******************* */
// }
// else if (insert_length > 2)
// {
// ////fprintf(stderr, "too long insertion\n");
// /*************************大于2个字符************************** */
// edgeID = get_insertion_Edges(g, edge, insert_length, insert);
// if (edgeID != -1)
// {
// ///这条路均只有一个出度
// edge->list[edgeID].weight++;
// }
// else
// {
// create_insertion_Edges(g, alignNodeID, insert_length, insert);
// }
// }
}
void addmatchedSeqToGraph_print(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end);
#endif
Process_Read.cpp
+81 -194
View File
@@ -21,30 +21,44 @@ pthread_cond_t i_readinputstallCond;
pthread_mutex_t i_doneMutex;
uint8_t seq_nt6_table[256] = {
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5
};
char bit_t_seq_table[256][4] = {0};
char bit_t_seq_table_rc[256][4] = {0};
char s_H[5] = {'A', 'C', 'G', 'T', 'N'};
char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'};
void init_All_reads(All_reads* r)
{
r->index_size = READ_INIT_NUMBER;
/**********should remove**********/
///r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size);
///r->index[0] = 0;
///r->read = NULL;
/**********should remove**********/
r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size);
r->read_sperate = NULL;
r->N_site = NULL;
r->total_reads_bases = 0;
r->name_index_size = READ_INIT_NUMBER;
r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size);
r->name_index[0] = 0;
r->name = NULL;
r->total_name_length = 0;
r->total_reads = 0;
}
void destory_All_reads(All_reads* r)
@@ -60,24 +74,19 @@ void destory_All_reads(All_reads* r)
}
free(r->N_site);
free(r->read_sperate);
///free(r->read);
free(r->name);
free(r->name_index);
free(r->read_length);
}
void write_All_reads(All_reads* r, char* read_file_name)
{
fprintf(stdout, "Writing reads to disk ...... \n");
fprintf(stderr, "Writing reads to disk... \n");
char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "w");
fwrite(&adapterLen, sizeof(adapterLen), 1, fp);
fwrite(&asm_opt.adapterLen, sizeof(asm_opt.adapterLen), 1, fp);
fwrite(&r->index_size, sizeof(r->index_size), 1, fp);
fwrite(&r->name_index_size, sizeof(r->name_index_size), 1, fp);
fwrite(&r->total_reads, sizeof(r->total_reads), 1, fp);
@@ -90,12 +99,10 @@ void write_All_reads(All_reads* r, char* read_file_name)
{
if (r->N_site[i] != NULL)
{
///这个实际上是N的个数
///number of Ns
fwrite(&r->N_site[i][0], sizeof(r->N_site[i][0]), 1, fp);
if (r->N_site[i][0])
{
///r->N_site[i]这实际是个长为r->N_site[i][0]+1
///这里从r->N_site[i] + 1写入了r->N_site[i][0]个元素
fwrite(r->N_site[i]+1, sizeof(r->N_site[i][0]), r->N_site[i][0], fp);
}
}
@@ -108,36 +115,25 @@ void write_All_reads(All_reads* r, char* read_file_name)
}
/**********should remove**********/
///fwrite(r->index, sizeof(uint64_t), r->index_size, fp);
/**********should remove**********/
fwrite(r->read_length, sizeof(uint64_t), r->total_reads, fp);
/**********should remove**********/
///fwrite(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp);
/**********should remove**********/
for (i = 0; i < r->total_reads; i++)
{
fwrite(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp);
}
fwrite(r->name, sizeof(char), r->total_name_length, fp);
fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
free(index_name);
fflush(fp);
fclose(fp);
fprintf(stdout, "Reads has been written.\n");
fprintf(stderr, "Reads has been written.\n");
}
int load_All_reads(All_reads* r, char* read_file_name)
{
fprintf(stdout, "Loading reads to disk ...... \n");
fprintf(stderr, "Loading reads from disk... \n");
char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "r");
@@ -145,20 +141,20 @@ int load_All_reads(All_reads* r, char* read_file_name)
{
return 0;
}
int local_adapterLen;
fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
if(local_adapterLen != adapterLen)
int f_flag;
f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
if(local_adapterLen != asm_opt.adapterLen)
{
fprintf(stdout, "the adapterLen of index is: %d, but the adapterLen set by user is: %d\n",
local_adapterLen, adapterLen);
fprintf(stderr, "the adapterLen of index is: %d, but the adapterLen set by user is: %d\n",
local_adapterLen, asm_opt.adapterLen);
exit(1);
}
fread(&r->index_size, sizeof(r->index_size), 1, fp);
fread(&r->name_index_size, sizeof(r->name_index_size), 1, fp);
fread(&r->total_reads, sizeof(r->total_reads), 1, fp);
fread(&r->total_reads_bases, sizeof(r->total_reads_bases), 1, fp);
fread(&r->total_name_length, sizeof(r->total_name_length), 1, fp);
f_flag += fread(&r->index_size, sizeof(r->index_size), 1, fp);
f_flag += fread(&r->name_index_size, sizeof(r->name_index_size), 1, fp);
f_flag += fread(&r->total_reads, sizeof(r->total_reads), 1, fp);
f_flag += fread(&r->total_reads_bases, sizeof(r->total_reads_bases), 1, fp);
f_flag += fread(&r->total_name_length, sizeof(r->total_name_length), 1, fp);
uint64_t i = 0;
uint64_t zero = 0;
@@ -166,7 +162,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
for (i = 0; i < r->total_reads; i++)
{
fread(&zero, sizeof(zero), 1, fp);
f_flag += fread(&zero, sizeof(zero), 1, fp);
if (zero)
{
@@ -175,9 +171,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->N_site[i][0] = zero;
if (r->N_site[i][0])
{
///r->N_site[i]这实际是个长为r->N_site[i][0]+1
///这里从r->N_site[i] + 1写入了r->N_site[i][0]个元素
fread(r->N_site[i]+1, sizeof(r->N_site[i][0]), r->N_site[i][0], fp);
f_flag += fread(r->N_site[i]+1, sizeof(r->N_site[i][0]), r->N_site[i][0], fp);
}
}
else
@@ -187,34 +181,25 @@ int load_All_reads(All_reads* r, char* read_file_name)
}
/**********should remove**********/
///r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size);
///fread(r->index, sizeof(uint64_t), r->index_size, fp);
/**********should remove**********/
r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
fread(r->read_length, sizeof(uint64_t), r->total_reads, fp);
f_flag += 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);
/**********should remove**********/
r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
{
r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
fread(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp);
f_flag += fread(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp);
}
r->name = (char*)malloc(sizeof(char)*r->total_name_length);
fread(r->name, sizeof(char), r->total_name_length, fp);
f_flag += fread(r->name, sizeof(char), r->total_name_length, fp);
r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size);
fread(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
f_flag += 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);
@@ -234,13 +219,9 @@ int load_All_reads(All_reads* r, char* read_file_name)
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
}
free(index_name);
fclose(fp);
fprintf(stdout, "Reads has been loaded.\n");
fprintf(stderr, "Reads has been loaded.\n");
return 1;
}
@@ -253,27 +234,17 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name)
r->total_reads_bases = r->total_reads_bases + read->l;
r->total_name_length = r->total_name_length + name->l;
///必须要+1
///must +1
if (r->index_size < r->total_reads + 2)
{
r->index_size = r->index_size * 2 + 2;
/**********should remove**********/
///r->index = (uint64_t*)realloc(r->index,sizeof(uint64_t)*(r->index_size));
/**********should remove**********/
r->read_length = (uint64_t*)realloc(r->read_length,sizeof(uint64_t)*(r->index_size));
r->name_index_size = r->name_index_size * 2 + 2;
r->name_index = (uint64_t*)realloc(r->name_index,sizeof(uint64_t)*(r->name_index_size));
}
/**********should remove**********/
///r->index[r->total_reads] = r->index[r->total_reads-1] + read->l;
/**********should remove**********/
r->read_length[r->total_reads - 1] = read->l;
//r->index[r->total_reads] = r->index[r->total_reads-1] + read->l/4 + 1;
r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l;
}
void malloc_All_reads(All_reads* r)
@@ -282,13 +253,9 @@ 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));
/**********should remove**********/
r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
long long i = 0;
for (i = 0; i < r->total_reads; i++)
for (i = 0; i < (long long)r->total_reads; i++)
{
r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
}
@@ -297,7 +264,7 @@ void malloc_All_reads(All_reads* r)
r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
for (i = 0; i < (long long)r->total_reads; i++)
{
r->second_round_cigar[i].size = r->cigars[i].size = 0;
r->second_round_cigar[i].length = r->cigars[i].length = 0;
@@ -310,15 +277,6 @@ void malloc_All_reads(All_reads* r)
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
}
r->name = (char*)malloc(sizeof(char)*r->total_name_length);
r->N_site = (uint64_t**)calloc(r->total_reads, sizeof(uint64_t*));
@@ -422,13 +380,13 @@ void recover_UC_Read_sub_region_begin_end
if (R_INF->N_site[ID])
{
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= start_pos && R_INF->N_site[ID][i] <= end_pos)
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)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)
else if((long long)R_INF->N_site[ID][i] > end_pos)
{
break;
}
@@ -468,14 +426,14 @@ void recover_UC_Read_sub_region_begin_end
{
long long offset = readLen - start_pos - 1;
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= end_pos && R_INF->N_site[ID][i] <= start_pos)
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)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)
else if((long long)R_INF->N_site[ID][i] > start_pos)
{
break;
}
@@ -529,13 +487,13 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
if (R_INF->N_site[ID])
{
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= start_pos && R_INF->N_site[ID][i] <= end_pos)
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)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)
else if((long long)R_INF->N_site[ID][i] > end_pos)
{
break;
}
@@ -575,14 +533,14 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
{
long long offset = readLen - start_pos - 1;
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if (R_INF->N_site[ID][i] >= end_pos && R_INF->N_site[ID][i] <= start_pos)
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)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)
else if((long long)R_INF->N_site[ID][i] > start_pos)
{
break;
}
@@ -609,7 +567,7 @@ void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID)
uint64_t i = 0;
while (i < r->length)
while ((long long)i < r->length)
{
memcpy(r->seq+i, bit_t_seq_table[src[i>>2]], 4);
i = i + 4;
@@ -660,7 +618,7 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
if (R_INF->N_site[ID])
{
for (i = 1; i <= R_INF->N_site[ID][0]; i++)
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N';
}
@@ -670,7 +628,7 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
#define COMPRESS_BASE {c = seq_nt6_table[src[i]];\
#define COMPRESS_BASE {c = seq_nt6_table[(uint8_t)src[i]];\
if (c >= 4)\
{\
c = 0;\
@@ -699,59 +657,31 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l
uint64_t dest_i = 0;
uint8_t tmp = 0;
uint8_t c = 0;
/**
fprintf(stderr, "src_l: %lld\n", src_l);
fflush(stderr);
**/
while (i + 4 <= src_l)
{
// fprintf(stderr, "0 i: %d, dest_i: %d, src_l: %d\n",
// i, dest_i, src_l);
// fflush(stderr);
tmp = 0;
COMPRESS_BASE;
tmp = tmp | (c<<6);
// fprintf(stderr, "*******1******1 i: %d, tmp: %d, c: %d\n",
// i, tmp, c);
// fflush(stderr);
COMPRESS_BASE;
tmp = tmp | (c<<4);
// fprintf(stderr, "*******2******1 i: %d, tmp: %d, c: %d\n",
// i, tmp, c);
// fflush(stderr);
COMPRESS_BASE;
tmp = tmp | (c<<2);
// fprintf(stderr, "*******3******1 i: %d, tmp: %d, c: %d\n",
// i, tmp, c);
// fflush(stderr);
COMPRESS_BASE;
tmp = tmp | c;
// fprintf(stderr, "*******4******1 i: %d, tmp: %d, c: %d\n",
// i, tmp, c);
// fflush(stderr);
dest[dest_i] = tmp;
// fprintf(stderr, "2 i: %d, dest_i: %d, src_l: %d\n",
// i, dest_i, src_l);
// fflush(stderr);
dest_i++;
}
//最多还剩3个字符
//at most 3 bases here
uint64_t shift = 6;
if (i < src_l)
{
@@ -773,6 +703,11 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l
void init_kseq(char* file)
{
fp = gzopen(file, "r");
if (fp == 0)
{
fprintf(stderr, "[ERROR] Cannot find the input file: %s\n", file);
exit(0);
}
seq = kseq_init(fp);
}
@@ -805,14 +740,14 @@ int get_read(kseq_t *s, int adapterLen)
if(adapterLen > 0)
{
if(s->seq.l <= adapterLen*2)
if((int)s->seq.l <= adapterLen*2)
{
s->seq.l = 0;
}
else
{
long long i;
for (i = 0; i < (s->seq.l - adapterLen*2); i++)
for (i = 0; i < ((int)s->seq.l - adapterLen*2); i++)
{
s->seq.s[i] = s->seq.s[i + adapterLen];
}
@@ -864,7 +799,6 @@ void init_R_buffer(int thread_num)
void destory_R_buffer_block(R_buffer_block* curr_sub_block)
{
kseq_destroy(curr_sub_block->read);
///free(curr_sub_block->read);
}
@@ -975,7 +909,6 @@ void* input_reads_muti_threads(void* arg)
total_reads = 0;
int i = 0;
int file_flag = 1;
R_buffer_block tmp_buf;
@@ -986,10 +919,7 @@ void* input_reads_muti_threads(void* arg)
while (1)
{
load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert, adapterLen);
load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert, asm_opt.adapterLen);
if (file_flag == 0)
{
@@ -1020,13 +950,11 @@ void* input_reads_muti_threads(void* arg)
destory_R_buffer_block(&tmp_buf);
fprintf(stdout, "total_reads: %llu\n",total_reads);
///fprintf(stdout, "R_INF.total_reads: %llu\n",R_INF.total_reads);
///fprintf(stdout, "R_INF.index[R_INF.total_reads]: %llu\n",R_INF.index[R_INF.total_reads]);
fprintf(stdout, "R_INF.total_reads_bases: %llu\n",R_INF.total_reads_bases);
///fprintf(stdout, "R_INF.name_index[R_INF.total_reads]: %llu\n",R_INF.name_index[R_INF.total_reads]);
fprintf(stdout, "R_INF.total_name_length: %llu\n",R_INF.total_name_length);
fprintf(stderr, "Reads #: %lu\n",total_reads);
fprintf(stderr, "Bases #: %lu\n",R_INF.total_reads_bases);
return NULL;
}
@@ -1078,7 +1006,7 @@ int get_reads_mul_thread(R_buffer_block* curr_sub_block)
void reverse_complement(char* pattern, uint64_t length)
{
int i = 0;
uint64_t i = 0;
uint64_t end = length / 2;
char k;
uint64_t index;
@@ -1099,44 +1027,3 @@ void reverse_complement(char* pattern, uint64_t length)
}
void Counting_block()
{
long long read_number = 0;
int i = 0;
int file_flag = 1;
R_buffer_block tmp_buf;
init_R_buffer_block(&tmp_buf);
while (1)
{
load_read_block(&tmp_buf, RDB.block_inner_size,
&file_flag, 0, adapterLen);
if (file_flag == 0)
{
break;
}
for (i = 0; i < tmp_buf.num; i++)
{
fprintf(stderr,"@%s\n", tmp_buf.read[i].name.s);
fprintf(stderr,"%s\n",tmp_buf.read[i].seq.s);
fprintf(stderr,"+\n");
fprintf(stderr,"%s\n",tmp_buf.read[i].qual.s);
read_number++;
}
}
fprintf(stdout, "read_number: %lld\n",read_number);
}
Process_Read.h
+6 -23
View File
@@ -30,29 +30,13 @@ KSEQ_INIT(gzFile, gzread)
static uint8_t seq_nt6_table[256] = {
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5
};
extern uint8_t seq_nt6_table[256];
extern char bit_t_seq_table[256][4];
extern char bit_t_seq_table_rc[256][4];
extern char s_H[5];
extern char rc_Table[5];
static char bit_t_seq_table[256][4] = {0};
static char bit_t_seq_table_rc[256][4] = {0};
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]]
@@ -204,7 +188,6 @@ void write_All_reads(All_reads* r, char* read_file_name);
int load_All_reads(All_reads* r, char* read_file_name);
void destory_All_reads(All_reads* r);
void Counting_block();
void destory_R_buffer_block(R_buffer_block* curr_sub_block);
void destory_R_buffer();
void clear_R_buffer();
debug_assembly.sh
-17
View File
@@ -1,17 +0,0 @@
#!/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
edlib.cpp
-1461
View File
File diff suppressed because it is too large Load Diff
edlib.h
-258
View File
@@ -1,258 +0,0 @@
#ifndef EDLIB_H
#define EDLIB_H
/**
* @file
* @author Martin Sosic
* @brief Main header file, containing all public functions and structures.
*/
#ifdef __cplusplus
extern "C" {
#endif
// Status codes
#define EDLIB_STATUS_OK 0
#define EDLIB_STATUS_ERROR 1
/**
* Alignment methods - how should Edlib treat gaps before and after query?
*/
typedef enum {
/**
* Global method. This is the standard method.
* Useful when you want to find out how similar is first sequence to second sequence.
*/
EDLIB_MODE_NW,
/**
* Prefix method. Similar to global method, but with a small twist - gap at query end is not penalized.
* What that means is that deleting elements from the end of second sequence is "free"!
* For example, if we had "AACT" and "AACTGGC", edit distance would be 0, because removing "GGC" from the end
* of second sequence is "free" and does not count into total edit distance. This method is appropriate
* when you want to find out how well first sequence fits at the beginning of second sequence.
*/
EDLIB_MODE_SHW,
/**
* Infix method. Similar as prefix method, but with one more twist - gaps at query end and start are
* not penalized. What that means is that deleting elements from the start and end of second sequence is "free"!
* For example, if we had ACT and CGACTGAC, edit distance would be 0, because removing CG from the start
* and GAC from the end of second sequence is "free" and does not count into total edit distance.
* This method is appropriate when you want to find out how well first sequence fits at any part of
* second sequence.
* For example, if your second sequence was a long text and your first sequence was a sentence from that text,
* but slightly scrambled, you could use this method to discover how scrambled it is and where it fits in
* that text. In bioinformatics, this method is appropriate for aligning read to a sequence.
*/
EDLIB_MODE_HW
} EdlibAlignMode;
/**
* Alignment tasks - what do you want Edlib to do?
*/
typedef enum {
EDLIB_TASK_DISTANCE, //!< Find edit distance and end locations.
EDLIB_TASK_LOC, //!< Find edit distance, end locations and start locations.
EDLIB_TASK_PATH //!< Find edit distance, end locations and start locations and alignment path.
} EdlibAlignTask;
/**
* Describes cigar format.
* @see http://samtools.github.io/hts-specs/SAMv1.pdf
* @see http://drive5.com/usearch/manual/cigar.html
*/
typedef enum {
EDLIB_CIGAR_STANDARD, //!< Match: 'M', Insertion: 'I', Deletion: 'D', Mismatch: 'M'.
EDLIB_CIGAR_EXTENDED //!< Match: '=', Insertion: 'I', Deletion: 'D', Mismatch: 'X'.
} EdlibCigarFormat;
// Edit operations.
#define EDLIB_EDOP_MATCH 0 //!< Match.
#define EDLIB_EDOP_INSERT 1 //!< Insertion to target = deletion from query.
#define EDLIB_EDOP_DELETE 2 //!< Deletion from target = insertion to query.
#define EDLIB_EDOP_MISMATCH 3 //!< Mismatch.
/**
* @brief Defines two given characters as equal.
*/
typedef struct {
char first;
char second;
} EdlibEqualityPair;
/**
* @brief Configuration object for edlibAlign() function.
*/
typedef struct {
/**
* Set k to non-negative value to tell edlib that edit distance is not larger than k.
* Smaller k can significantly improve speed of computation.
* If edit distance is larger than k, edlib will set edit distance to -1.
* Set k to negative value and edlib will internally auto-adjust k until score is found.
*/
int k;
/**
* Alignment method.
* EDLIB_MODE_NW: global (Needleman-Wunsch)
* EDLIB_MODE_SHW: prefix. Gap after query is not penalized.
* EDLIB_MODE_HW: infix. Gaps before and after query are not penalized.
*/
EdlibAlignMode mode;
/**
* Alignment task - tells Edlib what to calculate. Less to calculate, faster it is.
* EDLIB_TASK_DISTANCE - find edit distance and end locations of optimal alignment paths in target.
* EDLIB_TASK_LOC - find edit distance and start and end locations of optimal alignment paths in target.
* EDLIB_TASK_PATH - find edit distance, alignment path (and start and end locations of it in target).
*/
EdlibAlignTask task;
/**
* List of pairs of characters, where each pair defines two characters as equal.
* This way you can extend edlib's definition of equality (which is that each character is equal only
* to itself).
* This can be useful if you have some wildcard characters that should match multiple other characters,
* or e.g. if you want edlib to be case insensitive.
* Can be set to NULL if there are none.
*/
EdlibEqualityPair* additionalEqualities;
/**
* Number of additional equalities, which is non-negative number.
* 0 if there are none.
*/
int additionalEqualitiesLength;
} EdlibAlignConfig;
/**
* Helper method for easy construction of configuration object.
* @return Configuration object filled with given parameters.
*/
EdlibAlignConfig edlibNewAlignConfig(int k, EdlibAlignMode mode, EdlibAlignTask task,
EdlibEqualityPair* additionalEqualities,
int additionalEqualitiesLength);
/**
* @return Default configuration object, with following defaults:
* k = -1, mode = EDLIB_MODE_NW, task = EDLIB_TASK_DISTANCE, no additional equalities.
*/
EdlibAlignConfig edlibDefaultAlignConfig(void);
/**
* Container for results of alignment done by edlibAlign() function.
*/
typedef struct {
/**
* EDLIB_STATUS_OK or EDLIB_STATUS_ERROR. If error, all other fields will have undefined values.
*/
int status;
/**
* -1 if k is non-negative and edit distance is larger than k.
*/
int editDistance;
/**
* Array of zero-based positions in target where optimal alignment paths end.
* If gap after query is penalized, gap counts as part of query (NW), otherwise not.
* Set to NULL if edit distance is larger than k.
* If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free().
*/
int* endLocations;
/**
* Array of zero-based positions in target where optimal alignment paths start,
* they correspond to endLocations.
* If gap before query is penalized, gap counts as part of query (NW), otherwise not.
* Set to NULL if not calculated or if edit distance is larger than k.
* If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free().
*/
int* startLocations;
/**
* Number of end (and start) locations.
*/
int numLocations;
/**
* Alignment is found for first pair of start and end locations.
* Set to NULL if not calculated.
* Alignment is sequence of numbers: 0, 1, 2, 3.
* 0 stands for match.
* 1 stands for insertion to target.
* 2 stands for insertion to query.
* 3 stands for mismatch.
* Alignment aligns query to target from begining of query till end of query.
* If gaps are not penalized, they are not in alignment.
* If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free().
*/
unsigned char* alignment;
/**
* Length of alignment.
*/
int alignmentLength;
/**
* Number of different characters in query and target together.
*/
int alphabetLength;
} EdlibAlignResult;
/**
* Frees memory in EdlibAlignResult that was allocated by edlib.
* If you do not use it, make sure to free needed members manually using free().
*/
void edlibFreeAlignResult(EdlibAlignResult result);
/**
* Aligns two sequences (query and target) using edit distance (levenshtein distance).
* Through config parameter, this function supports different alignment methods (global, prefix, infix),
* as well as different modes of search (tasks).
* It always returns edit distance and end locations of optimal alignment in target.
* It optionally returns start locations of optimal alignment in target and alignment path,
* if you choose appropriate tasks.
* @param [in] query First sequence.
* @param [in] queryLength Number of characters in first sequence.
* @param [in] target Second sequence.
* @param [in] targetLength Number of characters in second sequence.
* @param [in] config Additional alignment parameters, like alignment method and wanted results.
* @return Result of alignment, which can contain edit distance, start and end locations and alignment path.
* Make sure to clean up the object using edlibFreeAlignResult() or by manually freeing needed members.
*/
EdlibAlignResult edlibAlign(const char* query, int queryLength,
const char* target, int targetLength,
const EdlibAlignConfig config);
/**
* Builds cigar string from given alignment sequence.
* @param [in] alignment Alignment sequence.
* 0 stands for match.
* 1 stands for insertion to target.
* 2 stands for insertion to query.
* 3 stands for mismatch.
* @param [in] alignmentLength
* @param [in] cigarFormat Cigar will be returned in specified format.
* @return Cigar string.
* I stands for insertion.
* D stands for deletion.
* X stands for mismatch. (used only in extended format)
* = stands for match. (used only in extended format)
* M stands for (mis)match. (used only in standard format)
* String is null terminated.
* Needed memory is allocated and given pointer is set to it.
* Do not forget to free it later using free()!
*/
char* edlibAlignmentToCigar(const unsigned char* alignment, int alignmentLength,
EdlibCigarFormat cigarFormat);
#ifdef __cplusplus
}
#endif
#endif // EDLIB_H
khash.h
+8 -7
View File
@@ -372,18 +372,19 @@ static const double __ac_HASH_UPPER = 0.77;
} \
SCOPE void kh_load_##name(kh_##name##_t *h, FILE* fp)\
{\
fread(&(h->n_buckets), sizeof(khint_t), 1, fp);\
fread(&(h->size), sizeof(khint_t), 1, fp);\
fread(&(h->n_occupied), sizeof(khint_t), 1, fp);\
fread(&(h->upper_bound), sizeof(khint_t), 1, fp);\
int f_flag;\
f_flag = fread(&(h->n_buckets), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->size), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->n_occupied), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->upper_bound), sizeof(khint_t), 1, fp);\
if (h->n_buckets)\
{\
h->flags = (khint32_t*)kmalloc(__ac_fsize(h->n_buckets) * sizeof(khint32_t));\
fread(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\
f_flag += fread(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\
h->keys = (khkey_t*)kmalloc(sizeof(khkey_t)*h->n_buckets);\
fread(h->keys, sizeof(khkey_t), h->n_buckets, fp);\
f_flag += fread(h->keys, sizeof(khkey_t), h->n_buckets, fp);\
h->vals = (khval_t*)kmalloc(sizeof(khval_t)*h->n_buckets);\
fread(h->vals, sizeof(khval_t), h->n_buckets, fp);\
f_flag += fread(h->vals, sizeof(khval_t), h->n_buckets, fp);\
}\
}
kmer.h
+2 -5
View File
@@ -15,11 +15,8 @@
typedef struct
{
///最大64-mer
///x[0]低位
///x[1]高位
//can represent at most 64-mer
uint64_t x[2];
} Hash_code;
typedef struct {
@@ -64,7 +61,7 @@ inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos)
for (; seq->i < seq->l; seq->i++)
{
///统计N的个数
///number of Ns
if (seq_nt6_table[(uint8_t)seq->str[seq->i]] >= 4)
{
seq->N_occ++;
ksw2.h
-177
View File
@@ -1,177 +0,0 @@
#ifndef KSW2_H_
#define KSW2_H_
#include <stdint.h>
#define KSW_NEG_INF -0x40000000
#define KSW_EZ_SCORE_ONLY 0x01 // don't record alignment path/cigar
#define KSW_EZ_RIGHT 0x02 // right-align gaps
#define KSW_EZ_GENERIC_SC 0x04 // without this flag: match/mismatch only; last symbol is a wildcard
#define KSW_EZ_APPROX_MAX 0x08 // approximate max; this is faster with sse
#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse
#define KSW_EZ_EXTZ_ONLY 0x40 // only perform extension
#define KSW_EZ_REV_CIGAR 0x80 // reverse CIGAR in the output
#define KSW_EZ_SPLICE_FOR 0x100
#define KSW_EZ_SPLICE_REV 0x200
#define KSW_EZ_SPLICE_FLANK 0x400
#ifdef __cplusplus
extern "C" {
#endif
typedef struct {
uint32_t max:31, zdropped:1;
int max_q, max_t; // max extension coordinate
int mqe, mqe_t; // max score when reaching the end of query
int mte, mte_q; // max score when reaching the end of target
int score; // max score reaching both ends; may be KSW_NEG_INF
int m_cigar, n_cigar;
int reach_end;
uint32_t *cigar;
} ksw_extz_t;
/**
* NW-like extension
*
* @param km memory pool, when used with kalloc
* @param qlen query length
* @param query query sequence with 0 <= query[i] < m
* @param tlen target length
* @param target target sequence with 0 <= target[i] < m
* @param m number of residue types
* @param mat m*m scoring mattrix in one-dimension array
* @param gapo gap open penalty; a gap of length l cost "-(gapo+l*gape)"
* @param gape gap extension penalty
* @param w band width (<0 to disable)
* @param zdrop off-diagonal drop-off to stop extension (positive; <0 to disable)
* @param flag flag (see KSW_EZ_* macros)
* @param ez (out) scores and cigar
*/
void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
/**
* Global alignment
*
* (first 10 parameters identical to ksw_extz_sse())
* @param m_cigar (modified) max CIGAR length; feed 0 if cigar==0
* @param n_cigar (out) number of CIGAR elements
* @param cigar (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, )
*
* @return score of the alignment
*/
int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat);
int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te);
#ifdef __cplusplus
}
#endif
/************************************
*** Private macros and functions ***
************************************/
#ifdef HAVE_KALLOC
#include "kalloc.h"
#else
#include <stdlib.h>
#define kmalloc(km, size) malloc((size))
#define kcalloc(km, count, size) calloc((count), (size))
#define krealloc(km, ptr, size) realloc((ptr), (size))
#define kfree(km, ptr) free((ptr))
#endif
static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t op, int len)
{
if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) {
if (*n_cigar == *m_cigar) {
*m_cigar = *m_cigar? (*m_cigar)<<1 : 4;
cigar = (uint32_t*)krealloc(km, cigar, (*m_cigar) << 2);
}
cigar[(*n_cigar)++] = len<<4 | op;
} else cigar[(*n_cigar)-1] += len<<4;
return cigar;
}
// In the backtrack matrix, value p[] has the following structure:
// bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
// bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
// bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
{ // p[] - lower 3 bits: which type gets the max; bit
int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
uint32_t *cigar = *cigar_, tmp;
while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check
int force_state = -1;
if (is_rot) {
r = i + j;
if (i < off[r]) force_state = 2;
if (off_end && i > off_end[r]) force_state = 1;
tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0;
} else {
if (j < off[i]) force_state = 2;
if (off_end && j > off_end[i]) force_state = 1;
tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0;
}
if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
if (force_state >= 0) state = force_state;
if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
}
if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
if (!is_rev)
for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;
*m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar;
}
static inline void ksw_reset_extz(ksw_extz_t *ez)
{
ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1;
ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF;
ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0;
}
static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e)
{
int r, t;
if (is_rot) r = a, t = b;
else r = a + b, t = a;
if (H > (int32_t)ez->max) {
ez->max = H, ez->max_t = t, ez->max_q = r - t;
} else if (t >= ez->max_t && r - t >= ez->max_q) {
int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l;
l = tl > ql? tl - ql : ql - tl;
if (zdrop >= 0 && ez->max - H > zdrop + l * e) {
ez->zdropped = 1;
return 1;
}
}
return 0;
}
#endif
ksw2_extz2_sse.cpp
-305
View File
@@ -1,305 +0,0 @@
#include <string.h>
#include <assert.h>
#include "ksw2.h"
#ifdef __SSE2__
#include <emmintrin.h>
#ifdef KSW_SSE2_ONLY
#undef __SSE4_1__
#endif
#ifdef __SSE4_1__
#include <smmintrin.h>
#endif
#ifdef KSW_CPU_DISPATCH
#ifdef __SSE4_1__
void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#else
void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#endif
#else
void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
#endif // ~KSW_CPU_DISPATCH
{
#define __dp_code_block1 \
z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \
xt1 = _mm_load_si128(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \
tmp = _mm_srli_si128(xt1, 15); /* tmp <- x[r-1][t+15] */ \
xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \
x1_ = tmp; \
vt1 = _mm_load_si128(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \
tmp = _mm_srli_si128(vt1, 15); /* tmp <- v[r-1][t+15] */ \
vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \
v1_ = tmp; \
a = _mm_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \
ut = _mm_load_si128(&u[t]); /* ut <- u[t..t+15] */ \
b = _mm_add_epi8(_mm_load_si128(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */
#define __dp_code_block2 \
z = _mm_max_epu8(z, b); /* z = max(z, b); this works because both are non-negative */ \
z = _mm_min_epu8(z, max_sc_); \
_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
_mm_store_si128(&v[t], _mm_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
z = _mm_sub_epi8(z, q_); \
a = _mm_sub_epi8(a, z); \
b = _mm_sub_epi8(b, z);
int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
int32_t *H = 0, H0 = 0, last_H0_t = 0;
uint8_t *qr, *sf, *mem, *mem2 = 0;
__m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_;
__m128i *u, *v, *x, *y, *s, *p = 0;
ksw_reset_extz(ez);
if (m <= 0 || qlen <= 0 || tlen <= 0) return;
zero_ = _mm_set1_epi8(0);
q_ = _mm_set1_epi8(q);
qe2_ = _mm_set1_epi8((q + e) * 2);
flag1_ = _mm_set1_epi8(1);
flag2_ = _mm_set1_epi8(2);
flag8_ = _mm_set1_epi8(0x08);
flag16_ = _mm_set1_epi8(0x10);
sc_mch_ = _mm_set1_epi8(mat[0]);
sc_mis_ = _mm_set1_epi8(mat[1]);
sc_N_ = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
m1_ = _mm_set1_epi8(m - 1); // wildcard
max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2);
if (w < 0) w = tlen > qlen? tlen : qlen;
wl = wr = w;
tlen_ = (tlen + 15) / 16;
n_col_ = qlen < tlen? qlen : tlen;
n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 15) / 16 + 1;
qlen_ = (qlen + 15) / 16;
for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) {
max_sc = max_sc > mat[t]? max_sc : mat[t];
min_sc = min_sc < mat[t]? min_sc : mat[t];
}
if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches
mem = (uint8_t*)kcalloc(km, tlen_ * 6 + qlen_ + 1, 16);
u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
if (!approx_max) {
H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
}
if (with_cigar) {
mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
off_end = off + qlen + tlen - 1;
}
for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
memcpy(sf, target, tlen);
for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
int st = 0, en = tlen - 1, st0, en0, st_, en_;
int8_t x1, v1;
uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
__m128i x1_, v1_;
// find the boundaries
if (st < r - qlen + 1) st = r - qlen + 1;
if (en > r) en = r;
if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil
if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor
if (st > en) {
ez->zdropped = 1;
break;
}
st0 = st, en0 = en;
st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1;
// set boundary conditions
if (st > 0) {
if (st - 1 >= last_st && st - 1 <= last_en)
x1 = ((uint8_t*)x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round
else x1 = v1 = 0; // not calculated; set to zeros
} else x1 = 0, v1 = r? q : 0;
if (en >= r) ((uint8_t*)y)[r] = 0, u8[r] = r? q : 0;
// loop fission: set scores first
if (!(flag & KSW_EZ_GENERIC_SC)) {
for (t = st0; t <= en0; t += 16) {
__m128i sq, st, tmp, mask;
sq = _mm_loadu_si128((__m128i*)&sf[t]);
st = _mm_loadu_si128((__m128i*)&qrr[t]);
mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_));
tmp = _mm_cmpeq_epi8(sq, st);
#ifdef __SSE4_1__
tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp);
tmp = _mm_blendv_epi8(tmp, sc_N_, mask);
#else
tmp = _mm_or_si128(_mm_andnot_si128(tmp, sc_mis_), _mm_and_si128(tmp, sc_mch_));
tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp), _mm_and_si128(mask, sc_N_));
#endif
_mm_storeu_si128((__m128i*)((uint8_t*)s + t), tmp);
}
} else {
for (t = st0; t <= en0; ++t)
((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]];
}
// core loop
x1_ = _mm_cvtsi32_si128(x1);
v1_ = _mm_cvtsi32_si128(v1);
st_ = st / 16, en_ = en / 16;
assert(en_ - st_ + 1 <= n_col_);
if (!with_cigar) { // score only
for (t = st_; t <= en_; ++t) {
__m128i z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
#endif
__dp_code_block2;
#ifdef __SSE4_1__
_mm_store_si128(&x[t], _mm_max_epi8(a, zero_));
_mm_store_si128(&y[t], _mm_max_epi8(b, zero_));
#else
tmp = _mm_cmpgt_epi8(a, zero_);
_mm_store_si128(&x[t], _mm_and_si128(a, tmp));
tmp = _mm_cmpgt_epi8(b, zero_);
_mm_store_si128(&y[t], _mm_and_si128(b, tmp));
#endif
}
} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
__m128i *pr = p + (size_t)r * n_col_ - st_;
off[r] = st, off_end[r] = en;
for (t = st_; t <= en_; ++t) {
__m128i d, z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
tmp = _mm_cmpgt_epi8(b, z);
d = _mm_blendv_epi8(d, flag2_, tmp); // d = b > z? 2 : d
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
tmp = _mm_cmpgt_epi8(b, z);
d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv
#endif
__dp_code_block2;
tmp = _mm_cmpgt_epi8(a, zero_);
_mm_store_si128(&x[t], _mm_and_si128(tmp, a));
d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_)); // d = a > 0? 0x08 : 0
tmp = _mm_cmpgt_epi8(b, zero_);
_mm_store_si128(&y[t], _mm_and_si128(tmp, b));
d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0
_mm_store_si128(&pr[t], d);
}
} else { // gap right-alignment
__m128i *pr = p + (size_t)r * n_col_ - st_;
off[r] = st, off_end[r] = en;
for (t = st_; t <= en_; ++t) {
__m128i d, z, a, b, xt1, vt1, ut, tmp;
__dp_code_block1;
d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1
#ifdef __SSE4_1__
z = _mm_max_epi8(z, a); // z = z > a? z : a (signed)
tmp = _mm_cmpgt_epi8(z, b);
d = _mm_blendv_epi8(flag2_, d, tmp); // d = z > b? d : 2
#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8()
z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0;
z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative
tmp = _mm_cmpgt_epi8(z, b);
d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv
#endif
__dp_code_block2;
tmp = _mm_cmpgt_epi8(zero_, a);
_mm_store_si128(&x[t], _mm_andnot_si128(tmp, a));
d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_)); // d = 0 > a? 0 : 0x08
tmp = _mm_cmpgt_epi8(zero_, b);
_mm_store_si128(&y[t], _mm_andnot_si128(tmp, b));
d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10
_mm_store_si128(&pr[t], d);
}
}
if (!approx_max) { // find the exact max with a 32-bit score array
int32_t max_H, max_t;
// compute H[], max_H and max_t
if (r > 0) {
int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i;
__m128i max_H_, max_t_, qe_;
max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element
max_t = en0;
max_H_ = _mm_set1_epi32(max_H);
max_t_ = _mm_set1_epi32(max_t);
qe_ = _mm_set1_epi32(q + e);
for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t;
__m128i H1, tmp, t_;
H1 = _mm_loadu_si128((__m128i*)&H[t]);
t_ = _mm_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3]);
H1 = _mm_add_epi32(H1, t_);
H1 = _mm_sub_epi32(H1, qe_);
_mm_storeu_si128((__m128i*)&H[t], H1);
t_ = _mm_set1_epi32(t);
tmp = _mm_cmpgt_epi32(H1, max_H_);
#ifdef __SSE4_1__
max_H_ = _mm_blendv_epi8(max_H_, H1, tmp);
max_t_ = _mm_blendv_epi8(max_t_, t_, tmp);
#else
max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_));
max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_));
#endif
}
_mm_storeu_si128((__m128i*)HH, max_H_);
_mm_storeu_si128((__m128i*)tt, max_t_);
for (i = 0; i < 4; ++i)
if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i;
for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE
H[t] += (int32_t)v8[t] - qe;
if (H[t] > max_H)
max_H = H[t], max_t = t;
}
} else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0
// update ez
if (en0 == tlen - 1 && H[en0] > ez->mte)
ez->mte = H[en0], ez->mte_q = r - en;
if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
ez->mqe = H[st0], ez->mqe_t = st0;
if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break;
if (r == qlen + tlen - 2 && en0 == tlen - 1)
ez->score = H[tlen - 1];
} else { // find approximate max; Z-drop might be inaccurate, too.
if (r > 0) {
if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) {
int32_t d0 = v8[last_H0_t] - qe;
int32_t d1 = u8[last_H0_t + 1] - qe;
if (d0 > d1) H0 += d0;
else H0 += d1, ++last_H0_t;
} else if (last_H0_t >= st0 && last_H0_t <= en0) {
H0 += v8[last_H0_t] - qe;
} else {
++last_H0_t, H0 += u8[last_H0_t] - qe;
}
if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) break;
} else H0 = v8[0] - qe - qe, last_H0_t = 0;
if (r == qlen + tlen - 2 && en0 == tlen - 1)
ez->score = H0;
}
last_st = st, last_en = en;
//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
}
kfree(km, mem);
if (!approx_max) kfree(km, H);
if (with_cigar) { // backtrack
int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
ez->reach_end = 1;
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
} else if (ez->max_t >= 0 && ez->max_q >= 0) {
ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
}
kfree(km, mem2); kfree(km, off);
}
}
#endif // __SSE2__
main.cpp
+3 -213
View File
@@ -4,224 +4,14 @@
#include "Process_Read.h"
#include "Assembly.h"
#include "Levenshtein_distance.h"
#include "edlib.h"
/********************************for debug***************************************/
///使用这个函数的时候，必须把Counting_multiple_thr()里的destory_Total_Count_Table(&TCB)注释掉
void debug_Counting()
{
init_kseq(read_file_name);
Verify_Counting();
fprintf(stderr, "debug over!\n");
destory_kseq();
}
int matrix[1000][1000] = {0};
///y_length > x_length
int edit_distance_normal(char* y, int y_length, char* x, int x_length)
{ memset(matrix, 0, sizeof(matrix));
int i, j;
for (i = 0; i <= x_length; i++)
{
matrix[i][0] = i;
}
int digonal, up, left, min;
///一列列算的
for (i = 0; i < x_length; i++)
{
for (j = 0; j < y_length; j++)
{
///matrix[i + 1][j + 1]
digonal = matrix[i][j] + (x[i] != y[j]);
up = matrix[i + 1][j] + 1;
left = matrix[i][j + 1] + 1;
min = digonal;
if (up < min)
{
min = up;
}
if (left< min)
{
min = left;
}
matrix[i + 1][j + 1] = min;
}
}
min = 999999;
for (j = 0; j <= y_length; j++)
{
if (matrix[i][j] < min)
{
min = matrix[i][j];
}
}
return min;
}
///y_length > x_length
int edit_distance_normal_banded(char* y, int y_length, char* x, int x_length, int error)
{ memset(matrix, 0, sizeof(matrix));
int i, j;
for (i = 0; i <= x_length; i++)
{
for (j = 0; j <= y_length; j++)
{
matrix[i][j] = 1000000;
}
}
for (i = 0; i <= x_length; i++)
{
matrix[i][0] = i;
}
for (i = 0; i <= y_length; i++)
{
matrix[0][i] = 0;
}
int banded_length = error*2 + 1;
int digonal, up, left, min;
for (i = 0; i < x_length; i++)
{
///for (j = 0; j < y_length; j++)
for (j = i; j < banded_length + i; j++)
{
///matrix[i + 1][j + 1]
digonal = matrix[i][j] + (x[i] != y[j]);
up = matrix[i + 1][j] + 1;
left = matrix[i][j + 1] + 1;
min = digonal;
if (up < min)
{
min = up;
}
if (left< min)
{
min = left;
}
matrix[i + 1][j + 1] = min;
}
}
min = 999999;
for (j = 0; j <= y_length; j++)
///for (j = i; j < banded_length + i; j++)
{
if (matrix[i][j] < min)
{
min = matrix[i][j];
}
}
return min;
}
void debug_edit_distance()
{
/**
char* x = "TTCCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCACTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCAGTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCGATTCCATTCAAATCAATTCCGTTCTATTCAATTTGATTCAT";
char* y = "CCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCATTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCACTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCCATTCCATTCAATTCAATTCCGTGCTATTCAATTTGATTCATTTCCATTTAATTCCATTCCATTAGATTCCATT";
**/
unsigned short toold = 15;
char* x
= "GAAAGAGAATCAAATGGAATTGAATCGAATGGAATCGAATGGATTGGAAAGGAATAGAATGGAATGGAATGGAATTGACTCAAATGGAATGGACTAGAATGGAATGGATTCGAATGGAAGGCAAAGGAATGGAATCTATCGGAATGGACTGTAATGGAATGGAATGGAAGGGATTGGAATGGATTCGAATGTAATGGACTGCAATAGAAAGGATTCGAATGGAATGAAAAAGAATTGAATGGAATAGAACAGAATGGAATCAAATCGAAGGAAATGGAATGGAATAGAAAGGAATGGAATGAAATGGAATGGAAAGGATTCGAATGGAATGCAATCGAATGGAATGGAATCGAACGGAATGGAATAAAATGGAAG";
char* y =
"GAAAGAGAATCAAATGGAATTGAATCGAATGGAATCGAATGGATTGGAAAGGAATAGAATGGAATGGAATGGAATGGACTCAAATGGAATGTACTAGAATGGAATGGATTCGAATGGAAGGCAAAGGAATGGAATCTATTGGAATGGACTGTAATGGAATGGAATGGAAGGGATTGGAATGGACTCGAATGGAATGGACTGCAATAGAAAGGATTCGAATGGAATGAAAAAGAATTGAATGGAATAGAACAGAATGGAATCAAATCGAATGAAATGGAATGGAATAGAAAGGAATGGAATGAAATGGAATGGAAAGGATTCGAATGGAATGCAATCGAATGGAATGGAATCGAACGGAATGGAATAAATTTTCTG";
fprintf(stderr, "x_length: %u\n", strlen(x));
fprintf(stderr, "y_length: %u\n", strlen(y));
EdlibAlignResult result = edlibAlign(x, strlen(x), y, strlen(y),
edlibNewAlignConfig(toold, EDLIB_MODE_HW, EDLIB_TASK_PATH, NULL, 0));
if (result.status == EDLIB_STATUS_OK) {
fprintf(stderr, "****\nedlib: %d, alignmentLength: %d, startLocations: %d, endLocations: %d\n",
result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0]);
char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD);
fprintf(stderr,"%s\n", cigar);
free(cigar);
}
edlibFreeAlignResult(result);
unsigned int error;
int end_site = Reserve_Banded_BPM(y, strlen(y), x, strlen(x), toold, &error);
fprintf(stderr, "BPM: error: %u, end_site: %u\n", error, end_site);
unsigned short band_length=(toold+1)*3-1-1-toold;
unsigned short band_down=toold-1;
unsigned short band_blew=2*(toold+1)-1-1;
int return_err = 99999;
///注意pattern/text和band_down/band_blew是反的
Reserve_Banded_BPM_new(y, strlen(y), x, strlen(x),
toold,band_blew,band_down,band_length, &return_err, 0);
fprintf(stderr, "new BPM: error: %u\n", return_err);
return_err = edit_distance_normal(y, strlen(y), x, strlen(x));
fprintf(stderr, "edit_distance_normal: error: %u\n", return_err);
return_err = edit_distance_normal_banded(y, strlen(y), x, strlen(x), toold);
fprintf(stderr, "edit_distance_normal_banded: error: %u\n", return_err);
end_site = Reserve_Banded_BPM_debug(y, strlen(y), x, strlen(x), toold, &error, matrix);
fprintf(stderr, "BPM debug: error: %u, end_site: %u\n", error, end_site);
}
int main(int argc, char *argv[])
{
init_opt(&asm_opt);
if (!CommandLine_process(argc, argv))
return 1;
if (!CommandLine_process(argc, argv, &asm_opt)) return 1;
fprintf(stdout, "Will perform %d round of error correction...\n", number_of_round);
Correct_Reads(asm_opt.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);
fprintf(stdout, "k-mer length: %d\n",k_mer_length);
fprintf(stdout, "coverage: %d\n",coverage);
fprintf(stdout, "read_graph: %d\n", read_graph);
fprintf(stdout, "adapterLen: %d\n", adapterLen);
fflush(stdout);
if(required_read_name)
{
fprintf(stdout, "required_read_name: %s\n", required_read_name);
}
Correct_Reads(number_of_round);
return 1;
}