?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
modified: nfig1.py
[GalaxyCodeBases.git] / BGI / SOAPdenovo2 / standardPregraph / lib.c
blob95e4eb753f8d49add71c47ec7e4ef617b9dc5ca6
1 /*
2 * lib.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 static char tabs[2][1024]; //for splitColumn()
30
31 /*************************************************
32 Function:
33 getMaxLongReadLen
34 Description:
35 Get the max length for long reads (asm_flags=4) in lib.
36 Input:
37 1. num_libs: number of lib
38 Output:
39 None.
40 Return:
41 Max read length.
42 *************************************************/
43 int getMaxLongReadLen ( int num_libs )
44 {
45 int i;
46 int maxLong = 0;
47 boolean Has = 0;
48
49 for ( i = 0; i < num_libs; i++ )
50 {
51 if ( lib_array[i].asm_flag != 4 )
52 {
53 continue;
54 }
55
56 Has = 1;
57 maxLong = maxLong < lib_array[i].rd_len_cutoff ? lib_array[i].rd_len_cutoff : maxLong;
58 }
59
60 if ( !Has )
61 {
62 return maxLong;
63 }
64 else
65 {
66 return maxLong > 0 ? maxLong : maxReadLen;
67 }
68 }
69
70 static boolean splitColumn ( char * line )
71 {
72 int len = strlen ( line );
73 int i = 0, j;
74 int tabs_n = 0;
75
76 while ( i < len )
77 {
78 if ( line[i] >= 32 && line[i] <= 126 && line[i] != '=' )
79 {
80 j = 0;
81
82 while ( i < len && line[i] >= 32 && line[i] <= 126 && line[i] != '=' )
83 {
84 tabs[tabs_n][j++] = line[i];
85 i++;
86 }
87
88 tabs[tabs_n][j] = '\0';
89 tabs_n++;
90
91 if ( tabs_n == 2 )
92 {
93 return 1;
94 }
95 }
96
97 i++;
98 }
99
100 if ( tabs_n == 2 )
101 {
102 return 1;
103 }
104 else
105 {
106 return 0;
107 }
108 }
109
110 static int cmp_lib ( const void * a, const void * b )
111 {
112 LIB_INFO * A, *B;
113 A = ( LIB_INFO * ) a;
114 B = ( LIB_INFO * ) b;
115
116 if ( A->avg_ins > B->avg_ins )
117 {
118 return 1;
119 }
120 else if ( A->avg_ins == B->avg_ins )
121 {
122 return 0;
123 }
124 else
125 {
126 return -1;
127 }
128 }
129
130 void scan_libInfo ( char * libfile )
131 {
132 FILE * fp;
133 char line[1024], ch;
134 int i, j, index;
135 int libCounter;
136 boolean flag;
137 boolean * pe;
138 fp = ckopen ( libfile, "r" );
139 num_libs = 0;
140
141 while ( fgets ( line, 1024, fp ) )
142 {
143 ch = line[5];
144 line[5] = '\0';
145
146 if ( strcmp ( line, "[LIB]" ) == 0 )
147 {
148 num_libs++;
149 }
150
151 if ( !num_libs )
152 {
153 line[5] = ch;
154 flag = splitColumn ( line );
155
156 if ( !flag )
157 {
158 continue;
159 }
160
161 if ( strcmp ( tabs[0], "max_rd_len" ) == 0 )
162 {
163 maxReadLen = atoi ( tabs[1] );
164 }
165 }
166 }
167
168 if ( num_libs == 0 )
169 {
170 fprintf ( stderr, "Config file error: no [LIB] in file\n" );
171 exit ( -1 );
172 }
173
174 //count file numbers of each type
175 lib_array = ( LIB_INFO * ) ckalloc ( num_libs * sizeof ( LIB_INFO ) );
176 pe = ( boolean * ) ckalloc ( num_libs * sizeof ( boolean ) );
177
178 for ( i = 0; i < num_libs; i++ )
179 {
180 lib_array[i].asm_flag = 3;
181 lib_array[i].rank = 0;
182 lib_array[i].pair_num_cut = 0;
183 lib_array[i].rd_len_cutoff = 0;
184 lib_array[i].map_len = 0;
185 lib_array[i].num_s_a_file = 0;
186 lib_array[i].num_s_q_file = 0;
187 lib_array[i].num_p_file = 0;
188 lib_array[i].num_a1_file = 0;
189 lib_array[i].num_a2_file = 0;
190 lib_array[i].num_q1_file = 0;
191 lib_array[i].num_q2_file = 0;
192 lib_array[i].num_b_file = 0; //init
193 pe[i] = false;
194 }
195
196 libCounter = -1;
197 rewind ( fp );
198 i = -1;
199
200 while ( fgets ( line, 1024, fp ) )
201 {
202 ch = line[5];
203 line[5] = '\0';
204
205 if ( strcmp ( line, "[LIB]" ) == 0 )
206 {
207 i++;
208 continue;
209 }
210
211 line[5] = ch;
212 flag = splitColumn ( line );
213
214 if ( !flag )
215 {
216 continue;
217 }
218
219 if ( strcmp ( tabs[0], "f1" ) == 0 )
220 {
221 lib_array[i].num_a1_file++;
222 pe[i] = true;
223 }
224 else if ( strcmp ( tabs[0], "q1" ) == 0 )
225 {
226 lib_array[i].num_q1_file++;
227 pe[i] = true;
228 }
229 else if ( strcmp ( tabs[0], "f2" ) == 0 )
230 {
231 lib_array[i].num_a2_file++;
232 pe[i] = true;
233 }
234 else if ( strcmp ( tabs[0], "q2" ) == 0 )
235 {
236 lib_array[i].num_q2_file++;
237 pe[i] = true;
238 }
239 else if ( strcmp ( tabs[0], "f" ) == 0 )
240 {
241 lib_array[i].num_s_a_file++;
242 }
243 else if ( strcmp ( tabs[0], "q" ) == 0 )
244 {
245 lib_array[i].num_s_q_file++;
246 }
247 else if ( strcmp ( tabs[0], "p" ) == 0 )
248 {
249 lib_array[i].num_p_file++;
250 pe[i] = true;
251 }
252 else if ( strcmp ( tabs[0], "b" ) == 0 ) // the bam file
253 {
254 lib_array[i].num_b_file++;
255 pe[i] = true;
256 }
257 }
258
259 //allocate memory for filenames
260 for ( i = 0; i < num_libs; i++ )
261 {
262 if ( lib_array[i].num_a2_file != lib_array[i].num_a1_file )
263 {
264 fprintf ( stderr, "Config file error: the number of mark \"f1\" is not the same as \"f2\"!\n" );
265 exit ( -1 );
266 }
267
268 if ( lib_array[i].num_q2_file != lib_array[i].num_q1_file )
269 {
270 fprintf ( stderr, "Config file error: the number of mark \"q1\" is not the same as \"q2\"!\n" );
271 exit ( -1 );
272 }
273
274 if ( lib_array[i].num_s_a_file )
275 {
276 lib_array[i].s_a_fname = ( char ** ) ckalloc ( lib_array[i].num_s_a_file * sizeof ( char * ) );
277
278 for ( j = 0; j < lib_array[i].num_s_a_file; j++ )
279 {
280 lib_array[i].s_a_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
281 }
282 }
283
284 if ( lib_array[i].num_s_q_file )
285 {
286 lib_array[i].s_q_fname = ( char ** ) ckalloc ( lib_array[i].num_s_q_file * sizeof ( char * ) );
287
288 for ( j = 0; j < lib_array[i].num_s_q_file; j++ )
289 {
290 lib_array[i].s_q_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
291 }
292 }
293
294 if ( lib_array[i].num_p_file )
295 {
296 lib_array[i].p_fname = ( char ** ) ckalloc ( lib_array[i].num_p_file * sizeof ( char * ) );
297
298 for ( j = 0; j < lib_array[i].num_p_file; j++ )
299 {
300 lib_array[i].p_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
301 }
302 }
303
304 if ( lib_array[i].num_a1_file )
305 {
306 lib_array[i].a1_fname = ( char ** ) ckalloc ( lib_array[i].num_a1_file * sizeof ( char * ) );
307
308 for ( j = 0; j < lib_array[i].num_a1_file; j++ )
309 {
310 lib_array[i].a1_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
311 }
312 }
313
314 if ( lib_array[i].num_a2_file )
315 {
316 lib_array[i].a2_fname = ( char ** ) ckalloc ( lib_array[i].num_a2_file * sizeof ( char * ) );
317
318 for ( j = 0; j < lib_array[i].num_a2_file; j++ )
319 {
320 lib_array[i].a2_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
321 }
322 }
323
324 if ( lib_array[i].num_q1_file )
325 {
326 lib_array[i].q1_fname = ( char ** ) ckalloc ( lib_array[i].num_q1_file * sizeof ( char * ) );
327
328 for ( j = 0; j < lib_array[i].num_q1_file; j++ )
329 {
330 lib_array[i].q1_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
331 }
332 }
333
334 if ( lib_array[i].num_q2_file )
335 {
336 lib_array[i].q2_fname = ( char ** ) ckalloc ( lib_array[i].num_q2_file * sizeof ( char * ) );
337
338 for ( j = 0; j < lib_array[i].num_q2_file; j++ )
339 {
340 lib_array[i].q2_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) );
341 }
342 }
343
344 if ( lib_array[i].num_b_file ) //allot memory for bam file name
345 {
346 lib_array[i].b_fname = ( char ** ) ckalloc ( lib_array[i].num_b_file * sizeof ( char * ) );
347
348 for ( j = 0; j < lib_array[i].num_b_file; j++ )
349 { lib_array[i].b_fname[j] = ( char * ) ckalloc ( 1024 * sizeof ( char ) ); }
350 }
351 }
352
353 // get file names
354 for ( i = 0; i < num_libs; i++ )
355 {
356 lib_array[i].curr_type = 1;
357 lib_array[i].curr_index = 0;
358 lib_array[i].fp1 = NULL;
359 lib_array[i].fp2 = NULL;
360 lib_array[i].num_s_a_file = 0;
361 lib_array[i].num_s_q_file = 0;
362 lib_array[i].num_p_file = 0;
363 lib_array[i].num_a1_file = 0;
364 lib_array[i].num_a2_file = 0;
365 lib_array[i].num_q1_file = 0;
366 lib_array[i].num_q2_file = 0;
367 lib_array[i].num_b_file = 0; //init
368 lib_array[i].fp3 = NULL;
369 }
370
371 libCounter = -1;
372 rewind ( fp );
373 i = -1;
374
375 while ( fgets ( line, 1024, fp ) )
376 {
377 ch = line[5];
378 line[5] = '\0';
379
380 if ( strcmp ( line, "[LIB]" ) == 0 )
381 {
382 i++;
383 continue;
384 }
385
386 line[5] = ch;
387 flag = splitColumn ( line );
388
389 if ( !flag )
390 {
391 continue;
392 }
393
394 if ( strcmp ( tabs[0], "f1" ) == 0 )
395 {
396 index = lib_array[i].num_a1_file++;
397 strcpy ( lib_array[i].a1_fname[index], tabs[1] );
398 }
399 else if ( strcmp ( tabs[0], "q1" ) == 0 )
400 {
401 index = lib_array[i].num_q1_file++;
402 strcpy ( lib_array[i].q1_fname[index], tabs[1] );
403 }
404 else if ( strcmp ( tabs[0], "f2" ) == 0 )
405 {
406 index = lib_array[i].num_a2_file++;
407 strcpy ( lib_array[i].a2_fname[index], tabs[1] );
408
409 if ( strcmp ( lib_array[i].a2_fname[index], lib_array[i].a1_fname[index] ) == 0 )
410 {
411 fprintf ( stderr, "Config file error: f2 file is the same as f1 file\n" );
412 fprintf ( stderr, "f1=%s\n", lib_array[i].a1_fname[index] );
413 fprintf ( stderr, "f2=%s\n", lib_array[i].a2_fname[index] );
414 exit ( -1 );
415 }
416 }
417 else if ( strcmp ( tabs[0], "q2" ) == 0 )
418 {
419 index = lib_array[i].num_q2_file++;
420 strcpy ( lib_array[i].q2_fname[index], tabs[1] );
421
422 if ( strcmp ( lib_array[i].q2_fname[index], lib_array[i].q1_fname[index] ) == 0 )
423 {
424 fprintf ( stderr, "Config file error: q2 file is the same as q1 file\n" );
425 fprintf ( stderr, "q1=%s\n", lib_array[i].q1_fname[index] );
426 fprintf ( stderr, "q2=%s\n", lib_array[i].q2_fname[index] );
427 exit ( -1 );
428 }
429 }
430 else if ( strcmp ( tabs[0], "f" ) == 0 )
431 {
432 index = lib_array[i].num_s_a_file++;
433 strcpy ( lib_array[i].s_a_fname[index], tabs[1] );
434 }
435 else if ( strcmp ( tabs[0], "q" ) == 0 )
436 {
437 index = lib_array[i].num_s_q_file++;
438 strcpy ( lib_array[i].s_q_fname[index], tabs[1] );
439 }
440 else if ( strcmp ( tabs[0], "p" ) == 0 )
441 {
442 index = lib_array[i].num_p_file++;
443 strcpy ( lib_array[i].p_fname[index], tabs[1] );
444 }
445 else if ( strcmp ( tabs[0], "b" ) == 0 )
446 {
447 //bam file
448 index = lib_array[i].num_b_file++;
449 strcpy ( lib_array[i].b_fname[index], tabs[1] );
450 }
451 else if ( strcmp ( tabs[0], "min_ins" ) == 0 )
452 {
453 lib_array[i].min_ins = atoi ( tabs[1] );
454 }
455 else if ( strcmp ( tabs[0], "max_ins" ) == 0 )
456 {
457 lib_array[i].max_ins = atoi ( tabs[1] );
458 }
459 else if ( strcmp ( tabs[0], "avg_ins" ) == 0 )
460 {
461 lib_array[i].avg_ins = atoi ( tabs[1] );
462 }
463 else if ( strcmp ( tabs[0], "rd_len_cutoff" ) == 0 )
464 {
465 lib_array[i].rd_len_cutoff = atoi ( tabs[1] );
466 }
467 else if ( strcmp ( tabs[0], "reverse_seq" ) == 0 )
468 {
469 lib_array[i].reverse = atoi ( tabs[1] );
470 }
471 else if ( strcmp ( tabs[0], "asm_flags" ) == 0 )
472 {
473 lib_array[i].asm_flag = atoi ( tabs[1] );
474 }
475 else if ( strcmp ( tabs[0], "rank" ) == 0 )
476 {
477 lib_array[i].rank = atoi ( tabs[1] );
478 }
479 else if ( strcmp ( tabs[0], "pair_num_cutoff" ) == 0 )
480 {
481 lib_array[i].pair_num_cut = atoi ( tabs[1] );
482 }
483 else if ( strcmp ( tabs[0], "rd_len_cutoff" ) == 0 )
484 {
485 lib_array[i].rd_len_cutoff = atoi ( tabs[1] );
486 }
487 else if ( strcmp ( tabs[0], "map_len" ) == 0 )
488 {
489 lib_array[i].map_len = atoi ( tabs[1] );
490 }
491 }
492
493 for ( i = 0; i < num_libs; i++ )
494 {
495 if ( pe[i] && lib_array[i].avg_ins == 0 )
496 {
497 fprintf ( stderr, "Config file error: PE reads need avg_ins in [LIB] %d\n", i + 1 );
498 exit ( -1 );
499 }
500 }
501
502 fclose ( fp );
503 qsort ( &lib_array[0], num_libs, sizeof ( LIB_INFO ), cmp_lib );
504 }
505
506 void free_libs ()
507 {
508 if ( !lib_array )
509 {
510 return;
511 }
512
513 int i, j;
514 fprintf ( stderr, "LIB(s) information:\n" );
515
516 for ( i = 0; i < num_libs; i++ )
517 {
518 fprintf ( stderr, " [LIB] %d, avg_ins %d, reverse %d.\n", i, lib_array[i].avg_ins, lib_array[i].reverse );
519
520 if ( lib_array[i].num_s_a_file )
521 {
522 //printf("%d single fasta files\n",lib_array[i].num_s_a_file);
523 for ( j = 0; j < lib_array[i].num_s_a_file; j++ )
524 {
525 free ( ( void * ) lib_array[i].s_a_fname[j] );
526 }
527
528 free ( ( void * ) lib_array[i].s_a_fname );
529 }
530
531 if ( lib_array[i].num_s_q_file )
532 {
533 //printf("%d single fastq files\n",lib_array[i].num_s_q_file);
534 for ( j = 0; j < lib_array[i].num_s_q_file; j++ )
535 {
536 free ( ( void * ) lib_array[i].s_q_fname[j] );
537 }
538
539 free ( ( void * ) lib_array[i].s_q_fname );
540 }
541
542 if ( lib_array[i].num_p_file )
543 {
544 //printf("%d paired fasta files\n",lib_array[i].num_p_file);
545 for ( j = 0; j < lib_array[i].num_p_file; j++ )
546 {
547 free ( ( void * ) lib_array[i].p_fname[j] );
548 }
549
550 free ( ( void * ) lib_array[i].p_fname );
551 }
552
553 if ( lib_array[i].num_a1_file )
554 {
555 //printf("%d read1 fasta files\n",lib_array[i].num_a1_file);
556 for ( j = 0; j < lib_array[i].num_a1_file; j++ )
557 {
558 free ( ( void * ) lib_array[i].a1_fname[j] );
559 }
560
561 free ( ( void * ) lib_array[i].a1_fname );
562 }
563
564 if ( lib_array[i].num_a2_file )
565 {
566 //printf("%d read2 fasta files\n",lib_array[i].num_a2_file);
567 for ( j = 0; j < lib_array[i].num_a2_file; j++ )
568 {
569 free ( ( void * ) lib_array[i].a2_fname[j] );
570 }
571
572 free ( ( void * ) lib_array[i].a2_fname );
573 }
574
575 if ( lib_array[i].num_q1_file )
576 {
577 //printf("%d read1 fastq files\n",lib_array[i].num_q1_file);
578 for ( j = 0; j < lib_array[i].num_q1_file; j++ )
579 {
580 free ( ( void * ) lib_array[i].q1_fname[j] );
581 }
582
583 free ( ( void * ) lib_array[i].q1_fname );
584 }
585
586 if ( lib_array[i].num_q2_file )
587 {
588 //printf("%d read2 fastq files\n",lib_array[i].num_q2_file);
589 for ( j = 0; j < lib_array[i].num_q2_file; j++ )
590 {
591 free ( ( void * ) lib_array[i].q2_fname[j] );
592 }
593
594 free ( ( void * ) lib_array[i].q2_fname );
595 }
596
597 if ( lib_array[i].num_b_file )
598 {
599 //free the bam file name
600 //printf("%d bam files\n",lib_array[i].num_b_file);
601 for ( j = 0; j < lib_array[i].num_b_file; j++ )
602 { free ( ( void * ) lib_array[i].b_fname[j] ); }
603
604 free ( ( void * ) lib_array[i].b_fname );
605 }
606 }
607
608 num_libs = 0;
609 free ( ( void * ) lib_array );
610 }
611
612 void alloc_pe_mem ( int gradsCounter )
613 {
614 if ( gradsCounter )
615 {
616 pes = ( PE_INFO * ) ckalloc ( gradsCounter * sizeof ( PE_INFO ) );
617 }
618 }
619
620 void free_pe_mem ()
621 {
622 if ( pes )
623 {
624 free ( ( void * ) pes );
625 pes = NULL;
626 }
627 }
A place to share codes with mates