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