diff options
author | Calvin Morrison <mutantturkey@gmail.com> | 2013-10-18 10:26:12 -0400 |
---|---|---|
committer | Calvin Morrison <mutantturkey@gmail.com> | 2013-10-18 10:26:12 -0400 |
commit | c1a34163771229aa269ada443c6baa38f3073c11 (patch) | |
tree | daa87dd67e45ad910f750587e4769acda68e9eec /src/c/quikr.c | |
parent | 2f62ad7c4ec7ed5bae74c4444442dd871222e6c6 (diff) |
add kmercounting code from dna-utils. update our make file to support debugging, and update quikr to count-kmer code internally
Diffstat (limited to 'src/c/quikr.c')
-rw-r--r-- | src/c/quikr.c | 35 |
1 files changed, 28 insertions, 7 deletions
diff --git a/src/c/quikr.c b/src/c/quikr.c index a954cce..8a238f6 100644 --- a/src/c/quikr.c +++ b/src/c/quikr.c @@ -3,12 +3,12 @@ #include <getopt.h> #include <math.h> #include <stdio.h> -#include <stdio.h> #include <stdlib.h> #include <string.h> #include <unistd.h> #include "nnls.h" +#include "kmer_utils.h" #include "quikr_functions.h" #define sensing_matrix(i,j) (sensing_matrix[width*i + j]) @@ -136,6 +136,11 @@ int main(int argc, char **argv) { exit(EXIT_FAILURE); } + if(kmer == 0) { + fprintf(stderr, "Error: zero is not a valid kmer\n"); + exit(EXIT_FAILURE); + } + // 4 "ACGT" ^ Kmer gives us the size of output rows width = pow(4, kmer); width = width + 1; @@ -149,24 +154,40 @@ int main(int argc, char **argv) { printf("width: %ld\nsequences %ld\n", width, sequences); } + + // load counts matrix + unsigned long long *integer_counts = get_kmer_counts_from_file(input_fasta_filename, kmer); + + double *count_matrix = malloc(sizeof(double) * width); + count_matrix[0] = 0; + + for(x = 1; x < width ; x++) + count_matrix[x] = (double)integer_counts[x-1]; + + free(integer_counts); + + // normalize our count_matrix + normalize_matrix(count_matrix, 1, width); + + for(x = 0; x < width; x++) + count_matrix[x] = count_matrix[x] * lambda; + + // load sensing matrix double *sensing_matrix = load_sensing_matrix(sensing_matrix_filename, sequences, width); - double *count_matrix = load_count_matrix(input_fasta_filename, width, kmer); - // multiply our matrix by lambda + // multiply our sensing matrix by lambda for(x = 1; x < sequences; x++) { for(y= 0; y < width - 1; y++) { sensing_matrix(x, y) = sensing_matrix(x, y) * lambda; } } + // set the first column of our sensing matrix to 0 for(x = 0; x < sequences; x++) { sensing_matrix(x, 0) = 1.0; } - // normalize our count_matrix - normalize_matrix(count_matrix, 1, width); - for(x = 0; x < width; x++) - count_matrix[x] = count_matrix[x] * lambda; + // run NNLS double *solution = nnls(sensing_matrix, count_matrix, sequences, width); // normalize our solution vector |