From d7ed90dc545e54a09ef206634bb6042e5893c04b Mon Sep 17 00:00:00 2001 From: Calvin Morrison Date: Thu, 7 Nov 2013 15:38:00 -0500 Subject: update to work with kmer rarity --- src/c/quikr.c | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 120 insertions(+), 11 deletions(-) (limited to 'src/c/quikr.c') diff --git a/src/c/quikr.c b/src/c/quikr.c index 9fc0759..1b8daa3 100644 --- a/src/c/quikr.c +++ b/src/c/quikr.c @@ -1,4 +1,3 @@ -#include #include #include #include @@ -24,6 +23,14 @@ int main(int argc, char **argv) { char *output_filename = NULL; unsigned long long x = 0; + unsigned long long y = 0; + unsigned long long z = 0; + + unsigned long long rare_value = 0; + unsigned long long rare_width = 0; + + double percentage = 1.0; + double rare_percent = 1.0; unsigned long long width = 0; @@ -39,6 +46,7 @@ int main(int argc, char **argv) { {"lambda", required_argument, 0, 'l'}, {"output", required_argument, 0, 'o'}, {"sensing-matrix", required_argument, 0, 's'}, + {"rare-percent", required_argument, 0, 'r'}, {"verbose", no_argument, 0, 'v'}, {"version", no_argument, 0, 'V'}, {"help", no_argument, 0, 'h'}, @@ -48,7 +56,7 @@ int main(int argc, char **argv) { int option_index = 0; - c = getopt_long (argc, argv, "k:l:s:i:o:hdvV", long_options, &option_index); + c = getopt_long (argc, argv, "k:l:s:r:i:o:r:hdvV", long_options, &option_index); if (c == -1) break; @@ -60,6 +68,9 @@ int main(int argc, char **argv) { case 'l': lambda = atoi(optarg); break; + case 'r': + rare_percent = atof(optarg); + break; case 's': sensing_matrix_filename = optarg; break; @@ -75,17 +86,14 @@ int main(int argc, char **argv) { case 'V': printf("%s\n", VERSION); exit(EXIT_SUCCESS); - break; case 'h': printf("%s\n", USAGE); exit(EXIT_SUCCESS); - break; default: break; } } - if(sensing_matrix_filename == NULL) { fprintf(stderr, "Error: sensing matrix filename (-s) must be specified\n\n"); fprintf(stderr, "%s\n", USAGE); @@ -102,7 +110,8 @@ int main(int argc, char **argv) { exit(EXIT_FAILURE); } - if(verbose) { + if(verbose) { + printf("rare width:%ld\n", rare_width); printf("kmer: %u\n", kmer); printf("lambda: %llu\n", lambda); printf("fasta: %s\n", input_fasta_filename); @@ -124,18 +133,109 @@ int main(int argc, char **argv) { exit(EXIT_FAILURE); } + // 4 "ACGT" ^ Kmer gives us the size of output rows - width = pow(4, kmer); - width = width + 1; + width = pow_four(kmer); // load counts matrix - double *count_matrix = setup_count_matrix(input_fasta_filename, kmer, lambda, width); + double *count_matrix = malloc(width * sizeof(double)); + check_malloc(count_matrix, NULL); + + { + unsigned long long *integer_counts = get_kmer_counts_from_file(input_fasta_filename, kmer); + + for(x = 0; x < width; x++) + count_matrix[x] = (double)integer_counts[x]; + + free(integer_counts); + } // load sensing matrix - struct matrix *sensing_matrix = setup_sensing_matrix(sensing_matrix_filename, kmer, lambda, width); + struct matrix *sensing_matrix = load_sensing_matrix(sensing_matrix_filename, kmer); + + + // get our "rare" counts + while(1) { + rare_width = 0; + for(x = 0; x < width; x++) { + if(count_matrix[x] < rare_value) { + rare_width++; + } + } + percentage = (float)rare_width / (float)width; + + if(percentage >= rare_percent) + break; + + rare_value++; + } + + if(verbose) + printf("there are %llu values less than %llu\n", rare_width, rare_value); + + // add a extra space for our zero's array + rare_width++; + + // store our count matrix + double *count_matrix_rare = calloc(rare_width, sizeof(double)); + check_malloc(count_matrix_rare, NULL); + + double *sensing_matrix_rare = malloc((rare_width) * sensing_matrix->sequences * sizeof(double)); + check_malloc(sensing_matrix_rare, NULL); + + // copy only kmers from our original counts that match our rareness percentage + // in both our count matrix and our sensing matrix + for(x = 0, y = 1; x < width; x++) { + if(count_matrix[x] < rare_value) { + count_matrix_rare[y] = count_matrix[x]; + + for(z = 0; z < sensing_matrix->sequences; z++) + sensing_matrix_rare[z*rare_width + y] = sensing_matrix->matrix[z*width + x]; + + y++; + } + } + + // normalize our kmer counts and our sensing_matrix + normalize_matrix(count_matrix_rare, 1, rare_width); + normalize_matrix(sensing_matrix_rare, sensing_matrix->sequences, rare_width); + + // multiply our kmer counts and sensing matrix by lambda + for(x = 1; x < rare_width; x++) + count_matrix_rare[x] *= lambda; + + for(x = 0; x < sensing_matrix->sequences; x++) { + for(y = 1; y < rare_width; y++) { + sensing_matrix_rare[rare_width*x + y] *= lambda; + } + } + + // count_matrix's first element should be zero + count_matrix_rare[0] = 0; + // stack one's on our first row of our sensing matrix + for(x = 0; x < sensing_matrix->sequences; x++) { + sensing_matrix_rare[x*rare_width] = 1.0; + } + + // DEBUG OUR MATRIX + if(verbose) { + printf("COUNT_MATRIX\n"); + for(x = 0; x < rare_width; x++) + printf("%lf\t", count_matrix_rare[x]); + printf("\n"); + + printf("SENSING_MATRIX\n"); + for(y = 0; y < sensing_matrix->sequences - 1; y++) { + for(x = 0; x < rare_width; x++) { + printf("%lf\t", sensing_matrix_rare[x*rare_width + y]); + } + printf("%lf\n", sensing_matrix_rare[x*rare_width + (rare_width - 1)]); + } + + } // run NNLS - double *solution = nnls(sensing_matrix->matrix, count_matrix, sensing_matrix->sequences, width); + double *solution = nnls(sensing_matrix_rare, count_matrix_rare, sensing_matrix->sequences, rare_width); // normalize our solution vector normalize_matrix(solution, 1, sensing_matrix->sequences); @@ -149,7 +249,16 @@ int main(int argc, char **argv) { for(x = 0; x < sensing_matrix->sequences; x++) { fprintf(output_fh, "%.10lf\n", solution[x]); } + fclose(output_fh); + free(solution); + + free(count_matrix); + free(sensing_matrix); + + free(count_matrix_rare); + free(sensing_matrix_rare); + return EXIT_SUCCESS; } -- cgit v1.2.3