Context Navigation

source: tags/arb-6.0/GDE/SUPPORT/CAP2.c

Visit:

Last change on this file was 6141, checked in by westram, 15 years ago
spellchecked all (phew)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 58.7 KB

Line
1	/* CONTIG ASSEMBLY PROGRAM (CAP)
2
3	copyright (c) 1991 Xiaoqiu Huang
4	The distribution of the program is granted provided no charge
5	is made and the copyright notice is included.
6
7	Proper attribution of the author as the source of the software
8	would be appreciated:
9	"A Contig Assembly Program Based on Sensitive Detection of
10	Fragment Overlaps" (submitted to Genomics, 1991)
11	Xiaoqiu Huang
12	Department of Computer Science
13	Michigan Technological University
14	Houghton, MI 49931
15	E-mail: huang@cs.mtu.edu
16
17	The CAP program uses a dynamic programming algorithm to compute
18	the maximal-scoring overlapping alignment between two fragments.
19	Fragments in random orientations are assembled into contigs by a
20	greedy approach in order of the overlap scores. CAP is efficient
21	in computer memory: a large number of arbitrarily long fragments
22	can be assembled. The time requirement is acceptable; for example,
23	CAP took 4 hours to assemble 1015 fragments of a total of 252 kb
24	nucleotides on a Sun SPARCstation SLC. The program is written in C
25	and runs on Sun workstations.
26
27	Below is a description of the parameters in the #define section of CAP.
28	Two specially chosen sets of substitution scores and indel penalties
29	are used by the dynamic programming algorithm: heavy set for regions
30	of low sequencing error rates and light set for fragment ends of high
31	sequencing error rates. (Use integers only.)
32	Heavy set: Light set:
33
34	MATCH = 2 MATCH = 2
35	MISMAT = -6 LTMISM = -3
36	EXTEND = 4 LTEXTEN = 2
37
38	In the initial assembly, any overlap must be of length at least OVERLEN,
39	and any overlap/containment must be of identity percentage at least
40	PERCENT. After the initial assembly, the program attempts to join
41	contigs together using weak overlaps. Two contigs are merged if the
42	score of the overlapping alignment is at least CUTOFF. The value for
43	CUTOFF is chosen according to the value for MATCH.
44
45	DELTA is a parameter in necessary conditions for overlap/containment.
46	Those conditions are used to quickly reject pairs of fragments that
47	could not possibly have an overlap/containment relationship.
48	The dynamic programming algorithm is only applied to pairs of fragments
49	that pass the screening. A large value for DELTA means stringent
50	conditions, where the value for DELTA is a real number at least 8.0.
51
52	POS5 and POS3 are fragment positions such that the 5' end between base 1
53	and base POS5, and the 3' end after base POS3 are of high sequencing
54	error rates, say more than 5%. For mismatches and indels occurring in
55	the two ends, light penalties are used.
56
57	A file of input fragments looks like:
58	>G019uabh
59	ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAA
60	GTCTTGCTTGAATTAAAGACTTGTTTAAACACAAAAATTTAGAGTTTTAC
61	TCAACAAAAGTGATTGATTGATTGATTGATTGATTGATGGTTTACAGTAG
62	GACTTCATTCTAGTCATTATAGCTGCTGGCAGTATAACTGGCCAGCCTTT
63	AATACATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTTGGTATGATTT
64	ATCTTTTTGGTCTTCTATAGCCTCCTTCCCCATCCCCATCAGTCTTAATC
65	AGTCTTGTTACGTTATGACTAATCTTTGGGGATTGTGCAGAATGTTATTT
66	TAGATAAGCAAAACGAGCAAAATGGGGAGTTACTTATATTTCTTTAAAGC
67	>G028uaah
68	CATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTT
69	TAAACACAAAATTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGAT
70	TGATTGATTGATGGTTTACAGTAGGACTTCATTCTAGTCATTATAGCTGC
71	TGGCAGTATAACTGGCCAGCCTTTAATACATTGCTGCTTAGAGTCAAAGC
72	ATGTACTTAGAGTTGGTATGATTTATCTTTTTGGTCTTCTATAGCCTCCT
73	TCCCCATCCCATCAGTCT
74	>G022uabh
75	TATTTTAGAGACCCAAGTTTTTGACCTTTTCCATGTTTACATCAATCCTG
76	TAGGTGATTGGGCAGCCATTTAAGTATTATTATAGACATTTTCACTATCC
77	CATTAAAACCCTTTATGCCCATACATCATAACACTACTTCCTACCCATAA
78	GCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTTTAAAC
79	ACAAAATTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGATTGATT
80	GATTGAT
81	>G023uabh
82	AATAAATACCAAAAAAATAGTATATCTACATAGAATTTCACATAAAATAA
83	ACTGTTTTCTATGTGAAAATTAACCTAAAAATATGCTTTGCTTATGTTTA
84	AGATGTCATGCTTTTTATCAGTTGAGGAGTTCAGCTTAATAATCCTCTAC
85	GATCTTAAACAAATAGGAAAAAAACTAAAAGTAGAAAATGGAAATAAAAT
86	GTCAAAGCATTTCTACCACTCAGAATTGATCTTATAACATGAAATGCTTT
87	TTAAAAGAAAATATTAAAGTTAAACTCCCCTATTTTGCTCGTTTTTGCTT
88	ATCTAAAATACATTCTGCACAATCCCCAAAGATTGATCATACGTTAC
89	>G006uaah
90	ACATAAAATAAACTGTTTTCTATGTGAAAATTAACCTANNATATGCTTTG
91	CTTATGTTTAAGATGTCATGCTTTTTATCAGTTGAGGAGTTCAGCTTAAT
92	AATCCTCTAAGATCTTAAACAAATAGGAAAAAAACTAAAAGTAGAAAATG
93	GAAATAAAATGTCAAAGCATTTCTACCACTCAGAATTGATCTTATAACAT
94	GAAATGCTTTTTAAAAGAAAATATTAAAGTTAAACTCCCC
95	A string after ">" is the name of the following fragment.
96	Only the five upper-case letters A, C, G, T and N are allowed
97	to appear in fragment data. No other characters are allowed.
98	A common mistake is the use of lower case letters in a fragment.
99
100	To run the program, type a command of form
101
102	cap file_of_fragments
103
104	The output goes to the terminal screen. So redirection of the
105	output into a file is necessary. The output consists of three parts:
106	overview of contigs at fragment level, detailed display of contigs
107	at nucleotide level, and consensus sequences.
108	'+' = direct orientation; '-' = reverse complement
109	The output of CAP on the sample input data looks like:
110
111	#Contig 1
112
113	#G022uabh+(0)
114	TATTTTAGAGACCCAAGTTTTTGACCTTTTCCATGTTTACATCAATCCTGTAGGTGATTG
115	GGCAGCCATTTAAGTATTATTATAGACATTTTCACTATCCCATTAAAACCCTTTATGCCC
116	ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTG
117	AATTAAAGACTTGTTTAAACACAAAA-TTTAGACTTTTACTCAACAAAAGTGATTGATTG
118	ATTGATTGATTGATTGAT
119	#G028uaah+(145)
120	CATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTTTAAACACAAA
121	A-TTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGATTGATTGATTGATGGTTTAC
122	AGTAGGACTTCATTCTAGTCATTATAGCTGCTGGCAGTATAACTGGCCAGCCTTTAATAC
123	ATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTTGGTATGATTTATCTTTTTGGTCTTC
124	TATAGCCTCCTTCCCCATCCC-ATCAGTCT
125	#G019uabh+(120)
126	ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTG
127	AATTAAAGACTTGTTTAAACACAAAAATTTAGAGTTTTACTCAACAAAAGTGATTGATTG
128	ATTGATTGATTGATTGATGGTTTACAGTAGGACTTCATTCTAGTCATTATAGCTGCTGGC
129	AGTATAACTGGCCAGCCTTTAATACATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTT
130	GGTATGATTTATCTTTTTGGTCTTCTATAGCCTCCTTCCCCATCCCCATCAGTCTTAATC
131	AGTCTTGTTACGTTATGACT-AATCTTTGGGGATTGTGCAGAATGTTATTTTAGATAAGC
132	AAAA-CGAGCAAAAT-GGGGAGTT-A-CTT-A-TATTT-CTTT-AAA--GC
133	#G023uabh-(426)
134	GTAACGT-ATGA-TCAATCTTTGGGGATTGTGCAGAATGT-ATTTTAGATAAGCAAAAAC
135	GAGCAAAATAGGGGAGTTTAACTTTAATATTTTCTTTTAAAAAGCATTTCATGTTATAAG
136	ATCAATTCTGAGTGGTAGAAATGCTTTGACATTTTATTTCCATTTTCTACTTTTAGTTTT
137	TTTCCTATTTGTTTAAGATCGTAGAGGATTATTAAGCTGAACTCCTCAACTGATAAAAAG
138	CATGACATCTTAAACATAAGCAAAGCATATTTTTAGGTTAATTTTCACATAGAAAACAGT
139	TTATTTTATGTGAAATTCTATGTAGATATACTATTTTTTTGGTATTTATT
140	#G006uaah-(496)
141	GGGGAGTTTAACTTTAATATTTTCTTTTAAAAAGCATTTCATGTTATAAGATCAATTCTG
142	AGTGGTAGAAATGCTTTGACATTTTATTTCCATTTTCTACTTTTAGTTTTTTTCCTATTT
143	GTTTAAGATCTTAGAGGATTATTAAGCTGAACTCCTCAACTGATAAAAAGCATGACATCT
144	TAAACATAAGCAAAGCATATNNT-AGGTTAATTTTCACATAGAAAACAGTTTATTTTATG
145	T
146
147
148
149	Slight modifications by S. Smith on Mon Feb 17 10:18:34 EST 1992.
150	These changes allow for command line arguments for several
151	of the hard coded parameters, as well as a slight modification to
152	the output routine to support GDE format. Changes are commented
153	as:
154
155	Mod by S.S.
156
157	*/
158
159	#include <stdio.h>
160
161	int OVERLEN; /* Minimum length of any overlap */
162	float PERCENT; /* Minimum identity percentage of any overlap */
163	#define CUTOFF 50 /* cutoff score for overlap or containment */
164	#define DELTA 9.0 /* Jump increment in check for overlap */
165	#define MATCH 2 /* score of a match */
166	#define MISMAT -6 /* score of a mismatch */
167	#define LTMISM -3 /* light score of a mismatch */
168	#define OPEN 0 /* gap open penalty */
169	#define EXTEND 4 /* gap extension penalty */
170	#define LTEXTEN 2 /* light gap extension penalty */
171	#define POS5 30 /* Sequencing errors often occur before base POS5 */
172	#define POS3 475 /* Sequencing errors often occur after base POS3 */
173	#define LINELEN 60 /* length of one printed line */
174	#define NAMELEN 20 /* length of printed fragment name */
175	#define TUPLELEN 4 /* Maximum length of words for lookup table */
176	#define SEQLEN 2000 /* initial size of an array for an output fragment */
177
178	static int over_len;
179	static float per_cent;
180	typedef struct OVERLAP /* of 5' and 3' segments */
181	{
182	int number; /* index of 3' segment */
183	int host; /* index of 5' segment */
184	int ind; /* used in reassembly */
185	int stari; /* start position of 5' suffix */
186	int endi; /* end position of 5' suffix */
187	int starj; /* start position of 3' prefix */
188	int endj; /* end position of 3' prefix */
189	short orienti; /* orientation of 5' segment: 0= rev. */
190	short orientj; /* orientation of 3' segment: 0= rev. */
191	int score; /* score of overlap alignment */
192	int length; /* length of alignment */
193	int match; /* number of matches in alignment */
194	short kind; /* 0 = containment; 1 = overlap */
195	int script; / script for constructing alignment */
196	struct OVERLAP *next;
197	} over, *overptr;
198	struct SEG
199	{
200	char name; / name string */
201	int len; /* length of segment name */
202	char seq; / segment sequence */
203	char rev; / reverse complement */
204	int length; /* length of sequence */
205	short kind; /* 0 = contain; 1 = overlap */
206	int lookup; / lookup table */
207	int pos; / location list */
208	overptr list; /* list of overlapping edges */
209	} segment; / array of segment records */
210	int seg_num; /* The number of segments */
211	overptr edge; / set of overlapping edges */
212	int edge_num; /* The number of overlaps */
213	struct CONS /* 1 = itself; 0 = reverse complement */
214	{
215	short isfive[2]; /* is 5' end free? */
216	short isthree[2]; /* is 3' end free? */
217	short orient[2]; /* orientation of 3' segment */
218	int group; /* contig serial number */
219	int next[2]; /* pointer to 3' adjacent segment */
220	int other; /* the other end of the contig */
221	int child; /* for the containment tree */
222	int brother;
223	int father;
224	overptr node[2]; /* pointers to overlapping edges */
225	} contigs; / list of contigs */
226	struct TTREE /* multiple alignment tree */
227	{
228	short head; /* head = 1 means the head of a contig */
229	short orient; /* orientation */
230	int begin; /* start position of previous segment */
231	int script; / alignment script */
232	int size; /* length of script */
233	int next; /* list of overlap segments */
234	int child; /* list of child segments */
235	int brother; /* list of sibling segments */
236	} *mtree;
237	int vert[128]; /* converted digits for 'A','C','G','T' */
238	int vertc[128]; /* for reverse complement */
239	int tuple; /* tuple = TUPLELEN - 1 */
240	int base[TUPLELEN]; /* locations of a lookup table */
241	int power[TUPLELEN]; /* powers of 4 */
242	typedef struct OUT
243	{
244	char line; / one print line */
245	char a; / pointer to slot in line */
246	char c; /* current char */
247	char seq; / pointer to sequence */
248	int length; /* length of segment */
249	int id; /* index of segment */
250	int s; / pointer to script vector */
251	int size; /* size of script vector */
252	int op; /* current operation */
253	char name[NAMELEN+2];/* name of segment */
254	short done; /* indicates if segment is done */
255	int loc; /* position of next segment symbol */
256	char kind; /* type of next symbol of segment */
257	char up; /* type of upper symbol of operation */
258	char dw; /* type of lower symbol of operation */
259	int offset; /* relative to consensus sequence */
260	int linesize; /* size of array line */
261	struct OUT child; / pointer to child subtree */
262	struct OUT brother; / pointer to brother node */
263	struct OUT link; / for operation linked list */
264	struct OUT father; / pointer to father node */
265	} row, rowptr; / node for segment */
266	rowptr work; / a set of working segments */
267	rowptr head, tail; /* first and last nodes of op list */
268	struct VX
269	{
270	int id; /* Segment index */
271	short kind; /* overlap or containment */
272	overptr list; /* list of overlapping edges */
273	} piece; / array of segment records */
274	char allconsen, allconpt; /* Storing consensus sequences */
275
276	main(argc, argv) int argc;
277	char *argv[];
278	{
279	int M; /* Sequence length */
280	int V[128][128], Q,R; /* Weights */
281	int V1[128][128], R1; /* Light weights */
282	int total; /* Total of segment lengths */
283	int number; /* The number of segments */
284	char sequence; / Storing all segments */
285	char reverse; / Storing all reverse segments */
286	int symbol, prev; /* The next character */
287	FILE Ap, ckopen(); /* For the sequence file */
288	char my_calloc(int); / space-allocating function */
289	register int i, j, k; /* index variables */
290	/* Mod by S.S. */
291	int jj;
292	short heading; /* 1: heading; 0: sequence */
293
294	/*
295	* Mod by S.S.
296	*
297	if ( argc != 2 )
298	fatalf("The proper form of command is: \n%s file_of_fragments",
299	argv[0]);
300	*/
301	if ( argc != 4 )
302	fatalf("usage: %s file_of_fragments MIN_OVERLAP PERCENT_MATCH",
303	argv[0]);
304	sscanf(argv[2],"%d",&OVERLEN);
305	sscanf(argv[3],"%d",&jj);
306	PERCENT = (float)jj/100.0;
307	if(PERCENT < 0.25) PERCENT = 0.25;
308	if(PERCENT > 1.0) PERCENT = 1.0;
309	if(OVERLEN < 1) OVERLEN = 1;
310	if(OVERLEN > 100) OVERLEN = 100;
311
312
313
314	/* determine number of segments and total lengths */
315	j = 0;
316
317	Ap = ckopen(argv[1], "r");
318	prev = '\n';
319	for (total = 3, number = 0; ( symbol = getc(Ap)) != EOF ; total++ )
320	{
321	if ( symbol == '>' && prev == '\n' )
322	number++;
323	prev = symbol;
324	}
325	(void) fclose(Ap);
326
327	total += number * 20;
328	/* allocate space for segments */
329	sequence = ( char * ) my_calloc( total * sizeof(char));
330	reverse = ( char * ) my_calloc( total * sizeof(char));
331	allconpt = allconsen = ( char * ) my_calloc( total * sizeof(char));
332	segment = ( struct SEG * ) my_calloc( number * sizeof(struct SEG));
333
334	/* read the segments into sequence */
335	M = 0;
336	Ap = ckopen(argv[1], "r");
337	number = -1;
338	heading = 0;
339	prev = '\n';
340	for ( i = 0, k = total ; ( symbol = getc(Ap)) != EOF ; )
341	{
342	if ( symbol != '\n' )
343	{
344	sequence[++i] = symbol;
345	switch ( symbol )
346	{
347	case 'A' :
348	reverse[--k] = 'T';
349	break;
350	case 'a' :
351	reverse[--k] = 't';
352	break;
353	case 'T' :
354	reverse[--k] = 'A';
355	break;
356	case 't' :
357	reverse[--k] = 'a';
358	break;
359	case 'C' :
360	reverse[--k] = 'G';
361	break;
362	case 'c' :
363	reverse[--k] = 'g';
364	break;
365	case 'G' :
366	reverse[--k] = 'C';
367	break;
368	case 'g' :
369	reverse[--k] = 'c';
370	break;
371	default :
372	reverse[--k] = symbol;
373	break;
374	}
375	}
376	if ( symbol == '>' && prev == '\n' )
377	{
378	heading = 1;
379	if ( number >= 0 )
380	{
381	segment[number].length = i - j - 1;
382	segment[number].rev = &(reverse[k]);
383	if ( i - j - 1 > M ) M = i - j -1;
384	}
385	number++;
386	j = i;
387	segment[number].name = &(sequence[i]);
388	segment[number].kind = 1;
389	segment[number].list = 0;
390	}
391	if ( heading && symbol == '\n' )
392	{
393	heading = 0;
394	segment[number].len = i - j;
395	segment[number].seq = &(sequence[i]);
396	j = i;
397	}
398
399	prev = symbol;
400	}
401	segment[number].length = i - j;
402	reverse[--k] = '>';
403	segment[number].rev = &(reverse[k]);
404	if ( i - j > M ) M = i - j;
405	seg_num = ++number;
406	(void) fclose(Ap);
407
408	Q = OPEN;
409	R = EXTEND;
410	R1 = LTEXTEN;
411	/* set match and mismatch weights */
412	for ( i = 0; i < 128 ; i++ )
413	for ( j = 0; j < 128 ; j++ )
414	if ((i\|32) == (j\|32) )
415	V[i][j] = V1[i][j] = MATCH;
416	else
417	{
418	V[i][j] = MISMAT;
419	V1[i][j] = LTMISM;
420	}
421	for ( i = 0; i < 128 ; i++ )
422	V['N'][i] = V[i]['N'] = MISMAT + 1;
423	V1['N']['N'] = MISMAT + 1;
424
425	over_len = OVERLEN;
426	per_cent = PERCENT;
427	edge_num = 0;
428	INIT(M);
429	MAKE();
430	PAIR(V,V1,Q,R,R1);
431	ASSEM();
432	REPAIR();
433	FORM_TREE();
434	/* GRAPH(); */
435	SHOW();
436	}
437
438	static int (v)[128]; / substitution scores */
439	static int q, r; /* gap penalties */
440	static int qr; /* qr = q + r */
441	static int (v1)[128]; / light substitution scores */
442	static int r1; /* light extension penalty */
443	static int qr1; /* qr1 = q + r1 */
444
445	static int SCORE;
446	static int STARI;
447	static int STARJ;
448	static int ENDI;
449	static int ENDJ;
450
451	static int CC, DD; /* saving matrix scores */
452	static int RR, SS; /* saving start-points */
453	static int S; / saving operations for diff */
454
455	/* The following definitions are for function diff() */
456
457	int diff();
458	static int zero = 0; /* int type zero */
459
460	#define gap(k) ((k) <= 0 ? 0 : q+r(k)) / k-symbol indel score */
461
462	static int sapp; / Current script append ptr */
463	static int last; /* Last script op appended */
464
465	static int no_mat; /* number of matches */
466
467	static int no_mis; /* number of mismatches */
468
469	static int al_len; /* length of alignment */
470	/* Append "Delete k" op */
471	#define DEL(k) \
472	{ al_len += k; \
473	if (last < 0) \
474	last = sapp[-1] -= (k); \
475	else \
476	last = *sapp++ = -(k); \
477	}
478	/* Append "Insert k" op */
479	#define INS(k) \
480	{ al_len += k; \
481	if (last < 0) \
482	{ sapp[-1] = (k); *sapp++ = last; } \
483	else \
484	last = *sapp++ = (k); \
485	}
486
487	/* Append "Replace" op */
488	#define REP \
489	{ last = *sapp++ = 0; \
490	al_len += 1; \
491	}
492
493	INIT(M) int M;
494	{
495	register int j; /* row and column indices */
496	char my_calloc(); / space-allocating function */
497
498	/* allocate space for all vectors */
499	j = (M + 1) * sizeof(int );
500	CC = ( int * ) my_calloc(j);
501	DD = ( int * ) my_calloc(j);
502	RR = ( int * ) my_calloc(j);
503	SS = ( int * ) my_calloc(j);
504	S = ( int * ) my_calloc(2 * j);
505	}
506
507	/* Make a lookup table for words of lengths up to TUPLELEN in each sequence.
508	The value of a word is used as an index to the table. */
509	MAKE()
510	{
511	int hash[TUPLELEN]; /* values of words of lengths up to TUPLELEN */
512	int table; / pointer to a lookup table */
513	int loc; / pointer to a table of sequence locations */
514	int size; /* size of a lookup table */
515	int limit, digit, step; /* temporary variables */
516	register int i, j, k, p, q; /* index varibles */
517	char my_calloc(); / space-allocating function */
518	char A; / pointer to a sequence */
519	int M; /* length of a sequence */
520
521	tuple = TUPLELEN - 1;
522	for ( i = 0; i < 128; i++ )
523	vert[i] = 4;
524	vert['A'] = vert['a'] = 0;
525	vert['C'] = vert['c'] = 1;
526	vert['G'] = vert['g'] = 2;
527	vert['T'] = vert['t'] = 3;
528	vertc['A'] = vertc['a'] = 3;
529	vertc['C'] = vertc['c'] = 2;
530	vertc['G'] = vertc['g'] = 1;
531	vertc['T'] = vertc['t'] = 0;
532	for ( i = 0, j = 1, size = 0; i <= tuple ; i++, j *= 4 )
533	{
534	base[i] = size;
535	power[i] = j;
536	size = ( size + 1 ) * 4;
537	}
538	for ( j = 0; j <= tuple; j++ )
539	hash[j] = 0;
540	/* make a lookup table for each sequence */
541	for ( i = 0; i < seg_num ; i++ )
542	{
543	A = segment[i].seq;
544	M = segment[i].length;
545	table = segment[i].lookup = (int * ) my_calloc(size * sizeof(int ));
546	loc = segment[i].pos = (int * ) my_calloc((M + 1) * sizeof(int ));
547	for ( j = 0; j < size; j++ )
548	table[j] = 0;
549	for ( k = 0, j = 1; j <= M; j++ )
550	if ( ( digit = vert[A[j]] ) != 4 )
551	{
552	for ( p = tuple; p > 0; p-- )
553	hash[p] = 4 * (hash[p-1] + 1) + digit;
554	hash[0] = digit;
555	step = j - k;
556	limit = tuple <= step ? tuple : step;
557	for ( p = 0; p < limit; p++ )
558	if ( ! table[(q = hash[p])] ) table[q] = 1;
559	if ( step > tuple )
560	{
561	loc[(p = j - tuple)] = table[(q = hash[tuple])];
562	table[q] = p;
563	}
564	}
565	else
566	k = j;
567	}
568	}
569
570	/* Perform pair-wise comparisons of sequences. The sequences not
571	satisfying the necessary condition for overlap are rejected quickly.
572	Those that do satisfy the condition are verified with a dynamic
573	programming algorithm to see if overlaps exist. */
574	PAIR(V,V1,Q,R,R1) int V[][128],V1[][128],Q,R,R1;
575	{
576	int endi, endj, stari, starj; /* endpoint and startpoint */
577
578	short orienti, orientj; /* orientation of segments */
579	short iscon; /* containment condition */
580	int score; /* the max score */
581	int count, limit; /* temporary variables */
582
583	register int i, j, d; /* row and column indices */
584	char my_calloc(); / space-allocating function */
585	int rl, cl;
586	char A, B;
587	int M, N;
588	overptr node1;
589	int total; /* total number of pairs */
590	int hit; /* number of pairs satisfying cond. */
591	int CHECK();
592
593	v = V;
594	q = Q;
595	r = R;
596	qr = q + r;
597	v1 = V1;
598	r1 = R1;
599	qr1 = q + r1;
600	total = hit = 0;
601	limit = 2 * ( seg_num - 1 );
602	for ( orienti = 0, d = 0; d < limit ; d++ )
603	{
604	i = d / 2;
605	orienti = 1 - orienti;
606	A = orienti ? segment[i].seq : segment[i].rev;
607	M = segment[i].length;
608	for ( j = i+1; j < seg_num ; j++ )
609	{
610	B = segment[j].seq;
611	orientj = 1;
612	N = segment[j].length;
613	total += 1;
614	if ( CHECK(orienti, i, j) )
615	{
616	hit += 1;
617	SCORE = 0;
618	big_pass(A,B,M,N,orienti,orientj);
619	if ( SCORE )
620	{
621	score = SCORE;
622	stari = ++STARI;
623	starj = ++STARJ;
624	endi = ENDI;
625	endj = ENDJ;
626	rl = endi - stari + 1;
627	cl = endj - starj + 1;
628	sapp = S;
629	last = 0;
630	al_len = no_mat = no_mis = 0;
631	(void) diff(&A[stari]-1, &B[starj]-1,rl,cl,q,q);
632	iscon = stari == 1 && endi == M \|\| starj == 1 && endj == N;
633	if ( no_mat >= al_len * per_cent &&
634	( al_len >= over_len \|\| iscon ) )
635	{
636	node1 = ( overptr ) my_calloc( (int ) sizeof(over));
637	if ( iscon )
638	node1->kind = 0; /* containment */
639	else
640	{
641	node1->kind = 1;
642	edge_num++;
643	} /* overlap */
644	if ( endi == M && ( endj != N \|\| starj != 1 ) ) /i is 5'/
645	{
646	node1->number = j;
647	node1->host = i;
648	node1->stari = stari;
649	node1->endi = endi;
650	node1->orienti = orienti;
651	node1->starj = starj;
652	node1->endj = endj;
653	node1->orientj = orientj;
654	}
655	else /* j is 5' */
656	{
657	node1->number = i;
658	node1->host = j;
659	node1->stari = starj;
660	node1->endi = endj;
661	node1->orienti = orientj;
662	node1->starj = stari;
663	node1->endj = endi;
664	node1->orientj = orienti;
665	}
666	node1->score = score;
667	node1->length = al_len;
668	node1->match = no_mat;
669	count = node1->number == i ? j : i;
670	node1->next = segment[count].list;
671	segment[count].list = node1;
672	if ( ! node1->kind )
673	segment[count].kind = 0;
674	}
675	}
676	}
677	}
678	}
679	}
680
681	/* Return 1 if two sequences satisfy the necessary condition for overlap,
682	and 0 otherwise. Parameters first and second are the indices of segments,
683	and parameter orient indicates the orientation of segment first. */
684	int CHECK(orient,first,second) short orient;
685	int first, second;
686	{
687	int limit, bound; /* maximum number of jumps */
688	int small; /* smaller of limit and bound */
689	float delta; /* cutoff factor */
690	float cut; /* cutoff score */
691	register int i; /* index variable */
692	int t, q; /* temporary variables */
693	int JUMP();
694	int RJUMP();
695	int JUMPC();
696	int RJUMPC();
697
698	delta = DELTA;
699	if ( orient )
700	limit = JUMP(CC, first, second, 1);
701	else
702	limit = JUMPC(CC, first, second);
703	bound = RJUMP(DD, second, first, orient);
704	small = limit <= bound ? limit : bound;
705	cut = 0;
706	for ( i = 1; i <= small; i++ )
707	{
708	if ( (t = DD[i] - 1) >= over_len && t >= cut
709	&& (q = CC[i] - 1) >= over_len && q >= cut )
710	return (1);
711	cut += delta;
712	}
713	if (DD[bound] >= delta(bound-1) \|\| CC[limit] >= delta(limit-1))
714	return (1);
715	limit = JUMP(CC, second, first, orient);
716	if ( orient )
717	bound = RJUMP(DD, first, second, 1);
718	else
719	bound = RJUMPC(DD, first, second);
720	small = limit <= bound ? limit : bound;
721	cut = 0;
722	for ( i = 1; i <= small; i++ )
723	{
724	if ( (t = DD[i] - 1) >= over_len && t >= cut
725	&& (q = CC[i] - 1) >= over_len && q >= cut )
726	return (1);
727	cut += delta;
728	}
729	return (0);
730	}
731
732	/* Compute a vector of lengths of jumps */
733	int JUMP(H,one,two,orient) int H[], one, two;
734	short orient;
735	{
736	char A, B; /* pointers to two sequences */
737	int M, N; /* lengths of two sequences */
738	int table; / pointer to a lookup table */
739	int loc; / pointer to a location table */
740	int value; /* value of a word */
741	int maxd; /* maximum length of an identical diagonal */
742	int d; /* length of current identical diagonal */
743	int s; /* length of jumps */
744	int k; /* number of jumps */
745
746	register int i, j, p; /* index variables */
747
748	A = segment[one].seq;
749	M = segment[one].length;
750	table = segment[one].lookup;
751	loc = segment[one].pos;
752	B = orient ? segment[two].seq : segment[two].rev;
753	N = segment[two].length;
754	for ( s = 1, k = 1; s <= N ; k++ )
755	{
756	maxd = 0;
757	for ( value = -1, d = 0, j = s; d <= tuple && j <= N; j++, d++)
758	{
759	if ( vert[B[j]] == 4 ) break;
760	value = 4 * (value + 1) + vert[B[j]];
761	if ( ! table[value] ) break;
762	}
763	if ( d > tuple )
764	{
765	for ( p = table[value]; p ; p = loc[p] )
766	{
767	d = tuple + 1;
768	for ( i = p+d, j = s+d; i <= M && j <= N; i++, j++, d++ )
769	if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
770	if ( maxd < d )
771	maxd = d;
772	if ( j > N ) break;
773	}
774	}
775	else
776	maxd = d;
777	s += maxd + 1;
778	H[k] = s;
779	}
780	return (k - 1);
781	}
782
783	/* Compute a vector of lengths of jumps for reverse complement of one */
784	int JUMPC(H,one,two) int H[], one, two;
785	{
786	char A, B; /* pointers to two sequences */
787	int M, N; /* lengths of two sequences */
788	int table; / pointer to a lookup table */
789	int loc; / pointer to a location table */
790	int value; /* value of a word */
791	int maxd; /* maximum length of an identical diagonal */
792	int d; /* length of current identical diagonal */
793	int s; /* length of jumps */
794	int k; /* number of jumps */
795
796	register int i, j, p; /* index variables */
797
798	A = segment[one].rev;
799	M = segment[one].length;
800	table = segment[one].lookup;
801	loc = segment[one].pos;
802	B = segment[two].seq;
803	N = segment[two].length;
804	for ( s = 1, k = 1; s <= N ; k++ )
805	{
806	maxd = 0;
807	for ( value = 0, d = 0, j = s; d <= tuple && j <= N; j++, d++)
808	{
809	if ( vert[B[j]] == 4 ) break;
810	value += vertc[B[j]] * power[d];
811	if ( ! table[value+base[d]] ) break;
812	}
813	if ( d > tuple )
814	{
815	for ( p = table[value+base[tuple]]; p ; p = loc[p] )
816	{
817	d = tuple + 1;
818	for ( i = M+2-p, j = s+d; i <= M && j <= N; i++, j++, d++ )
819	if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
820	if ( maxd < d )
821	maxd = d;
822	if ( j > N ) break;
823	}
824	}
825	else
826	maxd = d;
827	s += maxd + 1;
828	H[k] = s;
829	}
830	return (k - 1);
831	}
832
833	/* Compute a vector of lengths of reverse jumps */
834	int RJUMP(H,one,two,orient) int H[], one, two;
835	short orient;
836	{
837	char A, B; /* pointers to two sequences */
838	int N; /* length of a sequence */
839	int table; / pointer to a lookup table */
840	int loc; / pointer to a location table */
841	int value; /* value of a word */
842	int maxd; /* maximum length of an identical diagonal */
843	int d; /* length of current identical diagonal */
844	int s; /* length of jumps */
845	int k; /* number of jumps */
846
847	register int i, j, p; /* index variables */
848
849	A = segment[one].seq;
850	table = segment[one].lookup;
851	loc = segment[one].pos;
852	B = orient ? segment[two].seq : segment[two].rev;
853	N = segment[two].length;
854	for ( s = 1, k = 1; s <= N ; k++ )
855	{
856	maxd = 0;
857	for (value = 0, d = 0, j = N+1-s; d <= tuple && j >= 1; j--, d++)
858	{
859	if ( vert[B[j]] == 4 ) break;
860	value += vert[B[j]] * power[d];
861	if ( ! table[value+base[d]] ) break;
862	}
863	if ( d > tuple )
864	{
865	for ( p = table[value+base[tuple]]; p ; p = loc[p] )
866	{
867	d = tuple + 1;
868	for ( i = p-1, j = N+1-s-d; i >= 1 && j >= 1; i--, j--, d++ )
869	if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
870	if ( maxd < d )
871	maxd = d;
872	if ( j < 1 ) break;
873	}
874	}
875	else
876	maxd = d;
877	s += maxd + 1;
878	H[k] = s;
879	}
880	return (k - 1);
881	}
882
883	/* Compute a vector of lengths of reverse jumps for reverse complement */
884	int RJUMPC(H,one,two) int H[], one, two;
885	{
886	char A, B; /* pointers to two sequences */
887	int M, N; /* lengths of two sequences */
888	int table; / pointer to a lookup table */
889	int loc; / pointer to a location table */
890	int value; /* value of a word */
891	int maxd; /* maximum length of an identical diagonal */
892	int d; /* length of current identical diagonal */
893	int s; /* length of jumps */
894	int k; /* number of jumps */
895
896	register int i, j, p; /* index variables */
897
898	A = segment[one].rev;
899	M = segment[one].length;
900	table = segment[one].lookup;
901	loc = segment[one].pos;
902	B = segment[two].seq;
903	N = segment[two].length;
904	for ( s = 1, k = 1; s <= N ; k++ )
905	{
906	maxd = 0;
907	for (value = -1, d = 0, j = N+1-s; d <= tuple && j >= 1; j--, d++)
908	{
909	if ( vert[B[j]] == 4 ) break;
910	value = 4 * (value + 1) + vertc[B[j]];
911	if ( ! table[value] ) break;
912	}
913	if ( d > tuple )
914	{
915	for ( p = table[value]; p ; p = loc[p] )
916	{
917	d = tuple + 1;
918	i = M - p - tuple;
919	for ( j = N-s-tuple; i >= 1 && j >= 1; i--, j--, d++ )
920	if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
921	if ( maxd < d )
922	maxd = d;
923	if ( j < 1 ) break;
924	}
925	}
926	else
927	maxd = d;
928	s += maxd + 1;
929	H[k] = s;
930	}
931	return (k - 1);
932	}
933
934	/* Construct contigs */
935	ASSEM()
936	{
937	char my_calloc(); / space-allocating function */
938	register int i, j, k; /* index variables */
939	overptr node1, x, y; /* temporary pointer */
940	int five, three; /* indices of 5' and 3' segments */
941	short orienti; /* orientation of 5' segment */
942	short orientj; /* orientation of 3' segment */
943	short sorted; /* boolean variable */
944
945	contigs = ( struct CONS * ) my_calloc( seg_num * sizeof(struct CONS));
946	for ( i = 0; i < seg_num; i++ )
947	{
948	contigs[i].isfive[0] = contigs[i].isthree[0] = 1;
949	contigs[i].isfive[1] = contigs[i].isthree[1] = 1;
950	contigs[i].other = i;
951	contigs[i].group = contigs[i].child = -1;
952	contigs[i].brother = contigs[i].father = -1;
953	}
954	for ( i = 0; i < seg_num; i++ )
955	if ( ! segment[i].kind )
956	for ( ; ; )
957	{
958	for ( y = segment[i].list; y->kind; y = y->next )
959	;
960	for ( x = y->next; x != 0; x = x->next )
961	if ( ! x->kind && x->score > y->score )
962	y = x;
963	for ( j = y->number; (k = contigs[j].father) != -1; j = k )
964	;
965	if ( j != i )
966	{
967	contigs[i].father = j = y->number;
968	contigs[i].brother = contigs[j].child;
969	contigs[j].child = i;
970	contigs[i].node[1] = y;
971	break;
972	}
973	else
974	{
975	if ( segment[i].list->number == y->number )
976	segment[i].list = y->next;
977	else
978	{
979	for ( x = segment[i].list; x->next->number != y->number ; )
980	x = x->next;
981	x->next = y->next;
982	}
983	for ( x = segment[i].list; x != 0 && x->kind; x = x->next )
984	;
985	if ( x == 0 )
986	{
987	segment[i].kind = 1;
988	break;
989	}
990	}
991	}
992	edge = ( overptr * ) my_calloc( edge_num * sizeof(overptr) );
993	for ( j = 0, i = 0; i < seg_num; i++ )
994	if ( segment[i].kind )
995	for ( node1 = segment[i].list; node1 != 0; node1 = node1->next )
996	if ( segment[node1->number].kind )
997	edge[j++] = node1;
998	edge_num = j;
999	for ( i = edge_num - 1; i > 0; i-- )
1000	{
1001	sorted = 1;
1002	for ( j = 0; j < i; j++ )
1003	if ( edge[j]->score < edge[j+1]->score )
1004	{
1005	node1 = edge[j];
1006	edge[j] = edge[j+1];
1007	edge[j+1] = node1;
1008	sorted = 0;
1009	}
1010	if ( sorted )
1011	break;
1012	}
1013	for ( k = 0; k < edge_num; k++ )
1014	{
1015	five = edge[k]->host;
1016	three = edge[k]->number;
1017	orienti = edge[k]->orienti;
1018	orientj = edge[k]->orientj;
1019	if ( contigs[five].isthree[orienti] &&
1020	contigs[three].isfive[orientj] && contigs[five].other != three )
1021	{
1022	contigs[five].isthree[orienti] = 0;
1023	contigs[three].isfive[orientj] = 0;
1024	contigs[five].next[orienti] = three;
1025	contigs[five].orient[orienti] = orientj;
1026	contigs[five].node[orienti] = edge[k];
1027	contigs[three].isthree[(j = 1 - orientj)] = 0;
1028	contigs[five].isfive[(i = 1 - orienti)] = 0;
1029	contigs[three].next[j] = five;
1030	contigs[three].orient[j] = i;
1031	contigs[three].node[j] = edge[k];
1032	i = contigs[three].other;
1033	j = contigs[five].other;
1034	contigs[i].other = j;
1035	contigs[j].other = i;
1036	}
1037	}
1038	}
1039
1040	REPAIR()
1041	{
1042	int endi, endj, stari, starj; /* endpoint and startpoint */
1043
1044	short orienti, orientj; /* orientation of segments */
1045	short isconi, isconj; /* containment condition */
1046	int score; /* the max score */
1047	int i, j, f, d, e; /* row and column indices */
1048	char my_calloc(); / space-allocating function */
1049	char A, B;
1050	int M, N;
1051	overptr node1;
1052	int piece_num; /* The number of pieces */
1053	int count, limit;
1054	int number;
1055	int hit;
1056
1057	piece = ( struct VX * ) my_calloc( seg_num * sizeof(struct VX));
1058	for ( j = 0, i = 0; i < seg_num; i++ )
1059	if (segment[i].kind && (contigs[i].isfive[1] \|\| contigs[i].isfive[0]))
1060	piece[j++].id = i;
1061	piece_num = j;
1062	for ( i = 0; i < piece_num; i++ )
1063	{
1064	piece[i].kind = 1;
1065	piece[i].list = 0;
1066	}
1067	limit = 2 * ( piece_num - 1 );
1068	hit = number = 0;
1069	for ( orienti = 0, d = 0; d < limit ; d++ )
1070	{
1071	i = piece[(e = d / 2)].id;
1072	orienti = 1 - orienti;
1073	A = orienti ? segment[i].seq : segment[i].rev;
1074	M = segment[i].length;
1075	for ( f = e+1; f < piece_num ; f++ )
1076	{
1077	j = piece[f].id;
1078	B = segment[j].seq;
1079	orientj = 1;
1080	N = segment[j].length;
1081	SCORE = 0;
1082	hit++;
1083	big_pass(A,B,M,N,orienti,orientj);
1084	if ( SCORE > CUTOFF )
1085	{
1086	score = SCORE;
1087	stari = ++STARI;
1088	starj = ++STARJ;
1089	endi = ENDI;
1090	endj = ENDJ;
1091	isconi = stari == 1 && endi == M;
1092	isconj = starj == 1 && endj == N;
1093	node1 = ( overptr ) my_calloc( (int ) sizeof(over));
1094	if ( isconi \|\| isconj )
1095	node1->kind = 0; /* containment */
1096	else
1097	{
1098	node1->kind = 1;
1099	number++;
1100	} /* overlap */
1101	if ( endi == M && ! isconj ) /i is 5'/
1102	{
1103	node1->number = j;
1104	node1->host = i;
1105	node1->ind = f;
1106	node1->stari = stari;
1107	node1->endi = endi;
1108	node1->orienti = orienti;
1109	node1->starj = starj;
1110	node1->endj = endj;
1111	node1->orientj = orientj;
1112	}
1113	else /* j is 5' */
1114	{
1115	node1->number = i;
1116	node1->host = j;
1117	node1->ind = e;
1118	node1->stari = starj;
1119	node1->endi = endj;
1120	node1->orienti = orientj;
1121	node1->starj = stari;
1122	node1->endj = endi;
1123	node1->orientj = orienti;
1124	}
1125	node1->score = score;
1126	count = node1->number == i ? f : e;
1127	node1->next = piece[count].list;
1128	piece[count].list = node1;
1129	if ( ! node1->kind )
1130	piece[count].kind = 0;
1131	}
1132	}
1133	}
1134	REASSEM(piece_num, number);
1135	}
1136
1137	/* Construct contigs */
1138	REASSEM(piece_num, number) int piece_num, number;
1139	{
1140	char my_calloc(); / space-allocating function */
1141	int i, j, k, d; /* index variables */
1142	overptr node1, x, y; /* temporary pointer */
1143	int five, three; /* indices of 5' and 3' segments */
1144	short orienti; /* orientation of 5' segment */
1145	short orientj; /* orientation of 3' segment */
1146	short sorted; /* boolean variable */
1147
1148	for ( d = 0; d < piece_num; d++ )
1149	if ( ! piece[d].kind )
1150	for ( i = piece[d].id ; ; )
1151	{
1152	for ( y = piece[d].list; y->kind; y = y->next )
1153	;
1154	for ( x = y->next; x != 0; x = x->next )
1155	if ( ! x->kind && x->score > y->score )
1156	y = x;
1157	for ( j = y->number; (k = contigs[j].father) != -1; j = k )
1158	;
1159	if ( j != i && RECONCILE(y,&piece_num,&number) )
1160	{
1161	contigs[i].father = j = y->number;
1162	contigs[i].brother = contigs[j].child;
1163	contigs[j].child = i;
1164	contigs[i].node[1] = y;
1165	segment[i].kind = 0;
1166	break;
1167	}
1168	else
1169	{
1170	if ( piece[d].list->number == y->number )
1171	piece[d].list = y->next;
1172	else
1173	{
1174	for ( x = piece[d].list; x->next->number != y->number ; )
1175	x = x->next;
1176	x->next = y->next;
1177	}
1178	for ( x = piece[d].list; x != 0 && x->kind; x = x->next )
1179	;
1180	if ( x == 0 )
1181	{
1182	piece[d].kind = 1;
1183	break;
1184	}
1185	}
1186	}
1187	if ( number > edge_num )
1188	edge = ( overptr * ) my_calloc( number * sizeof(overptr) );
1189	for ( j = 0, d = 0; d < piece_num; d++ )
1190	if ( piece[d].kind )
1191	for ( node1 = piece[d].list; node1 != 0; node1 = node1->next )
1192	if ( piece[node1->ind].kind )
1193	edge[j++] = node1;
1194	edge_num = j;
1195	for ( i = edge_num - 1; i > 0; i-- )
1196	{
1197	sorted = 1;
1198	for ( j = 0; j < i; j++ )
1199	if ( edge[j]->score < edge[j+1]->score )
1200	{
1201	node1 = edge[j];
1202	edge[j] = edge[j+1];
1203	edge[j+1] = node1;
1204	sorted = 0;
1205	}
1206	if ( sorted )
1207	break;
1208	}
1209	for ( k = 0; k < edge_num; k++ )
1210	{
1211	five = edge[k]->host;
1212	three = edge[k]->number;
1213	orienti = edge[k]->orienti;
1214	orientj = edge[k]->orientj;
1215	if ( contigs[five].isthree[orienti] &&
1216	contigs[three].isfive[orientj] && contigs[five].other != three )
1217	{
1218	contigs[five].isthree[orienti] = 0;
1219	contigs[three].isfive[orientj] = 0;
1220	contigs[five].next[orienti] = three;
1221	contigs[five].orient[orienti] = orientj;
1222	contigs[five].node[orienti] = edge[k];
1223	contigs[three].isthree[(j = 1 - orientj)] = 0;
1224	contigs[five].isfive[(i = 1 - orienti)] = 0;
1225	contigs[three].next[j] = five;
1226	contigs[three].orient[j] = i;
1227	contigs[three].node[j] = edge[k];
1228	i = contigs[three].other;
1229	j = contigs[five].other;
1230	contigs[i].other = j;
1231	contigs[j].other = i;
1232	}
1233	}
1234	}
1235
1236	RECONCILE(y, pp,nn) overptr y;
1237	int pp,nn;
1238	{
1239	short orienti, orientj; /* orientation of segments */
1240	short orientk, orientd; /* orientation of segments */
1241	int i, j, k, d, f; /* row and column indices */
1242	char my_calloc(); / space-allocating function */
1243	char A, B;
1244	int M, N;
1245	overptr node1;
1246
1247	k = y->host;
1248	d = y->number;
1249	orientk = y->orienti;
1250	orientd = y->orientj;
1251	if ( ! contigs[k].isthree[orientk] )
1252	{
1253	if ( ! piece[y->ind].kind ) return (0);
1254	if ( contigs[d].isthree[orientd] )
1255	{
1256	orienti = orientd;
1257	i = d;
1258	orientj = contigs[k].orient[orientk];
1259	j = contigs[k].next[orientk];
1260	}
1261	else
1262	return (0);
1263	}
1264	else
1265	if ( ! contigs[k].isfive[orientk] )
1266	{
1267	if ( ! piece[y->ind].kind ) return (0);
1268	if ( contigs[d].isfive[orientd] )
1269	{
1270	orienti = contigs[k].orient[1-orientk];
1271	orienti = 1 - orienti;
1272	i = contigs[k].next[1-orientk];
1273	orientj = orientd;
1274	j = d;
1275	}
1276	else
1277	return (0);
1278	}
1279	else
1280	return (0);
1281	A = orienti ? segment[i].seq : segment[i].rev;
1282	M = segment[i].length;
1283	B = orientj ? segment[j].seq : segment[j].rev;
1284	N = segment[j].length;
1285	SCORE = 0;
1286	big_pass(A,B,M,N,orienti,orientj);
1287	if ( SCORE > CUTOFF && ENDI - STARI > over_len && ENDI == M && STARJ == 0 )
1288	{
1289	node1 = ( overptr ) my_calloc( (int ) sizeof(over));
1290	node1->kind = 1;
1291	node1->host = i;
1292	node1->number = j;
1293	node1->stari = ++STARI;
1294	node1->endi = ENDI;
1295	node1->orienti = orienti;
1296	node1->starj = ++STARJ;
1297	node1->endj = ENDJ;
1298	node1->orientj = orientj;
1299	node1->score = SCORE;
1300	piece[*pp].kind = 1;
1301	if ( i == d )
1302	{
1303	node1->ind = *pp;
1304	node1->next = piece[y->ind].list;
1305	piece[y->ind].list = node1;
1306	piece[*pp].id = j;
1307	piece[*pp].list = 0;
1308	}
1309	else
1310	{
1311	node1->ind = y->ind;
1312	piece[*pp].list = node1;
1313	node1->next = 0;
1314	piece[*pp].id = i;
1315	}
1316	(*nn)++;
1317	(*pp)++;
1318	f = contigs[k].other;
1319	if ( ! contigs[k].isthree[orientk] )
1320	{
1321	contigs[j].isfive[orientj] = 1;
1322	contigs[j].isthree[1 - orientj] = 1;
1323	contigs[k].isthree[orientk] = 1;
1324	contigs[k].isfive[1 - orientk] = 1;
1325	contigs[f].other = j;
1326	contigs[j].other = f;
1327	}
1328	else
1329	{
1330	contigs[i].isthree[orienti] = 1;
1331	contigs[i].isfive[1 - orienti] = 1;
1332	contigs[k].isfive[orientk] = 1;
1333	contigs[k].isthree[1 - orientk] = 1;
1334	contigs[f].other = i;
1335	contigs[i].other = f;
1336	}
1337	contigs[k].other = k;
1338	return (1);
1339	}
1340	return (0);
1341	}
1342
1343	/* Construct a tree of overlapping-containment segments */
1344	FORM_TREE()
1345	{
1346	register int i, j, k; /* index variables */
1347	char my_calloc(); / space-allocating function */
1348	overptr node1; /* temporary pointer */
1349	short orient; /* orientation of segment */
1350	int group; /* serial number of contigs */
1351	char A, B; /* pointers to segment sequences */
1352	int stari, endi, starj, endj;/* positions where alignment begins */
1353	int M, N; /* lengths of segment sequences */
1354	int count; /* temporary variables */
1355
1356	mtree = ( struct TTREE * ) my_calloc( seg_num * sizeof(struct TTREE));
1357	for ( i = 0; i < seg_num; i++ )
1358	{
1359	mtree[i].head = 0;
1360	mtree[i].next = mtree[i].child = mtree[i].brother = -1;
1361	}
1362	for ( group = 0, i = 0; i < seg_num; i++ )
1363	if ( segment[i].kind && contigs[i].group < 0 &&
1364	( contigs[i].isfive[1] \|\| contigs[i].isfive[0] ) )
1365	{
1366	orient = contigs[i].isfive[1] ? 1 : 0;
1367	mtree[i].head = 1;
1368	for ( j = i; ; )
1369	{
1370	contigs[j].group = group;
1371	mtree[j].orient = orient;
1372	SORT(j, orient);
1373	if ( contigs[j].isthree[orient] )
1374	break;
1375	else
1376	{
1377	k = contigs[j].next[orient];
1378	node1 = contigs[j].node[orient];
1379	if ( j == node1->host )
1380	{
1381	stari = node1->stari;
1382	endi = node1->endi;
1383	starj = node1->starj;
1384	endj = node1->endj;
1385	A = node1->orienti ? segment[j].seq : segment[j].rev;
1386	B = node1->orientj ? segment[k].seq : segment[k].rev;
1387	}
1388	else
1389	{
1390	M = segment[j].length;
1391	stari = M + 1 - node1->endj;
1392	endi = M + 1 - node1->starj;
1393	N = segment[k].length;
1394	starj = N + 1 - node1->endi;
1395	endj = N + 1 - node1->stari;
1396	A = node1->orientj ? segment[j].rev : segment[j].seq;
1397	B = node1->orienti ? segment[k].rev : segment[k].seq;
1398	}
1399	M = endi - stari + 1;
1400	N = endj - starj + 1;
1401	sapp = S;
1402	last = 0;
1403	al_len = no_mat = no_mis = 0;
1404	(void) diff(&A[stari]-1, &B[starj]-1,M,N,q,q);
1405	count = ( (N = sapp - S) + 1 ) * sizeof(int );
1406	mtree[k].script = ( int * ) my_calloc( count );
1407	for ( M = 0; M < N; M++)
1408	mtree[k].script[M] = S[M];
1409	mtree[k].size = N;
1410	mtree[k].begin = stari;
1411	mtree[j].next = k;
1412	orient = contigs[j].orient[orient];
1413	j = k;
1414	}
1415	}
1416	group++;
1417	}
1418	}
1419
1420	/* Sort the children of each node by the `begin' field */
1421	SORT(seg, ort) int seg;
1422	short ort;
1423	{
1424	register int i, j, k; /* index variables */
1425	char my_calloc(); / space-allocating function */
1426	overptr node1; /* temporary pointer */
1427	short orient; /* orientation of segment */
1428	char A, B; /* pointers to segment sequences */
1429	int stari, endi, starj, endj;/* positions where alignment begins */
1430	int M, N; /* lengths of segment sequences */
1431	int count; /* temporary variables */
1432
1433	for ( j = contigs[seg].child; j != -1; j = contigs[j].brother )
1434	{
1435	node1 = contigs[j].node[1];
1436	if ( ort == node1->orientj )
1437	{
1438	stari = node1->starj;
1439	endi = node1->endj;
1440	starj = node1->stari;
1441	endj = node1->endi;
1442	A = node1->orientj ? segment[seg].seq : segment[seg].rev;
1443	B = node1->orienti ? segment[j].seq : segment[j].rev;
1444	orient = node1->orienti;
1445	}
1446	else
1447	{
1448	M = segment[seg].length;
1449	stari = M + 1 - node1->endj;
1450	endi = M + 1 - node1->starj;
1451	N = segment[j].length;
1452	starj = N + 1 - node1->endi;
1453	endj = N + 1 - node1->stari;
1454	A = node1->orientj ? segment[seg].rev : segment[seg].seq;
1455	B = node1->orienti ? segment[j].rev : segment[j].seq;
1456	orient = 1 - node1->orienti;
1457	}
1458	M = endi - stari + 1;
1459	N = endj - starj + 1;
1460	sapp = S;
1461	last = 0;
1462	al_len = no_mat = no_mis = 0;
1463	(void) diff(&A[stari]-1, &B[starj]-1,M,N,q,q);
1464	count = ( (M = sapp - S ) + 1 ) * sizeof(int );
1465	mtree[j].script = ( int * ) my_calloc( count );
1466	for ( k = 0; k < M; k++)
1467	mtree[j].script[k] = S[k];
1468	mtree[j].size = M;
1469	mtree[j].begin = stari;
1470	mtree[j].orient = orient;
1471	if ( mtree[seg].child == -1 )
1472	mtree[seg].child = j;
1473	else
1474	{
1475	i = mtree[seg].child;
1476	if ( mtree[i].begin >= stari )
1477	{
1478	mtree[j].brother = i;
1479	mtree[seg].child = j;
1480	}
1481	else
1482	{
1483	M = mtree[i].brother;
1484	for ( ; M != -1; i = M, M = mtree[M].brother )
1485	if ( mtree[M].begin >= stari ) break;
1486	mtree[j].brother = M;
1487	mtree[i].brother = j;
1488	}
1489	}
1490	SORT(j, orient);
1491	}
1492	}
1493
1494	/* Display the alignments of segments */
1495	SHOW()
1496	{
1497	register int i, j, k; /* index variables */
1498	char my_calloc(); / space-allocating function */
1499	int n; /* number of working segments */
1500	int limit; /* number of slots in work */
1501	int col; /* number of output columns prepared */
1502	short done; /* tells if current group is done */
1503	rowptr root; /* pointer to root of op tree */
1504	int sym[6]; /* occurrence counts for six chars */
1505	char c; /* temp variable */
1506	rowptr t, w, yy; /* temp pointer */
1507	int x; /* temp variables */
1508	int group; /* Contigs number */
1509	char conlit[20], a; / String form of contig number */
1510	char spt; / pointer to the start of consensus */
1511
1512	work = ( rowptr * ) my_calloc( seg_num * sizeof(rowptr));
1513	group = 0;
1514	yy = 0;
1515	for ( j = 0; j < 6; j++ )
1516	sym[j] = 0;
1517	n = limit = col = 0;
1518	for ( i = 0; i < seg_num; i++ )
1519	if ( mtree[i].head )
1520	{
1521	(void) sprintf(conlit, ">Contig %d\n", group);
1522	for ( a = conlit; *a; )
1523	allconpt++ = a++;
1524	/* Mod by S.S.
1525	(void) printf("\n#Contig %d\n\n", group++);
1526	*/
1527	group++;
1528	done = 0;
1529	ENTER(&limit, &n, i, col, yy);
1530	root = work[0];
1531	spt = allconpt;
1532	while ( ! done )
1533	{
1534	for ( j = 0; j < n; j++ ) /* get segments into work */
1535	{
1536	t = work[j];
1537	k = t->id;
1538	if ((x = mtree[k].next) != -1 && mtree[x].begin == t->loc)
1539	{
1540	ENTER(&limit, &n, x, col, t);
1541	mtree[k].next = -1;
1542	}
1543	for ( x = mtree[k].child; x != -1; x = mtree[x].brother )
1544	if ( mtree[x].begin == t->loc )
1545	{
1546	ENTER(&limit, &n, x, col, t);
1547	mtree[k].child = mtree[x].brother;
1548	}
1549	else
1550	break;
1551	}
1552	COLUMN(root); /* determine next column */
1553	root->c = root->kind;
1554	for ( t = head; t != 0; t = t->link )
1555	t->c = t->kind;
1556	for ( j = 0; j < n; j++ )
1557	{
1558	t = work[j];
1559	if ( t->done )
1560	*t->a++ = ' ';
1561	else
1562	{
1563	if ( t->c == 'L' )
1564	{
1565	if ( t->loc == 1 )
1566	t->offset = allconpt - spt;
1567	c = *t->a++ = t->seq[t->loc++];
1568	}
1569	else
1570	if ( t->loc > 1 )
1571	c = *t->a++ = '-';
1572	else
1573	c = *t->a++ = ' ';
1574	if ( c != ' ' )
1575	if ( c == '-' )
1576	sym[5] += 1;
1577	else
1578	sym[vert[c]] += 1;
1579	t->c = ' ';
1580	}
1581	}
1582	/* determine consensus char */
1583	k = sym[0] + sym[1] + sym[2] + sym[3] + sym[4];
1584	if ( k < sym[5] )
1585	*allconpt++ = '-';
1586	else
1587	if ( sym[0] == sym[1] && sym[1] == sym[2] &&
1588	sym[2] == sym[3] )
1589	*allconpt++ = 'N';
1590	else
1591	{
1592	k = sym[0];
1593	c = 'A';
1594	if ( k < sym[1] ) {
1595	k = sym[1];
1596	c = 'C';
1597	}
1598	if ( k < sym[2] ) {
1599	k = sym[2];
1600	c = 'G';
1601	}
1602	if ( k < sym[3] ) c = 'T';
1603	*allconpt++ = c;
1604	}
1605	for ( j = 0; j < 6; j++ )
1606	sym[j] = 0;
1607	for ( t = head; t != 0; t = t->link )
1608	{
1609	NEXTOP(t);
1610	if ( t->done ) /* delete it from op tree */
1611	{
1612	w = t->father;
1613	if ( w->child->id == t->id )
1614	w->child = t->brother;
1615	else
1616	{
1617	w = w->child;
1618	for ( ; w->brother->id != t->id; w = w->brother )
1619	;
1620	w->brother = t->brother;
1621	}
1622	}
1623	}
1624	if ( root->loc > root->length ) /* check root node */
1625	{
1626	root->done = 1;
1627	if ( (w = root->child) != 0 )
1628	{
1629	w->father = 0;
1630	root = w;
1631	}
1632	else
1633	done = 1;
1634	}
1635	col++;
1636	if ( col == LINELEN \|\| done ) /* output */
1637	{
1638	col = 0;
1639	for ( j = 0; j < n; j++ )
1640	{
1641	t = work[j];
1642	if ( t->done )
1643	/*
1644	Mod by S.S.
1645	{ (void) printf("#");
1646	for ( a = t->name; *a; a++ )
1647	(void) printf("%c", *a);
1648	*/
1649	{
1650	int jj;
1651	(void) printf("{\nname ");
1652	for(jj=0;jj<strlen(t->name)-1;jj++)
1653	(void) printf("%c", t->name[jj]);
1654	printf("\nstrandedness %c\n",
1655	t->name[strlen(t->name)] == '+'? '1':'2');
1656
1657	printf("offset %d\ntype DNA\ngroup-ID %d\nsequence \"\n",
1658	t->offset,group);
1659	for ( k = 0, a = t->line ; a != t->a; a++ )
1660	if ( *a != ' ' )
1661	{
1662	k++;
1663	(void) printf("%c", *a);
1664	if ( k == LINELEN )
1665	{
1666	(void) printf("\n");
1667	k = 0;
1668	}
1669	}
1670	/*
1671	if ( k )
1672	*/
1673	(void) printf("\"\n}\n");
1674	}
1675	if ( t->linesize - (t->a - t->line) < LINELEN + 3 )
1676	ALOC_SEQ(t);
1677	}
1678	if ( !done )
1679	{
1680	for ( k = j = n - 1; j >= 0; j-- )
1681	if ( work[j]->done )
1682	{
1683	t = work[j];
1684	for ( x = j; x < k; x++ )
1685	work[x] = work[x+1];
1686	work[k--] = t;
1687	}
1688	n = k + 1;
1689	}
1690	else
1691	n = 0;
1692	}
1693	}
1694	}
1695	}
1696
1697	/* allocate more space for output fragment */
1698	ALOC_SEQ(t) rowptr t;
1699	{
1700	char my_calloc(); / space-allocating function */
1701	char start, end, *p;
1702	t->linesize *= 2;
1703	start = t->line;
1704	end = t->a;
1705	t->line = ( char * ) my_calloc( t->linesize * sizeof(char));
1706	for ( t->a = t->line, p = start ; p != end; )
1707	t->a++ = p++;
1708	free(start);
1709	return 0;
1710	}
1711
1712	/* enter a segment into working set */
1713	ENTER(b, d, id, pos, par) int b, d, id, pos;
1714	rowptr par;
1715	{
1716	int i;
1717	char my_calloc(); / space-allocating function */
1718	rowptr t;
1719
1720	if ( b <= d )
1721	{
1722	work[*b] = ( rowptr ) my_calloc( (int ) sizeof(row));
1723	work[b]->line = ( char ) my_calloc( SEQLEN * sizeof(char));
1724	work[*b]->linesize = SEQLEN;
1725	*b += 1;
1726	}
1727	t = work[*d];
1728	*d += 1;
1729	t->a = t->line;
1730	for ( i = 0; i < pos; i++ )
1731	*t->a++ = ' ';
1732	t->c = ' ';
1733	t->seq = mtree[id].orient ? segment[id].seq : segment[id].rev;
1734	t->length = segment[id].length;
1735	t->id = id;
1736	if ( par != 0 )
1737	{
1738	t->s = mtree[id].script;
1739	t->size = mtree[id].size;
1740	}
1741	t->op = 0;
1742	for ( i = 1; i <= segment[id].len && i <= NAMELEN; i++ )
1743	t->name[i-1] = segment[id].name[i];
1744	if ( mtree[id].orient )
1745	t->name[i-1] = '+';
1746	else
1747	t->name[i-1] = '-';
1748	t->name[i] = '\0';
1749	t->done = 0;
1750	t->loc = 1;
1751	t->child = 0;
1752	t->father = par;
1753	if ( par != 0 )
1754	{
1755	t->brother = par->child;
1756	par->child = t;
1757	NEXTOP(t);
1758	}
1759	}
1760
1761	/* get the next operation */
1762	NEXTOP(t) rowptr t;
1763	{
1764	if ( t->size \|\| t->op )
1765	if ( t->op == 0 && *t->s == 0 )
1766	{
1767	t->op = *t->s++;
1768	t->size--;
1769	t->up = 'L';
1770	t->dw = 'L';
1771	}
1772	else
1773	{
1774	if ( t->op == 0 )
1775	{
1776	t->op = *t->s++;
1777	t->size--;
1778	}
1779	if ( t->op > 0 )
1780	{
1781	t->up = '-';
1782	t->dw = 'L';
1783	t->op--;
1784	}
1785	else
1786	{
1787	t->up = 'L';
1788	t->dw = '-';
1789	t->op++;
1790	}
1791	}
1792	else
1793	if ( t->loc > t->length )
1794	t->done = 1;
1795	}
1796
1797	COLUMN(x) rowptr x;
1798	{
1799	rowptr y;
1800	rowptr start, end; /* first and last nodes for subtree */
1801
1802	if ( x->child == 0 )
1803	{
1804	head = tail = 0;
1805	x->kind = 'L';
1806	}
1807	else
1808	{
1809	start = end = 0;
1810	x->kind = 'L';
1811	for ( y = x->child; y != 0; y = y->brother )
1812	{
1813	COLUMN(y);
1814	if ( x->kind == y->up )
1815	if ( y->kind == y->dw )
1816	{
1817	if ( head == 0 )
1818	{
1819	y->link = 0;
1820	head = tail = y;
1821	}
1822	else
1823	{
1824	y->link = head;
1825	head = y;
1826	}
1827	if ( end == 0 )
1828	start = head;
1829	else
1830	end->link = head;
1831	end = tail;
1832	}
1833	else
1834	if ( y->kind == '-' )
1835	{
1836	start = head;
1837	end = tail;
1838	x->kind = '-';
1839	}
1840	else
1841	{
1842	y->link = 0;
1843	y->kind = '-';
1844	if ( end == 0 )
1845	start = end = y;
1846	else
1847	{
1848	end->link = y;
1849	end = y;
1850	}
1851	}
1852	else
1853	if ( y->kind == y->dw )
1854	if ( x->kind == '-' )
1855	;
1856	else
1857	{
1858	if ( head == 0 )
1859	{
1860	y->link = 0;
1861	head = tail = y;
1862	}
1863	else
1864	{
1865	y->link = head;
1866	head = y;
1867	}
1868	start = head;
1869	end = tail;
1870	x->kind = '-';
1871	}
1872	else
1873	if ( x->kind == '-' )
1874	if ( y->kind == '-' )
1875	{
1876	if ( end == 0 )
1877	{
1878	start = head;
1879	end = tail;
1880	}
1881	else
1882	if ( head == 0 )
1883	/* code folded from here */
1884	;
1885	/* unfolding */
1886	else
1887	{
1888	/* code folded from here */
1889	end->link = head;
1890	end = tail;
1891	/* unfolding */
1892	}
1893	}
1894	else
1895	;
1896	else
1897	{
1898	start = head;
1899	end = tail;
1900	x->kind = '-';
1901	}
1902	}
1903	head = start;
1904	tail = end;
1905	}
1906	}
1907
1908	/* Display a summary of contigs */
1909	GRAPH()
1910	{
1911	int i, j, k; /* index variables */
1912	int group; /* serial number of contigs */
1913	char name[NAMELEN+2]; /* name of segment */
1914	char t; / temp var */
1915	int length; /* length of name */
1916
1917	(void) printf("\nOVERLAPS CONTAINMENTS\n\n");
1918	group = 1;
1919	for ( i = 0; i < seg_num; i++ )
1920	if ( mtree[i].head )
1921	{
1922	(void) printf("***************** Contig %d ******************\n",
1923	group++ );
1924	for ( j = i; j != -1; j = mtree[j].next )
1925	{
1926	length = segment[j].len;
1927	t = segment[j].name + 1;
1928	for ( k = 0; k < length && k < NAMELEN; k++ )
1929	name[k] = *t++;
1930	if ( mtree[j].orient )
1931	name[k] = '+';
1932	else
1933	name[k] = '-';
1934	name[k+1] = '\0';
1935	(void) printf("%s\n", name);
1936	CONTAIN(mtree[j].child, name);
1937	}
1938	}
1939	}
1940
1941	CONTAIN(id, f) int id;
1942	char *f;
1943	{
1944	int k; /* index variable */
1945	char name[NAMELEN+2]; /* name of segment */
1946	char t; / temp var */
1947	int length; /* length of name */
1948
1949	if ( id != -1 )
1950	{
1951	length = segment[id].len;
1952	t = segment[id].name + 1;
1953	for ( k = 0; k < length && k < NAMELEN; k++ )
1954	name[k] = *t++;
1955	if ( mtree[id].orient )
1956	name[k] = '+';
1957	else
1958	name[k] = '-';
1959	name[k+1] = '\0';
1960	(void) printf(" %s is in %s\n", name,f);
1961	CONTAIN(mtree[id].child, name);
1962	CONTAIN(mtree[id].brother, f);
1963	}
1964	}
1965
1966	big_pass(A,B,M,N,orienti,orientj) char A[],B[];
1967	int M,N;
1968	short orienti, orientj;
1969	{
1970	register int i, j; /* row and column indices */
1971	register int c; /* best score at current point */
1972	register int f; /* best score ending with insertion */
1973	register int d; /* best score ending with deletion */
1974	register int p; /* best score at (i-1, j-1) */
1975	register int ci; /* end-point associated with c */
1976
1977	register int di; /* end-point associated with d */
1978	register int fi; /* end-point associated with f */
1979	register int pi; /* end-point associated with p */
1980	int va; / pointer to v(A[i], B[j]) */
1981	int x1, x2; /* regions of A before x1 or after x2 are lightly penalized */
1982	int y1, y2; /* regions of B before y1 or after y2 are lightly penalized */
1983	short heavy; /* 1 = heavy penalty */
1984	int ex, gx; /* current gap penalty scores */
1985
1986	/* determine x1, x2, y1, y2 */
1987	if ( POS5 >= POS3 )
1988	fatal("The value for POS5 must be less than the value for POS3");
1989	if ( orienti )
1990	{
1991	x1 = POS5 >= M ? 1 : POS5;
1992	x2 = POS3 >= M ? M : POS3;
1993	}
1994	else
1995	{
1996	x1 = POS3 >= M ? 1 : M - POS3 + 1;
1997	x2 = POS5 >= M ? M : M - POS5 + 1;
1998	}
1999	if ( orientj )
2000	{
2001	y1 = POS5 >= N ? 1 : POS5;
2002	y2 = POS3 >= N ? N : POS3;
2003	}
2004	else
2005	{
2006	y1 = POS3 >= N ? 1 : N - POS3 + 1;
2007	y2 = POS5 >= N ? N : N - POS5 + 1;
2008	}
2009	if ( x1 + 1 <= x2 ) x1++;
2010	if ( y1 + 1 <= y2 ) y1++;
2011	heavy = 0;
2012
2013	/* Compute the matrix.
2014	CC : the scores of the current row
2015	RR : the starting point that leads to score CC
2016	DD : the scores of the current row, ending with deletion
2017	SS : the starting point that leads to score DD */
2018	/* Initialize the 0 th row */
2019	for ( j = 1; j <= N ; j++ )
2020	{
2021	CC[j] = 0;
2022	DD[j] = - (q);
2023	RR[j] = SS[j] = -j;
2024	}
2025	for ( i = 1; i <= M; i++)
2026	{
2027	if ( i == x1 ) heavy = 1 - heavy;
2028	if ( i == x2 ) heavy = 1 - heavy;
2029	ex = r1;
2030	gx = qr1;
2031	va = v1[A[i]];
2032	c = 0; /* Initialize column 0 */
2033	f = - (q);
2034	ci = fi = i;
2035	p = 0;
2036	pi = i - 1;
2037	for ( j = 1 ; j <= N ; j++ )
2038	{
2039	if ( j == y1 )
2040	{
2041	if ( heavy )
2042	{
2043	ex = r;
2044	gx = qr;
2045	/*
2046	S.S.
2047	va = v[A[i]];
2048	*/
2049	va = v[A[i]];
2050	}
2051	}
2052	if ( j == y2 )
2053	{
2054	if ( heavy )
2055	{
2056	ex = r1;
2057	gx = qr1;
2058	va = v1[A[i]];
2059	}
2060	}
2061	if ( ( f = f - ex ) < ( c = c - gx ) )
2062	{
2063	f = c;
2064	fi = ci;
2065	}
2066	di = SS[j];
2067	if ( ( d = DD[j] - ex ) < ( c = CC[j] - gx ) )
2068	{
2069	d = c;
2070	di = RR[j];
2071	}
2072	c = p+va[B[j]]; /* diagonal */
2073	ci = pi;
2074	if ( c < d )
2075	{
2076	c = d;
2077	ci = di;
2078	}
2079	if ( c < f )
2080	{
2081	c = f;
2082	ci = fi;
2083	}
2084	p = CC[j];
2085	CC[j] = c;
2086	pi = RR[j];
2087	RR[j] = ci;
2088	DD[j] = d;
2089	SS[j] = di;
2090	if ( ( j == N \|\| i == M ) && c > SCORE )
2091	{
2092	SCORE = c;
2093	ENDI = i;
2094	ENDJ = j;
2095	STARI = ci;
2096	}
2097	}
2098	}
2099	if ( SCORE )
2100	if ( STARI < 0 )
2101	{
2102	STARJ = - STARI;
2103	STARI = 0;
2104	}
2105	else
2106	STARJ = 0;
2107	}
2108
2109	/* diff(A,B,M,N,tb,te) returns the score of an optimum conversion between
2110	A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
2111	and appends such a conversion to the current script. */
2112
2113	int diff(A,B,M,N,tb,te) char A, B;
2114	int M, N;
2115	int tb, te;
2116
2117	{
2118	int midi, midj, type; /* Midpoint, type, and cost */
2119	int midc;
2120
2121	{
2122	register int i, j;
2123	register int c, e, d, s;
2124	int t, *va;
2125	char *my_calloc();
2126
2127	/* Boundary cases: M <= 1 or N == 0 */
2128
2129	if (N <= 0)
2130	{
2131	if (M > 0) DEL(M)
2132	return - gap(M);
2133	}
2134	if (M <= 1)
2135	{
2136	if (M <= 0)
2137	{
2138	INS(N);
2139	return - gap(N);
2140	}
2141	if (tb > te) tb = te;
2142	midc = - (tb + r + gap(N) );
2143	midj = 0;
2144	va = v[A[1]];
2145	for (j = 1; j <= N; j++)
2146	{
2147	c = va[B[j]] - ( gap(j-1) + gap(N-j) );
2148	if (c > midc)
2149	{
2150	midc = c;
2151	midj = j;
2152	}
2153	}
2154	if (midj == 0)
2155	{
2156	INS(N) DEL(1) }
2157	else
2158	{
2159	if (midj > 1) INS(midj-1)
2160	REP
2161	if ( (A[1]\|32) == (B[midj]\|32) )
2162	no_mat += 1;
2163	else
2164	no_mis += 1;
2165	if (midj < N) INS(N-midj)
2166	}
2167	return midc;
2168	}
2169
2170	/* Divide: Find optimum midpoint (midi,midj) of cost midc */
2171
2172	midi = M/2; /* Forward phase: */
2173	CC[0] = 0; /* Compute C(M/2,k) & D(M/2,k) for all k */
2174	t = -q;
2175	for (j = 1; j <= N; j++)
2176	{
2177	CC[j] = t = t-r;
2178	DD[j] = t-q;
2179	}
2180	t = -tb;
2181	for (i = 1; i <= midi; i++)
2182	{
2183	s = CC[0];
2184	CC[0] = c = t = t-r;
2185	e = t-q;
2186	va = v[A[i]];
2187	for (j = 1; j <= N; j++)
2188	{
2189	if ((c = c - qr) > (e = e - r)) e = c;
2190	if ((c = CC[j] - qr) > (d = DD[j] - r)) d = c;
2191	c = s+va[B[j]];
2192	if (c < d) c = d;
2193	if (c < e) c = e;
2194	s = CC[j];
2195	CC[j] = c;
2196	DD[j] = d;
2197	}
2198	}
2199	DD[0] = CC[0];
2200
2201	RR[N] = 0; /* Reverse phase: */
2202	t = -q; /* Compute R(M/2,k) & S(M/2,k) for all k */
2203	for (j = N-1; j >= 0; j--)
2204	{
2205	RR[j] = t = t-r;
2206	SS[j] = t-q;
2207	}
2208	t = -te;
2209	for (i = M-1; i >= midi; i--)
2210	{
2211	s = RR[N];
2212	RR[N] = c = t = t-r;
2213	e = t-q;
2214	va = v[A[i+1]];
2215	for (j = N-1; j >= 0; j--)
2216	{
2217	if ((c = c - qr) > (e = e - r)) e = c;
2218	if ((c = RR[j] - qr) > (d = SS[j] - r)) d = c;
2219	c = s+va[B[j+1]];
2220	if (c < d) c = d;
2221	if (c < e) c = e;
2222	s = RR[j];
2223	RR[j] = c;
2224	SS[j] = d;
2225	}
2226	}
2227	SS[N] = RR[N];
2228
2229	midc = CC[0]+RR[0]; /* Find optimal midpoint */
2230	midj = 0;
2231	type = 1;
2232	for (j = 0; j <= N; j++)
2233	if ((c = CC[j] + RR[j]) >= midc)
2234	if (c > midc \|\| CC[j] != DD[j] && RR[j] == SS[j])
2235	{
2236	midc = c;
2237	midj = j;
2238	}
2239	for (j = N; j >= 0; j--)
2240	if ((c = DD[j] + SS[j] + q) > midc)
2241	{
2242	midc = c;
2243	midj = j;
2244	type = 2;
2245	}
2246	}
2247
2248	/* Conquer: recursively around midpoint */
2249
2250	if (type == 1)
2251	{
2252	(void) diff(A,B,midi,midj,tb,q);
2253	(void) diff(A+midi,B+midj,M-midi,N-midj,q,te);
2254	}
2255	else
2256	{
2257	(void) diff(A,B,midi-1,midj,tb,zero);
2258	DEL(2);
2259	(void) diff(A+midi+1,B+midj,M-midi-1,N-midj,zero,te);
2260	}
2261	return midc;
2262	}
2263
2264	/* lib.c - library of C procedures. */
2265
2266	/* fatal - print message and die */
2267	fatal(msg)
2268	char *msg;
2269	{
2270	(void) fprintf(stderr, "%s\n", msg);
2271	exit(1);
2272	}
2273
2274	/* fatalf - format message, print it, and die */
2275	fatalf(msg, val)
2276	char msg, val;
2277	{
2278	(void) fprintf(stderr, msg, val);
2279	(void) putc('\n', stderr);
2280	exit(1);
2281	}
2282
2283	/* ckopen - open file; check for success */
2284	FILE *ckopen(name, mode)
2285	char name, mode;
2286	{
2287	FILE fopen(), fp;
2288
2289	if ((fp = fopen(name, mode)) == NULL)
2290	fatalf("Cannot open %s.", name);
2291	return(fp);
2292	}
2293
2294	/* my_calloc - allocate space; check for success */
2295	char *my_calloc(amount)
2296	int amount;
2297	{
2298	char malloc(), p;
2299
2300	if ((p = malloc( (unsigned) amount)) == NULL)
2301	fatal("Ran out of memory.");
2302	return(p);
2303	}
2304

Note: See TracBrowser for help on using the repository browser.

Download in other formats: