Mercurial > repos > calkan > mrfast

comparison mrfast-2.1.0.5/HashTable.c @ 1:d4054b05b015 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Version update to 2.1.0.5
author calkan
date Fri, 09 Mar 2012 07:35:51 -0500
parents
children
comparison
equal deleted inserted replaced
0:7b3dc85dc7fd 1:d4054b05b015
1 /*
2 * Copyright (c) <2008 - 2012>, University of Washington, Simon Fraser University
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without modification,
6 * are permitted provided that the following conditions are met:
7 *
8 * Redistributions of source code must retain the above copyright notice, this list
9 * of conditions and the following disclaimer.
10 * - Redistributions in binary form must reproduce the above copyright notice, this
11 * list of conditions and the following disclaimer in the documentation and/or other
12 * materials provided with the distribution.
13 * - Neither the names of the University of Washington, Simon Fraser University,
14 * nor the names of its contributors may be
15 * used to endorse or promote products derived from this software without specific
16 * prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
25 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
26 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
27 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31 /*
32 Authors:
33 Farhad Hormozdiari
34 Faraz Hach
35 Can Alkan
36 Emails:
37 farhadh AT uw DOT edu
38 fhach AT cs DOT sfu DOT ca
39 calkan AT uw DOT edu
40 */
41
42
43
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <string.h>
47 #include <math.h>
48 #include "Common.h"
49 #include "RefGenome.h"
50 #include "HashTable.h"
51 /**********************************************/
52 FILE *_ih_fp = NULL;
53 IHashTable *_ih_hashTable = NULL;
54 int _ih_maxHashTableSize = 0;
55 char *_ih_refGen = NULL;
56 char *_ih_refGenName = NULL;
57 long long _ih_memUsage = 0;
58 int _ih_refGenOff = 0;
59 /**********************************************/
60
61 int hashVal(char *seq)
62 {
63 int i=0;
64 int val=0, numericVal=0;
65
66 while(i<WINDOW_SIZE)
67 {
68 switch (seq[i])
69 {
70 case 'A':
71 numericVal = 0;
72 break;
73 case 'C':
74 numericVal = 1;
75 break;
76 case 'G' :
77 numericVal = 2;
78 break;
79 case 'T':
80 numericVal = 3;
81 break;
82 default:
83 return -1;
84 break;
85 }
86 val = (val << 2)|numericVal;
87 i++;
88 }
89 return val;
90 }
91 /**********************************************/
92 void freeIHashTableContent(IHashTable *hashTable, unsigned int maxSize)
93 {
94 int i=0;
95 for (i=0; i<maxSize; i++)
96 {
97 if (hashTable[i].locs != NULL)
98 {
99 freeMem(hashTable[i].locs, (hashTable[i].locs[0]+1)*(sizeof(unsigned int)));
100 hashTable[i].locs = NULL;
101 }
102 }
103 }
104 /**********************************************/
105 void initSavingIHashTable(char *fileName)
106 {
107 int tmp;
108 _ih_fp = fileOpen(fileName, "w");
109 unsigned char bIndex = 0; // Bisulfite Index
110
111 // First Two bytes are indicating the type of the index & window size
112 tmp = fwrite(&bIndex, sizeof(bIndex), 1, _ih_fp);
113 tmp = fwrite(&WINDOW_SIZE, sizeof(WINDOW_SIZE), 1, _ih_fp);
114
115 if (tmp == 0)
116 fprintf(stderr, "Write error while initializing hash table.\n");
117
118 }
119 /**********************************************/
120 void finalizeSavingIHashTable()
121 {
122 fclose(_ih_fp);
123 }
124 /**********************************************/
125 void saveIHashTable(IHashTable *hashTable, unsigned int size, unsigned int maxSize, char *refGen, char *refGenName, int refGenOffset)
126 {
127 int tmp;
128
129 // Every Chunk starts with a byte indicating whether it has extra info;
130 unsigned char extraInfo = 0;
131 tmp = fwrite (&extraInfo, sizeof(extraInfo), 1, _ih_fp);
132
133 short len = strlen(refGenName);
134 tmp = fwrite(&len, sizeof(len), 1, _ih_fp);
135 tmp = fwrite(refGenName, sizeof(char), len, _ih_fp);
136
137 tmp = fwrite(&refGenOffset, sizeof(refGenOffset), 1, _ih_fp);
138
139 unsigned int refGenLength = strlen(refGen);
140 tmp = fwrite(&refGenLength, sizeof(refGenLength), 1, _ih_fp);
141
142 tmp = fwrite(refGen, sizeof(char), refGenLength, _ih_fp);
143 tmp = fwrite(&size, sizeof(size), 1, _ih_fp);
144
145 int i=0,j=0;
146 unsigned char cnt=0;
147 for (i=0; i<maxSize; i++)
148 {
149 if (hashTable[i].locs != NULL)
150 {
151 tmp = fwrite(&i, sizeof(i), 1, _ih_fp);
152
153 if (hashTable[i].locs[0] < 250)
154 {
155 cnt = hashTable[i].locs[0];
156 tmp = fwrite(&cnt, sizeof(cnt), 1, _ih_fp);
157 }
158 else
159 {
160 cnt =0;
161 tmp = fwrite (&cnt, sizeof(cnt), 1, _ih_fp);
162 tmp = fwrite (&(hashTable[i].locs[0]), sizeof(hashTable[i].locs[0]), 1, _ih_fp);
163 }
164
165 for (j=1; j<=hashTable[i].locs[0]; j++)
166 tmp = fwrite(&(hashTable[i].locs[j]), sizeof(hashTable[i].locs[j]), 1, _ih_fp);
167 }
168 }
169
170 if (tmp == 0)
171 fprintf(stderr, "Write error while saving hash table.\n");
172 }
173 /**********************************************/
174 unsigned int addIHashTableLocation(IHashTable *hashTable, int hv, int location)
175 {
176 unsigned int sizeInc = 0;
177
178 if (hashTable[hv].locs == NULL)
179 {
180 sizeInc = 1;
181 hashTable[hv].locs = getMem (sizeof(unsigned int)*2);
182 hashTable[hv].locs[0]=1;
183 hashTable[hv].locs[1]=location;
184 }
185 else
186 {
187 int size = hashTable[hv].locs[0];
188 int i;
189 unsigned int *tmp = getMem( (size + 2) * sizeof(unsigned int) );
190
191 for (i = 0; i <= size; i++)
192 {
193 tmp[i] = hashTable[hv].locs[i];
194 }
195 size++;
196 tmp[0] = size;
197 tmp[size] = location;
198
199 freeMem(hashTable[hv].locs, (hashTable[hv].locs[0]*(sizeof(unsigned int))));
200 hashTable[hv].locs = tmp;
201 }
202 return sizeInc;
203 }
204 /**********************************************/
205 void generateIHashTable(char *fileName, char *indexName)
206 {
207 double startTime = getTime();
208 unsigned int hashTableSize = 0;
209 unsigned int hashTableMaxSize = pow(4, WINDOW_SIZE);
210 IHashTable *hashTable = getMem(sizeof(IHashTable)*hashTableMaxSize);
211 char *refGenName;
212 char *refGen;
213 int refGenOff = 0;
214 int i, hv, l, flag;
215
216
217 for ( i = 0; i < hashTableMaxSize; i++)
218 {
219 hashTable[i].locs = NULL;
220 }
221
222 //Loading Fasta File
223 if (!initLoadingRefGenome(fileName))
224 return;
225 initSavingIHashTable(indexName);
226
227 fprintf(stdout, "Generating Index from %s", fileName);
228 fflush(stdout);
229
230 char *prev = getMem (CONTIG_NAME_SIZE);
231 prev[0]='\0';
232
233 do
234 {
235 flag = loadRefGenome (&refGen, &refGenName, &refGenOff);
236
237 if ( strcmp(prev, refGenName) != 0)
238 {
239 fprintf(stdout, "\n - %s ", refGenName);
240 fflush(stdout);
241 sprintf(prev, "%s", refGenName);
242 }
243 else
244 {
245 fprintf(stdout, ".");
246 fflush(stdout);
247 }
248
249 l = strlen(refGen) - WINDOW_SIZE;
250
251 for (i=0; i < l; i++)
252 {
253 hv = hashVal(refGen + i);
254 if (hv != -1)
255 {
256 hashTableSize += addIHashTableLocation (hashTable, hv, i+1);
257 }
258 }
259
260 saveIHashTable(hashTable, hashTableSize, hashTableMaxSize, refGen, refGenName, refGenOff);
261 freeIHashTableContent(hashTable, hashTableMaxSize);
262 hashTableSize = 0;
263 } while (flag);
264
265 freeMem(prev, CONTIG_NAME_SIZE);
266 freeMem(hashTable, sizeof(IHashTable)*hashTableMaxSize);
267
268 finalizeLoadingRefGenome();
269 finalizeSavingIHashTable();
270
271 fprintf(stdout, "\nDONE in %0.2fs!\n", (getTime()-startTime));
272 }
273
274 /**********************************************/
275 void finalizeLoadingIHashTable()
276 {
277 freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
278 freeMem(_ih_hashTable, sizeof(IHashTable)* _ih_maxHashTableSize);
279 freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
280 freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
281 fclose(_ih_fp);
282 }
283
284 /**********************************************/
285 int loadIHashTable(double *loadTime, int errThreshould)
286 {
287 double startTime = getTime();
288 unsigned char extraInfo = 0;
289 short len;
290 unsigned int refGenLength;
291 unsigned int hashTableSize;
292 unsigned int tmpSize;
293 int tmp;
294 int i=0,j=0;
295
296 if ( fread(&extraInfo, sizeof(extraInfo), 1, _ih_fp) != sizeof(extraInfo) )
297 return 0;
298
299 freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
300 freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
301 freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
302
303 // Reading Chr Name
304 tmp = fread(&len, sizeof(len), 1, _ih_fp);
305 if (tmp == 0){
306 fprintf(stderr, "Read error while loading hash table.\n");
307 exit(0);
308 }
309
310 _ih_refGenName = getMem(sizeof(char)* (len+1));
311 tmp = fread(_ih_refGenName, sizeof(char), len, _ih_fp);
312 _ih_refGenName [len] ='\0';
313
314 tmp = fread(&_ih_refGenOff, sizeof (_ih_refGenOff), 1, _ih_fp);
315
316 // Reading Size and Content of Ref Genome
317 tmp = fread(&refGenLength, sizeof(refGenLength), 1, _ih_fp);
318 _ih_refGen = getMem(sizeof(char)*(refGenLength+1));
319 tmp = fread(_ih_refGen, sizeof(char), refGenLength, _ih_fp);
320 _ih_refGen[refGenLength]='\0';
321
322 //Reading Hashtable Size and Content
323 tmp = fread(&hashTableSize, sizeof(hashTableSize), 1, _ih_fp);
324
325 unsigned int hv;
326 unsigned char cnt=0;
327 unsigned int index = 1;
328 unsigned int tmpValue = 0;
329 for (i=0; i<hashTableSize; i++)
330 {
331 index = 1;
332 tmp = fread(&hv, sizeof(hv), 1, _ih_fp);
333 tmp = fread(&cnt, sizeof(cnt), 1, _ih_fp);
334
335 if (cnt>0)
336 {
337 tmpSize = cnt;
338 }
339 else
340 {
341 tmp = fread(&tmpSize, sizeof(tmpSize), 1, _ih_fp);
342 }
343
344 _ih_hashTable[hv].locs = getMem( sizeof(unsigned int)* (tmpSize+1) );
345 _ih_hashTable[hv].locs[0] = tmpSize;
346
347 /* for (j=1; j<=tmpSize; j++)
348 tmp = fread(&(_ih_hashTable[hv].locs[j]), sizeof(_ih_hashTable[hv].locs[j]), 1, _ih_fp);
349 */
350
351 if(tmpSize != 0)
352 {
353 tmp = fread(&(_ih_hashTable[hv].locs[index]), sizeof(_ih_hashTable[hv].locs[index]), 1, _ih_fp);
354 index++;
355 }
356
357 for (j=2; j<=tmpSize; j++)
358 {
359 tmp = fread(&tmpValue, sizeof(_ih_hashTable[hv].locs[index-1]), 1, _ih_fp);
360 if(abs(tmpValue-_ih_hashTable[hv].locs[index-1]) > errThreshould)
361 {
362 _ih_hashTable[hv].locs[index] = tmpValue;
363 index++;
364 }
365 }
366 _ih_hashTable[hv].locs[0] = index-1;
367 }
368 *loadTime = getTime()-startTime;
369 return 1;
370 }
371 /**********************************************/
372 unsigned int *getIHashTableCandidates(int hv)
373 {
374 if ( hv != -1 )
375 return _ih_hashTable[hv].locs;
376 else
377 return NULL;
378 }
379 /**********************************************/
380 /**********************************************/
381 /**********************************************/
382 void configHashTable()
383 {
384 if (WINDOW_SIZE <= 15)
385 {
386 generateHashTable = &generateIHashTable;
387 loadHashTable = &loadIHashTable;
388 finalizeLoadingHashTable = &finalizeLoadingIHashTable;
389 getCandidates = &getIHashTableCandidates;
390 }
391 else
392 {
393
394
395 }
396 }
397 /**********************************************/
398 int initLoadingHashTable(char *fileName)
399 {
400 int i;
401 unsigned char bsIndex;
402 int tmp;
403
404 _ih_fp = fileOpen(fileName, "r");
405
406 if (_ih_fp == NULL)
407 return 0;
408
409 tmp = fread(&bsIndex, sizeof(bsIndex), 1, _ih_fp);
410
411 if (tmp == 0){
412 fprintf(stderr, "Read error while inilializing hash table loading.\n");
413 exit(0);
414 }
415
416 if (bsIndex)
417 {
418 fprintf(stdout, "Error: Wrong Type of Index indicated");
419 return 0;
420 }
421
422 tmp = fread(&WINDOW_SIZE, sizeof(WINDOW_SIZE), 1, _ih_fp);
423
424 configHashTable();
425
426 if (_ih_maxHashTableSize != pow(4, WINDOW_SIZE))
427 {
428
429 if (_ih_hashTable != NULL)
430 {
431 freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
432 freeMem(_ih_hashTable, sizeof(IHashTable)* _ih_maxHashTableSize);
433 freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
434 freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
435 }
436
437 _ih_maxHashTableSize = pow(4, WINDOW_SIZE);
438
439 _ih_hashTable = getMem (sizeof(IHashTable) * _ih_maxHashTableSize);
440 for (i=0; i<_ih_maxHashTableSize; i++)
441 _ih_hashTable[i].locs = NULL;
442 _ih_refGen = getMem(1);
443 _ih_refGen[0]='\0';
444 _ih_refGenName = getMem(1);
445 _ih_refGenName[0] = '\0';
446 }
447
448 return 1;
449 }
450 /**********************************************/
451 char *getRefGenome()
452 {
453 return _ih_refGen;
454
455 }
456 /**********************************************/
457 char *getRefGenomeName()
458 {
459 return _ih_refGenName;
460
461 }
462 /**********************************************/
463 int getRefGenomeOffset()
464 {
465 return _ih_refGenOff;
466
467 }
468 /**********************************************/
469 HashTable *getHashTable()
470 {
471 return NULL;
472 }