annotate IMSAME/src/combine_reads.c @ 5:353e5673bcb8 draft default tip

Uploaded
author bitlab
date Mon, 17 Dec 2018 12:20:41 -0500
parents 762009a91895
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
762009a91895 Uploaded
bitlab
parents:
diff changeset
1 #include <stdio.h>
762009a91895 Uploaded
bitlab
parents:
diff changeset
2 #include <stdlib.h>
762009a91895 Uploaded
bitlab
parents:
diff changeset
3 #include <string.h>
762009a91895 Uploaded
bitlab
parents:
diff changeset
4 #include <ctype.h>
762009a91895 Uploaded
bitlab
parents:
diff changeset
5 #include "structs.h"
762009a91895 Uploaded
bitlab
parents:
diff changeset
6 #include "commonFunctions.h"
762009a91895 Uploaded
bitlab
parents:
diff changeset
7
762009a91895 Uploaded
bitlab
parents:
diff changeset
8 #define MAX(x, y) (((x) > (y)) ? (x) : (y))
762009a91895 Uploaded
bitlab
parents:
diff changeset
9 #define MIN(x, y) (((x) <= (y)) ? (x) : (y))
762009a91895 Uploaded
bitlab
parents:
diff changeset
10 #define STARTING_SEQS 1000
762009a91895 Uploaded
bitlab
parents:
diff changeset
11 #define PIECE_OF_DB_REALLOC 3200000 //half a gigabyte if divided by 8 bytes
762009a91895 Uploaded
bitlab
parents:
diff changeset
12 #define MATVAL 101
762009a91895 Uploaded
bitlab
parents:
diff changeset
13
762009a91895 Uploaded
bitlab
parents:
diff changeset
14 // void terror(char *s) {
762009a91895 Uploaded
bitlab
parents:
diff changeset
15 // fprintf(stdout, "ERR**** %s ****\n", s);
762009a91895 Uploaded
bitlab
parents:
diff changeset
16 // exit(-1);
762009a91895 Uploaded
bitlab
parents:
diff changeset
17 // }
762009a91895 Uploaded
bitlab
parents:
diff changeset
18
762009a91895 Uploaded
bitlab
parents:
diff changeset
19
762009a91895 Uploaded
bitlab
parents:
diff changeset
20 int main(int argc, char ** av){
762009a91895 Uploaded
bitlab
parents:
diff changeset
21
762009a91895 Uploaded
bitlab
parents:
diff changeset
22 if(argc != 2) terror("USE:: combine_reads <file>");
762009a91895 Uploaded
bitlab
parents:
diff changeset
23
762009a91895 Uploaded
bitlab
parents:
diff changeset
24
762009a91895 Uploaded
bitlab
parents:
diff changeset
25 char names[10][100] = {"14np", "47np", "BR0716C", "BR11044P", "CA01044C", "CA01067C", "CA11149P", "CA11154P", "MA07066C", "MA09032P"};
762009a91895 Uploaded
bitlab
parents:
diff changeset
26 uint64_t llen[10] = {1758451, 959020, 1616163, 884193, 934256, 1007436, 1670521, 584151, 671010, 516632};
762009a91895 Uploaded
bitlab
parents:
diff changeset
27 uint64_t nseqs = 10;
762009a91895 Uploaded
bitlab
parents:
diff changeset
28
762009a91895 Uploaded
bitlab
parents:
diff changeset
29 FILE * results, * data;
762009a91895 Uploaded
bitlab
parents:
diff changeset
30 data = fopen(av[1], "rt");
762009a91895 Uploaded
bitlab
parents:
diff changeset
31 if(data == NULL) terror("Could not open input file");
762009a91895 Uploaded
bitlab
parents:
diff changeset
32
762009a91895 Uploaded
bitlab
parents:
diff changeset
33 uint64_t * mat = (uint64_t *) calloc(MATVAL, sizeof(uint64_t));
762009a91895 Uploaded
bitlab
parents:
diff changeset
34 if(mat == NULL) terror("Could not allocate matrix array");
762009a91895 Uploaded
bitlab
parents:
diff changeset
35 long double * mat_e = (long double *) calloc(MATVAL, sizeof(long double));
762009a91895 Uploaded
bitlab
parents:
diff changeset
36 if(mat_e == NULL) terror("Could not allocate float table");
762009a91895 Uploaded
bitlab
parents:
diff changeset
37
762009a91895 Uploaded
bitlab
parents:
diff changeset
38 char current_file[100];
762009a91895 Uploaded
bitlab
parents:
diff changeset
39 uint64_t i=5, idx=0;
762009a91895 Uploaded
bitlab
parents:
diff changeset
40 while(av[1][i] != '_' && av[1][i+1] != 'g'){
762009a91895 Uploaded
bitlab
parents:
diff changeset
41 current_file[idx++] = av[1][i];
762009a91895 Uploaded
bitlab
parents:
diff changeset
42 i++;
762009a91895 Uploaded
bitlab
parents:
diff changeset
43 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
44 current_file[idx] = '\0';
762009a91895 Uploaded
bitlab
parents:
diff changeset
45 fprintf(stdout, "Current file is %s\n", current_file);
762009a91895 Uploaded
bitlab
parents:
diff changeset
46
762009a91895 Uploaded
bitlab
parents:
diff changeset
47 char buffer[MAXLID];
762009a91895 Uploaded
bitlab
parents:
diff changeset
48 if ((results = fopen("accu.log", "r")) == NULL){
762009a91895 Uploaded
bitlab
parents:
diff changeset
49 results = fopen("accu.log", "wt");
762009a91895 Uploaded
bitlab
parents:
diff changeset
50 }else{
762009a91895 Uploaded
bitlab
parents:
diff changeset
51 // Load the matrix
762009a91895 Uploaded
bitlab
parents:
diff changeset
52
762009a91895 Uploaded
bitlab
parents:
diff changeset
53
762009a91895 Uploaded
bitlab
parents:
diff changeset
54 for(i=0;i<100;i++){
762009a91895 Uploaded
bitlab
parents:
diff changeset
55 if(0 == fgets(buffer, MAXLID, results)) terror("Missing number on load");
762009a91895 Uploaded
bitlab
parents:
diff changeset
56
762009a91895 Uploaded
bitlab
parents:
diff changeset
57 //fprintf(stdout, "Have %s\n", buffer);
762009a91895 Uploaded
bitlab
parents:
diff changeset
58 buffer[strlen(buffer)-1] = '\0';
762009a91895 Uploaded
bitlab
parents:
diff changeset
59 //mat[i] = asciiToUint64(buffer);
762009a91895 Uploaded
bitlab
parents:
diff changeset
60 mat_e[i] = (long double) atof(buffer);
762009a91895 Uploaded
bitlab
parents:
diff changeset
61 //fprintf(stdout, "%"PRIu64"\n", mat[i]);
762009a91895 Uploaded
bitlab
parents:
diff changeset
62 //getchar();
762009a91895 Uploaded
bitlab
parents:
diff changeset
63 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
64 fclose(results);
762009a91895 Uploaded
bitlab
parents:
diff changeset
65 results = fopen("accu.log", "wt"); // Re open
762009a91895 Uploaded
bitlab
parents:
diff changeset
66 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
67
762009a91895 Uploaded
bitlab
parents:
diff changeset
68 // Read file
762009a91895 Uploaded
bitlab
parents:
diff changeset
69 uint64_t read_id_1, read_id_2, coverage, identity, length, current, currmax, j, last = 100000, lastmax = 0xFFFFFFFFFFFFFFFF;
762009a91895 Uploaded
bitlab
parents:
diff changeset
70 while(!feof(data)){
762009a91895 Uploaded
bitlab
parents:
diff changeset
71 if(0 == fgets(buffer, MAXLID, data) && !feof(data)) terror("Missing values");
762009a91895 Uploaded
bitlab
parents:
diff changeset
72 // 2 77277 89 64 213
762009a91895 Uploaded
bitlab
parents:
diff changeset
73 sscanf(buffer, "%"PRIu64" %"PRIu64" %"PRIu64" %"PRIu64" %"PRIu64, &read_id_1, &read_id_2, &coverage, &identity, &length);
762009a91895 Uploaded
bitlab
parents:
diff changeset
74 //fprintf(stdout, "%"PRIu64" %"PRIu64" %"PRIu64" %"PRIu64" %"PRIu64"\n", read_id_1, read_id_2, coverage, identity, length);
762009a91895 Uploaded
bitlab
parents:
diff changeset
75 currmax = MIN(identity, coverage);
762009a91895 Uploaded
bitlab
parents:
diff changeset
76 //fprintf(stdout, "%"PRIu64"\n", currmax);
762009a91895 Uploaded
bitlab
parents:
diff changeset
77 current = read_id_1;
762009a91895 Uploaded
bitlab
parents:
diff changeset
78 /*
762009a91895 Uploaded
bitlab
parents:
diff changeset
79 for(j=currmax; j > 1; j--){
762009a91895 Uploaded
bitlab
parents:
diff changeset
80 if(current != lasts[j]){
762009a91895 Uploaded
bitlab
parents:
diff changeset
81 mat[j]++;
762009a91895 Uploaded
bitlab
parents:
diff changeset
82 lasts[j] = current;
762009a91895 Uploaded
bitlab
parents:
diff changeset
83 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
84 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
85 */
762009a91895 Uploaded
bitlab
parents:
diff changeset
86 if(lastmax == 0xFFFFFFFFFFFFFFFF){
762009a91895 Uploaded
bitlab
parents:
diff changeset
87 last = current;
762009a91895 Uploaded
bitlab
parents:
diff changeset
88 lastmax = 0;
762009a91895 Uploaded
bitlab
parents:
diff changeset
89 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
90 if(current == last){
762009a91895 Uploaded
bitlab
parents:
diff changeset
91 //mat[currmax] += 1;
762009a91895 Uploaded
bitlab
parents:
diff changeset
92 if(lastmax < currmax){
762009a91895 Uploaded
bitlab
parents:
diff changeset
93 lastmax = currmax;
762009a91895 Uploaded
bitlab
parents:
diff changeset
94 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
95 }else{
762009a91895 Uploaded
bitlab
parents:
diff changeset
96 mat[lastmax]++;
762009a91895 Uploaded
bitlab
parents:
diff changeset
97 last = current;
762009a91895 Uploaded
bitlab
parents:
diff changeset
98 lastmax = 0;
762009a91895 Uploaded
bitlab
parents:
diff changeset
99 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
100
762009a91895 Uploaded
bitlab
parents:
diff changeset
101 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
102
762009a91895 Uploaded
bitlab
parents:
diff changeset
103 uint64_t divider = 0;
762009a91895 Uploaded
bitlab
parents:
diff changeset
104 for(i=0; i < nseqs; i++){
762009a91895 Uploaded
bitlab
parents:
diff changeset
105
762009a91895 Uploaded
bitlab
parents:
diff changeset
106 if(strcmp(current_file, names[i]) == 0){
762009a91895 Uploaded
bitlab
parents:
diff changeset
107 divider = llen[i];
762009a91895 Uploaded
bitlab
parents:
diff changeset
108 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
109 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
110 if(divider == 0) terror("Could not find file for divider"); else fprintf(stdout, "The divider is %"PRIu64"\n", divider);
762009a91895 Uploaded
bitlab
parents:
diff changeset
111
762009a91895 Uploaded
bitlab
parents:
diff changeset
112 for(j=99; j>0; j--){
762009a91895 Uploaded
bitlab
parents:
diff changeset
113 mat[j] += mat[j+1];
762009a91895 Uploaded
bitlab
parents:
diff changeset
114 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
115 mat[0] = mat[1];
762009a91895 Uploaded
bitlab
parents:
diff changeset
116
762009a91895 Uploaded
bitlab
parents:
diff changeset
117 for(j=0; j<100; j++){
762009a91895 Uploaded
bitlab
parents:
diff changeset
118 fprintf(stdout, "%"PRIu64" --> %.5f\n", mat[j], (float)((100*(long double)mat[j])/(long double)divider));
762009a91895 Uploaded
bitlab
parents:
diff changeset
119
762009a91895 Uploaded
bitlab
parents:
diff changeset
120 fprintf(results, "%.5f\n", (float)(mat_e[j] + ((100*(long double)mat[j])/(long double)divider))/2);
762009a91895 Uploaded
bitlab
parents:
diff changeset
121 }
762009a91895 Uploaded
bitlab
parents:
diff changeset
122
762009a91895 Uploaded
bitlab
parents:
diff changeset
123
762009a91895 Uploaded
bitlab
parents:
diff changeset
124 fclose(results);
762009a91895 Uploaded
bitlab
parents:
diff changeset
125 fclose(data);
762009a91895 Uploaded
bitlab
parents:
diff changeset
126 free(mat);
762009a91895 Uploaded
bitlab
parents:
diff changeset
127 return 0;
762009a91895 Uploaded
bitlab
parents:
diff changeset
128 }