-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathmTR.h
180 lines (153 loc) · 7.59 KB
/
mTR.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
/*
Copyright (c) 2019, Shinichi Morishita
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of the FreeBSD Project.
*/
// Key default parameters
#define MAX_INPUT_LENGTH 1000000 // The maximum length of each read
#define MIN_MATCH_RATIO 0.6 // The minimum threshold of match ratio between the estimated repeat unit and the repeat in a given raw read
#define MIN_PERIOD 2 // Minimum period length
#define MAX_PERIOD 500 // Maximum period length
#define MIN_NUM_FREQ_UNIT 5 // The minimum threshold of number of units
#define ALIGNMENT_WIDTH_PRINTING 50
#define MAX_LEN_overlapping 10 // = MIN_PERIOD * MIN_NUM_FREQ_UNIT
// The following values are optimzed for a benchmark dataset.
#define MIN_WINDOW 5
#define MAX_WINDOW 10240
#define minKmer 5
#define maxKmer 15 // Increase this when no qualified repeats are found.
#define MAX_tiebreaks 1024
#define MIN_jaccard_index 0.98
#define BLK 4096 // Block size of input buffer.
#define WrapDPsize 200000000 // 200M = repeat_unit_size (200) x length_of_repeats (1,000,000)
// Choice of DeBruijn graph or progressive multiple alignment
#define ProgressiveMultipleAlignment 0
#define DeBruijnGraphSearch 1
#define count_maxKmer 6
#define PrimeMax 256019
// Global variables in the heap
int Manhattan_Distance; // The default setting is 1 in the main.
float min_match_ratio;
int *pow4; // A table of pow(4, k)
int *min_coverage_for_missing_bases;
int *orgInputString; // An input read string of length MAX_INPUT_LENGTH.
int *inputString; // 4 decimal encoding of the input read string of length MAX_INPUT_LENGTH.
int *inputString_w_rand; // 4 decimal encoding of the input read string of length MAX_INPUT_LENGTH.
int *count; // A table of size 4^k for counting sort.
int **freqNode; // A hash table for storing the frequency of each node (node, frequency)
int *sortedString; // Positions of 4-mers sorted wrt both 4-mers and their positions.
double *directional_index_tmp; // For storing all DI values temporarily for a given w
double *directional_index; // For storing a locally maximum DI
int *directional_index_end; // The end position of the local maximum
int *directional_index_w; // The window width w for the local maximum
int *vector0, *vector1, *vector2;
double *freq_interval_len; // For computing frequency distribution of interval lengths
// The length is MAX_INPUT_LENGTH
int *WrapDP; // 2D space for Wrap-around global alignment DP for handling tandem repeats
char *alignment_input;
char *alignment_symbols;
char *alignment_repeats;
// For printing the alignment of the predicted repeat unit with the input string
// The largest array, and the size is (MAX_PERIOD+1) * (MAX_INPUT_LENGTH+1)
int **consensus, **gaps; // Space for consensus
#define MAX_ID_LENGTH 1000
typedef struct{
int len;
int codedString[MAX_INPUT_LENGTH];
char ID[MAX_ID_LENGTH];
} Read;
int tmp_read_cnt;
char *nextReadID; // [MAX_ID_LENGTH];
typedef struct { // MAX_ID_LENGTH + MAX_EPRIOD + 28*4 = 612 bytes
int ID; // 0,1,2,...
char readID[BLK];
int inputLen;
int rep_start;
int rep_end;
int repeat_len;
int rep_period;
int Num_freq_unit;
int Num_matches;
int Num_mismatches;
int Num_insertions;
int Num_deletions;
int Kmer;
int match_gain;
int mismatch_penalty;
int indel_penalty;
char string[MAX_PERIOD];
int string_score[MAX_PERIOD];
int freq_2mer[16];
} repeat_in_read;
// External functions
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define DIFF(x, y) ((x) > (y) ? ((x) - (y)) : ((y) - (x)))
int handle_one_file(char *inputFile, int print_alignment);
void handle_one_read( char *readID, int inputLen, int read_cnt, int print_alignment);
void fill_directional_index_with_end(int DI_array_length, int inputLen, int random_string_length);
void init_inputString(int k, int query_start, int query_end, int inputLen);
int search_De_Bruijn_graph(int query_start, int query_end, repeat_in_read *rr);
void revise_representative_unit( repeat_in_read *rr);
void wrap_around_DP(int query_start, int query_end, repeat_in_read *rr);
void wrap_around_DP_sub( int query_start, int query_end, repeat_in_read *rr, int MATCH_GAIN, int MISMATCH_PENALTY, int INDEL_PENALTY );
void set_rr(repeat_in_read *rr_a, repeat_in_read *rr_b);
void clear_rr(repeat_in_read *rr_a);
void freq_2mer_array(int* val, int len, int *freq_2mer);
void print_one_repeat_in_read(repeat_in_read rr);
void print_freq(int rep_start, int rep_end, int rep_period, char* string, int inputLen, int k);
float time_all, time_memory, time_range, time_period, time_initialize_input_string, time_wrap_around_DP, time_count_table, time_chaining;
int query_counter;
// Interface between C and C++ functions
#ifndef __CSUB_H__
#define __CSUB_H__ 1
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
void insert_an_alignment_into_set(
char* readID,
int inputLen,
int rep_start,
int rep_end,
int repeat_len,
int rep_period,
int Num_freq_unit,
int Num_matches,
int Num_mismatches,
int Num_insertions,
int Num_deletions,
int Kmer,
int match_gain,
int mismatch_penalty,
int indel_penalty,
char* string,
int* string_score);
void chaining(int print_alignment);
extern void pretty_print_alignment(char *unit_string, int unit_len, int rep_start, int rep_end, int match_gain, int mismatch_penalty, int indel_penalty);
extern void print_freq(int rep_start, int rep_end, int rep_period, char* string, int inputLen, int k);
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif /* __CSUB_H__ */
//For debugging with #ifdef
//#define DEBUG
//#define DEBUG_match_gain