forked from refresh-bio/KMC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathkmc_query.cpp
138 lines (126 loc) · 3.6 KB
/
kmc_query.cpp
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
/*
* =====================================================================================
*
* Filename: kmer_query.cc
*
* Description:
*
* Version: 1.0
* Created: 04/27/2016 08:48:26 AM
* Revision: none
* Compiler: gcc
*
* Author: Prashant Pandey (), [email protected]
* Organization: Stony Brook University
*
* =====================================================================================
*/
#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>
#include <set>
#include <bitset>
#include <cassert>
#include <fstream>
#include <time.h>
#include <stdio.h>
#include <fcntl.h>
#include <unistd.h>
#include <sys/resource.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <sys/mman.h>
#include "kmc_api/kmc_file.h"
using namespace std;
/* Print elapsed time using the start and end timeval */
void print_time_elapsed(string desc, struct timeval* start, struct timeval* end)
{
struct timeval elapsed;
if (start->tv_usec > end->tv_usec) {
end->tv_usec += 1000000;
end->tv_sec--;
}
elapsed.tv_usec = end->tv_usec - start->tv_usec;
elapsed.tv_sec = end->tv_sec - start->tv_sec;
float time_elapsed = (elapsed.tv_sec * 1000000 + elapsed.tv_usec)/1000000.f;
cout << desc << "Total Time Elapsed: " << to_string(time_elapsed) << " seconds" << endl;
}
//void getRandomKmers(int n, vector<CKmerAPI>& kmers)
//{
//for (int j = 0; j < n; j++) {
//char kmer[28] = {0};
//for (int i = 0; i < 28; i++) {
//kmer[i] = bases[rand()/4];
//}
//CKmerAPI kmc_kmer(28);
//kmers.push_back(kmc_kmer.from_string(kmer));
//}
//}
/*
* === FUNCTION ======================================================================
* Name: main
* Description:
* =====================================================================================
*/
int main ( int argc, char *argv[] )
{
string filename = string(argv[1]);
int num_query = atoi(argv[2]);
int random = atoi(argv[3]);
int dump = atoi(argv[4]);
CKMCFile kmer_database_list;
CKMCFile kmer_database_rand;
vector<CKmerAPI> kmer_objects;
CKmerAPI kmer(28);
uint32_t counter;
struct timeval start, end;
struct timezone tzp;
char bases[] = {'C', 'A', 'T', 'G'}; // A=1, C=0, T=2, G=3
ofstream dumpfile;
dumpfile.open("RandomKmer.dump");
if (!kmer_database_list.OpenForListing(filename)) {
cout << "Can not open the database file" << endl;
return EXIT_FAILURE;
}
cout << "Reading kmers from the database list" << endl;
if (random) {
for (int j = 0; j < num_query; j++) {
string kmer;
for (int i = 0; i < 28; i++) {
kmer += bases[rand()%4];
}
if (dump)
dumpfile << kmer << endl;
CKmerAPI kmc_kmer(28);
kmer_objects.push_back(kmc_kmer.from_string(kmer));
}
} else {
uint64_t i = 0;
while (kmer_database_list.ReadNextKmer(kmer, counter)) {
i++;
kmer_objects.push_back(kmer);
}
cout << "Total kmers: " << i << endl;
kmer_database_list.Close();
}
if (!kmer_database_rand.OpenForRA(filename)) {
cout << "Can not open the database file" << endl;
return EXIT_FAILURE;
}
srand(time(NULL));
cout << "Querying kmers in the KMC database rand" << endl;
gettimeofday(&start, &tzp);
for (int i = 0; i < num_query; i++) {
uint64_t id = rand()%kmer_objects.size();
/*cout << "index: " << id << endl;*/
if (!kmer_database_rand.CheckKmer(kmer_objects[id], counter)) {
cout << "Can not find the kmer: " << kmer_objects[id].to_string() << endl;
abort();
}
}
gettimeofday(&end, &tzp);
print_time_elapsed("", &start, &end);
kmer_database_rand.Close();
return EXIT_SUCCESS;
} /* ---------- end of function main ---------- */