1
|
1 /*****************************************************************
|
|
2 * SQUID - a library of functions for biological sequence analysis
|
|
3 * Copyright (C) 1992-2002 Washington University School of Medicine
|
|
4 *
|
|
5 * This source code is freely distributed under the terms of the
|
|
6 * GNU General Public License. See the files COPYRIGHT and LICENSE
|
|
7 * for details.
|
|
8 *****************************************************************/
|
|
9
|
|
10 /* File: sqio.c
|
|
11 * From: ureadseq.c in Don Gilbert's sequence i/o package
|
|
12 *
|
|
13 * Reads and writes nucleic/protein sequence in various
|
|
14 * formats. Data files may have multiple sequences.
|
|
15 *
|
|
16 * Heavily modified from READSEQ package
|
|
17 * Copyright (C) 1990 by D.G. Gilbert
|
|
18 * Biology Dept., Indiana University, Bloomington, IN 47405
|
|
19 * email: gilbertd@bio.indiana.edu
|
|
20 * Thanks Don!
|
|
21 *
|
|
22 * SRE: Modifications as noted. Fri Jul 3 09:44:54 1992
|
|
23 * Packaged for squid, Thu Oct 1 10:07:11 1992
|
|
24 * ANSI conversion in full swing, Mon Jul 12 12:22:21 1993
|
|
25 *
|
|
26 * CVS $Id: sqio.c,v 1.29 2002/08/26 23:10:52 eddy Exp)
|
|
27 *
|
|
28 *****************************************************************
|
|
29 * Basic API for single sequence reading:
|
|
30 *
|
|
31 * SQFILE *sqfp;
|
|
32 * char *seqfile;
|
|
33 * int format; - see squid.h for formats; example: SQFILE_FASTA
|
|
34 * char *seq;
|
|
35 * SQINFO sqinfo;
|
|
36 *
|
|
37 * if ((sqfp = SeqfileOpen(seqfile, format, "BLASTDB")) == NULL)
|
|
38 * Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
|
|
39 * while (ReadSeq(sqfp, sqfp->format, &seq, &sqinfo)) {
|
|
40 * do_stuff;
|
|
41 * FreeSequence(seq, &sqinfo);
|
|
42 * }
|
|
43 * SeqfileClose(sqfp);
|
|
44 *
|
|
45 *****************************************************************
|
|
46 */
|
|
47
|
|
48 #include <stdio.h>
|
|
49 #include <stdlib.h>
|
|
50 #include <string.h>
|
|
51 #include <ctype.h>
|
|
52
|
|
53 #ifndef SEEK_SET
|
|
54 #include <unistd.h>
|
|
55 #endif
|
|
56
|
|
57 #include "squid.h"
|
|
58 #include "msa.h"
|
|
59 #include "ssi.h"
|
|
60
|
|
61 static void SeqfileGetLine(SQFILE *V);
|
|
62
|
|
63 #define kStartLength 500
|
|
64
|
|
65 static char *aminos = "ABCDEFGHIKLMNPQRSTVWXYZ*";
|
|
66 static char *primenuc = "ACGTUN";
|
|
67 static char *protonly = "EFIPQZ";
|
|
68
|
|
69 static SQFILE *seqfile_open(char *filename, int format, char *env, int ssimode);
|
|
70
|
|
71 /* Function: SeqfileOpen()
|
|
72 *
|
|
73 * Purpose : Open a sequence database file and prepare for reading
|
|
74 * sequentially.
|
|
75 *
|
|
76 * Args: filename - name of file to open
|
|
77 * format - format of file
|
|
78 * env - environment variable for path (e.g. BLASTDB)
|
|
79 * ssimode - -1, SSI_OFFSET_I32, or SSI_OFFSET_I64
|
|
80 *
|
|
81 * Returns opened SQFILE ptr, or NULL on failure.
|
|
82 */
|
|
83 SQFILE *
|
|
84 SeqfileOpen(char *filename, int format, char *env)
|
|
85 {
|
|
86 return seqfile_open(filename, format, env, -1);
|
|
87 }
|
|
88 SQFILE *
|
|
89 SeqfileOpenForIndexing(char *filename, int format, char *env, int ssimode)
|
|
90 {
|
|
91 return seqfile_open(filename, format, env, ssimode);
|
|
92 }
|
|
93 static SQFILE *
|
|
94 seqfile_open(char *filename, int format, char *env, int ssimode)
|
|
95 {
|
|
96 SQFILE *dbfp;
|
|
97
|
|
98 dbfp = (SQFILE *) MallocOrDie (sizeof(SQFILE));
|
|
99
|
|
100 dbfp->ssimode = ssimode;
|
|
101 dbfp->rpl = -1; /* flag meaning "unset" */
|
|
102 dbfp->lastrpl = 0;
|
|
103 dbfp->maxrpl = 0;
|
|
104 dbfp->bpl = -1; /* flag meaning "unset" */
|
|
105 dbfp->lastbpl = 0;
|
|
106 dbfp->maxbpl = 0;
|
|
107
|
|
108 /* Open our file handle.
|
|
109 * Three possibilities:
|
|
110 * 1. normal file open
|
|
111 * 2. filename = "-"; read from stdin
|
|
112 * 3. filename = "*.gz"; read thru pipe from gzip
|
|
113 * If we're reading from stdin or a pipe, we can't reliably
|
|
114 * back up, so we can't do two-pass parsers like the interleaved alignment
|
|
115 * formats.
|
|
116 */
|
|
117 if (strcmp(filename, "-") == 0)
|
|
118 {
|
|
119 dbfp->f = stdin;
|
|
120 dbfp->do_stdin = TRUE;
|
|
121 dbfp->do_gzip = FALSE;
|
|
122 dbfp->fname = sre_strdup("[STDIN]", -1);
|
|
123 }
|
|
124 #ifndef SRE_STRICT_ANSI
|
|
125 /* popen(), pclose() aren't portable to non-POSIX systems; disable */
|
|
126 else if (Strparse("^.*\\.gz$", filename, 0))
|
|
127 {
|
|
128 char cmd[256];
|
|
129
|
|
130 /* Note that popen() will return "successfully"
|
|
131 * if file doesn't exist, because gzip works fine
|
|
132 * and prints an error! So we have to check for
|
|
133 * existence of file ourself.
|
|
134 */
|
|
135 if (! FileExists(filename))
|
|
136 Die("%s: file does not exist", filename);
|
|
137
|
|
138 if (strlen(filename) + strlen("gzip -dc ") >= 256)
|
|
139 Die("filename > 255 char in SeqfileOpen()");
|
|
140 sprintf(cmd, "gzip -dc %s", filename);
|
|
141 if ((dbfp->f = popen(cmd, "r")) == NULL)
|
|
142 return NULL;
|
|
143
|
|
144 dbfp->do_stdin = FALSE;
|
|
145 dbfp->do_gzip = TRUE;
|
|
146 dbfp->fname = sre_strdup(filename, -1);
|
|
147 }
|
|
148 #endif /*SRE_STRICT_ANSI*/
|
|
149 else
|
|
150 {
|
|
151 if ((dbfp->f = fopen(filename, "r")) == NULL &&
|
|
152 (dbfp->f = EnvFileOpen(filename, env, NULL)) == NULL)
|
|
153 return NULL;
|
|
154
|
|
155 dbfp->do_stdin = FALSE;
|
|
156 dbfp->do_gzip = FALSE;
|
|
157 dbfp->fname = sre_strdup(filename, -1);
|
|
158 }
|
|
159
|
|
160
|
|
161 /* Invoke autodetection if we haven't already been told what
|
|
162 * to expect.
|
|
163 */
|
|
164 if (format == SQFILE_UNKNOWN)
|
|
165 {
|
|
166 if (dbfp->do_stdin == TRUE || dbfp->do_gzip)
|
|
167 Die("Can't autodetect sequence file format from a stdin or gzip pipe");
|
|
168 format = SeqfileFormat(dbfp->f);
|
|
169 if (format == SQFILE_UNKNOWN)
|
|
170 Die("Can't determine format of sequence file %s", dbfp->fname);
|
|
171 }
|
|
172
|
|
173 /* The hack for sequential access of an interleaved alignment file:
|
|
174 * read the alignment in, we'll copy sequences out one at a time.
|
|
175 */
|
|
176 dbfp->msa = NULL;
|
|
177 dbfp->afp = NULL;
|
|
178 dbfp->format = format;
|
|
179 dbfp->linenumber = 0;
|
|
180 dbfp->buf = NULL;
|
|
181 dbfp->buflen = 0;
|
|
182 if (IsAlignmentFormat(format))
|
|
183 {
|
|
184 /* We'll be reading from the MSA interface. Copy our data
|
|
185 * to the MSA afp's structure.
|
|
186 */
|
|
187 dbfp->afp = MallocOrDie(sizeof(MSAFILE));
|
|
188 dbfp->afp->f = dbfp->f; /* just a ptr, don't close */
|
|
189 dbfp->afp->do_stdin = dbfp->do_stdin;
|
|
190 dbfp->afp->do_gzip = dbfp->do_gzip;
|
|
191 dbfp->afp->fname = dbfp->fname; /* just a ptr, don't free */
|
|
192 dbfp->afp->format = dbfp->format; /* e.g. format */
|
|
193 dbfp->afp->linenumber = dbfp->linenumber; /* e.g. 0 */
|
|
194 dbfp->afp->buf = NULL;
|
|
195 dbfp->afp->buflen = 0;
|
|
196
|
|
197 if ((dbfp->msa = MSAFileRead(dbfp->afp)) == NULL)
|
|
198 Die("Failed to read any alignment data from file %s", dbfp->fname);
|
|
199 /* hack: overload/reuse msa->lastidx; indicates
|
|
200 next seq to return upon a ReadSeq() call */
|
|
201 dbfp->msa->lastidx = 0;
|
|
202
|
|
203 return dbfp;
|
|
204 }
|
|
205
|
|
206 /* Load the first line.
|
|
207 */
|
|
208 SeqfileGetLine(dbfp);
|
|
209 return dbfp;
|
|
210 }
|
|
211
|
|
212 /* Function: SeqfilePosition()
|
|
213 *
|
|
214 * Purpose: Move to a particular offset in a seqfile.
|
|
215 * Will not work on alignment files.
|
|
216 */
|
|
217 void
|
|
218 SeqfilePosition(SQFILE *sqfp, SSIOFFSET *offset)
|
|
219 {
|
|
220 if (sqfp->do_stdin || sqfp->do_gzip || IsAlignmentFormat(sqfp->format))
|
|
221 Die("SeqfilePosition() failed: in a nonrewindable data file or stream");
|
|
222
|
|
223 if (SSISetFilePosition(sqfp->f, offset) != 0)
|
|
224 Die("SSISetFilePosition failed, but that shouldn't happen.");
|
|
225 SeqfileGetLine(sqfp);
|
|
226 }
|
|
227
|
|
228
|
|
229 /* Function: SeqfileRewind()
|
|
230 *
|
|
231 * Purpose: Set a sequence file back to the first sequence.
|
|
232 *
|
|
233 * Won't work on alignment files. Although it would
|
|
234 * seem that it could (just set msa->lastidx back to 0),
|
|
235 * that'll fail on "multiple multiple" alignment file formats
|
|
236 * (e.g. Stockholm).
|
|
237 */
|
|
238 void
|
|
239 SeqfileRewind(SQFILE *sqfp)
|
|
240 {
|
|
241 if (sqfp->do_stdin || sqfp->do_gzip)
|
|
242 Die("SeqfileRewind() failed: in a nonrewindable data file or stream");
|
|
243
|
|
244 rewind(sqfp->f);
|
|
245 SeqfileGetLine(sqfp);
|
|
246 }
|
|
247
|
|
248 /* Function: SeqfileLineParameters()
|
|
249 * Date: SRE, Thu Feb 15 17:00:41 2001 [St. Louis]
|
|
250 *
|
|
251 * Purpose: After all the sequences have been read from the file,
|
|
252 * but before closing it, retrieve overall bytes-per-line and
|
|
253 * residues-per-line info. If non-zero, these mean that
|
|
254 * the file contains homogeneous sequence line lengths (except
|
|
255 * the last line in each record).
|
|
256 *
|
|
257 * If either of bpl or rpl is determined to be inhomogeneous,
|
|
258 * both are returned as 0.
|
|
259 *
|
|
260 * Args: *sqfp - an open but fully read sequence file
|
|
261 * ret_bpl - RETURN: bytes per line, or 0 if inhomogeneous
|
|
262 * ret_rpl - RETURN: residues per line, or 0 if inhomogenous.
|
|
263 *
|
|
264 * Returns: void
|
|
265 */
|
|
266 void
|
|
267 SeqfileLineParameters(SQFILE *V, int *ret_bpl, int *ret_rpl)
|
|
268 {
|
|
269 if (V->rpl > 0 && V->maxrpl == V->rpl &&
|
|
270 V->bpl > 0 && V->maxbpl == V->bpl) {
|
|
271 *ret_bpl = V->bpl;
|
|
272 *ret_rpl = V->rpl;
|
|
273 } else {
|
|
274 *ret_bpl = 0;
|
|
275 *ret_rpl = 0;
|
|
276 }
|
|
277 }
|
|
278
|
|
279
|
|
280 void
|
|
281 SeqfileClose(SQFILE *sqfp)
|
|
282 {
|
|
283 /* note: don't test for sqfp->msa being NULL. Now that
|
|
284 * we're holding afp open and allowing access to multi-MSA
|
|
285 * databases (e.g. Stockholm format, Pfam), msa ends
|
|
286 * up being NULL when we run out of alignments.
|
|
287 */
|
|
288 if (sqfp->afp != NULL) {
|
|
289 if (sqfp->msa != NULL) MSAFree(sqfp->msa);
|
|
290 if (sqfp->afp->buf != NULL) free(sqfp->afp->buf);
|
|
291 free(sqfp->afp);
|
|
292 }
|
|
293 #ifndef SRE_STRICT_ANSI /* gunzip functionality only on POSIX systems */
|
|
294 if (sqfp->do_gzip) pclose(sqfp->f);
|
|
295 #endif
|
|
296 else if (! sqfp->do_stdin) fclose(sqfp->f);
|
|
297 if (sqfp->buf != NULL) free(sqfp->buf);
|
|
298 if (sqfp->fname != NULL) free(sqfp->fname);
|
|
299 free(sqfp);
|
|
300 }
|
|
301
|
|
302
|
|
303 /* Function: SeqfileGetLine()
|
|
304 * Date: SRE, Tue Jun 22 09:15:49 1999 [Sanger Centre]
|
|
305 *
|
|
306 * Purpose: read a line from a sequence file into V->buf
|
|
307 * If the fgets() is NULL, sets V->buf[0] to '\0'.
|
|
308 *
|
|
309 * Args: V
|
|
310 *
|
|
311 * Returns: void
|
|
312 */
|
|
313 static void
|
|
314 SeqfileGetLine(SQFILE *V)
|
|
315 {
|
|
316 if (V->ssimode >= 0)
|
|
317 if (0 != SSIGetFilePosition(V->f, V->ssimode, &(V->ssioffset)))
|
|
318 Die("SSIGetFilePosition() failed");
|
|
319 if (sre_fgets(&(V->buf), &(V->buflen), V->f) == NULL)
|
|
320 *(V->buf) = '\0';
|
|
321 V->linenumber++;
|
|
322 }
|
|
323
|
|
324
|
|
325 void
|
|
326 FreeSequence(char *seq, SQINFO *sqinfo)
|
|
327 {
|
|
328 if (seq != NULL) free(seq); /* FS, r244, here is potential problem in profile/profile */
|
|
329 if (sqinfo->flags & SQINFO_SS){
|
|
330 if (NULL != sqinfo->ss){ /* FS, r244 -> r245 */
|
|
331 free(sqinfo->ss);
|
|
332 }
|
|
333 }
|
|
334 if (sqinfo->flags & SQINFO_SA){
|
|
335 if (NULL != sqinfo->sa){ /* FS, r244 -> r245 */
|
|
336 free(sqinfo->sa);
|
|
337 }
|
|
338 }
|
|
339 }
|
|
340
|
|
341 int
|
|
342 SetSeqinfoString(SQINFO *sqinfo, char *sptr, int flag)
|
|
343 {
|
|
344 int len;
|
|
345 int pos;
|
|
346
|
|
347 /* silently ignore NULL. */
|
|
348 if (sptr == NULL) return 1;
|
|
349
|
|
350 while (*sptr == ' ') sptr++; /* ignore leading whitespace */
|
|
351 for (pos = strlen(sptr)-1; pos >= 0; pos--)
|
|
352 if (! isspace((int) sptr[pos])) break;
|
|
353 sptr[pos+1] = '\0'; /* ignore trailing whitespace */
|
|
354
|
|
355 switch (flag) {
|
|
356 case SQINFO_NAME:
|
|
357 if (*sptr != '-')
|
|
358 {
|
|
359 strncpy(sqinfo->name, sptr, SQINFO_NAMELEN-1);
|
|
360 sqinfo->name[SQINFO_NAMELEN-1] = '\0';
|
|
361 sqinfo->flags |= SQINFO_NAME;
|
|
362 }
|
|
363 break;
|
|
364
|
|
365 case SQINFO_ID:
|
|
366 if (*sptr != '-')
|
|
367 {
|
|
368 strncpy(sqinfo->id, sptr, SQINFO_NAMELEN-1);
|
|
369 sqinfo->id[SQINFO_NAMELEN-1] = '\0';
|
|
370 sqinfo->flags |= SQINFO_ID;
|
|
371 }
|
|
372 break;
|
|
373
|
|
374 case SQINFO_ACC:
|
|
375 if (*sptr != '-')
|
|
376 {
|
|
377 strncpy(sqinfo->acc, sptr, SQINFO_NAMELEN-1);
|
|
378 sqinfo->acc[SQINFO_NAMELEN-1] = '\0';
|
|
379 sqinfo->flags |= SQINFO_ACC;
|
|
380 }
|
|
381 break;
|
|
382
|
|
383 case SQINFO_DESC:
|
|
384 if (*sptr != '-')
|
|
385 {
|
|
386 if (sqinfo->flags & SQINFO_DESC) /* append? */
|
|
387 {
|
|
388 len = strlen(sqinfo->desc);
|
|
389 if (len < SQINFO_DESCLEN-2) /* is there room? */
|
|
390 {
|
|
391 strncat(sqinfo->desc, " ", SQINFO_DESCLEN-1-len); len++;
|
|
392 strncat(sqinfo->desc, sptr, SQINFO_DESCLEN-1-len);
|
|
393 }
|
|
394 }
|
|
395 else /* else copy */
|
|
396 strncpy(sqinfo->desc, sptr, SQINFO_DESCLEN-1);
|
|
397 sqinfo->desc[SQINFO_DESCLEN-1] = '\0';
|
|
398 sqinfo->flags |= SQINFO_DESC;
|
|
399 }
|
|
400 break;
|
|
401
|
|
402 case SQINFO_START:
|
|
403 if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
|
|
404 sqinfo->start = atoi(sptr);
|
|
405 if (sqinfo->start != 0) sqinfo->flags |= SQINFO_START;
|
|
406 break;
|
|
407
|
|
408 case SQINFO_STOP:
|
|
409 if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
|
|
410 sqinfo->stop = atoi(sptr);
|
|
411 if (sqinfo->stop != 0) sqinfo->flags |= SQINFO_STOP;
|
|
412 break;
|
|
413
|
|
414 case SQINFO_OLEN:
|
|
415 if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
|
|
416 sqinfo->olen = atoi(sptr);
|
|
417 if (sqinfo->olen != 0) sqinfo->flags |= SQINFO_OLEN;
|
|
418 break;
|
|
419
|
|
420 default:
|
|
421 Die("Invalid flag %d to SetSeqinfoString()", flag);
|
|
422 }
|
|
423 return 1;
|
|
424 }
|
|
425
|
|
426 void
|
|
427 SeqinfoCopy(SQINFO *sq1, SQINFO *sq2)
|
|
428 {
|
|
429 sq1->flags = sq2->flags;
|
|
430 if (sq2->flags & SQINFO_NAME) strcpy(sq1->name, sq2->name);
|
|
431 if (sq2->flags & SQINFO_ID) strcpy(sq1->id, sq2->id);
|
|
432 if (sq2->flags & SQINFO_ACC) strcpy(sq1->acc, sq2->acc);
|
|
433 if (sq2->flags & SQINFO_DESC) strcpy(sq1->desc, sq2->desc);
|
|
434 if (sq2->flags & SQINFO_LEN) sq1->len = sq2->len;
|
|
435 if (sq2->flags & SQINFO_START) sq1->start = sq2->start;
|
|
436 if (sq2->flags & SQINFO_STOP) sq1->stop = sq2->stop;
|
|
437 if (sq2->flags & SQINFO_OLEN) sq1->olen = sq2->olen;
|
|
438 if (sq2->flags & SQINFO_TYPE) sq1->type = sq2->type;
|
|
439 if (sq2->flags & SQINFO_SS) sq1->ss = Strdup(sq2->ss);
|
|
440 if (sq2->flags & SQINFO_SA) sq1->sa = Strdup(sq2->sa);
|
|
441 }
|
|
442
|
|
443 /* Function: ToDNA()
|
|
444 *
|
|
445 * Purpose: Convert a sequence to DNA.
|
|
446 * U --> T
|
|
447 */
|
|
448 void
|
|
449 ToDNA(char *seq)
|
|
450 {
|
|
451 for (; *seq != '\0'; seq++)
|
|
452 {
|
|
453 if (*seq == 'U') *seq = 'T';
|
|
454 else if (*seq == 'u') *seq = 't';
|
|
455 }
|
|
456 }
|
|
457
|
|
458 /* Function: ToRNA()
|
|
459 *
|
|
460 * Purpose: Convert a sequence to RNA.
|
|
461 * T --> U
|
|
462 */
|
|
463 void
|
|
464 ToRNA(char *seq)
|
|
465 {
|
|
466 for (; *seq != '\0'; seq++)
|
|
467 {
|
|
468 if (*seq == 'T') *seq = 'U';
|
|
469 else if (*seq == 't') *seq = 'u';
|
|
470 }
|
|
471 }
|
|
472
|
|
473
|
|
474 /* Function: ToIUPAC()
|
|
475 *
|
|
476 * Purpose: Convert X's, o's, other junk in a nucleic acid sequence to N's,
|
|
477 * to comply with IUPAC code. If is_aseq is TRUE, will allow gap
|
|
478 * characters though, so we can call ToIUPAC() on aligned seqs.
|
|
479 *
|
|
480 * NUCLEOTIDES is defined in squid.h as:
|
|
481 * "ACGTUNRYMKSWHBVDacgtunrymkswhbvd"
|
|
482 * gap chars allowed by isgap() are defined in squid.h as:
|
|
483 * " ._-~"
|
|
484 *
|
|
485 * WU-BLAST's pressdb will
|
|
486 * choke on X's, for instance, necessitating conversion
|
|
487 * of certain genome centers' data.
|
|
488 */
|
|
489 void
|
|
490 ToIUPAC(char *seq, int is_aseq)
|
|
491 {
|
|
492 if (is_aseq) {
|
|
493 for (; *seq != '\0'; seq++)
|
|
494 if (strchr(NUCLEOTIDES, *seq) == NULL && ! isgap(*seq)) *seq = 'N';
|
|
495 } else {
|
|
496 for (; *seq != '\0'; seq++)
|
|
497 if (strchr(NUCLEOTIDES, *seq) == NULL) *seq = 'N';
|
|
498 }
|
|
499 }
|
|
500
|
|
501
|
|
502 /* Function: addseq()
|
|
503 *
|
|
504 * Purpose: Add a line of sequence to the growing string in V.
|
|
505 *
|
|
506 * In the seven supported unaligned formats, all sequence
|
|
507 * lines may contain whitespace that must be filtered out;
|
|
508 * four formats (PIR, EMBL, Genbank, GCG) include coordinates
|
|
509 * that must be filtered out. Thus an (!isdigit && !isspace)
|
|
510 * test on each character before we accept it.
|
|
511 */
|
|
512 static void
|
|
513 addseq(char *s, struct ReadSeqVars *V)
|
|
514 {
|
|
515 char *s0;
|
|
516 char *sq;
|
|
517 int rpl; /* valid residues per line */
|
|
518 int bpl; /* characters per line */
|
|
519
|
|
520 if (V->ssimode == -1)
|
|
521 { /* Normal mode: keeping the seq */
|
|
522 /* Make sure we have enough room. We know that s is <= buflen,
|
|
523 * so just make sure we've got room for a whole new buflen worth
|
|
524 * of sequence.
|
|
525 */
|
|
526 if (V->seqlen + V->buflen > V->maxseq) {
|
|
527 V->maxseq += MAX(V->buflen, kStartLength);
|
|
528 V->seq = ReallocOrDie (V->seq, V->maxseq+1);
|
|
529 }
|
|
530
|
|
531 sq = V->seq + V->seqlen;
|
|
532 while (*s != 0) {
|
|
533 #ifdef CLUSTALO
|
|
534 if (! isdigit((int) *s) && ! isspace((int) *s) && isprint((int) *s)) {
|
|
535 #else
|
|
536 if (! isdigit((int) *s) && ! isspace((int) *s)) {
|
|
537 #endif
|
|
538 *sq = *s;
|
|
539 sq++;
|
|
540 }
|
|
541 s++;
|
|
542 }
|
|
543 V->seqlen = sq - V->seq;
|
|
544 }
|
|
545 else /* else: indexing mode, discard the seq */
|
|
546 {
|
|
547 s0 = s;
|
|
548 rpl = 0;
|
|
549 while (*s != 0) {
|
|
550 if (! isdigit((int) *s) && ! isspace((int) *s)) {
|
|
551 rpl++;
|
|
552 }
|
|
553 s++;
|
|
554 }
|
|
555 V->seqlen += rpl;
|
|
556 bpl = s - s0;
|
|
557
|
|
558 /* Keep track of the global rpl, bpl for the file.
|
|
559 * This is overly complicated because we have to
|
|
560 * allow the last line of each record (e.g. the last addseq() call
|
|
561 * on each sequence) to have a different length - and sometimes
|
|
562 * we'll have one-line sequence records, too. Thus we only
|
|
563 * do something with the global V->rpl when we have *passed over*
|
|
564 * a line - we keep the last line's rpl in last_rpl. And because
|
|
565 * a file might consist entirely of single-line records, we keep
|
|
566 * a third guy, maxrpl, that tells us the maximum rpl of any line
|
|
567 * in the file. If we reach the end of file and rpl is still unset,
|
|
568 * we'll set it to maxrpl. If we reach eof and rpl is set, but is
|
|
569 * less than maxrpl, that's a weird case where a last line in some
|
|
570 * record is longer than every other line.
|
|
571 */
|
|
572 if (V->rpl != 0) { /* 0 means we already know rpl is invalid */
|
|
573 if (V->lastrpl > 0) { /* we're on something that's not the first line */
|
|
574 if (V->rpl > 0 && V->lastrpl != V->rpl) V->rpl = 0;
|
|
575 else if (V->rpl == -1) V->rpl = V->lastrpl;
|
|
576 }
|
|
577 V->lastrpl = rpl;
|
|
578 if (rpl > V->maxrpl) V->maxrpl = rpl; /* make sure we check max length of final lines */
|
|
579 }
|
|
580 if (V->bpl != 0) { /* 0 means we already know bpl is invalid */
|
|
581 if (V->lastbpl > 0) { /* we're on something that's not the first line */
|
|
582 if (V->bpl > 0 && V->lastbpl != V->bpl) V->bpl = 0;
|
|
583 else if (V->bpl == -1) V->bpl = V->lastbpl;
|
|
584 }
|
|
585 V->lastbpl = bpl;
|
|
586 if (bpl > V->maxbpl) V->maxbpl = bpl; /* make sure we check max length of final lines */
|
|
587 }
|
|
588 } /* end of indexing mode of addseq(). */
|
|
589
|
|
590 }
|
|
591
|
|
592 static void
|
|
593 readLoop(int addfirst, int (*endTest)(char *,int *), struct ReadSeqVars *V)
|
|
594 {
|
|
595 int addend = 0;
|
|
596 int done = 0;
|
|
597
|
|
598 V->seqlen = 0;
|
|
599 V->lastrpl = V->lastbpl = 0;
|
|
600 if (addfirst) {
|
|
601 if (V->ssimode >= 0) V->d_off = V->ssioffset;
|
|
602 addseq(V->buf, V);
|
|
603 } else if (V->ssimode >= 0)
|
|
604 if (0 != SSIGetFilePosition(V->f, V->ssimode, &(V->d_off)))
|
|
605 Die("SSIGetFilePosition() failed");
|
|
606
|
|
607 do {
|
|
608 SeqfileGetLine(V);
|
|
609 /* feof() alone is a bug; files not necessarily \n terminated */
|
|
610 if (*(V->buf) == '\0' && feof(V->f))
|
|
611 done = TRUE;
|
|
612 done |= (*endTest)(V->buf, &addend);
|
|
613 if (addend || !done)
|
|
614 addseq(V->buf, V);
|
|
615 } while (!done);
|
|
616 }
|
|
617
|
|
618
|
|
619 static int
|
|
620 endPIR(char *s, int *addend)
|
|
621 {
|
|
622 *addend = 0;
|
|
623 if ((strncmp(s, "///", 3) == 0) ||
|
|
624 (strncmp(s, "ENTRY", 5) == 0))
|
|
625 return 1;
|
|
626 else
|
|
627 return 0;
|
|
628 }
|
|
629
|
|
630 static void
|
|
631 readPIR(struct ReadSeqVars *V)
|
|
632 {
|
|
633 char *sptr;
|
|
634 /* load first line of entry */
|
|
635 while (!feof(V->f) && strncmp(V->buf, "ENTRY", 5) != 0) {
|
|
636 SeqfileGetLine(V);
|
|
637 }
|
|
638 if (feof(V->f)) return;
|
|
639 if (V->ssimode >= 0) V->r_off = V->ssioffset;
|
|
640
|
|
641 if ((sptr = strtok(V->buf + 15, "\n\t ")) != NULL)
|
|
642 {
|
|
643 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
644 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
|
|
645 }
|
|
646 do {
|
|
647 SeqfileGetLine(V);
|
|
648 if (!feof(V->f) && strncmp(V->buf, "TITLE", 5) == 0)
|
|
649 SetSeqinfoString(V->sqinfo, V->buf+15, SQINFO_DESC);
|
|
650 else if (!feof(V->f) && strncmp(V->buf, "ACCESSION", 9) == 0)
|
|
651 {
|
|
652 if ((sptr = strtok(V->buf+15, " \t\n")) != NULL)
|
|
653 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
|
|
654 }
|
|
655 } while (! feof(V->f) && (strncmp(V->buf,"SEQUENCE", 8) != 0));
|
|
656 SeqfileGetLine(V); /* skip next line, coords */
|
|
657
|
|
658 readLoop(0, endPIR, V);
|
|
659
|
|
660 /* reading a real PIR-CODATA database file, we keep the source coords
|
|
661 */
|
|
662 V->sqinfo->start = 1;
|
|
663 V->sqinfo->stop = V->seqlen;
|
|
664 V->sqinfo->olen = V->seqlen;
|
|
665 V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;
|
|
666
|
|
667 /* get next line
|
|
668 */
|
|
669 while (!feof(V->f) && strncmp(V->buf, "ENTRY", 5) != 0) {
|
|
670 SeqfileGetLine(V);
|
|
671 }
|
|
672 }
|
|
673
|
|
674
|
|
675
|
|
676 static int
|
|
677 endIG(char *s, int *addend)
|
|
678 {
|
|
679 *addend = 1; /* 1 or 2 occur in line w/ bases */
|
|
680 return((strchr(s,'1')!=NULL) || (strchr(s,'2')!=NULL));
|
|
681 }
|
|
682
|
|
683 static void
|
|
684 readIG(struct ReadSeqVars *V)
|
|
685 {
|
|
686 char *nm;
|
|
687 /* position past ';' comments */
|
|
688 do {
|
|
689 SeqfileGetLine(V);
|
|
690 } while (! (feof(V->f) || ((*V->buf != 0) && (*V->buf != ';')) ));
|
|
691
|
|
692 if (!feof(V->f))
|
|
693 {
|
|
694 if ((nm = strtok(V->buf, "\n\t ")) != NULL)
|
|
695 SetSeqinfoString(V->sqinfo, nm, SQINFO_NAME);
|
|
696
|
|
697 readLoop(0, endIG, V);
|
|
698 }
|
|
699
|
|
700 while (!(feof(V->f) || ((*V->buf != '\0') && (*V->buf == ';'))))
|
|
701 SeqfileGetLine(V);
|
|
702 }
|
|
703
|
|
704 static int
|
|
705 endStrider(char *s, int *addend)
|
|
706 {
|
|
707 *addend = 0;
|
|
708 return (strstr( s, "//") != NULL);
|
|
709 }
|
|
710
|
|
711 static void
|
|
712 readStrider(struct ReadSeqVars *V)
|
|
713 {
|
|
714 char *nm;
|
|
715
|
|
716 while ((!feof(V->f)) && (*V->buf == ';'))
|
|
717 {
|
|
718 if (strncmp(V->buf,"; DNA sequence", 14) == 0)
|
|
719 {
|
|
720 if ((nm = strtok(V->buf+16, ",\n\t ")) != NULL)
|
|
721 SetSeqinfoString(V->sqinfo, nm, SQINFO_NAME);
|
|
722 }
|
|
723 SeqfileGetLine(V);
|
|
724 }
|
|
725
|
|
726 if (! feof(V->f))
|
|
727 readLoop(1, endStrider, V);
|
|
728
|
|
729 /* load next line
|
|
730 */
|
|
731 while ((!feof(V->f)) && (*V->buf != ';'))
|
|
732 SeqfileGetLine(V);
|
|
733 }
|
|
734
|
|
735
|
|
736 static int
|
|
737 endGB(char *s, int *addend)
|
|
738 {
|
|
739 *addend = 0;
|
|
740 return ((strstr(s,"//") != NULL) || (strstr(s,"LOCUS") == s));
|
|
741 }
|
|
742
|
|
743 static void
|
|
744 readGenBank(struct ReadSeqVars *V)
|
|
745 {
|
|
746 char *sptr;
|
|
747 int in_definition;
|
|
748
|
|
749 /* We'll map three genbank identifiers onto names:
|
|
750 * LOCUS -> sqinfo.name
|
|
751 * ACCESSION -> sqinfo.acc [primary accession only]
|
|
752 * VERSION -> sqinfo.id
|
|
753 * We don't currently store the GI number, or secondary accessions.
|
|
754 */
|
|
755 while (strncmp(V->buf, "LOCUS", 5) != 0) {
|
|
756 SeqfileGetLine(V);
|
|
757 }
|
|
758 if (V->ssimode >= 0) V->r_off = V->ssioffset;
|
|
759
|
|
760 if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
|
|
761 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
762
|
|
763 in_definition = FALSE;
|
|
764 while (! feof(V->f))
|
|
765 {
|
|
766 SeqfileGetLine(V);
|
|
767 if (! feof(V->f) && strstr(V->buf, "DEFINITION") == V->buf)
|
|
768 {
|
|
769 if ((sptr = strtok(V->buf+12, "\n")) != NULL)
|
|
770 SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
|
|
771 in_definition = TRUE;
|
|
772 }
|
|
773 else if (! feof(V->f) && strstr(V->buf, "ACCESSION") == V->buf)
|
|
774 {
|
|
775 if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
|
|
776 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
|
|
777 in_definition = FALSE;
|
|
778 }
|
|
779 else if (! feof(V->f) && strstr(V->buf, "VERSION") == V->buf)
|
|
780 {
|
|
781 if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
|
|
782 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
|
|
783 in_definition = FALSE;
|
|
784 }
|
|
785 else if (strncmp(V->buf,"ORIGIN", 6) != 0)
|
|
786 {
|
|
787 if (in_definition)
|
|
788 SetSeqinfoString(V->sqinfo, V->buf, SQINFO_DESC);
|
|
789 }
|
|
790 else
|
|
791 break;
|
|
792 }
|
|
793
|
|
794 readLoop(0, endGB, V);
|
|
795
|
|
796 /* reading a real GenBank database file, we keep the source coords
|
|
797 */
|
|
798 V->sqinfo->start = 1;
|
|
799 V->sqinfo->stop = V->seqlen;
|
|
800 V->sqinfo->olen = V->seqlen;
|
|
801 V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;
|
|
802
|
|
803
|
|
804 while (!(feof(V->f) || ((*V->buf!=0) && (strstr(V->buf,"LOCUS") == V->buf))))
|
|
805 SeqfileGetLine(V);
|
|
806 /* SRE: V->s now holds "//", so sequential
|
|
807 reads are wedged: fixed Tue Jul 13 1993 */
|
|
808 while (!feof(V->f) && strstr(V->buf, "LOCUS ") != V->buf)
|
|
809 SeqfileGetLine(V);
|
|
810 }
|
|
811
|
|
812 static int
|
|
813 endGCGdata(char *s, int *addend)
|
|
814 {
|
|
815 *addend = 0;
|
|
816 return (*s == '>');
|
|
817 }
|
|
818
|
|
819 static void
|
|
820 readGCGdata(struct ReadSeqVars *V)
|
|
821 {
|
|
822 int binary = FALSE; /* whether data are binary or not */
|
|
823 int blen = 0; /* length of binary sequence */
|
|
824
|
|
825 /* first line contains ">>>>" followed by name */
|
|
826 if (Strparse(">>>>([^ ]+) .+2BIT +Len: ([0-9]+)", V->buf, 2))
|
|
827 {
|
|
828 binary = TRUE;
|
|
829 SetSeqinfoString(V->sqinfo, sqd_parse[1], SQINFO_NAME);
|
|
830 blen = atoi(sqd_parse[2]);
|
|
831 }
|
|
832 else if (Strparse(">>>>([^ ]+) .+ASCII +Len: [0-9]+", V->buf, 1))
|
|
833 SetSeqinfoString(V->sqinfo, sqd_parse[1], SQINFO_NAME);
|
|
834 else
|
|
835 Die("bogus GCGdata format? %s", V->buf);
|
|
836
|
|
837 /* second line contains free text description */
|
|
838 SeqfileGetLine(V);
|
|
839 SetSeqinfoString(V->sqinfo, V->buf, SQINFO_DESC);
|
|
840
|
|
841 if (binary) {
|
|
842 /* allocate for blen characters +3... (allow for 3 bytes of slop) */
|
|
843 if (blen >= V->maxseq) {
|
|
844 V->maxseq = blen;
|
|
845 if ((V->seq = (char *) realloc (V->seq, sizeof(char)*(V->maxseq+4)))==NULL)
|
|
846 Die("malloc failed");
|
|
847 }
|
|
848 /* read (blen+3)/4 bytes from file */
|
|
849 if (fread(V->seq, sizeof(char), (blen+3)/4, V->f) < (size_t) ((blen+3)/4))
|
|
850 Die("fread failed");
|
|
851 V->seqlen = blen;
|
|
852 /* convert binary code to seq */
|
|
853 GCGBinaryToSequence(V->seq, blen);
|
|
854 }
|
|
855 else readLoop(0, endGCGdata, V);
|
|
856
|
|
857 while (!(feof(V->f) || ((*V->buf != 0) && (*V->buf == '>'))))
|
|
858 SeqfileGetLine(V);
|
|
859 }
|
|
860
|
|
861 static int
|
|
862 endPearson(char *s, int *addend)
|
|
863 {
|
|
864 *addend = 0;
|
|
865 return(*s == '>');
|
|
866 }
|
|
867
|
|
868 static void
|
|
869 readPearson(struct ReadSeqVars *V)
|
|
870 {
|
|
871 char *sptr;
|
|
872
|
|
873 if (V->ssimode >= 0) V->r_off = V->ssioffset;
|
|
874
|
|
875 if (*V->buf != '>')
|
|
876 Die("\
|
|
877 File %s does not appear to be in FASTA format at line %d.\n\
|
|
878 You may want to specify the file format on the command line.\n\
|
|
879 Usually this is done with an option --informat <fmt>.\n",
|
|
880 V->fname, V->linenumber);
|
|
881
|
|
882 if ((sptr = strtok(V->buf+1, "\n\t ")) != NULL)
|
|
883 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
884 if ((sptr = strtok(NULL, "\n")) != NULL)
|
|
885 SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
|
|
886
|
|
887 readLoop(0, endPearson, V);
|
|
888
|
|
889 while (!(feof(V->f) || ((*V->buf != 0) && (*V->buf == '>')))) {
|
|
890 SeqfileGetLine(V);
|
|
891 }
|
|
892 }
|
|
893
|
|
894
|
|
895 static int
|
|
896 endEMBL(char *s, int *addend)
|
|
897 {
|
|
898 *addend = 0;
|
|
899 /* Some people (Berlin 5S rRNA database, f'r instance) use
|
|
900 * an extended EMBL format that attaches extra data after
|
|
901 * the sequence -- watch out for that. We use the fact that
|
|
902 * real EMBL sequence lines begin with five spaces.
|
|
903 *
|
|
904 * We can use this as the sole end test because readEMBL() will
|
|
905 * advance to the next ID line before starting to read again.
|
|
906 */
|
|
907 return (strncmp(s," ",5) != 0);
|
|
908 /* return ((strstr(s,"//") != NULL) || (strstr(s,"ID ") == s)); */
|
|
909 }
|
|
910
|
|
911 static void
|
|
912 readEMBL(struct ReadSeqVars *V)
|
|
913 {
|
|
914 char *sptr;
|
|
915
|
|
916 /* make sure we have first line */
|
|
917 while (!feof(V->f) && strncmp(V->buf, "ID ", 4) != 0) {
|
|
918 SeqfileGetLine(V);
|
|
919 }
|
|
920 if (V->ssimode >= 0) V->r_off = V->ssioffset;
|
|
921
|
|
922 if ((sptr = strtok(V->buf+5, "\n\t ")) != NULL)
|
|
923 {
|
|
924 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
925 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
|
|
926 }
|
|
927
|
|
928 do {
|
|
929 SeqfileGetLine(V);
|
|
930 if (!feof(V->f) && strstr(V->buf, "AC ") == V->buf)
|
|
931 {
|
|
932 if ((sptr = strtok(V->buf+5, "; \t\n")) != NULL)
|
|
933 SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
|
|
934 }
|
|
935 else if (!feof(V->f) && strstr(V->buf, "DE ") == V->buf)
|
|
936 {
|
|
937 if ((sptr = strtok(V->buf+5, "\n")) != NULL)
|
|
938 SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
|
|
939 }
|
|
940 } while (! feof(V->f) && strncmp(V->buf,"SQ",2) != 0);
|
|
941
|
|
942 readLoop(0, endEMBL, V);
|
|
943
|
|
944 /* Hack for Staden experiment files: convert - to N
|
|
945 */
|
|
946 if (V->ssimode == -1) /* if we're in ssi mode, we're not keeping the seq */
|
|
947 for (sptr = V->seq; *sptr != '\0'; sptr++)
|
|
948 if (*sptr == '-') *sptr = 'N';
|
|
949
|
|
950 /* reading a real EMBL database file, we keep the source coords
|
|
951 */
|
|
952 V->sqinfo->start = 1;
|
|
953 V->sqinfo->stop = V->seqlen;
|
|
954 V->sqinfo->olen = V->seqlen;
|
|
955 V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;
|
|
956
|
|
957 /* load next record's ID line */
|
|
958 while (!feof(V->f) && strncmp(V->buf, "ID ", 4) != 0) {
|
|
959 SeqfileGetLine(V);
|
|
960 }
|
|
961
|
|
962 }
|
|
963
|
|
964
|
|
965 static int
|
|
966 endZuker(char *s, int *addend)
|
|
967 {
|
|
968 *addend = 0;
|
|
969 return( *s == '(' );
|
|
970 }
|
|
971
|
|
972 static void
|
|
973 readZuker(struct ReadSeqVars *V)
|
|
974 {
|
|
975 char *sptr;
|
|
976
|
|
977 SeqfileGetLine(V); /*s == "seqLen seqid string..."*/
|
|
978
|
|
979 if ((sptr = strtok(V->buf+6, " \t\n")) != NULL)
|
|
980 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
981
|
|
982 if ((sptr = strtok(NULL, "\n")) != NULL)
|
|
983 SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
|
|
984
|
|
985 readLoop(0, endZuker, V);
|
|
986
|
|
987 while (!(feof(V->f) | ((*V->buf != '\0') & (*V->buf == '('))))
|
|
988 SeqfileGetLine(V);
|
|
989 }
|
|
990
|
|
991 static void
|
|
992 readUWGCG(struct ReadSeqVars *V)
|
|
993 {
|
|
994 char *si;
|
|
995 char *sptr;
|
|
996 int done;
|
|
997
|
|
998 V->seqlen = 0;
|
|
999
|
|
1000 /*writeseq: " %s Length: %d (today) Check: %d ..\n" */
|
|
1001 /*drop above or ".." from id*/
|
|
1002 if ((si = strstr(V->buf," Length: ")) != NULL) *si = 0;
|
|
1003 else if ((si = strstr(V->buf,"..")) != NULL) *si = 0;
|
|
1004
|
|
1005 if ((sptr = strtok(V->buf, "\n\t ")) != NULL)
|
|
1006 SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
|
|
1007
|
|
1008 do {
|
|
1009 done = feof(V->f);
|
|
1010 SeqfileGetLine(V);
|
|
1011 if (! done) addseq(V->buf, V);
|
|
1012 } while (!done);
|
|
1013 }
|
|
1014
|
|
1015
|
|
1016 /* Function: ReadSeq()
|
|
1017 *
|
|
1018 * Purpose: Read next sequence from an open database file.
|
|
1019 * Return the sequence and associated info.
|
|
1020 *
|
|
1021 * Args: fp - open sequence database file pointer
|
|
1022 * format - format of the file (previously determined
|
|
1023 * by call to SeqfileFormat()).
|
|
1024 * Currently unused, since we carry it in V.
|
|
1025 * ret_seq - RETURN: sequence
|
|
1026 * sqinfo - RETURN: filled in w/ other information
|
|
1027 *
|
|
1028 * Limitations: uses squid_errno, so it's not threadsafe.
|
|
1029 *
|
|
1030 * Return: 1 on success, 0 on failure.
|
|
1031 * ret_seq and some field of sqinfo are allocated here,
|
|
1032 * The preferred call mechanism to properly free the memory is:
|
|
1033 *
|
|
1034 * SQINFO sqinfo;
|
|
1035 * char *seq;
|
|
1036 *
|
|
1037 * ReadSeq(fp, format, &seq, &sqinfo);
|
|
1038 * ... do something...
|
|
1039 * FreeSequence(seq, &sqinfo);
|
|
1040 */
|
|
1041 int
|
|
1042 ReadSeq(SQFILE *V, int format, char **ret_seq, SQINFO *sqinfo)
|
|
1043 {
|
|
1044 int gotuw;
|
|
1045
|
|
1046 squid_errno = SQERR_OK;
|
|
1047
|
|
1048 /* Here's the hack for sequential access of sequences from
|
|
1049 * the multiple sequence alignment formats
|
|
1050 */
|
|
1051 if (IsAlignmentFormat(V->format))
|
|
1052 {
|
|
1053 if (V->msa->lastidx >= V->msa->nseq)
|
|
1054 { /* out of data. try to read another alignment */
|
|
1055 MSAFree(V->msa);
|
|
1056 if ((V->msa = MSAFileRead(V->afp)) == NULL)
|
|
1057 return 0;
|
|
1058 V->msa->lastidx = 0;
|
|
1059 }
|
|
1060 /* copy and dealign the appropriate aligned seq */
|
|
1061 /* AW: stopping squid from dealigning sequences and corresponding info */
|
|
1062 #ifdef CLUSTALO
|
|
1063 V->seq = sre_strdup(V->msa->aseq[V->msa->lastidx], V->msa->alen);
|
|
1064 #else
|
|
1065 MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen,
|
|
1066 V->msa->aseq[V->msa->lastidx], &(V->seq));
|
|
1067 #endif
|
|
1068 V->seqlen = strlen(V->seq);
|
|
1069
|
|
1070 /* Extract sqinfo stuff for this sequence from the msa.
|
|
1071 * Tedious; code that should be cleaned.
|
|
1072 */
|
|
1073 sqinfo->flags = 0;
|
|
1074 if (V->msa->sqname[V->msa->lastidx] != NULL)
|
|
1075 SetSeqinfoString(sqinfo, V->msa->sqname[V->msa->lastidx], SQINFO_NAME);
|
|
1076 if (V->msa->sqacc != NULL && V->msa->sqacc[V->msa->lastidx] != NULL)
|
|
1077 SetSeqinfoString(sqinfo, V->msa->sqacc[V->msa->lastidx], SQINFO_ACC);
|
|
1078 if (V->msa->sqdesc != NULL && V->msa->sqdesc[V->msa->lastidx] != NULL)
|
|
1079 SetSeqinfoString(sqinfo, V->msa->sqdesc[V->msa->lastidx], SQINFO_DESC);
|
|
1080 if (V->msa->ss != NULL && V->msa->ss[V->msa->lastidx] != NULL) {
|
|
1081 /* AW: stopping squid from dealigning sequences and corresponding info */
|
|
1082 #ifdef CLUSTALO
|
|
1083 sqinfo->ss = sre_strdup(V->msa->ss[V->msa->lastidx], V->msa->alen);
|
|
1084 #else
|
|
1085 MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen,
|
|
1086 V->msa->ss[V->msa->lastidx], &(sqinfo->ss));
|
|
1087 #endif
|
|
1088 sqinfo->flags |= SQINFO_SS;
|
|
1089 }
|
|
1090 if (V->msa->sa != NULL && V->msa->sa[V->msa->lastidx] != NULL) {
|
|
1091 /* AW: stopping squid from dealigning sequences and corresponding info */
|
|
1092 #ifdef CLUSTALO
|
|
1093 sqinfo->sa = sre_strdup(V->msa->sa[V->msa->lastidx], V->msa->alen);
|
|
1094 #else
|
|
1095 MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen,
|
|
1096 V->msa->sa[V->msa->lastidx], &(sqinfo->sa));
|
|
1097 #endif
|
|
1098 sqinfo->flags |= SQINFO_SA;
|
|
1099 }
|
|
1100 V->msa->lastidx++;
|
|
1101 }
|
|
1102 else {
|
|
1103 if (feof(V->f)) return 0;
|
|
1104
|
|
1105 if (V->ssimode == -1) { /* normal mode */
|
|
1106 V->seq = (char*) calloc (kStartLength+1, sizeof(char));
|
|
1107 V->maxseq = kStartLength;
|
|
1108 } else { /* index mode: discarding seq */
|
|
1109 V->seq = NULL;
|
|
1110 V->maxseq = 0;
|
|
1111 }
|
|
1112 V->seqlen = 0;
|
|
1113 V->sqinfo = sqinfo;
|
|
1114 V->sqinfo->flags = 0;
|
|
1115
|
|
1116 switch (V->format) {
|
|
1117 case SQFILE_IG : readIG(V); break;
|
|
1118 case SQFILE_STRIDER : readStrider(V); break;
|
|
1119 case SQFILE_GENBANK : readGenBank(V); break;
|
|
1120 case SQFILE_FASTA : readPearson(V); break;
|
|
1121 case SQFILE_EMBL : readEMBL(V); break;
|
|
1122 case SQFILE_ZUKER : readZuker(V); break;
|
|
1123 case SQFILE_PIR : readPIR(V); break;
|
|
1124 case SQFILE_GCGDATA : readGCGdata(V); break;
|
|
1125
|
|
1126 case SQFILE_GCG :
|
|
1127 do { /* skip leading comments on GCG file */
|
|
1128 gotuw = (strstr(V->buf,"..") != NULL);
|
|
1129 if (gotuw) readUWGCG(V);
|
|
1130 SeqfileGetLine(V);
|
|
1131 } while (! feof(V->f));
|
|
1132 break;
|
|
1133
|
|
1134 case SQFILE_IDRAW: /* SRE: no attempt to read idraw postscript */
|
|
1135 default:
|
|
1136 squid_errno = SQERR_FORMAT;
|
|
1137 free(V->seq);
|
|
1138 return 0;
|
|
1139 }
|
|
1140 if (V->seq != NULL) /* (it can be NULL in indexing mode) */
|
|
1141 V->seq[V->seqlen] = 0; /* stick a string terminator on it */
|
|
1142 }
|
|
1143
|
|
1144 /* Cleanup
|
|
1145 */
|
|
1146 sqinfo->len = V->seqlen;
|
|
1147 sqinfo->flags |= SQINFO_LEN;
|
|
1148 *ret_seq = V->seq;
|
|
1149 if (squid_errno == SQERR_OK) return 1; else return 0;
|
|
1150 }
|
|
1151
|
|
1152 /* Function: SeqfileFormat()
|
|
1153 * Date: SRE, Tue Jun 22 10:58:58 1999 [Sanger Centre]
|
|
1154 *
|
|
1155 * Purpose: Determine format of an open file.
|
|
1156 * Returns format code.
|
|
1157 * Rewinds the file.
|
|
1158 *
|
|
1159 * Autodetects the following unaligned formats:
|
|
1160 * SQFILE_FASTA
|
|
1161 * SQFILE_GENBANK
|
|
1162 * SQFILE_EMBL
|
|
1163 * SQFILE_GCG
|
|
1164 * SQFILE_GCGDATA
|
|
1165 * SQFILE_PIR
|
|
1166 * Also autodetects the following alignment formats:
|
|
1167 * MSAFILE_STOCKHOLM
|
|
1168 * MSAFILE_MSF
|
|
1169 * MSAFILE_CLUSTAL
|
|
1170 * MSAFILE_SELEX
|
|
1171 * MSAFILE_PHYLIP
|
|
1172 *
|
|
1173 * Can't autodetect MSAFILE_A2M, calls it SQFILE_FASTA.
|
|
1174 * MSAFileFormat() does the opposite.
|
|
1175 *
|
|
1176 * Args: sfp - open SQFILE
|
|
1177 *
|
|
1178 * Return: format code, or SQFILE_UNKNOWN if unrecognized
|
|
1179 */
|
|
1180 int
|
|
1181 SeqfileFormat(FILE *fp)
|
|
1182 {
|
|
1183 char *buf;
|
|
1184 int len;
|
|
1185 int fmt = SQFILE_UNKNOWN;
|
|
1186 int ndataline;
|
|
1187 char *bufcpy, *s, *s1, *s2;
|
|
1188 int has_junk;
|
|
1189
|
|
1190 buf = NULL;
|
|
1191 len = 0;
|
|
1192 ndataline = 0;
|
|
1193 has_junk = FALSE;
|
|
1194 while (sre_fgets(&buf, &len, fp) != NULL)
|
|
1195 {
|
|
1196 if (IsBlankline(buf)) continue;
|
|
1197
|
|
1198 /* Well-behaved formats identify themselves in first nonblank line.
|
|
1199 */
|
|
1200 if (ndataline == 0)
|
|
1201 {
|
|
1202 if (strncmp(buf, ">>>>", 4) == 0 && strstr(buf, "Len: "))
|
|
1203 { fmt = SQFILE_GCGDATA; goto DONE; }
|
|
1204
|
|
1205 if (buf[0] == '>')
|
|
1206 { fmt = SQFILE_FASTA; goto DONE; }
|
|
1207
|
|
1208 if (strncmp(buf, "!!AA_SEQUENCE", 13) == 0 ||
|
|
1209 strncmp(buf, "!!NA_SEQUENCE", 13) == 0)
|
|
1210 { fmt = SQFILE_GCG; goto DONE; }
|
|
1211
|
|
1212 if (strncmp(buf, "# STOCKHOLM 1.", 14) == 0)
|
|
1213 { fmt = MSAFILE_STOCKHOLM; goto DONE; }
|
|
1214
|
|
1215 if (strncmp(buf, "CLUSTAL", 7) == 0 &&
|
|
1216 strstr(buf, "multiple sequence alignment") != NULL)
|
|
1217 { fmt = MSAFILE_CLUSTAL; goto DONE; }
|
|
1218
|
|
1219 if (strncmp(buf, "!!AA_MULTIPLE_ALIGNMENT", 23) == 0 ||
|
|
1220 strncmp(buf, "!!NA_MULTIPLE_ALIGNMENT", 23) == 0)
|
|
1221 { fmt = MSAFILE_MSF; goto DONE; }
|
|
1222
|
|
1223 /* PHYLIP id: also just a good bet */
|
|
1224 bufcpy = sre_strdup(buf, -1);
|
|
1225 s = bufcpy;
|
|
1226 if ((s1 = sre_strtok(&s, WHITESPACE, NULL)) != NULL &&
|
|
1227 (s2 = sre_strtok(&s, WHITESPACE, NULL)) != NULL &&
|
|
1228 IsInt(s1) &&
|
|
1229 IsInt(s2))
|
|
1230 { free(bufcpy); fmt = MSAFILE_PHYLIP; goto DONE; }
|
|
1231 free(bufcpy);
|
|
1232 }
|
|
1233
|
|
1234 /* We trust that other formats identify themselves soon.
|
|
1235 */
|
|
1236 /* dead giveaways for extended SELEX */
|
|
1237 if (strncmp(buf, "#=AU", 4) == 0 ||
|
|
1238 strncmp(buf, "#=ID", 4) == 0 ||
|
|
1239 strncmp(buf, "#=AC", 4) == 0 ||
|
|
1240 strncmp(buf, "#=DE", 4) == 0 ||
|
|
1241 strncmp(buf, "#=GA", 4) == 0 ||
|
|
1242 strncmp(buf, "#=TC", 4) == 0 ||
|
|
1243 strncmp(buf, "#=NC", 4) == 0 ||
|
|
1244 strncmp(buf, "#=SQ", 4) == 0 ||
|
|
1245 strncmp(buf, "#=SS", 4) == 0 ||
|
|
1246 strncmp(buf, "#=CS", 4) == 0 ||
|
|
1247 strncmp(buf, "#=RF", 4) == 0)
|
|
1248 { fmt = MSAFILE_SELEX; goto DONE; }
|
|
1249
|
|
1250 if (strncmp(buf, "///", 3) == 0 || strncmp(buf, "ENTRY ", 6) == 0)
|
|
1251 { fmt = SQFILE_PIR; goto DONE; }
|
|
1252
|
|
1253 /* a ha, diagnostic of an (old) MSF file */
|
|
1254 if ((strstr(buf, "..") != NULL) &&
|
|
1255 (strstr(buf, "MSF:") != NULL) &&
|
|
1256 (strstr(buf, "Check:")!= NULL))
|
|
1257 { fmt = MSAFILE_MSF; goto DONE; }
|
|
1258
|
|
1259 /* unaligned GCG (must follow MSF test!) */
|
|
1260 if (strstr(buf, " Check: ") != NULL && strstr(buf, "..") != NULL)
|
|
1261 { fmt = SQFILE_GCG; goto DONE; }
|
|
1262
|
|
1263 if (strncmp(buf,"LOCUS ",6) == 0 || strncmp(buf,"ORIGIN ",6) == 0)
|
|
1264 { fmt = SQFILE_GENBANK; goto DONE; }
|
|
1265
|
|
1266 if (strncmp(buf,"ID ",5) == 0 || strncmp(buf,"SQ ",5) == 0)
|
|
1267 { fmt = SQFILE_EMBL; goto DONE; }
|
|
1268
|
|
1269 /* But past here, we're being desperate. A simple SELEX file is
|
|
1270 * very difficult to detect; we can only try to disprove it.
|
|
1271 */
|
|
1272 s = buf;
|
|
1273 if ((s1 = sre_strtok(&s, WHITESPACE, NULL)) == NULL) continue; /* skip blank lines */
|
|
1274 if (strchr("#%", *s1) != NULL) continue; /* skip comment lines */
|
|
1275
|
|
1276 /* Disproof 1. Noncomment, nonblank lines in a SELEX file
|
|
1277 * must have at least two space-delimited fields (name/seq)
|
|
1278 */
|
|
1279 if ((s2 = sre_strtok(&s, WHITESPACE, NULL)) == NULL)
|
|
1280 has_junk = TRUE;
|
|
1281
|
|
1282 /* Disproof 2.
|
|
1283 * The sequence field should look like a sequence.
|
|
1284 */
|
|
1285 if (s2 != NULL && Seqtype(s2) == kOtherSeq)
|
|
1286 has_junk = TRUE;
|
|
1287
|
|
1288 ndataline++;
|
|
1289 if (ndataline == 300) break; /* only look at first 300 lines */
|
|
1290 }
|
|
1291
|
|
1292 if (ndataline == 0)
|
|
1293 Die("Sequence file contains no data");
|
|
1294
|
|
1295 /* If we've made it this far, we've run out of data, but there
|
|
1296 * was at least one line of it; check if we've
|
|
1297 * disproven SELEX. If not, cross our fingers, pray, and guess SELEX.
|
|
1298 */
|
|
1299 if (has_junk == TRUE) fmt = SQFILE_UNKNOWN;
|
|
1300 else fmt = MSAFILE_SELEX;
|
|
1301
|
|
1302 DONE:
|
|
1303 if (buf != NULL) free(buf);
|
|
1304 rewind(fp);
|
|
1305 return fmt;
|
|
1306 }
|
|
1307
|
|
1308 /* Function: GCGBinaryToSequence()
|
|
1309 *
|
|
1310 * Purpose: Convert a GCG 2BIT binary string to DNA sequence.
|
|
1311 * 0 = C 1 = T 2 = A 3 = G
|
|
1312 * 4 nts/byte
|
|
1313 *
|
|
1314 * Args: seq - binary sequence. Converted in place to DNA.
|
|
1315 * len - length of DNA. binary is (len+3)/4 bytes
|
|
1316 */
|
|
1317 int
|
|
1318 GCGBinaryToSequence(char *seq, int len)
|
|
1319 {
|
|
1320 int bpos; /* position in binary */
|
|
1321 int spos; /* position in sequence */
|
|
1322 char twobit;
|
|
1323 int i;
|
|
1324
|
|
1325 for (bpos = (len-1)/4; bpos >= 0; bpos--)
|
|
1326 {
|
|
1327 twobit = seq[bpos];
|
|
1328 spos = bpos*4;
|
|
1329
|
|
1330 for (i = 3; i >= 0; i--)
|
|
1331 {
|
|
1332 switch (twobit & 0x3) {
|
|
1333 case 0: seq[spos+i] = 'C'; break;
|
|
1334 case 1: seq[spos+i] = 'T'; break;
|
|
1335 case 2: seq[spos+i] = 'A'; break;
|
|
1336 case 3: seq[spos+i] = 'G'; break;
|
|
1337 }
|
|
1338 twobit = twobit >> 2;
|
|
1339 }
|
|
1340 }
|
|
1341 seq[len] = '\0';
|
|
1342 return 1;
|
|
1343 }
|
|
1344
|
|
1345
|
|
1346 /* Function: GCGchecksum()
|
|
1347 * Date: SRE, Mon May 31 11:13:21 1999 [St. Louis]
|
|
1348 *
|
|
1349 * Purpose: Calculate a GCG checksum for a sequence.
|
|
1350 * Code provided by Steve Smith of Genetics
|
|
1351 * Computer Group.
|
|
1352 *
|
|
1353 * Args: seq - sequence to calculate checksum for.
|
|
1354 * may contain gap symbols.
|
|
1355 * len - length of sequence (usually known,
|
|
1356 * so save a strlen() call)
|
|
1357 *
|
|
1358 * Returns: GCG checksum.
|
|
1359 */
|
|
1360 int
|
|
1361 GCGchecksum(char *seq, int len)
|
|
1362 {
|
|
1363 int i; /* position in sequence */
|
|
1364 int chk = 0; /* calculated checksum */
|
|
1365
|
|
1366 for (i = 0; i < len; i++)
|
|
1367 chk = (chk + (i % 57 + 1) * (sre_toupper((int) seq[i]))) % 10000;
|
|
1368 return chk;
|
|
1369 }
|
|
1370
|
|
1371
|
|
1372 /* Function: GCGMultchecksum()
|
|
1373 *
|
|
1374 * Purpose: GCG checksum for a multiple alignment: sum of
|
|
1375 * individual sequence checksums (including their
|
|
1376 * gap characters) modulo 10000.
|
|
1377 *
|
|
1378 * Implemented using spec provided by Steve Smith of
|
|
1379 * Genetics Computer Group.
|
|
1380 *
|
|
1381 * Args: seqs - sequences to be checksummed; aligned or not
|
|
1382 * nseq - number of sequences
|
|
1383 *
|
|
1384 * Return: the checksum, a number between 0 and 9999
|
|
1385 */
|
|
1386 int
|
|
1387 GCGMultchecksum(char **seqs, int nseq)
|
|
1388 {
|
|
1389 int chk = 0;
|
|
1390 int idx;
|
|
1391
|
|
1392 for (idx = 0; idx < nseq; idx++)
|
|
1393 chk = (chk + GCGchecksum(seqs[idx], strlen(seqs[idx]))) % 10000;
|
|
1394 return chk;
|
|
1395 }
|
|
1396
|
|
1397
|
|
1398
|
|
1399
|
|
1400 /* Function: Seqtype()
|
|
1401 *
|
|
1402 * Purpose: Returns a (very good) guess about type of sequence:
|
|
1403 * kDNA, kRNA, kAmino, or kOtherSeq.
|
|
1404 *
|
|
1405 * Modified from, and replaces, Gilbert getseqtype().
|
|
1406 */
|
|
1407 int
|
|
1408 Seqtype(char *seq)
|
|
1409 {
|
|
1410 int saw; /* how many non-gap characters I saw */
|
|
1411 char c;
|
|
1412 int po = 0; /* count of protein-only */
|
|
1413 int nt = 0; /* count of t's */
|
|
1414 int nu = 0; /* count of u's */
|
|
1415 int na = 0; /* count of nucleotides */
|
|
1416 int aa = 0; /* count of amino acids */
|
|
1417 int no = 0; /* count of others */
|
|
1418
|
|
1419 /* Look at the first 300 non-gap characters
|
|
1420 */
|
|
1421
|
|
1422 #ifdef CLUSTALO
|
|
1423 /* VGGNGDDYLSGGTGNDTL is recognized as unknown using squid's default
|
|
1424 * approach.
|
|
1425 * We change it to the following:
|
|
1426
|
|
1427 * 1. counting: ignore gaps and not alpha characters. if protein-only then
|
|
1428 * count as such (po). otherwise decide if amino-acid (aa) or nucleic-acid
|
|
1429 * (na) or unknown (no)
|
|
1430 *
|
|
1431 * 2. determine type: if we saw more unknown than aa or na, return unknown.
|
|
1432 * if encountered protein-only return protein-only. otherwise decide based
|
|
1433 * on majority. (if aa==na return na)
|
|
1434 */
|
|
1435 for (saw = 0; *seq != '\0' && saw < 300; seq++) {
|
|
1436 c = sre_toupper((int) *seq);
|
|
1437 int unknown = 1;
|
|
1438
|
|
1439 if (isgap(c) || ! isalpha((int) c)) {
|
|
1440 continue;
|
|
1441 }
|
|
1442
|
|
1443 if (strchr(protonly, c)) {
|
|
1444 po++;
|
|
1445 unknown = 0;
|
|
1446 }
|
|
1447
|
|
1448 if (strchr(aminos,c)) {
|
|
1449 aa++;
|
|
1450 unknown = 0;
|
|
1451 }
|
|
1452
|
|
1453 if (strchr(primenuc,c)) {
|
|
1454 na++;
|
|
1455 unknown = 0;
|
|
1456
|
|
1457 if (c == 'T')
|
|
1458 nt++;
|
|
1459 else if (c == 'U')
|
|
1460 nu++;
|
|
1461 }
|
|
1462
|
|
1463 if (unknown) {
|
|
1464 no ++;
|
|
1465 }
|
|
1466
|
|
1467 saw++;
|
|
1468 }
|
|
1469
|
|
1470 if (no > aa && no > na)
|
|
1471 return kOtherSeq;
|
|
1472
|
|
1473 if (po > 0 || aa>na)
|
|
1474 return kAmino;
|
|
1475
|
|
1476 if (na >= aa) {
|
|
1477 if (nu > nt)
|
|
1478 return kRNA;
|
|
1479 else
|
|
1480 return kDNA;
|
|
1481 }
|
|
1482
|
|
1483 return kOtherSeq;
|
|
1484
|
|
1485
|
|
1486 #else
|
|
1487 for (saw = 0; *seq != '\0' && saw < 300; seq++)
|
|
1488 {
|
|
1489 c = sre_toupper((int) *seq);
|
|
1490 if (! isgap(c))
|
|
1491 {
|
|
1492 if (strchr(protonly, c)) po++;
|
|
1493 else if (strchr(primenuc,c)) {
|
|
1494 na++;
|
|
1495 if (c == 'T') nt++;
|
|
1496 else if (c == 'U') nu++;
|
|
1497 }
|
|
1498 else if (strchr(aminos,c)) aa++;
|
|
1499 else if (isalpha((int) c)) no++;
|
|
1500 saw++;
|
|
1501 }
|
|
1502 }
|
|
1503
|
|
1504 if (no > 0) return kOtherSeq;
|
|
1505 else if (po > 0) return kAmino;
|
|
1506 else if (na > aa) {
|
|
1507 if (nu > nt) return kRNA;
|
|
1508 else return kDNA;
|
|
1509 }
|
|
1510 else return kAmino; /* ooooh. risky. */
|
|
1511 #endif
|
|
1512
|
|
1513 }
|
|
1514
|
|
1515
|
|
1516 /* Function: GuessAlignmentSeqtype()
|
|
1517 * Date: SRE, Wed Jul 7 09:42:34 1999 [St. Louis]
|
|
1518 *
|
|
1519 * Purpose: Try to guess whether an alignment is protein
|
|
1520 * or nucleic acid; return a code for the
|
|
1521 * type (kRNA, kDNA, or kAmino).
|
|
1522 *
|
|
1523 * Args: aseq - array of aligned sequences. (Could also
|
|
1524 * be an rseq unaligned sequence array)
|
|
1525 * nseq - number of aseqs
|
|
1526 *
|
|
1527 * Returns: kRNA, kDNA, kAmino;
|
|
1528 * kOtherSeq if inconsistency is detected.
|
|
1529 */
|
|
1530 int
|
|
1531 GuessAlignmentSeqtype(char **aseq, int nseq)
|
|
1532 {
|
|
1533 int idx;
|
|
1534 int nrna = 0;
|
|
1535 int ndna = 0;
|
|
1536 int namino = 0;
|
|
1537 int nother = 0;
|
|
1538
|
|
1539 for (idx = 0; idx < nseq; idx++)
|
|
1540 switch (Seqtype(aseq[idx])) {
|
|
1541 case kRNA: nrna++; break;
|
|
1542 case kDNA: ndna++; break;
|
|
1543 case kAmino: namino++; break;
|
|
1544 default: nother++;
|
|
1545 }
|
|
1546
|
|
1547 /* Unambiguous decisions:
|
|
1548 */
|
|
1549 if (nother) return kOtherSeq;
|
|
1550 if (namino == nseq) return kAmino;
|
|
1551 if (ndna == nseq) return kDNA;
|
|
1552 if (nrna == nseq) return kRNA;
|
|
1553
|
|
1554 /* Ambiguous decisions:
|
|
1555 */
|
|
1556 if (namino == 0) return kRNA; /* it's nucleic acid, but seems mixed RNA/DNA */
|
|
1557 return kAmino; /* some amino acid seen; others probably short seqs, some
|
|
1558 of which may be entirely ACGT (ala,cys,gly,thr). We
|
|
1559 could be a little more sophisticated: U would be a giveaway
|
|
1560 that we're not in protein seqs */
|
|
1561 }
|
|
1562
|
|
1563 /* Function: WriteSimpleFASTA()
|
|
1564 * Date: SRE, Tue Nov 16 18:06:00 1999 [St. Louis]
|
|
1565 *
|
|
1566 * Purpose: Just write a FASTA format sequence to a file;
|
|
1567 * minimal interface, mostly for quick and dirty programs.
|
|
1568 *
|
|
1569 * Args: fp - open file handle (stdout, possibly)
|
|
1570 * seq - sequence to output
|
|
1571 * name - name for the sequence
|
|
1572 * desc - optional description line, or NULL.
|
|
1573 *
|
|
1574 * Returns: void
|
|
1575 */
|
|
1576 void
|
|
1577 WriteSimpleFASTA(FILE *fp, char *seq, char *name, char *desc)
|
|
1578 {
|
|
1579 char buf[61];
|
|
1580 int len;
|
|
1581 int pos;
|
|
1582
|
|
1583 len = strlen(seq);
|
|
1584 buf[60] = '\0';
|
|
1585 fprintf(fp, ">%s %s\n", name, desc != NULL ? desc : "");
|
|
1586 for (pos = 0; pos < len; pos += 60)
|
|
1587 {
|
|
1588 strncpy(buf, seq+pos, 60);
|
|
1589 fprintf(fp, "%s\n", buf);
|
|
1590 }
|
|
1591 }
|
|
1592
|
|
1593 int
|
|
1594 WriteSeq(FILE *outf, int outform, char *seq, SQINFO *sqinfo)
|
|
1595 {
|
|
1596 int numline = 0;
|
|
1597 int lines = 0, spacer = 0, width = 50, tab = 0;
|
|
1598 int i, j, l, l1, ibase;
|
|
1599 char endstr[10];
|
|
1600 char s[100]; /* buffer for sequence */
|
|
1601 char ss[100]; /* buffer for structure */
|
|
1602 int checksum = 0;
|
|
1603 int seqlen;
|
|
1604 int which_case; /* 0 = do nothing. 1 = upper case. 2 = lower case */
|
|
1605 int dostruc; /* TRUE to print structure lines*/
|
|
1606
|
|
1607 which_case = 0;
|
|
1608 dostruc = FALSE;
|
|
1609 seqlen = (sqinfo->flags & SQINFO_LEN) ? sqinfo->len : strlen(seq);
|
|
1610
|
|
1611 if (IsAlignmentFormat(outform))
|
|
1612 Die("Tried to write an aligned format with WriteSeq() -- bad, bad.");
|
|
1613
|
|
1614
|
|
1615 strcpy( endstr,"");
|
|
1616 l1 = 0;
|
|
1617 checksum = GCGchecksum(seq, seqlen);
|
|
1618
|
|
1619 switch (outform) {
|
|
1620 case SQFILE_UNKNOWN: /* no header, just sequence */
|
|
1621 strcpy(endstr,"\n"); /* end w/ extra blank line */
|
|
1622 break;
|
|
1623
|
|
1624 case SQFILE_GENBANK:
|
|
1625 fprintf(outf,"LOCUS %s %d bp\n",
|
|
1626 sqinfo->name, seqlen);
|
|
1627 fprintf(outf,"ACCESSION %s\n",
|
|
1628 (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : ".");
|
|
1629 fprintf(outf,"DEFINITION %s\n",
|
|
1630 (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : ".");
|
|
1631 fprintf(outf,"VERSION %s\n",
|
|
1632 (sqinfo->flags & SQINFO_ID) ? sqinfo->id : ".");
|
|
1633 fprintf(outf,"ORIGIN \n");
|
|
1634 spacer = 11;
|
|
1635 numline = 1;
|
|
1636 strcpy(endstr, "\n//");
|
|
1637 break;
|
|
1638
|
|
1639 case SQFILE_GCGDATA:
|
|
1640 fprintf(outf, ">>>>%s 9/95 ASCII Len: %d\n", sqinfo->name, seqlen);
|
|
1641 fprintf(outf, "%s\n", (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
|
|
1642 break;
|
|
1643
|
|
1644 case SQFILE_PIR:
|
|
1645 fprintf(outf, "ENTRY %s\n",
|
|
1646 (sqinfo->flags & SQINFO_ID) ? sqinfo->id : sqinfo->name);
|
|
1647 fprintf(outf, "TITLE %s\n",
|
|
1648 (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
|
|
1649 fprintf(outf, "ACCESSION %s\n",
|
|
1650 (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : "-");
|
|
1651 fprintf(outf, "SUMMARY #Length %d #Checksum %d\n",
|
|
1652 sqinfo->len, checksum);
|
|
1653 fprintf(outf, "SEQUENCE\n");
|
|
1654 fprintf(outf, " 5 10 15 20 25 30\n");
|
|
1655 spacer = 2; /* spaces after every residue */
|
|
1656 numline = 1; /* number lines w/ coords */
|
|
1657 width = 30; /* 30 aa per line */
|
|
1658 strcpy(endstr, "\n///");
|
|
1659 break;
|
|
1660
|
|
1661 case SQFILE_SQUID:
|
|
1662 fprintf(outf, "NAM %s\n", sqinfo->name);
|
|
1663 if (sqinfo->flags & (SQINFO_ID | SQINFO_ACC | SQINFO_START | SQINFO_STOP | SQINFO_OLEN))
|
|
1664 fprintf(outf, "SRC %s %s %d..%d::%d\n",
|
|
1665 (sqinfo->flags & SQINFO_ID) ? sqinfo->id : "-",
|
|
1666 (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : "-",
|
|
1667 (sqinfo->flags & SQINFO_START) ? sqinfo->start : 0,
|
|
1668 (sqinfo->flags & SQINFO_STOP) ? sqinfo->stop : 0,
|
|
1669 (sqinfo->flags & SQINFO_OLEN) ? sqinfo->olen : 0);
|
|
1670 if (sqinfo->flags & SQINFO_DESC)
|
|
1671 fprintf(outf, "DES %s\n", sqinfo->desc);
|
|
1672 if (sqinfo->flags & SQINFO_SS)
|
|
1673 {
|
|
1674 fprintf(outf, "SEQ +SS\n");
|
|
1675 dostruc = TRUE; /* print structure lines too */
|
|
1676 }
|
|
1677 else
|
|
1678 fprintf(outf, "SEQ\n");
|
|
1679 numline = 1; /* number seq lines w/ coords */
|
|
1680 strcpy(endstr, "\n++");
|
|
1681 break;
|
|
1682
|
|
1683 case SQFILE_EMBL:
|
|
1684 fprintf(outf,"ID %s\n",
|
|
1685 (sqinfo->flags & SQINFO_ID) ? sqinfo->id : sqinfo->name);
|
|
1686 fprintf(outf,"AC %s\n",
|
|
1687 (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : "-");
|
|
1688 fprintf(outf,"DE %s\n",
|
|
1689 (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
|
|
1690 fprintf(outf,"SQ %d BP\n", seqlen);
|
|
1691 strcpy(endstr, "\n//"); /* 11Oct90: bug fix*/
|
|
1692 tab = 5; /** added 31jan91 */
|
|
1693 spacer = 11; /** added 31jan91 */
|
|
1694 break;
|
|
1695
|
|
1696 case SQFILE_GCG:
|
|
1697 fprintf(outf,"%s\n", sqinfo->name);
|
|
1698 if (sqinfo->flags & SQINFO_ACC)
|
|
1699 fprintf(outf,"ACCESSION %s\n", sqinfo->acc);
|
|
1700 if (sqinfo->flags & SQINFO_DESC)
|
|
1701 fprintf(outf,"DEFINITION %s\n", sqinfo->desc);
|
|
1702 fprintf(outf," %s Length: %d (today) Check: %d ..\n",
|
|
1703 sqinfo->name, seqlen, checksum);
|
|
1704 spacer = 11;
|
|
1705 numline = 1;
|
|
1706 strcpy(endstr, "\n"); /* this is insurance to help prevent misreads at eof */
|
|
1707 break;
|
|
1708
|
|
1709 case SQFILE_STRIDER: /* ?? map ?*/
|
|
1710 fprintf(outf,"; ### from DNA Strider ;-)\n");
|
|
1711 fprintf(outf,"; DNA sequence %s, %d bases, %d checksum.\n;\n",
|
|
1712 sqinfo->name, seqlen, checksum);
|
|
1713 strcpy(endstr, "\n//");
|
|
1714 break;
|
|
1715
|
|
1716 /* SRE: Don had Zuker default to Pearson, which is not
|
|
1717 intuitive or helpful, since Zuker's MFOLD can't read
|
|
1718 Pearson format. More useful to use kIG */
|
|
1719 case SQFILE_ZUKER:
|
|
1720 which_case = 1; /* MFOLD requires upper case. */
|
|
1721 /*FALLTHRU*/
|
|
1722 case SQFILE_IG:
|
|
1723 fprintf(outf,";%s %s\n",
|
|
1724 sqinfo->name,
|
|
1725 (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "");
|
|
1726 fprintf(outf,"%s\n", sqinfo->name);
|
|
1727 strcpy(endstr,"1"); /* == linear dna */
|
|
1728 break;
|
|
1729
|
|
1730 case SQFILE_RAW: /* Raw: no header at all. */
|
|
1731 break;
|
|
1732
|
|
1733 default :
|
|
1734 case SQFILE_FASTA:
|
|
1735 fprintf(outf,">%s %s\n", sqinfo->name,
|
|
1736 (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "");
|
|
1737 break;
|
|
1738 }
|
|
1739
|
|
1740 if (which_case == 1) s2upper(seq);
|
|
1741 if (which_case == 2) s2lower(seq);
|
|
1742
|
|
1743
|
|
1744 width = MIN(width,100);
|
|
1745 for (i=0, l=0, ibase = 1, lines = 0; i < seqlen; ) {
|
|
1746 if (l1 < 0) l1 = 0;
|
|
1747 else if (l1 == 0) {
|
|
1748 if (numline) fprintf(outf,"%8d ",ibase);
|
|
1749 for (j=0; j<tab; j++) fputc(' ',outf);
|
|
1750 }
|
|
1751 if ((spacer != 0) && ((l+1) % spacer == 1))
|
|
1752 { s[l] = ' '; ss[l] = ' '; l++; }
|
|
1753 s[l] = seq[i];
|
|
1754 ss[l] = (sqinfo->flags & SQINFO_SS) ? sqinfo->ss[i] : '.';
|
|
1755 l++; i++;
|
|
1756 l1++; /* don't count spaces for width*/
|
|
1757 if (l1 == width || i == seqlen) {
|
|
1758 s[l] = ss[l] = '\0';
|
|
1759 l = 0; l1 = 0;
|
|
1760 if (dostruc)
|
|
1761 {
|
|
1762 fprintf(outf, "%s\n", s);
|
|
1763 if (numline) fprintf(outf," ");
|
|
1764 for (j=0; j<tab; j++) fputc(' ',outf);
|
|
1765 if (i == seqlen) fprintf(outf,"%s%s\n",ss,endstr);
|
|
1766 else fprintf(outf,"%s\n",ss);
|
|
1767 }
|
|
1768 else
|
|
1769 {
|
|
1770 if (i == seqlen) fprintf(outf,"%s%s\n",s,endstr);
|
|
1771 else fprintf(outf,"%s\n",s);
|
|
1772 }
|
|
1773 lines++;
|
|
1774 ibase = i+1;
|
|
1775 }
|
|
1776 }
|
|
1777 return lines;
|
|
1778 }
|
|
1779
|
|
1780
|
|
1781 /* Function: ReadMultipleRseqs()
|
|
1782 *
|
|
1783 * Purpose: Open a data file and
|
|
1784 * parse it into an array of rseqs (raw, unaligned
|
|
1785 * sequences).
|
|
1786 *
|
|
1787 * Caller is responsible for free'ing memory allocated
|
|
1788 * to ret_rseqs, ret_weights, and ret_names.
|
|
1789 *
|
|
1790 * Weights are currently only supported for MSF format.
|
|
1791 * Sequences read from all other formats will be assigned
|
|
1792 * weights of 1.0. If the caller isn't interested in
|
|
1793 * weights, it passes NULL as ret_weights.
|
|
1794 *
|
|
1795 * Returns 1 on success. Returns 0 on failure and sets
|
|
1796 * squid_errno to indicate the cause.
|
|
1797 */
|
|
1798 int
|
|
1799 ReadMultipleRseqs(char *seqfile,
|
|
1800 int fformat,
|
|
1801 char ***ret_rseqs,
|
|
1802 SQINFO **ret_sqinfo,
|
|
1803 int *ret_num)
|
|
1804 {
|
|
1805 SQINFO *sqinfo; /* array of sequence optional info */
|
|
1806 SQFILE *dbfp; /* open ptr for sequential access of file */
|
|
1807 char **rseqs; /* sequence array */
|
|
1808 int numalloced; /* num of seqs currently alloced for */
|
|
1809 int num;
|
|
1810
|
|
1811
|
|
1812 num = 0;
|
|
1813 numalloced = 16;
|
|
1814 rseqs = (char **) MallocOrDie (numalloced * sizeof(char *));
|
|
1815 sqinfo = (SQINFO *) MallocOrDie (numalloced * sizeof(SQINFO));
|
|
1816 if ((dbfp = SeqfileOpen(seqfile, fformat, NULL)) == NULL) return 0;
|
|
1817
|
|
1818 while (ReadSeq(dbfp, dbfp->format, &rseqs[num], &(sqinfo[num])))
|
|
1819 {
|
|
1820 num++;
|
|
1821 if (num == numalloced) /* more seqs coming, alloc more room */
|
|
1822 {
|
|
1823 numalloced += 16;
|
|
1824 rseqs = (char **) ReallocOrDie (rseqs, numalloced*sizeof(char *));
|
|
1825 sqinfo = (SQINFO *) ReallocOrDie (sqinfo, numalloced * sizeof(SQINFO));
|
|
1826 }
|
|
1827 }
|
|
1828 SeqfileClose(dbfp);
|
|
1829
|
|
1830 *ret_rseqs = rseqs;
|
|
1831 *ret_sqinfo = sqinfo;
|
|
1832 *ret_num = num;
|
|
1833 return 1;
|
|
1834 }
|
|
1835
|
|
1836
|
|
1837 /* Function: String2SeqfileFormat()
|
|
1838 * Date: SRE, Sun Jun 27 15:25:54 1999 [TW 723 over Canadian Shield]
|
|
1839 *
|
|
1840 * Purpose: Convert a string (e.g. from command line option arg)
|
|
1841 * to a format code. Case insensitive. Return
|
|
1842 * MSAFILE_UNKNOWN/SQFILE_UNKNOWN if string is bad.
|
|
1843 * Uses codes defined in squid.h (unaligned formats) and
|
|
1844 * msa.h (aligned formats).
|
|
1845 *
|
|
1846 * Args: s - string to convert; e.g. "stockholm"
|
|
1847 *
|
|
1848 * Returns: format code; e.g. MSAFILE_STOCKHOLM
|
|
1849 */
|
|
1850 int
|
|
1851 String2SeqfileFormat(char *s)
|
|
1852 {
|
|
1853 char *s2;
|
|
1854 int code = SQFILE_UNKNOWN;
|
|
1855
|
|
1856 if (s == NULL) return SQFILE_UNKNOWN;
|
|
1857 s2 = sre_strdup(s, -1);
|
|
1858 s2upper(s2);
|
|
1859
|
|
1860 if (strcmp(s2, "FASTA") == 0) code = SQFILE_FASTA;
|
|
1861 #ifdef CLUSTALO
|
|
1862 if (strcmp(s2, "FA") == 0) code = SQFILE_FASTA;
|
|
1863 else if (strcmp(s2, "VIENNA") == 0) code = SQFILE_VIENNA;
|
|
1864 else if (strcmp(s2, "VIE") == 0) code = SQFILE_VIENNA;
|
|
1865 #endif
|
|
1866 else if (strcmp(s2, "GENBANK") == 0) code = SQFILE_GENBANK;
|
|
1867 #ifdef CLUSTALO
|
|
1868 else if (strcmp(s2, "GB") == 0) code = SQFILE_GENBANK;
|
|
1869 #endif
|
|
1870 else if (strcmp(s2, "EMBL") == 0) code = SQFILE_EMBL;
|
|
1871 else if (strcmp(s2, "GCG") == 0) code = SQFILE_GCG;
|
|
1872 else if (strcmp(s2, "GCGDATA") == 0) code = SQFILE_GCGDATA;
|
|
1873 else if (strcmp(s2, "RAW") == 0) code = SQFILE_RAW;
|
|
1874 else if (strcmp(s2, "IG") == 0) code = SQFILE_IG;
|
|
1875 else if (strcmp(s2, "STRIDER") == 0) code = SQFILE_STRIDER;
|
|
1876 else if (strcmp(s2, "IDRAW") == 0) code = SQFILE_IDRAW;
|
|
1877 else if (strcmp(s2, "ZUKER") == 0) code = SQFILE_ZUKER;
|
|
1878 else if (strcmp(s2, "PIR") == 0) code = SQFILE_PIR;
|
|
1879 else if (strcmp(s2, "SQUID") == 0) code = SQFILE_SQUID;
|
|
1880 else if (strcmp(s2, "STOCKHOLM") == 0) code = MSAFILE_STOCKHOLM;
|
|
1881 #ifdef CLUSTALO
|
|
1882 else if (strcmp(s2, "ST") == 0) code = MSAFILE_STOCKHOLM;
|
|
1883 else if (strcmp(s2, "STK") == 0) code = MSAFILE_STOCKHOLM;
|
|
1884 #endif
|
|
1885 else if (strcmp(s2, "SELEX") == 0) code = MSAFILE_SELEX;
|
|
1886 else if (strcmp(s2, "MSF") == 0) code = MSAFILE_MSF;
|
|
1887 else if (strcmp(s2, "CLUSTAL") == 0) code = MSAFILE_CLUSTAL;
|
|
1888 #ifdef CLUSTALO
|
|
1889 else if (strcmp(s2, "CLU") == 0) code = MSAFILE_CLUSTAL;
|
|
1890 #endif
|
|
1891 else if (strcmp(s2, "A2M") == 0) code = MSAFILE_A2M;
|
|
1892 else if (strcmp(s2, "PHYLIP") == 0) code = MSAFILE_PHYLIP;
|
|
1893 #ifdef CLUSTALO
|
|
1894 else if (strcmp(s2, "PHY") == 0) code = MSAFILE_PHYLIP;
|
|
1895 #endif
|
|
1896 else if (strcmp(s2, "EPS") == 0) code = MSAFILE_EPS;
|
|
1897 #ifdef CLUSTALO
|
|
1898 else code = SQFILE_UNKNOWN;
|
|
1899 #endif
|
|
1900 free(s2);
|
|
1901 return code;
|
|
1902 }
|
|
1903 char *
|
|
1904 SeqfileFormat2String(int code)
|
|
1905 {
|
|
1906 switch (code) {
|
|
1907 case SQFILE_UNKNOWN: return "unknown";
|
|
1908 case SQFILE_FASTA: return "FASTA";
|
|
1909 #ifdef CLUSTALO
|
|
1910 case SQFILE_VIENNA: return "Vienna";
|
|
1911 #endif
|
|
1912 case SQFILE_GENBANK: return "Genbank";
|
|
1913 case SQFILE_EMBL: return "EMBL";
|
|
1914 case SQFILE_GCG: return "GCG";
|
|
1915 case SQFILE_GCGDATA: return "GCG data library";
|
|
1916 case SQFILE_RAW: return "raw";
|
|
1917 case SQFILE_IG: return "Intelligenetics";
|
|
1918 case SQFILE_STRIDER: return "MacStrider";
|
|
1919 case SQFILE_IDRAW: return "Idraw Postscript";
|
|
1920 case SQFILE_ZUKER: return "Zuker";
|
|
1921 case SQFILE_PIR: return "PIR";
|
|
1922 case SQFILE_SQUID: return "SQUID";
|
|
1923 case MSAFILE_STOCKHOLM: return "Stockholm";
|
|
1924 case MSAFILE_SELEX: return "SELEX";
|
|
1925 case MSAFILE_MSF: return "MSF";
|
|
1926 case MSAFILE_CLUSTAL: return "Clustal";
|
|
1927 case MSAFILE_A2M: return "a2m";
|
|
1928 case MSAFILE_PHYLIP: return "Phylip";
|
|
1929 case MSAFILE_EPS: return "EPS";
|
|
1930 default:
|
|
1931 Die("Bad code passed to MSAFormat2String()");
|
|
1932 }
|
|
1933 /*NOTREACHED*/
|
|
1934 return NULL;
|
|
1935 }
|
|
1936
|
|
1937
|
|
1938 /* Function: MSAToSqinfo()
|
|
1939 * Date: SRE, Tue Jul 20 14:36:56 1999 [St. Louis]
|
|
1940 *
|
|
1941 * Purpose: Take an MSA and generate a SQINFO array suitable
|
|
1942 * for use in annotating the unaligned sequences.
|
|
1943 * Return the array.
|
|
1944 *
|
|
1945 * Permanent temporary code. sqinfo was poorly designed.
|
|
1946 * it must eventually be replaced, but the odds
|
|
1947 * of this happening soon are nil, so I have to deal.
|
|
1948 *
|
|
1949 * Args: msa - the alignment
|
|
1950 *
|
|
1951 * Returns: ptr to allocated sqinfo array.
|
|
1952 * Freeing is ghastly: free in each individual sqinfo[i]
|
|
1953 * with FreeSequence(NULL, &(sqinfo[i])), then
|
|
1954 * free(sqinfo).
|
|
1955 */
|
|
1956 SQINFO *
|
|
1957 MSAToSqinfo(MSA *msa)
|
|
1958 {
|
|
1959 int idx;
|
|
1960 SQINFO *sqinfo;
|
|
1961
|
|
1962 sqinfo = MallocOrDie(sizeof(SQINFO) * msa->nseq);
|
|
1963
|
|
1964 for (idx = 0; idx < msa->nseq; idx++)
|
|
1965 {
|
|
1966 sqinfo[idx].flags = 0;
|
|
1967 SetSeqinfoString(&(sqinfo[idx]),
|
|
1968 msa->sqname[idx], SQINFO_NAME);
|
|
1969 SetSeqinfoString(&(sqinfo[idx]),
|
|
1970 MSAGetSeqAccession(msa, idx), SQINFO_ACC);
|
|
1971 SetSeqinfoString(&(sqinfo[idx]),
|
|
1972 MSAGetSeqDescription(msa, idx), SQINFO_DESC);
|
|
1973
|
|
1974 if (msa->ss != NULL && msa->ss[idx] != NULL) {
|
|
1975 MakeDealignedString(msa->aseq[idx], msa->alen,
|
|
1976 msa->ss[idx], &(sqinfo[idx].ss));
|
|
1977 sqinfo[idx].flags |= SQINFO_SS;
|
|
1978 }
|
|
1979
|
|
1980 if (msa->sa != NULL && msa->sa[idx] != NULL) {
|
|
1981 MakeDealignedString(msa->aseq[idx], msa->alen,
|
|
1982 msa->sa[idx], &(sqinfo[idx].sa));
|
|
1983 sqinfo[idx].flags |= SQINFO_SA;
|
|
1984 }
|
|
1985
|
|
1986 sqinfo[idx].len = DealignedLength(msa->aseq[idx]);
|
|
1987 sqinfo[idx].flags |= SQINFO_LEN;
|
|
1988 }
|
|
1989 return sqinfo;
|
|
1990 }
|
|
1991
|
|
1992
|
|
1993
|
|
1994 /* cc -o sqio_test -DA_QUIET_DAY -L. sqio.c -lsquid */
|
|
1995 #ifdef A_QUIET_DAY
|
|
1996 #include "ssi.h"
|
|
1997 int
|
|
1998 main(int argc, char **argv)
|
|
1999 {
|
|
2000 FILE *fp;
|
|
2001 char *filename;
|
|
2002 char *buf;
|
|
2003 int len;
|
|
2004 int mode = 3;
|
|
2005 SSIOFFSET off;
|
|
2006
|
|
2007 filename = argv[1];
|
|
2008
|
|
2009 if (mode == 1) {
|
|
2010 buf = malloc(sizeof(char) * 256);
|
|
2011 if ((fp = fopen(filename, "r")) == NULL)
|
|
2012 Die("open of %s failed", filename);
|
|
2013 while (fgets(buf, 255, fp) != NULL)
|
|
2014 ;
|
|
2015 fclose(fp);
|
|
2016 free(buf);
|
|
2017 } else if (mode == 2) {
|
|
2018 if ((fp = fopen(filename, "r")) == NULL)
|
|
2019 Die("open of %s failed", filename);
|
|
2020 buf = NULL; len = 0;
|
|
2021 while (sre_fgets(&buf, &len, fp) != NULL)
|
|
2022 SSIGetFilePosition(fp, SSI_OFFSET_I32, &off);
|
|
2023 fclose(fp);
|
|
2024 free(buf);
|
|
2025 } else if (mode == 3) {
|
|
2026 SQFILE *dbfp;
|
|
2027 SQINFO info;
|
|
2028
|
|
2029 if ((dbfp = SeqfileOpen(filename, SQFILE_FASTA, NULL)) == NULL)
|
|
2030 Die("open of %s failed", filename);
|
|
2031 while (ReadSeq(dbfp, dbfp->format, &buf, &info)) {
|
|
2032 SSIGetFilePosition(dbfp->f, SSI_OFFSET_I32, &off);
|
|
2033 FreeSequence(buf, &info);
|
|
2034 }
|
|
2035 SeqfileClose(dbfp);
|
|
2036 }
|
|
2037
|
|
2038 }
|
|
2039
|
|
2040
|
|
2041 #endif
|