forked from weizhongli/cdhit
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cdhit-est-2d.c++
115 lines (93 loc) · 3.79 KB
/
cdhit-est-2d.c++
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
// =============================================================================
// CD-HI-EST
// http://cd-hit.org/
// Cluster Database at High Identity (EST version)
// modified from CD-HI
//
// program written by
// Weizhong Li
// UCSD, San Diego Supercomputer Center
// La Jolla, CA, 92093
// Email liwz@sdsc.edu
// at
// Adam Godzik's lab
// The Burnham Institute
// La Jolla, CA, 92037
// Email adam@burnham-inst.org
//
// Modified by:
// Limin Fu
// Center for Research in Biological Systems (CRBS), UCSD
// La Jolla, CA, 92093
// Email: l2fu@ucsd.edu, fu@daovm.net
// =============================================================================
#include "cdhit-common.h"
#include "cdhit-utility.h"
//over-write some defs in cd-hi.h for est version
#undef MAX_UAA
#define MAX_UAA 4
//over-write some defs in cd-hi-init.h for est version
void setaa_to_na();
void make_comp_iseq(int len, char *iseq_comp, char *iseq);
void make_comp_short_word_index(int NAA, int *NAAN_array, Vector<int> & Comp_AAN_idx);
Options options;
SequenceDB seq_db;
SequenceDB seq_db2;
// next two control how if seqs in db2 is longer than reps in db1
// by deault, only seqs in db2 that are shorter than rep in db1
// are clustered to the rep in db1
//////////////////////////////////// MAIN /////////////////////////////////////
int main(int argc, char **argv)
{
string db_in;
string db_in2;
string db_out;
string db_in_pe;
string db_in2_pe;
string db_out_pe;
options.cluster_thd = 0.95;
options.NAA = 10;
options.NAAN = NAA8;
seq_db.NAAN = NAA8;
options.NAA_top_limit = 12;
setaa_to_na();
mat.set_to_na(); //mat.set_gap(-6,-1);
float begin_time = current_time();
float end_time;
// *********************************** parse command line and open file
if (argc < 7) print_usage_est_2d(argv[0]);
if (options.SetOptions( argc, argv, true, true ) == 0) print_usage_est_2d(argv[0]);
options.Validate();
db_in = options.input;
db_in_pe = options.input_pe;
db_in2 = options.input2;
db_in2_pe = options.input2_pe;
db_out = options.output;
db_out_pe = options.output_pe;
InitNAA( MAX_UAA );
options.NAAN = NAAN_array[options.NAA];
seq_db.NAAN = NAAN_array[options.NAA];
seq_db2.NAAN = NAAN_array[options.NAA];
if ( options.option_r ) {
Comp_AAN_idx.resize( seq_db.NAAN );
make_comp_short_word_index(options.NAA, NAAN_array, Comp_AAN_idx);
}
if ( options.PE_mode ) {seq_db.Read( db_in.c_str(), db_in_pe.c_str(), options );}
else {seq_db.Read( db_in.c_str(), options );}
cout << "total seq in db1: " << seq_db.sequences.size() << endl;
if ( options.PE_mode ) { seq_db2.Read( db_in2.c_str(), db_in2_pe.c_str(), options );}
else { seq_db2.Read( db_in2.c_str(), options );}
cout << "total seq in db2: " << seq_db2.sequences.size() << endl;
seq_db.SortDivide( options );
seq_db2.SortDivide( options, false );
seq_db2.ClusterTo( seq_db, options );
cout << "writing non-redundant sequences from db2" << endl;
seq_db2.WriteClusters( db_in2.c_str(), db_out.c_str(), options );
if ( options.PE_mode ) { seq_db2.WriteClusters( db_in2.c_str(), db_in2_pe.c_str(), db_out.c_str(), db_out_pe.c_str(), options ); }
else { seq_db2.WriteClusters( db_in2.c_str(), db_out.c_str(), options ); }
seq_db2.WriteExtra2D( seq_db, options );
cout << "program completed !" << endl << endl;
end_time = current_time();
printf( "Total CPU time %.2f\n", end_time - begin_time );
return 0;
} // END int main