?  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 / prlRead2path.c
blobf584e2a53150859278bcd233d88f626d639a4308
1 /*
2 * prlRead2path.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 #define preARCBLOCKSIZE 100000
30 #ifdef MER127
31 static const Kmer kmerZero = { 0, 0, 0, 0 };
32 #else
33 static const Kmer kmerZero = { 0, 0 };
34 #endif
35
36 static unsigned int * arcCounters;
37 static int buffer_size = 100000000;
38 static long long markCounter = 0;
39 static unsigned int * fwriteBuf;
40 static unsigned char * markerOnEdge; //edge occured times for each edge
41
42 //buffer related varibles for chop kmer
43 static int read_c;
44 static char ** rcSeq;
45 static char ** seqBuffer;
46 static int * lenBuffer;
47
48 //edge and (K+1)mer related variables
49 static preARC ** preArc_array;
50 static Kmer * mixBuffer; //kmer buffer; after searching, mixBuffer[j].low = node->l_links; 'node->l_links'
51 static boolean * flagArray; //indicate each item in mixBuffer whether it's a (K+1)mer
52
53 // kmer related variables
54 static char ** flags;
55 static int kmer_c;
56 static Kmer * kmerBuffer;
57 static ubyte8 * hashBanBuffer;
58 static kmer_t ** nodeBuffer; //kmer_t buffer related to 'kmerBuffer'
59 static boolean * smallerBuffer;
60 static int * indexArray;
61
62 static int * deletion; //read deletion number for each thread
63
64 static struct aiocb aio1;
65 static struct aiocb aio2;
66 static char * aioBuffer1;
67 static char * aioBuffer2;
68 static char * readBuffer1;
69 static char * readBuffer2;
70
71 static void parse1read ( int t, int threadID );
72 static void search1kmerPlus ( int j, unsigned char thrdID );
73 static void threadRoutine ( void * thrdID );
74 static void searchKmer ( int t, KmerSet * kset );
75 static void chopKmer4read ( int t, int threadID );
76 static void thread_wait ( pthread_t * threads );
77 static void thread_add1preArc ( unsigned int from_ed, unsigned int to_ed, unsigned int thrdID );
78
79 static void creatThrds ( pthread_t * threads, PARAMETER * paras )
80 {
81 unsigned char i;
82 int temp;
83
84 for ( i = 0; i < thrd_num; i++ )
85 {
86 //printf("to create %dth thread\n",(*(char *)&(threadID[i])));
87 if ( ( temp = pthread_create ( &threads[i], NULL, ( void * ) threadRoutine, & ( paras[i] ) ) ) != 0 )
88 {
89 fprintf ( stderr, "Create threads failed.\n" );
90 exit ( 1 );
91 }
92 }
93
94 fprintf ( stderr, "%d thread(s) initialized.\n", thrd_num );
95 }
96
97 static void threadRoutine ( void * para )
98 {
99 PARAMETER * prm;
100 int i, t, j, start, finish;
101 unsigned char id;
102 prm = ( PARAMETER * ) para;
103 id = prm->threadID;
104
105 //printf("%dth thread with task %d, hash_table %p\n",id,prm.task,prm.hash_table);
106 while ( 1 )
107 {
108 if ( * ( prm->selfSignal ) == 1 )
109 {
110 for ( i = 0; i < kmer_c; i++ )
111 {
112 //if((hashBanBuffer[i]&taskMask)!=prm.threadID)
113 if ( ( hashBanBuffer[i] % thrd_num ) != id )
114 {
115 continue;
116 }
117
118 searchKmer ( i, KmerSets[id] );
119 }
120
121 * ( prm->selfSignal ) = 0;
122 }
123 else if ( * ( prm->selfSignal ) == 2 )
124 {
125 for ( i = 0; i < read_c; i++ )
126 {
127 if ( i % thrd_num != id )
128 {
129 continue;
130 }
131
132 chopKmer4read ( i, id + 1 );
133 }
134
135 * ( prm->selfSignal ) = 0;
136 }
137 else if ( * ( prm->selfSignal ) == 3 )
138 {
139 // parse reads
140 for ( t = 0; t < read_c; t++ )
141 {
142 if ( t % thrd_num != id )
143 {
144 continue;
145 }
146
147 parse1read ( t, id + 1 );
148 }
149
150 * ( prm->selfSignal ) = 0;
151 }
152 else if ( * ( prm->selfSignal ) == 4 )
153 {
154 //printf("thread %d, reads %d splay kmerplus\n",id,read_c);
155 for ( t = 0; t < read_c; t++ )
156 {
157 start = indexArray[t];
158 finish = indexArray[t + 1];
159 #ifdef MER127
160
161 for ( j = start; j < finish; j++ )
162 {
163 if ( flagArray[j] == 0 )
164 {
165 if ( mixBuffer[j].low2 == 0 )
166 {
167 break;
168 }
169 }
170 else if ( hashBanBuffer[j] % thrd_num == id )
171 {
172 //fprintf(stderr,"thread %d search for ban %lld\n",id,hashBanBuffer[j]);
173 search1kmerPlus ( j, id );
174 }
175
176 /*
177 if(flagArray[j]==0&&mixBuffer[j]==0)
178 break;
179 if(!flagArray[j]||(hashBanBuffer[j]%thrd_num)!=id)
180 continue;
181 search1kmerPlus(j,id);
182 */
183 }
184
185 #else
186
187 for ( j = start; j < finish; j++ )
188 {
189 if ( flagArray[j] == 0 )
190 {
191 if ( mixBuffer[j].low == 0 )
192 {
193 break;
194 }
195 }
196 else if ( hashBanBuffer[j] % thrd_num == id )
197 {
198 //fprintf(stderr,"thread %d search for ban %lld\n",id,hashBanBuffer[j]);
199 search1kmerPlus ( j, id );
200 }
201
202 /*
203 if(flagArray[j]==0&&mixBuffer[j]==0)
204 break;
205 if(!flagArray[j]||(hashBanBuffer[j]%thrd_num)!=id)
206 continue;
207 search1kmerPlus(j,id);
208 */
209 }
210
211 #endif
212 }
213
214 * ( prm->selfSignal ) = 0;
215 }
216 else if ( * ( prm->selfSignal ) == 6 )
217 {
218 for ( t = 0; t < read_c; t++ )
219 {
220 start = indexArray[t];
221 finish = indexArray[t + 1];
222 #ifdef MER127
223
224 for ( j = start; j < finish - 1; j++ )
225 {
226 if ( mixBuffer[j].low2 == 0 || mixBuffer[j + 1].low2 == 0 )
227 {
228 break;
229 }
230
231 if ( mixBuffer[j].low2 % thrd_num != id )
232 {
233 continue;
234 }
235
236 thread_add1preArc ( mixBuffer[j].low2, mixBuffer[j + 1].low2, id );
237 }
238
239 #else
240
241 for ( j = start; j < finish - 1; j++ )
242 {
243 if ( mixBuffer[j].low == 0 || mixBuffer[j + 1].low == 0 )
244 {
245 break;
246 }
247
248 if ( mixBuffer[j].low % thrd_num != id )
249 {
250 continue;
251 }
252
253 thread_add1preArc ( mixBuffer[j].low, mixBuffer[j + 1].low, id );
254 }
255
256 #endif
257 }
258
259 * ( prm->selfSignal ) = 0;
260 }
261 else if ( * ( prm->selfSignal ) == 5 )
262 {
263 * ( prm->selfSignal ) = 0;
264 break;
265 }
266
267 usleep ( 1 );
268 }
269 }
270
271 static void chopKmer4read ( int t, int threadID )
272 {
273 char * src_seq = seqBuffer[t];
274 char * bal_seq = rcSeq[threadID];
275 int len_seq = lenBuffer[t];
276 int j, bal_j;
277 ubyte8 hash_ban, bal_hash_ban;
278 Kmer word, bal_word;
279 int index;
280 #ifdef MER127
281 word = kmerZero;
282
283 for ( index = 0; index < overlaplen; index++ )
284 {
285 word = KmerLeftBitMoveBy2 ( word );
286 word.low2 |= src_seq[index];
287 }
288
289 #else
290 word = kmerZero;
291
292 for ( index = 0; index < overlaplen; index++ )
293 {
294 word = KmerLeftBitMoveBy2 ( word );
295 word.low |= src_seq[index];
296 }
297
298 #endif
299 reverseComplementSeq ( src_seq, len_seq, bal_seq );
300 // complementary node
301 bal_word = reverseComplement ( word, overlaplen );
302 bal_j = len_seq - 0 - overlaplen; // 0;
303 index = indexArray[t];
304
305 if ( KmerSmaller ( word, bal_word ) )
306 {
307 hash_ban = hash_kmer ( word );
308 kmerBuffer[index] = word;
309 smallerBuffer[index] = 1;
310 hashBanBuffer[index++] = hash_ban;
311 }
312 else
313 {
314 bal_hash_ban = hash_kmer ( bal_word );
315 kmerBuffer[index] = bal_word;
316 smallerBuffer[index] = 0;
317 hashBanBuffer[index++] = bal_hash_ban;
318 }
319
320 //printf("%dth: %p with %p\n",kmer_c-1,bal_word,bal_hash_ban);
321 for ( j = 1; j <= len_seq - overlaplen; j++ )
322 {
323 word = nextKmer ( word, src_seq[j - 1 + overlaplen] );
324 bal_j = len_seq - j - overlaplen; // j;
325 bal_word = prevKmer ( bal_word, bal_seq[bal_j] );
326
327 if ( KmerSmaller ( word, bal_word ) )
328 {
329 hash_ban = hash_kmer ( word );
330 kmerBuffer[index] = word;
331 smallerBuffer[index] = 1;
332 hashBanBuffer[index++] = hash_ban;
333 //printf("%dth: %p with %p\n",kmer_c-1,word,hashBanBuffer[kmer_c-1]);
334 }
335 else
336 {
337 // complementary node
338 bal_hash_ban = hash_kmer ( bal_word );
339 kmerBuffer[index] = bal_word;
340 smallerBuffer[index] = 0;
341 hashBanBuffer[index++] = bal_hash_ban;
342 //printf("%dth: %p with %p\n",kmer_c-1,bal_word,hashBanBuffer[kmer_c-1]);
343 }
344 }
345 }
346
347 //splay for one kmer in buffer and save the node to nodeBuffer
348 static void searchKmer ( int t, KmerSet * kset )
349 {
350 kmer_t * node;
351 boolean found = search_kmerset ( kset, kmerBuffer[t], &node );
352
353 if ( !found )
354 {
355 fprintf ( stderr, "SearchKmer: kmer " );
356 PrintKmer ( stderr, kmerBuffer[t] );
357 fprintf ( stderr, " is not found.\n" );
358 /*
359 #ifdef MER127
360 fprintf (stderr,"searchKmer: kmer %llx %llx %llx %llx is not found\n", kmerBuffer[t].high1, kmerBuffer[t].low1, kmerBuffer[t].high2, kmerBuffer[t].low2);
361 #else
362 fprintf (stderr,"searchKmer: kmer %llx %llx is not found\n", kmerBuffer[t].high, kmerBuffer[t].low);
363 #endif
364 */
365 }
366
367 nodeBuffer[t] = node;
368 }
369
370 static preARC * getPreArcBetween ( unsigned int from_ed, unsigned int to_ed )
371 {
372 preARC * parc;
373 parc = preArc_array[from_ed];
374
375 while ( parc )
376 {
377 if ( parc->to_ed == to_ed )
378 {
379 return parc;
380 }
381
382 parc = parc->next;
383 }
384
385 return parc;
386 }
387
388 static void thread_add1preArc ( unsigned int from_ed, unsigned int to_ed, unsigned int thrdID )
389 {
390 preARC * parc = getPreArcBetween ( from_ed, to_ed );
391
392 if ( parc )
393 {
394 parc->multiplicity++;
395 }
396 else
397 {
398 parc = prlAllocatePreArc ( to_ed, preArc_mem_managers[thrdID] );
399 arcCounters[thrdID]++;
400 parc->next = preArc_array[from_ed];
401 preArc_array[from_ed] = parc;
402 }
403 }
404
405 static void memoAlloc4preArc ()
406 {
407 unsigned int i;
408 preArc_array = ( preARC ** ) ckalloc ( ( num_ed + 1 ) * sizeof ( preARC * ) );
409
410 for ( i = 0; i <= num_ed; i++ )
411 {
412 preArc_array[i] = NULL;
413 }
414 }
415
416 static void memoFree4preArc ()
417 {
418 prlDestroyPreArcMem ();
419
420 if ( preArc_array )
421 {
422 free ( ( void * ) preArc_array );
423 }
424 }
425
426 static void output_arcs ( char * outfile )
427 {
428 unsigned int i;
429 char name[256];
430 FILE * outfp, *outfp2 = NULL;
431 preARC * parc;
432 sprintf ( name, "%s.preArc", outfile );
433 outfp = ckopen ( name, "w" );
434
435 if ( repsTie )
436 {
437 sprintf ( name, "%s.markOnEdge", outfile );
438 outfp2 = ckopen ( name, "w" );
439 }
440
441 markCounter = 0;
442
443 for ( i = 1; i <= num_ed; i++ )
444 {
445 if ( repsTie )
446 {
447 markCounter += markerOnEdge[i];
448 fprintf ( outfp2, "%d\n", markerOnEdge[i] );
449 }
450
451 parc = preArc_array[i];
452
453 if ( !parc )
454 {
455 continue;
456 }
457
458 fprintf ( outfp, "%u", i );
459
460 while ( parc )
461 {
462 fprintf ( outfp, " %u %u", parc->to_ed, parc->multiplicity );
463 parc = parc->next;
464 }
465
466 fprintf ( outfp, "\n" );
467 }
468
469 fclose ( outfp );
470
471 if ( repsTie )
472 {
473 fclose ( outfp2 );
474 // fprintf (stderr,"%lld marker(s) counted.\n", markCounter);
475 }
476 }
477
478 static void recordPathBin ( FILE * outfp )
479 {
480 int t, j, start, finish;
481 unsigned char counter;
482
483 for ( t = 0; t < read_c; t++ )
484 {
485 start = indexArray[t];
486 finish = indexArray[t + 1];
487 #ifdef MER127
488
489 if ( finish - start < 3 || mixBuffer[start].low2 == 0 || mixBuffer[start + 1].low2 == 0 || mixBuffer[start + 2].low2 == 0 )
490 {
491 continue;
492 }
493
494 counter = 0;
495
496 for ( j = start; j < finish; j++ )
497 {
498 if ( mixBuffer[j].low2 == 0 )
499 {
500 break;
501 }
502
503 fwriteBuf[counter++] = ( unsigned int ) mixBuffer[j].low2;
504
505 if ( markerOnEdge[mixBuffer[j].low2] < 255 )
506 {
507 markerOnEdge[mixBuffer[j].low2]++;
508 }
509
510 markCounter++;
511 }
512
513 #else
514
515 if ( finish - start < 3 || mixBuffer[start].low == 0 || mixBuffer[start + 1].low == 0 || mixBuffer[start + 2].low == 0 )
516 {
517 continue;
518 }
519
520 counter = 0;
521
522 for ( j = start; j < finish; j++ )
523 {
524 if ( mixBuffer[j].low == 0 )
525 {
526 break;
527 }
528
529 fwriteBuf[counter++] = ( unsigned int ) mixBuffer[j].low;
530
531 if ( markerOnEdge[mixBuffer[j].low] < 255 )
532 {
533 markerOnEdge[mixBuffer[j].low]++;
534 }
535
536 markCounter++;
537 }
538
539 #endif
540 fwrite ( &counter, sizeof ( char ), 1, outfp );
541 fwrite ( fwriteBuf, sizeof ( unsigned int ), ( int ) counter, outfp );
542 }
543 }
544
545 /*************************************************
546 Function:
547 search1kmerPlus
548 Description:
549 Searchs (k+1) kmer in hashset.
550 Input:
551 1. j: (k+1) mer's index
552 2. thrdID: thread id
553 Output:
554 None.
555 Return:
556 None.
557 *************************************************/
558 static void search1kmerPlus ( int j, unsigned char thrdID )
559 {
560 kmer_t * node;
561 boolean found = search_kmerset ( KmerSetsPatch[thrdID], mixBuffer[j], &node );
562
563 if ( !found )
564 {
565 /*
566 fprintf(stderr,"kmerPlus %llx %llx (hashban %lld) not found, flag %d!\n",
567 mixBuffer[j].high,mixBuffer[j].low,hashBanBuffer[j],flagArray[j]);
568 */
569 mixBuffer[j] = kmerZero;
570 return;
571 } //else fprintf(stderr,"kmerPlus found\n");
572
573 #ifdef MER127
574
575 if ( smallerBuffer[j] )
576 {
577 mixBuffer[j].low2 = node->l_links;
578 }
579 else
580 {
581 mixBuffer[j].low2 = node->l_links + node->twin - 1;
582 }
583
584 #else
585
586 if ( smallerBuffer[j] )
587 {
588 mixBuffer[j].low = node->l_links;
589 }
590 else
591 {
592 mixBuffer[j].low = node->l_links + node->twin - 1;
593 }
594
595 #endif
596 }
597
598 static void parse1read ( int t, int threadID )
599 {
600 unsigned int j, retain = 0;
601 unsigned int edge_index = 0;
602 kmer_t * node;
603 boolean isSmaller;
604 Kmer wordplus, bal_wordplus;
605 unsigned int start, finish, pos;
606 Kmer prevKmer, currentKmer;
607 boolean IsPrevKmer = 0;
608 start = indexArray[t];
609 finish = indexArray[t + 1];
610 pos = start;
611
612 for ( j = start; j < finish; j++ )
613 {
614 node = nodeBuffer[j];
615
616 //extract edges or keep kmers
617 if ( ( node->deleted ) || ( node->linear && !node->inEdge ) ) // deleted or in a floating loop
618 {
619 if ( retain < 2 )
620 {
621 retain = 0;
622 pos = start;
623 }
624 else
625 {
626 break;
627 }
628
629 continue;
630 }
631
632 isSmaller = smallerBuffer[j];
633
634 if ( node->linear )
635 {
636 if ( isSmaller )
637 {
638 edge_index = node->l_links;
639 }
640 else
641 {
642 edge_index = node->l_links + node->twin - 1;
643 }
644
645 #ifdef MER127
646
647 if ( retain == 0 || IsPrevKmer )
648 {
649 retain++;
650 mixBuffer[pos].low2 = edge_index;
651 flagArray[pos++] = 0;
652 IsPrevKmer = 0;
653 }
654 else if ( edge_index != mixBuffer[pos - 1].low2 )
655 {
656 retain++;
657 mixBuffer[pos].low2 = edge_index;
658 flagArray[pos++] = 0;
659 }
660
661 #else
662
663 if ( retain == 0 || IsPrevKmer )
664 {
665 retain++;
666 mixBuffer[pos].low = edge_index;
667 flagArray[pos++] = 0;
668 IsPrevKmer = 0;
669 }
670 else if ( edge_index != mixBuffer[pos - 1].low )
671 {
672 retain++;
673 mixBuffer[pos].low = edge_index;
674 flagArray[pos++] = 0;
675 }
676
677 #endif
678 }
679 else
680 {
681 if ( isSmaller )
682 {
683 currentKmer = node->seq;
684 }
685 else
686 {
687 currentKmer = reverseComplement ( node->seq, overlaplen );
688 }
689
690 if ( IsPrevKmer )
691 {
692 retain++;
693 wordplus = KmerPlus ( prevKmer, lastCharInKmer ( currentKmer ) );
694 bal_wordplus = reverseComplement ( wordplus, overlaplen + 1 );
695
696 if ( KmerSmaller ( wordplus, bal_wordplus ) )
697 {
698 smallerBuffer[pos] = 1;
699 hashBanBuffer[pos] = hash_kmer ( wordplus );
700 mixBuffer[pos] = wordplus;
701 }
702 else
703 {
704 smallerBuffer[pos] = 0;
705 hashBanBuffer[pos] = hash_kmer ( bal_wordplus );
706 mixBuffer[pos] = bal_wordplus;
707 }
708
709 // fprintf(stderr,"%lld\n",hashBanBuffer[pos]);
710 flagArray[pos++] = 1;
711 }
712
713 IsPrevKmer = 1;
714 prevKmer = currentKmer;
715 }
716 }
717
718 /*
719 for(j=start;j<pos;j++)
720 fprintf(stderr,"%d ",flagArray[j]);
721 fprintf(stderr,"\n");
722 */
723 if ( retain < 1 )
724 {
725 deletion[threadID]++;
726 }
727
728 if ( retain < 2 )
729 {
730 flagArray[start] = 0;
731 mixBuffer[start] = kmerZero;
732 return;
733 }
734
735 if ( ( pos - start ) != retain )
736 {
737 fprintf ( stderr, "Read %d, %d vs %d.\n", t, retain, edge_index - start );
738 }
739
740 if ( pos < finish )
741 {
742 flagArray[pos] = 0;
743 mixBuffer[pos] = kmerZero;
744 }
745 }
746
747 static void sendWorkSignal ( unsigned char SIG, unsigned char * thrdSignals )
748 {
749 int t;
750
751 for ( t = 0; t < thrd_num; t++ )
752 {
753 thrdSignals[t + 1] = SIG;
754 }
755
756 while ( 1 )
757 {
758 usleep ( 10 );
759
760 for ( t = 0; t < thrd_num; t++ )
761 if ( thrdSignals[t + 1] )
762 {
763 break;
764 }
765
766 if ( t == thrd_num )
767 {
768 break;
769 }
770 }
771 }
772
773 /*************************************************
774 Function:
775 prlRead2edge
776 Description:
777 Maps the reads to edges and builds pre-arcs between edges.
778 Input:
779 1. libfile: the reads config file
780 2. outfile: the output file prefix
781 Output:
782 None.
783 Return:
784 None.
785 *************************************************/
786 void prlRead2edge ( char * libfile, char * outfile )
787 {
788 char * cach1;
789 char * cach2;
790 unsigned char asm_ctg = 1;
791 long long i;
792 char name[256], *src_name, *next_name;
793 FILE * outfp = NULL;
794 int maxReadNum, libNo;
795 boolean flag, pairs = 0;
796 pthread_t threads[thrd_num];
797 unsigned char thrdSignal[thrd_num + 1];
798 PARAMETER paras[thrd_num];
799 //init
800 maxReadLen = 0;
801 maxNameLen = 256;
802 scan_libInfo ( libfile );
803 alloc_pe_mem ( num_libs );
804
805 if ( !maxReadLen )
806 {
807 maxReadLen = 100;
808 }
809
810 maxReadLen4all = maxReadLen;
811 fprintf ( stderr, "In file: %s, max seq len %d, max name len %d.\n", libfile, maxReadLen, maxNameLen );
812
813 if ( repsTie )
814 {
815 sprintf ( name, "%s.path", outfile );
816 outfp = ckopen ( name, "wb" );
817 }
818
819 src_name = ( char * ) ckalloc ( ( maxNameLen + 1 ) * sizeof ( char ) );
820 next_name = ( char * ) ckalloc ( ( maxNameLen + 1 ) * sizeof ( char ) );
821 kmerBuffer = ( Kmer * ) ckalloc ( buffer_size * sizeof ( Kmer ) );
822 mixBuffer = ( Kmer * ) ckalloc ( buffer_size * sizeof ( Kmer ) );
823 hashBanBuffer = ( ubyte8 * ) ckalloc ( buffer_size * sizeof ( ubyte8 ) );
824 nodeBuffer = ( kmer_t ** ) ckalloc ( buffer_size * sizeof ( kmer_t * ) );
825 smallerBuffer = ( boolean * ) ckalloc ( buffer_size * sizeof ( boolean ) );
826 flagArray = ( boolean * ) ckalloc ( buffer_size * sizeof ( boolean ) );
827 maxReadNum = buffer_size / ( maxReadLen - overlaplen + 1 );
828 //printf("buffer for at most %d reads\n",maxReadNum);
829 int maxAIOSize = 32768;
830 aioBuffer1 = ( char * ) ckalloc ( ( maxAIOSize ) * sizeof ( char ) );
831 aioBuffer2 = ( char * ) ckalloc ( ( maxAIOSize ) * sizeof ( char ) );
832 readBuffer1 = ( char * ) ckalloc ( ( maxAIOSize + ( maxReadLen * 4 + 1024 ) ) * sizeof ( char ) ); //(char *)ckalloc(maxAIOSize*sizeof(char)); //1024
833 readBuffer2 = ( char * ) ckalloc ( ( maxAIOSize + ( maxReadLen * 4 + 1024 ) ) * sizeof ( char ) ); //1024
834 cach1 = ( char * ) ckalloc ( ( maxReadLen * 4 + 1024 ) * sizeof ( char ) ); //1024
835 cach2 = ( char * ) ckalloc ( ( maxReadLen * 4 + 1024 ) * sizeof ( char ) ); //1024
836 memset ( cach1, '\0', ( maxReadLen * 4 + 1024 ) ); //1024
837 memset ( cach2, '\0', ( maxReadLen * 4 + 1024 ) ); //1024
838 seqBuffer = ( char ** ) ckalloc ( maxReadNum * sizeof ( char * ) );
839 lenBuffer = ( int * ) ckalloc ( maxReadNum * sizeof ( int ) );
840 indexArray = ( int * ) ckalloc ( ( maxReadNum + 1 ) * sizeof ( int ) );
841
842 for ( i = 0; i < maxReadNum; i++ )
843 {
844 seqBuffer[i] = ( char * ) ckalloc ( maxReadLen * sizeof ( char ) );
845 }
846
847 memoAlloc4preArc ();
848 flags = ( char ** ) ckalloc ( ( thrd_num + 1 ) * sizeof ( char * ) );
849 deletion = ( int * ) ckalloc ( ( thrd_num + 1 ) * sizeof ( int ) );
850 rcSeq = ( char ** ) ckalloc ( ( thrd_num + 1 ) * sizeof ( char * ) );
851
852 if ( repsTie )
853 {
854 markerOnEdge = ( unsigned char * ) ckalloc ( ( num_ed + 1 ) * sizeof ( unsigned char ) );
855
856 for ( i = 1; i <= num_ed; i++ )
857 {
858 markerOnEdge[i] = 0;
859 }
860
861 fwriteBuf = ( unsigned int * ) ckalloc ( ( maxReadLen - overlaplen + 1 ) * sizeof ( unsigned int ) );
862 }
863
864 thrdSignal[0] = 0;
865
866 if ( 1 )
867 {
868 preArc_mem_managers = ( MEM_MANAGER ** ) ckalloc ( thrd_num * sizeof ( MEM_MANAGER * ) );
869 arcCounters = ( unsigned int * ) ckalloc ( thrd_num * sizeof ( unsigned int ) );
870
871 for ( i = 0; i < thrd_num; i++ )
872 {
873 arcCounters[i] = 0;
874 preArc_mem_managers[i] = createMem_manager ( preARCBLOCKSIZE, sizeof ( preARC ) );
875 deletion[i + 1] = 0;
876 flags[i + 1] = ( char * ) ckalloc ( 2 * maxReadLen * sizeof ( char ) );
877 rcSeq[i + 1] = ( char * ) ckalloc ( maxReadLen * sizeof ( char ) );
878 thrdSignal[i + 1] = 0;
879 paras[i].threadID = i;
880 paras[i].mainSignal = &thrdSignal[0];
881 paras[i].selfSignal = &thrdSignal[i + 1];
882 }
883
884 creatThrds ( threads, paras );
885 }
886
887 if ( 1 )
888 {
889 deletion[0] = 0;
890 flags[0] = ( char * ) ckalloc ( 2 * maxReadLen * sizeof ( char ) );
891 rcSeq[0] = ( char * ) ckalloc ( maxReadLen * sizeof ( char ) );
892 }
893
894 kmer_c = n_solexa = read_c = i = libNo = readNumBack = gradsCounter = 0;
895 int t0, t1, t2, t3, t4, t5, t6;
896 t0 = t1 = t2 = t3 = t4 = t5 = t6 = 0;
897 time_t read_start, read_end, time_bef, time_aft;
898 time ( &read_start );
899
900 while ( openNextFile ( &libNo, pairs, asm_ctg ) )
901 {
902 if ( lib_array[libNo].curr_type == 4 )
903 {
904 int type = 0; //deside the PE reads is good or bad
905
906 while ( ( flag = read1seqInLibBam ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), &libNo, pairs, 1, &type ) ) != 0 )
907 {
908 if ( type == -1 ) //if the reads is bad, go back.
909 {
910 i--;
911
912 if ( lenBuffer[read_c - 1] >= overlaplen + 1 )
913 {
914 kmer_c -= lenBuffer[read_c - 1] - overlaplen + 1;
915 read_c--;
916 }
917
918 n_solexa -= 2;
919 continue;
920 }
921
922 if ( ( ++i ) % 100000000 == 0 )
923 {
924 fprintf ( stderr, "--- %lldth reads.\n", i );
925 }
926
927 if ( lenBuffer[read_c] < overlaplen + 1 )
928 {
929 continue;
930 }
931
932 //if(lenBuffer[read_c]>70)
933 // lenBuffer[read_c] = 70;
934 //else if(lenBuffer[read_c]>40)
935 // lenBuffer[read_c] = 40;
936 indexArray[read_c] = kmer_c;
937 kmer_c += lenBuffer[read_c] - overlaplen + 1;
938 read_c++;
939
940 if ( read_c == maxReadNum )
941 {
942 indexArray[read_c] = kmer_c;
943 time ( &read_end );
944 t0 += read_end - read_start;
945 time ( &time_bef );
946 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
947 time ( &time_aft );
948 t1 += time_aft - time_bef;
949 time ( &time_bef );
950 sendWorkSignal ( 1, thrdSignal ); //searchKmer
951 time ( &time_aft );
952 t2 += time_aft - time_bef;
953 time ( &time_bef );
954 sendWorkSignal ( 3, thrdSignal ); //parse1read
955 time ( &time_aft );
956 t3 += time_aft - time_bef;
957 time ( &time_bef );
958 sendWorkSignal ( 4, thrdSignal ); //search1kmerPlus
959 time ( &time_aft );
960 t4 += time_aft - time_bef;
961 time ( &time_bef );
962 sendWorkSignal ( 6, thrdSignal ); //thread_add1preArc
963 time ( &time_aft );
964 t5 += time_aft - time_bef;
965 time ( &time_bef );
966
967 //recordPreArc();
968 if ( repsTie )
969 {
970 recordPathBin ( outfp );
971 }
972
973 time ( &time_aft );
974 t6 += time_aft - time_bef;
975 //output_path(read_c,edge_no,flags,outfp);
976 kmer_c = 0;
977 read_c = 0;
978 time ( &read_start );
979 }
980 }
981 }
982 else if ( lib_array[libNo].curr_type == 1 || lib_array[libNo].curr_type == 2 )
983 {
984 initAIO ( &aio1, aioBuffer1, fileno ( lib_array[libNo].fp1 ), maxAIOSize );
985 initAIO ( &aio2, aioBuffer2, fileno ( lib_array[libNo].fp2 ), maxAIOSize );
986 int offset1, offset2, flag1, flag2, rt1, rt2;
987 offset1 = offset2 = 0;
988 rt1 = aio_read ( &aio1 );
989 rt2 = aio_read ( &aio2 );
990 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
991 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
992
993 if ( flag1 && flag2 )
994 {
995 int start1, start2, turn;
996 start1 = start2 = 0;
997 turn = 1;
998
999 while ( start1 < offset1 || start2 < offset2 )
1000 {
1001 if ( turn == 1 )
1002 {
1003 turn = 2;
1004 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer1, &start1, offset1, libNo );
1005
1006 if ( ( ++i ) % 100000000 == 0 )
1007 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
1008
1009 if ( lenBuffer[read_c] < overlaplen + 1 )
1010 {
1011 if ( start1 >= offset1 )
1012 {
1013 start1 = 0;
1014 offset1 = 0;
1015 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
1016 }
1017
1018 continue;
1019 }
1020
1021 indexArray[read_c] = kmer_c;
1022 kmer_c += lenBuffer[read_c] - overlaplen + 1;
1023 read_c++;
1024
1025 if ( start1 >= offset1 )
1026 {
1027 start1 = 0;
1028 offset1 = 0;
1029 flag1 = AIORead ( &aio1, &offset1, readBuffer1, cach1, &rt1, lib_array[libNo].curr_type );
1030 }
1031
1032 if ( read_c == maxReadNum )
1033 {
1034 indexArray[read_c] = kmer_c;
1035 time ( &read_end );
1036 t0 += read_end - read_start;
1037 time ( &time_bef );
1038 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
1039 time ( &time_aft );
1040 t1 += time_aft - time_bef;
1041 time ( &time_bef );
1042 sendWorkSignal ( 1, thrdSignal ); //searchKmer
1043 time ( &time_aft );
1044 t2 += time_aft - time_bef;
1045 time ( &time_bef );
1046 sendWorkSignal ( 3, thrdSignal ); //parse1read
1047 time ( &time_aft );
1048 t3 += time_aft - time_bef;
1049 time ( &time_bef );
1050 sendWorkSignal ( 4, thrdSignal ); //search1kmerPlus
1051 time ( &time_aft );
1052 t4 += time_aft - time_bef;
1053 time ( &time_bef );
1054 sendWorkSignal ( 6, thrdSignal ); //thread_add1preArc
1055 time ( &time_aft );
1056 t5 += time_aft - time_bef;
1057 time ( &time_bef );
1058
1059 //recordPreArc();
1060 if ( repsTie )
1061 { recordPathBin ( outfp ); }
1062
1063 time ( &time_aft );
1064 t6 += time_aft - time_bef;
1065 //output_path(read_c,edge_no,flags,outfp);
1066 kmer_c = 0;
1067 read_c = 0;
1068 time ( &read_start );
1069 }
1070
1071 continue;
1072 }
1073
1074 if ( turn == 2 )
1075 {
1076 turn = 1;
1077 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer2, &start2, offset2, libNo );
1078
1079 if ( ( ++i ) % 100000000 == 0 )
1080 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
1081
1082 if ( lenBuffer[read_c] < overlaplen + 1 )
1083 {
1084 if ( ( flag2 == 2 ) && ( start2 >= offset2 ) )
1085 { break; }
1086
1087 if ( start2 >= offset2 )
1088 {
1089 start2 = 0;
1090 offset2 = 0;
1091 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
1092 }
1093
1094 continue;
1095 }
1096
1097 indexArray[read_c] = kmer_c;
1098 kmer_c += lenBuffer[read_c] - overlaplen + 1;
1099 read_c++;
1100
1101 if ( ( flag2 == 2 ) && ( start2 >= offset2 ) )
1102 { break; }
1103
1104 if ( start2 >= offset2 )
1105 {
1106 start2 = 0;
1107 offset2 = 0;
1108 flag2 = AIORead ( &aio2, &offset2, readBuffer2, cach2, &rt2, lib_array[libNo].curr_type );
1109 }
1110
1111 if ( read_c == maxReadNum )
1112 {
1113 indexArray[read_c] = kmer_c;
1114 time ( &read_end );
1115 t0 += read_end - read_start;
1116 time ( &time_bef );
1117 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
1118 time ( &time_aft );
1119 t1 += time_aft - time_bef;
1120 time ( &time_bef );
1121 sendWorkSignal ( 1, thrdSignal ); //searchKmer
1122 time ( &time_aft );
1123 t2 += time_aft - time_bef;
1124 time ( &time_bef );
1125 sendWorkSignal ( 3, thrdSignal ); //parse1read
1126 time ( &time_aft );
1127 t3 += time_aft - time_bef;
1128 time ( &time_bef );
1129 sendWorkSignal ( 4, thrdSignal ); //search1kmerPlus
1130 time ( &time_aft );
1131 t4 += time_aft - time_bef;
1132 time ( &time_bef );
1133 sendWorkSignal ( 6, thrdSignal ); //thread_add1preArc
1134 time ( &time_aft );
1135 t5 += time_aft - time_bef;
1136 time ( &time_bef );
1137
1138 //recordPreArc();
1139 if ( repsTie )
1140 { recordPathBin ( outfp ); }
1141
1142 time ( &time_aft );
1143 t6 += time_aft - time_bef;
1144 //output_path(read_c,edge_no,flags,outfp);
1145 kmer_c = 0;
1146 read_c = 0;
1147 time ( &read_start );
1148 }
1149
1150 continue;
1151 }
1152 }
1153 }
1154 else
1155 {
1156 fprintf(stderr, "Error: aio_read error.\n");
1157 }
1158 }
1159 else
1160 {
1161 initAIO ( &aio1, aioBuffer1, fileno ( lib_array[libNo].fp1 ), maxAIOSize );
1162 int offset, flag1, rt;
1163 offset = 0;
1164 rt = aio_read ( &aio1 );
1165
1166 while ( ( flag1 = AIORead ( &aio1, &offset, readBuffer1, cach1, &rt, lib_array[libNo].curr_type ) ) )
1167 {
1168 int start = 0;
1169
1170 while ( start < offset )
1171 {
1172 readseqInLib ( seqBuffer[read_c], next_name, & ( lenBuffer[read_c] ), readBuffer1, &start, offset, libNo );
1173
1174 if ( ( ++i ) % 100000000 == 0 )
1175 { fprintf ( stderr, "--- %lldth reads.\n", i ); }
1176
1177 if ( lenBuffer[read_c] < overlaplen + 1 )
1178 { continue; }
1179
1180 indexArray[read_c] = kmer_c;
1181 kmer_c += lenBuffer[read_c] - overlaplen + 1;
1182 read_c++;
1183 }
1184
1185 if ( read_c > maxReadNum - 1024 )
1186 {
1187 indexArray[read_c] = kmer_c;
1188 time ( &read_end );
1189 t0 += read_end - read_start;
1190 time ( &time_bef );
1191 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
1192 time ( &time_aft );
1193 t1 += time_aft - time_bef;
1194 time ( &time_bef );
1195 sendWorkSignal ( 1, thrdSignal ); //searchKmer
1196 time ( &time_aft );
1197 t2 += time_aft - time_bef;
1198 time ( &time_bef );
1199 sendWorkSignal ( 3, thrdSignal ); //parse1read
1200 time ( &time_aft );
1201 t3 += time_aft - time_bef;
1202 time ( &time_bef );
1203 sendWorkSignal ( 4, thrdSignal ); //search1kmerPlus
1204 time ( &time_aft );
1205 t4 += time_aft - time_bef;
1206 time ( &time_bef );
1207 sendWorkSignal ( 6, thrdSignal ); //thread_add1preArc
1208 time ( &time_aft );
1209 t5 += time_aft - time_bef;
1210 time ( &time_bef );
1211
1212 //recordPreArc();
1213 if ( repsTie )
1214 { recordPathBin ( outfp ); }
1215
1216 time ( &time_aft );
1217 t6 += time_aft - time_bef;
1218 //output_path(read_c,edge_no,flags,outfp);
1219 kmer_c = 0;
1220 read_c = 0;
1221 time ( &read_start );
1222 }
1223
1224 if ( flag1 == 2 )
1225 { break; }
1226 }
1227 }
1228 }
1229
1230 fprintf ( stderr, "%lld read(s) processed.\n", i );
1231 // fprintf (stderr,"Time ReadingReads: %d,chopKmer4read: %d,searchKmer: %d,parse1read: %d,search1kmerPlus: %d,thread_add1preArc: %d,recordPathBin: %d\n", t0, t1, t2, t3, t4, t5, t6);
1232 fprintf ( stderr, "Time spent on:\n" );
1233 fprintf ( stderr, " importing reads: %ds,\n", t0 );
1234 fprintf ( stderr, " chopping reads to kmers: %ds,\n", t1 );
1235 fprintf ( stderr, " searching kmers: %ds,\n", t2 );
1236 fprintf ( stderr, " aligning reads to edges: %ds,\n", t3 );
1237 fprintf ( stderr, " searching (K+1)mers: %ds,\n", t4 );
1238 fprintf ( stderr, " adding pre-arcs: %ds,\n", t5 );
1239 fprintf ( stderr, " recording read paths: %ds.\n", t6 );
1240
1241 if ( read_c )
1242 {
1243 indexArray[read_c] = kmer_c;
1244 sendWorkSignal ( 2, thrdSignal ); //chopKmer4read
1245 sendWorkSignal ( 1, thrdSignal ); //searchKmer
1246 sendWorkSignal ( 3, thrdSignal ); //parse1read
1247 sendWorkSignal ( 4, thrdSignal ); //search1kmerPlus
1248 sendWorkSignal ( 6, thrdSignal ); //thread_add1preArc
1249
1250 //recordPreArc();
1251 if ( repsTie )
1252 {
1253 recordPathBin ( outfp );
1254 }
1255 }
1256
1257 fprintf ( stderr, "%lld marker(s) output.\n", markCounter );
1258 sendWorkSignal ( 5, thrdSignal ); //over
1259 thread_wait ( threads );
1260 output_arcs ( outfile );
1261 memoFree4preArc ();
1262
1263 if ( 1 ) // multi-threads
1264 {
1265 arcCounter = 0;
1266
1267 for ( i = 0; i < thrd_num; i++ )
1268 {
1269 arcCounter += arcCounters[i];
1270 free ( ( void * ) flags[i + 1] );
1271 deletion[0] += deletion[i + 1];
1272 free ( ( void * ) rcSeq[i + 1] );
1273 }
1274 }
1275
1276 if ( 1 )
1277 {
1278 free ( ( void * ) flags[0] );
1279 free ( ( void * ) rcSeq[0] );
1280 }
1281
1282 fprintf ( stderr, "Reads alignment done, %d read(s) deleted, %lld pre-arc(s) added.\n", deletion[0], arcCounter );
1283
1284 if ( repsTie )
1285 {
1286 free ( ( void * ) markerOnEdge );
1287 free ( ( void * ) fwriteBuf );
1288 }
1289
1290 free ( ( void * ) arcCounters );
1291 free ( ( void * ) rcSeq );
1292
1293 for ( i = 0; i < maxReadNum; i++ )
1294 {
1295 free ( ( void * ) seqBuffer[i] );
1296 }
1297
1298 free ( ( void * ) seqBuffer );
1299 free ( ( void * ) lenBuffer );
1300 free ( ( void * ) indexArray );
1301 free ( ( void * ) flags );
1302 free ( ( void * ) deletion );
1303 free ( ( void * ) kmerBuffer );
1304 free ( ( void * ) mixBuffer );
1305 free ( ( void * ) smallerBuffer );
1306 free ( ( void * ) flagArray );
1307 free ( ( void * ) hashBanBuffer );
1308 free ( ( void * ) nodeBuffer );
1309 free ( ( void * ) src_name );
1310 free ( ( void * ) next_name );
1311 free ( ( void * ) aioBuffer1 );
1312 free ( ( void * ) aioBuffer2 );
1313 free ( ( void * ) readBuffer1 );
1314 free ( ( void * ) readBuffer2 );
1315 free ( ( void * ) cach1 );
1316 free ( ( void * ) cach2 );
1317
1318 if ( gLineLen < maxReadLen )
1319 {
1320 free ( ( void * ) gStr );
1321 gStr = NULL;
1322 }
1323
1324 if ( repsTie )
1325 {
1326 fclose ( outfp );
1327 }
1328
1329 free_pe_mem ();
1330 free_libs ();
1331 }
1332
1333 static void thread_wait ( pthread_t * threads )
1334 {
1335 int i;
1336
1337 for ( i = 0; i < thrd_num; i++ )
1338 if ( threads[i] != 0 )
1339 {
1340 pthread_join ( threads[i], NULL );
1341 }
1342 }
A place to share codes with mates