[Bio] / FigKernelScripts / fasta_cksum.c Repository:
ViewVC logotype

Annotation of /FigKernelScripts/fasta_cksum.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (view) (download) (as text)

1 : olson 1.4 /*
2 :     * Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     * for Interpretations of Genomes. All Rights Reserved.
4 :     *
5 :     * This file is part of the SEED Toolkit.
6 :     *
7 :     * The SEED Toolkit is free software. You can redistribute
8 :     * it and/or modify it under the terms of the SEED Toolkit
9 :     * Public License.
10 :     *
11 :     * You should have received a copy of the SEED Toolkit Public License
12 :     * along with this program; if not write to the University of Chicago
13 :     * at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     * Genomes at veronika@thefig.info or download a copy from
15 :     * http://www.theseed.org/LICENSE.TXT.
16 :     */
17 :    
18 : olson 1.3
19 : olson 1.2 /* fasta_cksum.c
20 :     *
21 :     * Usage: fasta_cksum -t < fasta_file > nseq_nresidues_cksum
22 :     * or fasta_cksum < fasta_file > id_nresidues_cksum...
23 :     *
24 :     /* These include files are appropriate for Machintosh OS X */
25 :    
26 :     #include <stdio.h>
27 :     #include <ctype.h> /* isspace() */
28 :     #include <stdlib.h> /* exit() */
29 :     #include <unistd.h> /* read() */
30 :    
31 :     #define BUFLEN (256*1024)
32 :     #define INPLEN ( 64*1024)
33 :     #define IDLEN ( 16*1024)
34 :     #define DFLT_INDEX_INTERVAL 10000
35 :    
36 :     #define fillbuf(buf, len) read( (int) 0, (void *) buf, (size_t) len )
37 :    
38 :     /* From the discussion of cksum:
39 :     *
40 :     * The Open Group Base Specifications Issue 6
41 :     * IEEE Std 1003.1, 2004 Edition
42 :     * Copyright 2001-2004 The IEEE and The Open Group, All Rights reserved.
43 :     *
44 :     * http://www.opengroup.org/onlinepubs/009695399/utilities/cksum.html
45 :     */
46 :    
47 :     static unsigned long crctab[] = {
48 :     0x00000000,
49 :     0x04c11db7, 0x09823b6e, 0x0d4326d9, 0x130476dc, 0x17c56b6b,
50 :     0x1a864db2, 0x1e475005, 0x2608edb8, 0x22c9f00f, 0x2f8ad6d6,
51 :     0x2b4bcb61, 0x350c9b64, 0x31cd86d3, 0x3c8ea00a, 0x384fbdbd,
52 :     0x4c11db70, 0x48d0c6c7, 0x4593e01e, 0x4152fda9, 0x5f15adac,
53 :     0x5bd4b01b, 0x569796c2, 0x52568b75, 0x6a1936c8, 0x6ed82b7f,
54 :     0x639b0da6, 0x675a1011, 0x791d4014, 0x7ddc5da3, 0x709f7b7a,
55 :     0x745e66cd, 0x9823b6e0, 0x9ce2ab57, 0x91a18d8e, 0x95609039,
56 :     0x8b27c03c, 0x8fe6dd8b, 0x82a5fb52, 0x8664e6e5, 0xbe2b5b58,
57 :     0xbaea46ef, 0xb7a96036, 0xb3687d81, 0xad2f2d84, 0xa9ee3033,
58 :     0xa4ad16ea, 0xa06c0b5d, 0xd4326d90, 0xd0f37027, 0xddb056fe,
59 :     0xd9714b49, 0xc7361b4c, 0xc3f706fb, 0xceb42022, 0xca753d95,
60 :     0xf23a8028, 0xf6fb9d9f, 0xfbb8bb46, 0xff79a6f1, 0xe13ef6f4,
61 :     0xe5ffeb43, 0xe8bccd9a, 0xec7dd02d, 0x34867077, 0x30476dc0,
62 :     0x3d044b19, 0x39c556ae, 0x278206ab, 0x23431b1c, 0x2e003dc5,
63 :     0x2ac12072, 0x128e9dcf, 0x164f8078, 0x1b0ca6a1, 0x1fcdbb16,
64 :     0x018aeb13, 0x054bf6a4, 0x0808d07d, 0x0cc9cdca, 0x7897ab07,
65 :     0x7c56b6b0, 0x71159069, 0x75d48dde, 0x6b93dddb, 0x6f52c06c,
66 :     0x6211e6b5, 0x66d0fb02, 0x5e9f46bf, 0x5a5e5b08, 0x571d7dd1,
67 :     0x53dc6066, 0x4d9b3063, 0x495a2dd4, 0x44190b0d, 0x40d816ba,
68 :     0xaca5c697, 0xa864db20, 0xa527fdf9, 0xa1e6e04e, 0xbfa1b04b,
69 :     0xbb60adfc, 0xb6238b25, 0xb2e29692, 0x8aad2b2f, 0x8e6c3698,
70 :     0x832f1041, 0x87ee0df6, 0x99a95df3, 0x9d684044, 0x902b669d,
71 :     0x94ea7b2a, 0xe0b41de7, 0xe4750050, 0xe9362689, 0xedf73b3e,
72 :     0xf3b06b3b, 0xf771768c, 0xfa325055, 0xfef34de2, 0xc6bcf05f,
73 :     0xc27dede8, 0xcf3ecb31, 0xcbffd686, 0xd5b88683, 0xd1799b34,
74 :     0xdc3abded, 0xd8fba05a, 0x690ce0ee, 0x6dcdfd59, 0x608edb80,
75 :     0x644fc637, 0x7a089632, 0x7ec98b85, 0x738aad5c, 0x774bb0eb,
76 :     0x4f040d56, 0x4bc510e1, 0x46863638, 0x42472b8f, 0x5c007b8a,
77 :     0x58c1663d, 0x558240e4, 0x51435d53, 0x251d3b9e, 0x21dc2629,
78 :     0x2c9f00f0, 0x285e1d47, 0x36194d42, 0x32d850f5, 0x3f9b762c,
79 :     0x3b5a6b9b, 0x0315d626, 0x07d4cb91, 0x0a97ed48, 0x0e56f0ff,
80 :     0x1011a0fa, 0x14d0bd4d, 0x19939b94, 0x1d528623, 0xf12f560e,
81 :     0xf5ee4bb9, 0xf8ad6d60, 0xfc6c70d7, 0xe22b20d2, 0xe6ea3d65,
82 :     0xeba91bbc, 0xef68060b, 0xd727bbb6, 0xd3e6a601, 0xdea580d8,
83 :     0xda649d6f, 0xc423cd6a, 0xc0e2d0dd, 0xcda1f604, 0xc960ebb3,
84 :     0xbd3e8d7e, 0xb9ff90c9, 0xb4bcb610, 0xb07daba7, 0xae3afba2,
85 :     0xaafbe615, 0xa7b8c0cc, 0xa379dd7b, 0x9b3660c6, 0x9ff77d71,
86 :     0x92b45ba8, 0x9675461f, 0x8832161a, 0x8cf30bad, 0x81b02d74,
87 :     0x857130c3, 0x5d8a9099, 0x594b8d2e, 0x5408abf7, 0x50c9b640,
88 :     0x4e8ee645, 0x4a4ffbf2, 0x470cdd2b, 0x43cdc09c, 0x7b827d21,
89 :     0x7f436096, 0x7200464f, 0x76c15bf8, 0x68860bfd, 0x6c47164a,
90 :     0x61043093, 0x65c52d24, 0x119b4be9, 0x155a565e, 0x18197087,
91 :     0x1cd86d30, 0x029f3d35, 0x065e2082, 0x0b1d065b, 0x0fdc1bec,
92 :     0x3793a651, 0x3352bbe6, 0x3e119d3f, 0x3ad08088, 0x2497d08d,
93 :     0x2056cd3a, 0x2d15ebe3, 0x29d4f654, 0xc5a92679, 0xc1683bce,
94 :     0xcc2b1d17, 0xc8ea00a0, 0xd6ad50a5, 0xd26c4d12, 0xdf2f6bcb,
95 :     0xdbee767c, 0xe3a1cbc1, 0xe760d676, 0xea23f0af, 0xeee2ed18,
96 :     0xf0a5bd1d, 0xf464a0aa, 0xf9278673, 0xfde69bc4, 0x89b8fd09,
97 :     0x8d79e0be, 0x803ac667, 0x84fbdbd0, 0x9abc8bd5, 0x9e7d9662,
98 :     0x933eb0bb, 0x97ffad0c, 0xafb010b1, 0xab710d06, 0xa6322bdf,
99 :     0xa2f33668, 0xbcb4666d, 0xb8757bda, 0xb5365d03, 0xb1f740b4
100 :     };
101 :    
102 :    
103 :     /* Function prototypes: */
104 :    
105 :     void report_seq( char * id, unsigned long seqlen, unsigned long crc );
106 :    
107 :     void report_ttl( int n_seq, unsigned long long ttllen, unsigned long ttlcrc );
108 :    
109 :     unsigned long finish_cksum( unsigned long crc, unsigned long seqlen );
110 :    
111 :     void usage( char *prog );
112 :    
113 :    
114 :     unsigned char buffer[BUFLEN];
115 :     char idbuf[IDLEN+1];
116 :    
117 :    
118 :     int main ( int argc, char **argv ) {
119 :     unsigned long long ttllen;
120 :     unsigned char *bptr;
121 :     unsigned long c, seqlen, crc, ttlcrc;
122 :     int totalonly, n_seq, idlen, ntogo;
123 :    
124 :     /* -t flag returns only the total */
125 :    
126 :     totalonly = 0;
127 :     if ( ( argc > 1 ) && ( argv[1][0] == '-' )
128 :     && ( argv[1][1] == 't' )
129 :     && ( argv[1][2] == '\0' )
130 :     ) {
131 :     totalonly = 1;
132 :     argc--;
133 :     }
134 :     else if (argc > 1 ) { usage( argv[0] ); }
135 :    
136 :     idbuf[0] = '\0'; /* initialize to empty string */
137 :     bptr = buffer; /* pointer to next character in buffer */
138 :     ntogo = 0; /* unused characters in input buffer */
139 :     n_seq = 0;
140 :     seqlen = 0;
141 :     ttllen = 0;
142 :     crc = ~0; /* this will get "finished" to 0 */
143 :     ttlcrc = 0;
144 :    
145 :     /* Process one line of input */
146 :    
147 :     while ( 1 ) {
148 :     if ( ntogo <= 0 ) {
149 :     if ( ( ntogo = fillbuf( buffer, BUFLEN ) ) <= 0 ) {
150 :     if ( totalonly ) {
151 :     ttllen += seqlen;
152 :     ttlcrc ^= finish_cksum( crc, seqlen );
153 :     report_ttl( n_seq, ttllen, ttlcrc );
154 :     }
155 :     else report_seq( idbuf, seqlen, crc );
156 :     exit( ntogo );
157 :     }
158 :     bptr= buffer;
159 :     }
160 :     c = *bptr++; ntogo--;
161 :    
162 :     /* Line could start with >, or be sequence data */
163 :    
164 :     if ( c == '>' ) {
165 :    
166 :     /* New sequence. Is there an previous sequence to report? */
167 :    
168 :     if ( ! totalonly ) report_seq( idbuf, seqlen, crc );
169 :    
170 :     /* Adjust cumulative values and reset sequence values */
171 :    
172 :     ttllen += seqlen;
173 :     ttlcrc ^= finish_cksum( crc, seqlen );
174 :     seqlen = crc = 0;
175 :    
176 :     if ( ntogo <= 0 ) {
177 :     if ( ( ntogo = fillbuf( buffer, BUFLEN ) ) <= 0 ) {
178 :     if ( totalonly ) report_ttl( n_seq, ttllen, ttlcrc );
179 :     exit( ntogo );
180 :     }
181 :     bptr = buffer;
182 :     }
183 :     c = *bptr++; ntogo--;
184 :    
185 :     /* Make a copy of the new id */
186 :    
187 :     idlen = 0;
188 :     while ( ( ! isspace(c) ) && ( idlen < IDLEN ) ) {
189 :     idbuf[ idlen++ ] = c;
190 :     if ( ntogo <= 0 ) {
191 :     if ( ( ntogo = fillbuf( buffer, BUFLEN ) ) <= 0 ) {
192 :     if ( totalonly ) report_ttl( n_seq, ttllen, ttlcrc );
193 :     exit(ntogo);
194 :     }
195 :     bptr = buffer;
196 :     }
197 :     c = *bptr++; ntogo--;
198 :     }
199 :     idbuf[ idlen ] = '\0';
200 :    
201 :     /* report truncated id */
202 :    
203 :     if ( ! isspace(c) ) {
204 :     fprintf( stderr, "Sequence id truncated to %d characters:\n", (int) IDLEN );
205 :     fprintf( stderr, ">%s\n", idbuf );
206 :     }
207 :    
208 :     /* Flush the rest of the input line */
209 :    
210 :     while ( c != '\n' ) {
211 :     if ( ntogo <= 0 ) {
212 :     if ( ( ntogo = fillbuf( buffer, BUFLEN ) ) <= 0 ) {
213 :     if ( totalonly ) report_ttl( n_seq, ttllen, ttlcrc );
214 :     exit(ntogo);
215 :     }
216 :     bptr = buffer;
217 :     }
218 :     c = *bptr++; ntogo--;
219 :     }
220 :    
221 :     n_seq++; /* First data for a new sequence */
222 :     }
223 :    
224 :     /* Not an id line, so it's data: */
225 :    
226 :     else {
227 :     while ( c != '\n' ) { /* finish the line */
228 :     if ( ! isspace( c ) ) {
229 :     seqlen++;
230 :     crc = (crc << 8) ^ crctab[ (crc >> 24) ^ c ];
231 :     }
232 :    
233 :     /* Next character */
234 :    
235 :     if ( ntogo <= 0 ) {
236 :     if ( ( ntogo = fillbuf( buffer, BUFLEN ) ) <= 0 ) {
237 :     if ( totalonly ) {
238 :     ttllen += seqlen;
239 :     ttlcrc ^= finish_cksum( crc, seqlen );
240 :     report_ttl( n_seq, ttllen, ttlcrc );
241 :     }
242 :     else report_seq( idbuf, seqlen, crc );
243 :     exit( ntogo );
244 :     }
245 :     bptr = buffer;
246 :     }
247 :     c = *bptr++; ntogo--;
248 :     }
249 :     }
250 :     /* Go back to top to process next line */
251 :     }
252 :    
253 :     exit( 0 ); /* never get here */
254 :     }
255 :    
256 :    
257 :     void report_seq( char * id, unsigned long seqlen, unsigned long crc ) {
258 :     if ( id && id[0] ) {
259 :     printf( "%s\t%lu\t%lu\n", id, seqlen, finish_cksum( crc, seqlen ) );
260 :     }
261 :     }
262 :    
263 :    
264 :     /* Finish the crc calculation with length bytes, and bitwise complement */
265 :    
266 :     unsigned long finish_cksum( unsigned long crc, unsigned long seqlen ) {
267 :     while ( seqlen != 0 ) {
268 :     crc = (crc << 8) ^ crctab[ (crc >> 24) ^ (seqlen & 0xFF) ];
269 :     seqlen >>= 8;
270 :     }
271 :     return ~crc;
272 :     }
273 :    
274 :    
275 :     void report_ttl( int n_seq, unsigned long long ttllen, unsigned long ttlcrc ) {
276 :     printf( "%d\t%llu\t%lu\n", n_seq, ttllen, ttlcrc );
277 :     }
278 :    
279 :    
280 :     void usage( char *prog ) {
281 :     fprintf( stderr,
282 :     "Usage: %s -t < fasta_file > nseq \\t nresidues \\t cksum\n"
283 :     "or %s < fasta_file > id \\t nresidues \\t cksum \\n ...\n",
284 :     prog, prog
285 :     );
286 :     exit(0);
287 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3