36 #include <cuda_runtime_api.h>
37 #include <vector_types.h>
38 #include <vector_functions.h>
47 using namespace nvbio;
48 using namespace alndiff;
51 int main(
int argc,
char* argv[])
57 cudaGetDeviceCount(&device_count);
58 log_verbose(stderr,
" cuda devices : %d\n", device_count);
60 const char* report_name = NULL;
62 const char* filter_name = NULL;
65 int32 filter_delta = 5;
72 if (strcmp( argv[arg],
"-device" ) == 0)
74 cuda_device = atoi(argv[++arg]);
77 else if (strcmp( argv[arg],
"-paired" ) == 0)
82 else if (strcmp( argv[arg],
"-report" ) == 0)
84 report_name = argv[++arg];
87 else if (strcmp( argv[arg],
"-no-ids" ) == 0)
92 else if (strcmp( argv[arg],
"-filter" ) == 0)
94 filter_name = argv[++arg];
96 const std::string flags_string( argv[++arg] );
98 filter_delta = atoi(argv[++arg]);
104 const char*
begin = flags_string.c_str();
105 const char* end =
begin;
110 while (*end !=
'|' && *end !=
'\0')
112 temp[end -
begin] = *end;
116 temp[end -
begin] =
'\0';
118 if (strcmp( temp,
"distant" ) == 0)
120 else if (strcmp( temp,
"discordant" ) == 0)
122 else if (strcmp( temp,
"diff-ref" ) == 0)
132 begin = stats_string.c_str();
138 while (*end !=
'|' && *end !=
'\0')
140 temp[end -
begin] = *end;
144 temp[end -
begin] =
'\0';
146 if (strcmp( temp,
"ed" ) == 0)
148 else if (strcmp( temp,
"mapQ" ) == 0)
150 else if (strcmp( temp,
"mms" ) == 0)
152 else if (strcmp( temp,
"ins" ) == 0)
154 else if (strcmp( temp,
"dels" ) == 0)
170 if (cuda_device == -1)
173 cudaDeviceProp best_device_prop;
174 cudaGetDeviceProperties( &best_device_prop, best_device );
176 for (
int device = 0; device < device_count; ++device)
178 cudaDeviceProp device_prop;
179 cudaGetDeviceProperties( &device_prop, device );
180 log_verbose(stderr,
" device %d has compute capability %d.%d\n", device, device_prop.major, device_prop.minor);
181 log_verbose(stderr,
" SM count : %u\n", device_prop.multiProcessorCount);
182 log_verbose(stderr,
" SM clock rate : %u Mhz\n", device_prop.clockRate / 1000);
183 log_verbose(stderr,
" memory clock rate : %.1f Ghz\n",
float(device_prop.memoryClockRate) * 1.0e-6f);
185 if (device_prop.major >= best_device_prop.major &&
186 device_prop.minor >= best_device_prop.minor)
188 best_device_prop = device_prop;
189 best_device = device;
192 cuda_device = best_device;
194 log_verbose(stderr,
" chosen device %d\n", cuda_device);
196 cudaDeviceProp device_prop;
197 cudaGetDeviceProperties( &device_prop, cuda_device );
198 log_verbose(stderr,
" device name : %s\n", device_prop.name);
199 log_verbose(stderr,
" compute capability : %d.%d\n", device_prop.major, device_prop.minor);
201 cudaSetDevice( cuda_device );
206 log_info(stderr,
"nvbio-aln-diff [OPTIONS] <file1> <file2>\n");
208 log_info(stderr,
" -paired # paired-end input\n" );
209 log_info(stderr,
" -no-ids # do not perform id checks\n" );
210 log_info(stderr,
" -report <file-name> # HTML report\n" );
211 log_info(stderr,
" -filter <file-name>\n" );
212 log_info(stderr,
" <flags={distant|discordant|diff-ref}>\n" );
213 log_info(stderr,
" <stats={ed|mapQ|mms|ins|dels}>\n" );
220 std::string filter_flags_string;
222 filter_flags_string += std::string(filter_flags_string.length() ?
"|" :
"") +
"distant";
224 filter_flags_string += std::string(filter_flags_string.length() ?
"|" :
"") +
"discordant";
226 filter_flags_string += std::string(filter_flags_string.length() ?
"|" :
"") +
"diff-ref";
228 std::string filter_stats_string;
230 filter_stats_string += std::string(filter_stats_string.length() ?
"|" :
"") +
"ed";
232 filter_stats_string += std::string(filter_stats_string.length() ?
"|" :
"") +
"mapQ";
234 filter_stats_string += std::string(filter_stats_string.length() ?
"|" :
"") +
"mms";
236 filter_stats_string += std::string(filter_stats_string.length() ?
"|" :
"") +
"ins";
238 filter_stats_string += std::string(filter_stats_string.length() ?
"|" :
"") +
"dels";
240 log_info(stderr,
"filter file : %s\n", filter_name);
241 log_info(stderr,
" flags : %s\n", filter_flags_string.c_str());
242 log_info(stderr,
" stats : %s\n", filter_stats_string.c_str());
243 log_info(stderr,
" delta : %d\n", filter_delta);
248 const char *aln_file_nameL = argv[arg];
249 const char *aln_file_nameR = argv[arg+1];
253 const char *aln_file_nameL = argv[arg];
254 const char *aln_file_nameR = argv[arg+1];
259 if (aln_streamL == NULL || aln_streamL->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameL); exit(1); }
260 if (aln_streamR == NULL || aln_streamR->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameR); exit(1); }
262 const uint32 BATCH_SIZE = 500000;
263 std::vector<Alignment> batchL( BATCH_SIZE );
264 std::vector<Alignment> batchR( BATCH_SIZE );
266 Filter filter( filter_name, filter_flags, filter_stats, filter_delta );
272 const uint32 batch_sizeL = aln_streamL->
next_batch( BATCH_SIZE, &batchL[0] );
273 const uint32 batch_sizeR = aln_streamR->
next_batch( BATCH_SIZE, &batchR[0] );
278 if (min_batch_size != max_batch_size)
279 log_warning(stderr,
"alignment files have different size\n");
281 log_info(stderr,
"analizing batch[%u]: %u alignments (%.1f M)\n", n_batch, min_batch_size,
float(n_batch * BATCH_SIZE + min_batch_size)*1.0e-6f);
283 for (
uint32 i = 0; i < min_batch_size; i += 2)
292 log_error(stderr,
"alignments %u and %u in \"%s\" refer to the same mate, must come from different reads\n", i, i+1, aln_file_nameL);
297 log_error(stderr,
"alignments %u and %u in \"%s\" refer to the same mate, must come from different reads\n", i, i+1, aln_file_nameR);
310 analyzer->
generate_report( aln_file_nameL, aln_file_nameR, report_name );
330 if (min_batch_size < BATCH_SIZE)
341 if (aln_streamL == NULL || aln_streamL->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameL); exit(1); }
342 if (aln_streamR == NULL || aln_streamR->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameR); exit(1); }
344 const uint32 BATCH_SIZE = 500000;
345 std::vector<Alignment> batchL( BATCH_SIZE );
346 std::vector<Alignment> batchR( BATCH_SIZE );
348 Filter filter( filter_name, filter_flags, filter_stats, filter_delta );
354 const uint32 batch_sizeL = aln_streamL->
next_batch( BATCH_SIZE, &batchL[0] );
355 const uint32 batch_sizeR = aln_streamR->
next_batch( BATCH_SIZE, &batchR[0] );
357 if (batch_sizeL != batch_sizeR)
358 log_warning(stderr,
"alignment files have different size\n");
362 log_info(stderr,
"analizing batch[%u]: %u alignments (%.1f M)\n", n_batch, min_batch_size,
float(n_batch * BATCH_SIZE + min_batch_size)*1.0e-6f);
364 for (
uint32 i = 0; i < min_batch_size; ++i)
365 analyzer->
push( batchL[i], batchR[i] );
368 analyzer->
generate_report( aln_file_nameL, aln_file_nameR, report_name );
383 if (min_batch_size < BATCH_SIZE)
390 else if (argc == arg + 4)
392 const char *aln_file_nameL1 = argv[arg];
393 const char *aln_file_nameL2 = argv[arg+1];
394 const char *aln_file_nameR1 = argv[arg+2];
395 const char *aln_file_nameR2 = argv[arg+3];
397 const char* aln_file_nameL = aln_file_nameL1;
398 const char* aln_file_nameR = aln_file_nameR1;
405 if (aln_streamL1 == NULL || aln_streamL1->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameL1); exit(1); }
406 if (aln_streamL2 == NULL || aln_streamL2->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameL2); exit(1); }
408 if (aln_streamR2 == NULL || aln_streamR1->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameR1); exit(1); }
409 if (aln_streamR2 == NULL || aln_streamR2->
is_ok() ==
false) {
log_error(stderr,
"failed opening \"%s\"\n", aln_file_nameR2); exit(1); }
411 const uint32 BATCH_SIZE = 500000;
412 std::vector<Alignment> batchL1( BATCH_SIZE );
413 std::vector<Alignment> batchL2( BATCH_SIZE );
414 std::vector<Alignment> batchR1( BATCH_SIZE );
415 std::vector<Alignment> batchR2( BATCH_SIZE );
417 Filter filter( filter_name, filter_flags, filter_stats, filter_delta );
423 const uint32 batch_sizeL1 = aln_streamL1->
next_batch( BATCH_SIZE, &batchL1[0] );
424 const uint32 batch_sizeL2 = aln_streamL2->
next_batch( BATCH_SIZE, &batchL2[0] );
425 const uint32 batch_sizeR1 = aln_streamR1->
next_batch( BATCH_SIZE, &batchR1[0] );
426 const uint32 batch_sizeR2 = aln_streamR2->
next_batch( BATCH_SIZE, &batchR2[0] );
436 if (min_batch_size != max_batch_size)
437 log_warning(stderr,
"alignment files have different size\n");
439 log_info(stderr,
"analizing batch[%u]: %u alignments (%.1f M)\n", n_batch, min_batch_size,
float(n_batch * BATCH_SIZE + min_batch_size)*1.0e-6f);
441 for (
uint32 i = 0; i < min_batch_size; ++i)
447 analyzer->
generate_report( aln_file_nameL, aln_file_nameR, report_name );
467 if (min_batch_size < BATCH_SIZE)