source: tags/ms_r16q3/ARB_GDE/GDE_Genbank.cxx

Last change on this file was 15176, checked in by westram, 8 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 KB
Line 
1#include "GDE_extglob.h"
2#include <ctime>
3#include <algorithm>
4
5/*
6  Copyright (c) 1989-1990, University of Illinois board of trustees.  All
7  rights reserved.  Written by Steven Smith at the Center for Prokaryote Genome
8  Analysis.  Design and implementation guidance by Dr. Gary Olsen and Dr.
9  Carl Woese.
10
11  Copyright (c) 1990,1991,1992 Steven Smith at the Harvard Genome Laboratory.
12  all rights reserved.
13
14  Copyright (c) 1993, Steven Smith, all rights reserved.
15
16*/
17
18static bool CheckType(char *seq, int len) {
19    /*   CheckType:  Check base composition to see if the sequence
20     *   appears to be an amino acid sequence.
21     */
22    int j, count1 = 0, count2 = 0;
23
24    for (j=0; j<len; j++)
25        if (((seq[j]|32) < 'z') && ((seq[j]|32) > 'a'))
26        {
27            count1++;
28            if (strchr("ACGTUNacgtun", seq[j]) == NULL)
29                count2++;
30        }
31
32    return ((count2 > count1/4) ? true : false);
33}
34
35// ARB
36struct ARB_TIME {
37    int yy;
38    int mm;
39    int dd;
40    int hr;
41    int mn;
42    int sc;
43};
44
45static void AsciiTime(void *b, char *asciitime)
46{
47    ARB_TIME *a=(ARB_TIME*)b;
48    int j;
49    char temp[GBUFSIZ];
50
51    a->dd = 0;
52    a->yy = 0;
53    a->mm = 0;
54    sscanf(asciitime, "%d%5c%d", &(a->dd), temp, &(a->yy));
55    temp[5] = '\0';
56    for (j=0; j<12; j++)
57        if (strcmp(temp, GDEmonth[j]) == 0)
58            a->mm = j+1;
59    if (a->dd <0 || a->dd > 31 || a->yy < 0 || a->mm > 11)
60        SetTime(a);
61    return;
62}
63// ENDARB
64
65GB_ERROR ReadGen(char *filename, NA_Alignment& dataset) {
66    GB_ERROR  error = NULL;
67    FILE     *file  = fopen(filename, "r");
68    if (!file) {
69        error = GB_IO_error("reading", filename);
70    }
71    else {
72        bool    done         = false;
73        size_t  len          = 0;
74        bool    IS_REALLY_AA = false;
75        char    in_line[GBUFSIZ], c;
76        char   *buffer       = 0, *gencomments = NULL, fields[8][GBUFSIZ];
77        size_t  buflen       = 0;
78        int     genclen      = 0, curelem = 0, n = 0;
79        int     start_col    = -1;
80
81        NA_Sequence *this_elem = 0;
82
83        for (; fgets(in_line, GBUFSIZ, file) != 0;)
84        {
85            if (in_line[strlen(in_line)-1] == '\n')
86                in_line[strlen(in_line)-1] = '\0';
87            if (Find(in_line, "LOCUS"))
88            {
89                curelem   = Arbdb_get_curelem(dataset);
90                this_elem = &(dataset.element[curelem]);
91                n         = sscanf(in_line, "%s %s %s %s %s %s %s %s",
92                                   fields[0], fields[1], fields[2], fields[3], fields[4],
93                                   fields[5], fields[6], fields[7]);
94
95                if (IS_REALLY_AA)
96                {
97                    InitNASeq(this_elem, PROTEIN);
98                }
99                else if (Find(in_line, "DNA"))
100                {
101                    InitNASeq(this_elem, DNA);
102                }
103                else if (Find(in_line, "RNA"))
104                {
105                    InitNASeq(this_elem, RNA);
106                }
107                else if (Find(in_line, "MASK"))
108                {
109                    InitNASeq(this_elem, MASK);
110                }
111                else if (Find(in_line, "TEXT"))
112                {
113                    InitNASeq(this_elem, TEXT);
114                }
115                else if (Find(in_line, "PROT"))
116                {
117                    InitNASeq(this_elem, PROTEIN);
118                }
119                else
120                    InitNASeq(this_elem, DNA);
121
122                strncpy_terminate(this_elem->short_name, fields[1], SIZE_SHORT_NAME);
123                AsciiTime(&(this_elem->t_stamp.origin), fields[n-1]);
124                this_elem->attr = DEFAULT_X_ATTR;
125
126                if (Find(in_line, "Circular"))
127                    this_elem->attr |= IS_CIRCULAR;
128
129                gencomments = NULL;
130                genclen = 0;
131            }
132            else if (Find(in_line, "DEFINITION"))
133                strncpy_terminate(this_elem->description, in_line+12, SIZE_DESCRIPTION);
134
135            else if (Find(in_line, "AUTHOR"))
136                strncpy_terminate(this_elem->authority, in_line+12, SIZE_AUTHORITY);
137
138            else if (Find(in_line, "  ORGANISM"))
139                strncpy_terminate(this_elem->seq_name, in_line+12, SIZE_SEQ_NAME);
140
141            else if (Find(in_line, "ACCESSION"))
142                strncpy_terminate(this_elem->id, in_line+12, SIZE_ID);
143
144            else if (Find(in_line, "ORIGIN"))
145            {
146                done = false;
147                len = 0;
148                for (; !done && fgets(in_line, GBUFSIZ, file) != 0;)
149                {
150                    if (in_line[0] != '/')
151                    {
152                        if (buflen == 0) {
153                            buflen = GBUFSIZ;
154                            ARB_calloc(buffer, buflen);
155                        }
156                        else if (len+strlen(in_line) >= buflen) {
157                            size_t new_buflen = buflen+GBUFSIZ;
158                            ARB_recalloc(buffer, buflen, new_buflen);
159                            buflen            = new_buflen;
160                        }
161                        // Search for the fist column of data (whitespace-number-whitespace)data
162                        if (start_col == -1)
163                        {
164                            for (start_col=0; in_line[start_col] == ' ' || in_line[start_col] == '\t'; start_col++) ;
165                            for (start_col++; strchr("1234567890", in_line[start_col]) != NULL; start_col++) ;
166                            for (start_col++; in_line[start_col] == ' ' || in_line[start_col] == '\t'; start_col++) ;
167                        }
168                        for (int j=start_col; (c = in_line[j]) != '\0'; j++)
169                        {
170                            if ((c != '\n') && ((j-start_col + 1) % 11 != 0)) {
171                                buffer[len++] = c;
172                            }
173                        }
174                    }
175                    else
176                    {
177                        AppendNA((NA_Base*)buffer, len, &(dataset.element[curelem]));
178                        for (size_t j=0; j<len; j++) buffer[j] = '\0';
179                        len = 0;
180                        done = true;
181                        dataset.element[curelem].comments = gencomments;
182                        dataset.element[curelem].comments_len= genclen - 1;
183                        dataset.element[curelem].comments_maxlen = genclen;
184
185                        gencomments = NULL;
186                        genclen = 0;
187                    }
188                }
189                /* Test if sequence should be converted by the translation table
190                 * If it looks like a protein...
191                 */
192                if (dataset.element[curelem].rmatrix && !IS_REALLY_AA) {
193                    IS_REALLY_AA = CheckType((char*)dataset.element[curelem]. sequence, dataset.element[curelem].seqlen);
194
195                    if (!IS_REALLY_AA)
196                        Ascii2NA((char*)dataset.element[curelem].sequence,
197                                 dataset.element[curelem].seqlen,
198                                 dataset.element[curelem].rmatrix);
199                    else {
200                        // Force the sequence to be AA
201                        dataset.element[curelem].elementtype = PROTEIN;
202                        dataset.element[curelem].rmatrix = NULL;
203                        dataset.element[curelem].tmatrix = NULL;
204                        dataset.element[curelem].col_lut = Default_PROColor_LKUP;
205                    }
206                }
207            }
208            else if (Find(in_line, "ZZZZZ"))
209            {
210                free(gencomments);
211                genclen = 0;
212            }
213            else
214            {
215                if (gencomments == NULL)
216                {
217                    gencomments = ARB_strdup(in_line);
218                    genclen = strlen(gencomments)+1;
219                }
220                else
221                {
222                    genclen += strlen(in_line)+1;
223                    ARB_realloc(gencomments, genclen);
224                    strncat(gencomments, in_line, GBUFSIZ);
225                    strncat(gencomments, "\n", GBUFSIZ);
226                }
227            }
228        }
229        free(buffer);
230        fclose(file);
231    }
232    for (size_t j=0; j<dataset.numelements; j++) {
233        dataset.maxlen = std::max(dataset.maxlen,
234                                  dataset.element[j].seqlen+dataset.element[j].offset);
235    }
236    return error;
237}
238
239int WriteGen(NA_Alignment& aln, char *filename, int method)
240{
241    int i;
242    size_t j;
243    int k;
244    NA_Sequence *this_elem;
245
246    FILE *file = fopen(filename, "w");
247    if (file == NULL)
248    {
249        Warning("Cannot open file for output");
250        return (1);
251    }
252
253    for (j=0; j<aln.numelements; j++)
254    {
255        this_elem = &(aln.element[j]);
256        if (method == ALL)
257        {
258            fprintf(file,
259                    "LOCUS       %10s%8d bp    %4s  %10s   %2d%5s%4d\n",
260                    this_elem->short_name, this_elem->seqlen+this_elem->offset,
261                    (this_elem->elementtype == DNA) ? "DNA" :
262                    (this_elem->elementtype == RNA) ? "RNA" :
263                    (this_elem->elementtype == MASK) ? "MASK" :
264                    (this_elem->elementtype == PROTEIN) ? "PROT" : "TEXT",
265                    this_elem->attr & IS_CIRCULAR ? "Circular" : "",
266                    this_elem->t_stamp.origin.dd,
267                    GDEmonth[this_elem->t_stamp.origin.mm-1],
268                    (this_elem->t_stamp.origin.yy>1900) ? this_elem->t_stamp.origin.yy :
269                    this_elem->t_stamp.origin.yy+1900);
270
271            if (this_elem->description[0])
272                fprintf(file, "DEFINITION  %s\n", this_elem->description);
273            if (this_elem->seq_name[0])
274                fprintf(file, "  ORGANISM  %s\n", this_elem->seq_name);
275            if (this_elem->id[0])
276                fprintf(file, " ACCESSION  %s\n", this_elem->id);
277            if (this_elem->authority[0])
278                fprintf(file, "   AUTHORS  %s\n", this_elem->authority);
279            if (this_elem->comments)
280                fprintf(file, "%s\n", this_elem->comments);
281            fprintf(file, "ORIGIN");
282            if (this_elem->tmatrix)
283            {
284                for (i=0, k=0; k<this_elem->seqlen+this_elem->offset; k++)
285                {
286                    if (i%60 == 0)
287                        fprintf(file, "\n%9d", i+1);
288                    if (i%10 == 0)
289                        fprintf(file, " ");
290                    fprintf(file, "%c", this_elem->tmatrix
291                            [getelem(this_elem, k)]);
292                    i++;
293                }
294            }
295            else
296            {
297                for (i=0, k=0; k<this_elem->seqlen+this_elem->offset; k++)
298                {
299                    if (i%60 == 0)
300                        fprintf(file, "\n%9d", i+1);
301                    if (i%10 == 0)
302                        fprintf(file, " ");
303                    fprintf(file, "%c", getelem(this_elem, k));
304                    i++;
305                }
306            }
307            fprintf(file, "\n//\n");
308        }
309    }
310    fclose(file);
311    return (0);
312}
313
314
315void SetTime(void *b)
316{
317    ARB_TIME *a=(ARB_TIME*)b;
318    struct tm *tim;
319    time_t clock;
320
321    clock = time(0);
322    tim = localtime(&clock);
323
324    a->yy = tim->tm_year;
325    a->mm = tim->tm_mon+1;
326    a->dd = tim->tm_mday;
327    a->hr = tim->tm_hour;
328    a->mn = tim->tm_min;
329    a->sc = tim->tm_sec;
330    return;
331}
332
Note: See TracBrowser for help on using the repository browser.