?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
limit fstBC to 30bp in Python3 ver.
[GalaxyCodeBases.git] / BGI / SOAPdenovo2 / standardPregraph / prlHashReads.c
blob6a9158e21082f8f3c4204953a0ba0521e3556062
1 /*
2 * prlHashReads.c
3 *
4 * Copyright (c) 2008-2012 BGI-Shenzhen <soap at genomics dot org dot cn>.
5 *
6 * This file is part of SOAPdenovo.
7 *
8 * SOAPdenovo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * SOAPdenovo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with SOAPdenovo. If not, see <http://www.gnu.org/licenses/>.
20 *
21 */
22
23 #include "stdinc.h"
24 #include "newhash.h"
25 #include "kmerhash.h"
26 #include "extfunc.h"
27 #include "extvab.h"
28
29 //debugging variables
30 static long long * tips;
31 static long long * kmerCounter; //kmer number for each KmerSet[thread id]
32
33 static long long ** kmerFreq;
34
35 //buffer related varibles for chop kmer
36 static int read_c;
37 static char ** rcSeq; //sequence pointer for each thread
38 static char ** seqBuffer; //read buffer
39 static int * lenBuffer; //read length buffer
40 static int * indexArray; // the read's begin kmer 's kmer_c
41
42 //buffer related varibles for splay tree
43 static int buffer_size = 100000000; //buffer size for kmerBuffer...
44 static volatile int kmer_c; // kmer count number
45 static Kmer * kmerBuffer; //kmer buffer array
46 static ubyte8 * hashBanBuffer; //the buffered hash value for 'kmerBuffer'
47 static char * nextcBuffer, *prevcBuffer; //next char buffer , previous char buffer for 'kmerBuffer'
48
49 static struct aiocb aio1;
50 static struct aiocb aio2;
51 static char * aioBuffer1;
52 static char * aioBuffer2;
53 static char * readBuffer1;
54 static char * readBuffer2;
55
56 static void thread_mark ( KmerSet * set, unsigned char thrdID );
57 static void Mark1in1outNode ( unsigned char * thrdSignal );
58 static void thread_delow ( KmerSet * set, unsigned char thrdID );
59 static void deLowCov ( unsigned char * thrdSignal );
60
61 static void singleKmer ( int t, KmerSet * kset );
62 static void chopKmer4read ( int t, int threadID );
63
64 static void freqStat ( char * outfile );
65
66 static void threadRoutine ( void * para )
67 {
68 PARAMETER * prm;
69 int i;
70 unsigned char id;
71 prm = ( PARAMETER * ) para;
72 id = prm->threadID;
73
74 //printf("%dth thread with threadID %d, hash_table %p\n",id,prm.threadID,prm.hash_table);
75 while ( 1 )
76 {
77 if ( * ( prm->selfSignal ) == 1 )
78 {
79 for ( i = 0; i < kmer_c; i++ )
80 {
81 //if((unsigned char)(magic_seq(hashBanBuffer[i])%thrd_num)!=id)
82 //if((kmerBuffer[i]%thrd_num)!=id)
83 if ( ( hashBanBuffer[i] % thrd_num ) != id )
84 {
85 continue;
86 }
87
88 kmerCounter[id + 1]++;
89 singleKmer ( i, KmerSets[id] );
90 }
91
92 * ( prm->selfSignal ) = 0;
93 }
94 else if ( * ( prm->selfSignal ) == 2 )
95 {
96 for ( i = 0; i < read_c; i++ )
97 {
98 if ( i % thrd_num != id )
99 {
100 continue;
101 }
102
103 chopKmer4read ( i, id + 1 );
104 }
105
106 * ( prm->selfSignal ) = 0;
107 }
108 else if ( * ( prm->selfSignal ) == 3 )
109 {
110 * ( prm->selfSignal ) = 0;
111 break;
112 }
113 else if ( * ( prm->selfSignal ) == 4 )
114 {
115 thread_mark ( KmerSets[id], id );
116 * ( prm->selfSignal ) = 0;
117 }
118 else if ( * ( prm->selfSignal ) == 5 )
119 {
120 thread_delow ( KmerSets[id], id );
121 * ( prm->selfSignal ) = 0;
122 }
123
124 usleep ( 1 );
125 }
126 }
127
128 static void singleKmer ( int t, KmerSet * kset )
129 {
130 kmer_t * pos;
131 put_kmerset ( kset, kmerBuffer[t], prevcBuffer[t], nextcBuffer[t], &pos );
132 }
133
134 static void creatThrds ( pthread_t * threads, PARAMETER * paras )
135 {
136 unsigned char i;
137 int temp;
138
139 for ( i = 0; i < thrd_num; i++ )
140 {
141 //printf("to create %dth thread\n",(*(char *)&(threadID[i])));
142 if ( ( temp = pthread_create ( &threads[i], NULL, ( void * ) threadRoutine, & ( paras[i] ) ) ) != 0 )
143 {
144 fprintf ( stderr, "Create threads failed.\n" );
145 exit ( 1 );
146 }
147 }
148
149 fprintf ( stderr, "%d thread(s) initialized.\n", thrd_num );
150 }
151
152 static void thread_wait ( pthread_t * threads )
153 {
154 int i;
155
156 for ( i = 0; i < thrd_num; i++ )
157 if ( threads[i] != 0 )
158 {
159 pthread_join ( threads[i], NULL );
160 }
161 }
162
163 static void chopKmer4read ( int t, int threadID )
164 {
165 char * src_seq = seqBuffer[t];
166 char * bal_seq = rcSeq[threadID];
167 int len_seq = lenBuffer[t];
168 int j, bal_j;
169 ubyte8 hash_ban, bal_hash_ban;
170 Kmer word, bal_word;
171 int index;
172 char InvalidCh = 4;
173 #ifdef MER127
174 word.high1 = word.low1 = word.high2 = word.low2 = 0;
175
176 for ( index = 0; index < overlaplen; index++ )
177 {
178 word = KmerLeftBitMoveBy2 ( word );
179 word.low2 |= src_seq[index];
180 }
181
182 #else
183 word.high = word.low = 0;
184
185 for ( index = 0; index < overlaplen; index++ )
186 {
187 word = KmerLeftBitMoveBy2 ( word );
188 word.low |= src_seq[index];
189 }
190
191 #endif
192 reverseComplementSeq ( src_seq, len_seq, bal_seq );
193 // complementary node
194 bal_word = reverseComplement ( word, overlaplen );
195 bal_j = len_seq - overlaplen;
196 index = indexArray[t];
197
198 if ( KmerSmaller ( word, bal_word ) )
199 {
200 hash_ban = hash_kmer ( word );
201 hashBanBuffer[index] = hash_ban;
202 kmerBuffer[index] = word;
203 prevcBuffer[index] = InvalidCh;
204 nextcBuffer[index++] = src_seq[0 + overlaplen];
205 }
206 else
207 {
208 bal_hash_ban = hash_kmer ( bal_word );
209 hashBanBuffer[index] = bal_hash_ban;
210 kmerBuffer[index] = bal_word;
211 prevcBuffer[index] = bal_seq[bal_j - 1];
212 nextcBuffer[index++] = InvalidCh;
213 }
214
215 for ( j = 1; j <= len_seq - overlaplen; j++ )
216 {
217 word = nextKmer ( word, src_seq[j - 1 + overlaplen] );
218 bal_j = len_seq - j - overlaplen;
219 bal_word = prevKmer ( bal_word, bal_seq[bal_j] );
220
221 if ( KmerSmaller ( word, bal_word ) )
222 {
223 hash_ban = hash_kmer ( word );
224 hashBanBuffer[index] = hash_ban;
225 kmerBuffer[index] = word;
226 prevcBuffer[index] = src_seq[j - 1];
227
228 if ( j < len_seq - overlaplen )
229 {
230 nextcBuffer[index++] = src_seq[j + overlaplen];
231 }
232 else
233 {
234 nextcBuffer[index++] = InvalidCh;
235 }
236
237 //printf("%dth: %p with %p\n",kmer_c-1,word,hashBanBuffer[kmer_c-1]);
238 }
239 else
240 {
241 // complementary node
242 bal_hash_ban = hash_kmer ( bal_word );
243 hashBanBuffer[index] = bal_hash_ban;
244 kmerBuffer[index] = bal_word;
245
246 if ( bal_j > 0 )
247 {
248 prevcBuffer[index] = bal_seq[bal_j - 1];
249 }
250 else
251 {
252 prevcBuffer[index] = InvalidCh;
253 }
254
255 nextcBuffer[index++] = bal_seq[bal_j + overlaplen];
256 //printf("%dth: %p with %p\n",kmer_c-1,bal_word,hashBanBuffer[kmer_c-1]);
257 }
258 }
259 }
260
261 static void sendWorkSignal ( unsigned char SIG, unsigned char * thrdSignals )
262 {
263 int t;
264
265 for ( t = 0; t < thrd_num; t++ )
266 {
267 thrdSignals[t + 1] = SIG;
268 }
269
270 while ( 1 )
271 {
272 usleep ( 10 );
273
274 for ( t = 0; t < thrd_num; t++ )
275 if ( thrdSignals[t + 1] )
276 {
277 break;
278 }
279
280 if ( t == thrd_num )
281 {
282 break;
283 }
284 }
285 }
286
287 /*************************************************
288 Function:
289 prlRead2HashTable
290 Description:
291 1. Imports the reads from the lib file one by one.
292 2. Chops the reads into kmers and store them in KmerSets.
293 3. Removes the kmers with low coverage.
294 4. Marks the linear kmers.
295 5. Counts the kmer frequences.
296 Input:
297 1. libfile : the reads config file
298 2. outfile : the output file prefix
299 Output:
300 None.
301 Return:
302 1 if exits normally.
303 *************************************************/
304 boolean prlRead2HashTable ( char * libfile, char * outfile )
305 {
306 char * cach1;
307 char * cach2;
308 unsigned char asm_ctg = 1;
309 long long i;
310 char * next_name, name[256];
311 FILE * fo;
312 time_t start_t, stop_t;
313 int maxReadNum;
314 int libNo;
315 pthread_t threads[thrd_num];
316 unsigned char thrdSignal[thrd_num + 1];
317 PARAMETER paras[thrd_num];
318 boolean flag, pairs = 0;
319 WORDFILTER = createFilter ( overlaplen );
320 maxReadLen = 0;
321 maxNameLen = 256;
322 scan_libInfo ( libfile );
323 alloc_pe_mem ( num_libs );
324
325 if ( !maxReadLen )
326 {
327 maxReadLen = 100;
328 }
329
330 if ( gLineLen < maxReadLen )
331 {
332 gStr = ( char * ) ckalloc ( ( maxReadLen + 1 ) * sizeof ( char ) );
333 }
334
335 //init
336 maxReadLen4all = maxReadLen;
337 fprintf ( stderr, "In %s, %d lib(s), maximum read length %d, maximum name length %d.\n\n", libfile, num_libs, maxReadLen, maxNameLen );
338 next_name = ( char * ) ckalloc ( ( maxNameLen + 1 ) * sizeof ( char ) );
339 kmerBuffer = ( Kmer * ) ckalloc ( buffer_size * sizeof ( Kmer ) );
340 hashBanBuffer = ( ubyte8 * ) ckalloc ( buffer_size * sizeof ( ubyte8 ) );
341 prevcBuffer = ( char * ) ckalloc ( buffer_size * sizeof ( char ) );
342 nextcBuffer = ( char * ) ckalloc ( buffer_size * sizeof ( char ) );
343 maxReadNum = buffer_size / ( maxReadLen - overlaplen + 1 );
344 //printf("buffer size %d, max read len %d, max read num %d\n",buffer_size,maxReadLen,maxReadNum);
345 int maxAIOSize = 32768;
346 aioBuffer1 = ( char * ) ckalloc ( ( maxAIOSize ) * sizeof ( char ) );
347 aioBuffer2 = ( char * ) ckalloc ( ( maxAIOSize ) * sizeof ( char ) );
348 readBuffer1 = ( char * ) ckalloc ( ( maxAIOSize + ( maxReadLen * 4 + 1024 ) ) * sizeof ( char ) ); //(char *)ckalloc(maxAIOSize*sizeof(char)); //1024
349 readBuffer2 = ( char * ) ckalloc ( ( maxAIOSize + ( maxReadLen * 4 + 1024 ) ) * sizeof ( char ) ); //1024
350 cach1 = ( char * ) ckalloc ( ( maxReadLen * 4 + 1024 ) * sizeof ( char ) ); //1024
351 cach2 = ( char * ) ckalloc ( ( maxReadLen * 4 + 1024 ) * sizeof ( char ) ); //1024
352 memset ( cach1, '\0', ( maxReadLen * 4 + 1024 ) ); //1024
353 memset ( cach2, '\0', ( maxReadLen * 4 + 1024 ) ); //1024
354 seqBuffer = ( char ** ) ckalloc ( maxReadNum * sizeof ( char * ) );
355 lenBuffer = ( int * ) ckalloc ( maxReadNum * sizeof ( int ) );
356 indexArray = ( int * ) ckalloc ( maxReadNum * sizeof ( int ) );
357
358 for ( i = 0; i < maxReadNum; i++ )
359 {
360 seqBuffer[i] = ( char * ) ckalloc ( maxReadLen * sizeof ( char ) );
361 }
362
363 rcSeq = ( char ** ) ckalloc ( ( thrd_num + 1 ) * sizeof ( char * ) );
364
365 if ( 1 )
366 {
367 kmerCounter = ( long long * ) ckalloc ( ( thrd_num + 1 ) * sizeof ( long long ) );
368 KmerSets = ( KmerSet ** ) ckalloc ( thrd_num * sizeof ( KmerSet * ) );
369 ubyte8 init_size = 1024;
370 ubyte8 k = 0;
371
372 if ( initKmerSetSize )
373 {
374 #ifdef MER127
375 init_size = ( ubyte8 ) ( ( double ) initKmerSetSize * 1024.0f * 1024.0f * 1024.0f / ( double ) thrd_num / 40 );
376 #else
377 init_size = ( ubyte8 ) ( ( double ) initKmerSetSize * 1024.0f * 1024.0f * 1024.0f / ( double ) thrd_num / 24 ); //is it true?
378 #endif
379
380 do
381 {
382 ++k;
383 }
384 while ( k * 0xFFFFFFLLU < init_size );
385 }
386
387 for ( i = 0; i < thrd_num; i++ )
388 {
389 //KmerSets[i] = init_kmerset(1024,0.77f);
390 KmerSets[i] = init_kmerset ( ( ( initKmerSetSize ) ? ( k * 0xFFFFFFLLU ) : ( init_size ) ), 0.77f );
391 thrdSignal[i + 1] = 0;
392 paras[i].threadID = i;
393 paras[i].mainSignal = &thrdSignal[0];
394 paras[i].selfSignal = &thrdSignal[i + 1];
395 kmerCounter[i + 1] = 0;
396 rcSeq[i + 1] = ( char * ) ckalloc ( maxReadLen * sizeof ( char ) );
397 }
398
399 creatThrds ( threads, paras );
400 }
401
402 thrdSignal[0] = kmerCounter[0] = 0;
403 time ( &start_t );
404 kmer_c = n_solexa = read_c = i = libNo = readNumBack = gradsCounter = 0;
405
406 while ( openNextFile ( &libNo, pairs, asm_ctg ) )
407 {
408 //read bam file
409 if ( lib_array[libNo].curr_type == 4 )
410 {
411 int type = 0; //deside the PE reads is good or bad
412
413 while ( ( flag = read1seqInLibBam ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), &libNo, pairs, 1, &type ) ) != 0 )
414 {
415 if ( type == -1 ) //if the reads is bad, go back.
416 {
417 i--;
418
419 if ( lenBuffer[read_c - 1] >= overlaplen + 1 )
420 {
421 kmer_c -= lenBuffer[read_c - 1] - overlaplen + 1;
422 read_c--;
423 }
424
425 n_solexa -= 2;
426 continue;
427 }
428
429 if ( ( ++i ) % 100000000 == 0 )
430 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
431
432 if ( lenBuffer[read_c] < 0 )
433 { fprintf ( stderr, "Read len %d.\n", lenBuffer[read_c] ); }
434
435 if ( lenBuffer[read_c] < overlaplen + 1 )
436 { continue; }
437
438 /*
439 if(lenBuffer[read_c]>70)
440 lenBuffer[read_c] = 50;
441 else if(lenBuffer[read_c]>40)
442 lenBuffer[read_c] = 40;
443 */
444 indexArray[read_c] = kmer_c;
445 kmer_c += lenBuffer[read_c] - overlaplen + 1;
446 read_c++;
447
448 if ( read_c == maxReadNum )
449 {
450 kmerCounter[0] += kmer_c;
451 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
452 sendWorkSignal ( 1, thrdSignal ); //singleKmer
453 kmer_c = read_c = 0;
454 }
455 }
456 }
457 //read PE fasta or fastq
458 else if ( lib_array[libNo].curr_type == 1 || lib_array[libNo].curr_type == 2 )
459 {
460 initAIO ( &aio1, aioBuffer1, fileno ( lib_array[libNo].fp1 ), maxAIOSize );
461 initAIO ( &aio2, aioBuffer2, fileno ( lib_array[libNo].fp2 ), maxAIOSize );
462 int offset1, offset2, flag1, flag2, rt1, rt2;
463 offset1 = offset2 = 0;
464 rt1 = aio_read ( &aio1 );
465 rt2 = aio_read ( &aio2 );
466 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
467 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
468
469 if ( flag1 && flag2 )
470 {
471 int start1, start2, turn;
472 start1 = start2 = 0;
473 turn = 1;
474
475 while ( start1 < offset1 || start2 < offset2 )
476 {
477 if ( turn == 1 )
478 {
479 turn = 2;
480 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer1, &start1, offset1, libNo );
481
482 if ( ( ++i ) % 100000000 == 0 )
483 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
484
485 if ( lenBuffer[read_c] < 0 )
486 { fprintf ( stderr, "Read len %d.\n", lenBuffer[read_c] ); }
487
488 if ( lenBuffer[read_c] < overlaplen + 1 )
489 {
490 if ( start1 >= offset1 )
491 {
492 start1 = 0;
493 offset1 = 0;
494 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
495 }
496
497 continue;
498 }
499
500 indexArray[read_c] = kmer_c;
501 kmer_c += lenBuffer[read_c] - overlaplen + 1;
502 read_c++;
503
504 if ( start1 >= offset1 )
505 {
506 start1 = 0;
507 offset1 = 0;
508 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
509 }
510
511 if ( read_c == maxReadNum )
512 {
513 kmerCounter[0] += kmer_c;
514 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
515 sendWorkSignal ( 1, thrdSignal ); //singleKmer
516 kmer_c = read_c = 0;
517 }
518
519 continue;
520 }
521
522 if ( turn == 2 )
523 {
524 turn = 1;
525 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer2, &start2, offset2, libNo );
526
527 if ( ( ++i ) % 100000000 == 0 )
528 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
529
530 if ( lenBuffer[read_c] < 0 )
531 { fprintf ( stderr, "Read len %d.\n", lenBuffer[read_c] ); }
532
533 if ( lenBuffer[read_c] < overlaplen + 1 )
534 {
535 if ( ( flag2 == 2 ) && ( start2 >= offset2 ) )
536 { break; }
537
538 if ( start2 >= offset2 )
539 {
540 start2 = 0;
541 offset2 = 0;
542 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
543 }
544
545 continue;
546 }
547
548 indexArray[read_c] = kmer_c;
549 kmer_c += lenBuffer[read_c] - overlaplen + 1;
550 read_c++;
551
552 if ( ( flag2 == 2 ) && ( start2 >= offset2 ) )
553 { break; }
554
555 if ( start2 >= offset2 )
556 {
557 start2 = 0;
558 offset2 = 0;
559 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
560 }
561
562 if ( read_c == maxReadNum )
563 {
564 kmerCounter[0] += kmer_c;
565 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
566 sendWorkSignal ( 1, thrdSignal ); //singleKmer
567 kmer_c = read_c = 0;
568 }
569
570 continue;
571 }
572 }
573 }
574 else
575 {
576 fprintf(stderr, "Error: aio_read error.\n");
577 }
578 }
579 //read single fasta, single fastq and PE fasta in one file
580 else
581 {
582 initAIO ( &aio1, aioBuffer1, fileno ( lib_array[libNo].fp1 ), maxAIOSize );
583 int offset, flag1, rt;
584 offset = 0;
585 rt = aio_read ( &aio1 );
586
587 while ( ( flag1 = AIORead ( &aio1, &offset, readBuffer1, cach1, &rt, lib_array[libNo].curr_type ) ) )
588 {
589 int start = 0;
590
591 while ( start < offset )
592 {
593 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer1, &start, offset, libNo );
594
595 if ( ( ++i ) % 100000000 == 0 )
596 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
597
598 if ( lenBuffer[read_c] < 0 )
599 { fprintf ( stderr, "Read len %d.\n", lenBuffer[read_c] ); }
600
601 if ( lenBuffer[read_c] < overlaplen + 1 )
602 { continue; }
603
604 indexArray[read_c] = kmer_c;
605 kmer_c += lenBuffer[read_c] - overlaplen + 1;
606 read_c++;
607 }
608
609 if ( read_c > maxReadNum - 1024 )
610 {
611 kmerCounter[0] += kmer_c;
612 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
613 sendWorkSignal ( 1, thrdSignal ); //singleKmer
614 kmer_c = read_c = 0;
615 }
616
617 if ( flag1 == 2 )
618 { break; }
619 }
620 }
621 }
622
623 if ( read_c )
624 {
625 kmerCounter[0] += kmer_c;
626 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
627 sendWorkSignal ( 1, thrdSignal ); //singleKmer
628 }
629
630 time ( &stop_t );
631 fprintf ( stderr, "Time spent on hashing reads: %ds, %lld read(s) processed.\n", ( int ) ( stop_t - start_t ), i );
632
633 //record insert size info
634 if ( pairs )
635 {
636 if ( gradsCounter )
637 { fprintf ( stderr, "%d pe insert size, the largest boundary is %lld.\n\n", gradsCounter, pes[gradsCounter - 1].PE_bound ); }
638 else
639 {
640 fprintf ( stderr, "No paired reads found.\n" );
641 }
642
643 sprintf ( name, "%s.peGrads", outfile );
644 fo = ckopen ( name, "w" );
645 fprintf ( fo, "grads&num: %d\t%lld\n", gradsCounter, n_solexa );
646
647 for ( i = 0; i < gradsCounter; i++ )
648 {
649 fprintf ( fo, "%d\t%lld\t%d\n", pes[i].insertS, pes[i].PE_bound, pes[i].rank );
650 }
651
652 fclose ( fo );
653 }
654
655 free_pe_mem ();
656 free_libs ();
657
658 if ( 1 )
659 {
660 unsigned long long alloCounter = 0;
661 unsigned long long allKmerCounter = 0;
662
663 for ( i = 0; i < thrd_num; i++ )
664 {
665 alloCounter += count_kmerset ( ( KmerSets[i] ) );
666 allKmerCounter += kmerCounter[i + 1];
667 free ( ( void * ) rcSeq[i + 1] );
668 }
669
670 fprintf ( stderr, "%lli node(s) allocated, %lli kmer(s) in reads, %lli kmer(s) processed.\n", alloCounter, kmerCounter[0], allKmerCounter );
671 }
672
673 free ( ( void * ) rcSeq );
674 free ( ( void * ) kmerCounter );
675
676 for ( i = 0; i < maxReadNum; i++ )
677 {
678 free ( ( void * ) seqBuffer[i] );
679 }
680
681 free ( ( void * ) seqBuffer );
682 free ( ( void * ) lenBuffer );
683 free ( ( void * ) indexArray );
684 free ( ( void * ) kmerBuffer );
685 free ( ( void * ) hashBanBuffer );
686 free ( ( void * ) nextcBuffer );
687 free ( ( void * ) prevcBuffer );
688 free ( ( void * ) next_name );
689 free ( ( void * ) aioBuffer1 );
690 free ( ( void * ) aioBuffer2 );
691 free ( ( void * ) readBuffer1 );
692 free ( ( void * ) readBuffer2 );
693 free ( ( void * ) cach1 );
694 free ( ( void * ) cach2 );
695 fprintf ( stderr, "done hashing nodes\n" );
696
697 if ( deLowKmer )
698 {
699 time ( &start_t );
700 deLowCov ( thrdSignal );
701 time ( &stop_t );
702 fprintf ( stderr, "Time spent on delowcvgNode: %ds.\n", ( int ) ( stop_t - start_t ) );
703 }
704
705 time ( &start_t );
706 Mark1in1outNode ( thrdSignal );
707 freqStat ( outfile );
708 time ( &stop_t );
709 fprintf ( stderr, "Time spent on marking linear nodes: %ds.\n", ( int ) ( stop_t - start_t ) );
710 sendWorkSignal ( 3, thrdSignal ); //exit
711 thread_wait ( threads );
712 return 1;
713 }
714
715 void initAIO ( struct aiocb * aio, char * buf, int fd, int size )
716 {
717 bzero ( aio, sizeof ( struct aiocb ) );
718 aio->aio_buf = ( void * ) buf;
719 aio->aio_fildes = fd;
720 aio->aio_nbytes = size;
721 aio->aio_offset = 0;
722 }
723
724 int AIORead ( struct aiocb * mycb, int * offset, char * buf, char * cach, int * rt, int curr_type )
725 {
726 int i, i2, i3, j;
727 int num;
728 size_t mode, get, max_list;
729
730 // rt = aio_read(mycb);
731 if ( *rt == 0 )
732 {
733 struct aiocb * aiocb_list[1];
734 aiocb_list[0] = mycb;
735 max_list = 1;
736
737 while ( 1 )
738 {
739 mode = aio_suspend ( ( const struct aiocb * const * ) aiocb_list, max_list, NULL );
740
741 if ( mode == -1 )
742 {
743 if ( errno != EAGAIN && errno != EINTR )
744 {
745 fprintf ( stderr, "Error:%s.\n", errno );
746 return 0;
747 }
748 else
749 { continue; }
750 }
751 else
752 {
753 //while(aio_error(mycb) == EINPROGRESS);
754 get = aio_return ( mycb );
755 j = strlen ( cach );
756
757 if ( get > 0 )
758 {
759 char * temp = ( char * ) ( ( *mycb ).aio_buf );
760
761 if ( ( get % 32768 ) != 0 )
762 {
763 strcpy ( buf, cach );
764 memcpy ( &buf[j], temp, get );
765 memset ( cach, '\0', j );
766 //printf("%s",buf);
767 *offset = j + get;
768 return 2;
769 }
770
771 if ( ( curr_type == 2 ) || ( curr_type == 6 ) )
772 {
773 num = 0;
774
775 for ( i = get - 1; ( temp[i] != '@' ) || ( temp[i - 1] != '\n' ); i-- )
776 {
777 if ( temp[i] == '\n' ) {num++;}
778 }
779
780 if ( num <= 1 )
781 {
782 for ( i2 = i - 2; temp[i2] != '\n'; i2-- ) { ; }
783
784 if ( temp[i2 + 1] == '+' )
785 {
786 for ( i2 = i2 - 1; temp[i2] != '\n'; i2-- ) { ; }
787
788 if ( temp[i2 + 1] != '+' ) {for ( i = i2 - 1; ( temp[i] != '@' ) || ( temp[i - 1] != '\n' ); i-- ) { ; }}
789 }
790 }
791 }
792 else if ( ( curr_type == 1 ) || ( curr_type == 3 ) || ( curr_type == 5 ) )
793 for ( i = get - 1; temp[i] != '>'; i-- ) { ; }
794
795 //for (i = get - 1; temp[i] != '>' && temp[i] != '@'; i--) ;
796 strcpy ( buf, cach );
797 memcpy ( &buf[j], temp, i );
798 //printf("%s",buf);
799 *offset = i + j;
800 memset ( cach, '\0', j );
801 memcpy ( cach, &temp[i], get - i );
802 ( *mycb ).aio_offset += get;
803 *rt = aio_read ( mycb );
804 return 1;
805 }
806 else
807 {
808 fprintf(stderr, "Error: aio_return error.\n");
809 }
810 /*else
811 {
812 char *temp = (char *)((*mycb).aio_buf);
813 strcpy(buf,cach);
814 strcpy(&buf[j],temp);
815 *offset = j + get;
816 return 2;
817 } */
818 }
819 }
820 }
821 else
822 {
823 fprintf(stderr, "Error: (*rt != 0) in AIORead.\n");
824 }
825
826 return 0;
827 }
828
829 boolean openNextFile ( int * libNo, boolean pairs, unsigned char asm_ctg )
830 {
831 int i = *libNo;
832 int prevLib = i;
833
834 if ( lib_array[i].fp1 )
835 { closeFp1InLab ( i ); }
836
837 if ( lib_array[i].fp2 )
838 { closeFp2InLab ( i ); }
839
840 *libNo = nextValidIndex ( i, pairs, asm_ctg );
841 i = *libNo;
842
843 if ( lib_array[i].rd_len_cutoff > 0 )
844 { maxReadLen = lib_array[i].rd_len_cutoff < maxReadLen4all ? lib_array[i].rd_len_cutoff : maxReadLen4all; }
845 else
846 { maxReadLen = maxReadLen4all; }
847
848 //record insert size info
849 //printf("from lib %d to %d, read %lld to %ld\n",prevLib,i,readNumBack,n_solexa);
850 if ( pairs && i != prevLib )
851 {
852 if ( readNumBack < n_solexa )
853 {
854 pes[gradsCounter].PE_bound = n_solexa;
855 pes[gradsCounter].rank = lib_array[prevLib].rank;
856 pes[gradsCounter].pair_num_cut = lib_array[prevLib].pair_num_cut;
857 pes[gradsCounter++].insertS = lib_array[prevLib].avg_ins;
858 readNumBack = n_solexa;
859 }
860 }
861
862 if ( i >= num_libs )
863 { return 0; }
864
865 openFileInLib ( i );
866 return 1;
867 }
868
869 static void thread_delow ( KmerSet * set, unsigned char thrdID )
870 {
871 int i, in_num, out_num, cvgSingle;
872 int l_cvg, r_cvg;
873 kmer_t * rs;
874 set->iter_ptr = 0;
875
876 while ( set->iter_ptr < set->size )
877 {
878 if ( !is_kmer_entity_null ( set->flags, set->iter_ptr ) )
879 {
880 in_num = out_num = l_cvg = r_cvg = 0;
881 rs = set->array + set->iter_ptr;
882
883 for ( i = 0; i < 4; i++ )
884 {
885 cvgSingle = get_kmer_left_cov ( *rs, i );
886
887 if ( cvgSingle > 0 && cvgSingle <= deLowKmer )
888 {
889 set_kmer_left_cov ( *rs, i, 0 );
890 }
891
892 cvgSingle = get_kmer_right_cov ( *rs, i );
893
894 if ( cvgSingle > 0 && cvgSingle <= deLowKmer )
895 {
896 set_kmer_right_cov ( *rs, i, 0 );
897 }
898 }
899
900 if ( rs->l_links == 0 && rs->r_links == 0 )
901 {
902 rs->deleted = 1;
903 tips[thrdID]++;
904 }
905 }
906
907 set->iter_ptr++;
908 }
909
910 //printf("%lld single nodes, %lld linear\n",counter,tips[thrdID]);
911 }
912
913 static void deLowCov ( unsigned char * thrdSignal )
914 {
915 int i;
916 long long counter = 0;
917 tips = ( long long * ) ckalloc ( thrd_num * sizeof ( long long ) );
918
919 for ( i = 0; i < thrd_num; i++ )
920 {
921 tips[i] = 0;
922 }
923
924 sendWorkSignal ( 5, thrdSignal ); //mark linear nodes
925
926 for ( i = 0; i < thrd_num; i++ )
927 {
928 counter += tips[i];
929 }
930
931 free ( ( void * ) tips );
932 fprintf ( stderr, "%lld kmer(s) removed.\n", counter );
933 }
934
935 static void thread_mark ( KmerSet * set, unsigned char thrdID )
936 {
937 int i, in_num, out_num, cvgSingle;
938 int l_cvg, r_cvg;
939 kmer_t * rs;
940 long long counter = 0;
941 set->iter_ptr = 0;
942
943 while ( set->iter_ptr < set->size )
944 {
945 if ( !is_kmer_entity_null ( set->flags, set->iter_ptr ) )
946 {
947 in_num = out_num = l_cvg = r_cvg = 0;
948 rs = set->array + set->iter_ptr;
949
950 for ( i = 0; i < 4; i++ )
951 {
952 cvgSingle = get_kmer_left_cov ( *rs, i );
953
954 if ( cvgSingle > 0 )
955 {
956 in_num++;
957 l_cvg += cvgSingle;
958 }
959
960 cvgSingle = get_kmer_right_cov ( *rs, i );
961
962 if ( cvgSingle > 0 )
963 {
964 out_num++;
965 r_cvg += cvgSingle;
966 }
967 }
968
969 if ( rs->single )
970 {
971 kmerFreq[thrdID][1]++;
972 counter++;
973 }
974 else
975 {
976 kmerFreq[thrdID][ ( l_cvg > r_cvg ? l_cvg : r_cvg )]++;
977 }
978
979 if ( in_num == 1 && out_num == 1 )
980 {
981 rs->linear = 1;
982 tips[thrdID]++;
983 }
984 }
985
986 set->iter_ptr++;
987 }
988
989 //printf("%lld single nodes, %lld linear\n",counter,tips[thrdID]);
990 }
991
992 static void Mark1in1outNode ( unsigned char * thrdSignal )
993 {
994 int i;
995 long long counter = 0;
996 tips = ( long long * ) ckalloc ( thrd_num * sizeof ( long long ) );
997 kmerFreq = ( long long ** ) ckalloc ( thrd_num * sizeof ( long long * ) );
998
999 for ( i = 0; i < thrd_num; i++ )
1000 {
1001 kmerFreq[i] = ( long long * ) ckalloc ( 257 * sizeof ( long long ) );
1002 memset ( kmerFreq[i], 0, 257 * sizeof ( long long ) );
1003 tips[i] = 0;
1004 }
1005
1006 sendWorkSignal ( 4, thrdSignal ); //mark linear nodes
1007
1008 for ( i = 0; i < thrd_num; i++ )
1009 {
1010 counter += tips[i];
1011 }
1012
1013 free ( ( void * ) tips );
1014 fprintf ( stderr, "%lld linear node(s) marked.\n", counter );
1015 }
1016
1017 static void freqStat ( char * outfile )
1018 {
1019 FILE * fo;
1020 char name[256];
1021 int i, j;
1022 long long sum;
1023 sprintf ( name, "%s.kmerFreq", outfile );
1024 fo = ckopen ( name, "w" );
1025
1026 for ( i = 1; i < 256; i++ )
1027 {
1028 sum = 0;
1029
1030 for ( j = 0; j < thrd_num; j++ )
1031 {
1032 sum += kmerFreq[j][i];
1033 }
1034
1035 fprintf ( fo, "%lld\n", sum );
1036 }
1037
1038 for ( i = 0; i < thrd_num; i++ )
1039 {
1040 free ( ( void * ) kmerFreq[i] );
1041 }
1042
1043 free ( ( void * ) kmerFreq );
1044 fclose ( fo );
1045 }
A place to share codes with mates