source: tags/ms_r16q2/ARB_GDE/GDE_Genbank.cxx

Last change on this file was 12833, checked in by westram, 11 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.8 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 = dataset.numelements++;
90                if (curelem == 0)
91                {
92                    dataset.element=(NA_Sequence*) Calloc(5, sizeof(NA_Sequence));
93                    dataset.maxnumelements = 5;
94                }
95                else if (curelem==dataset.maxnumelements)
96                {
97                    dataset.maxnumelements *= 2;
98                    dataset.element = (NA_Sequence*) Realloc((char*)dataset.element, dataset.maxnumelements * sizeof(NA_Sequence));
99                }
100                this_elem = &(dataset.element[curelem]);
101                n = sscanf(in_line, "%s %s %s %s %s %s %s %s",
102                           fields[0], fields[1], fields[2], fields[3], fields[4],
103                           fields[5], fields[6], fields[7]);
104                if (IS_REALLY_AA)
105                {
106                    InitNASeq(this_elem, PROTEIN);
107                }
108                else if (Find(in_line, "DNA"))
109                {
110                    InitNASeq(this_elem, DNA);
111                }
112                else if (Find(in_line, "RNA"))
113                {
114                    InitNASeq(this_elem, RNA);
115                }
116                else if (Find(in_line, "MASK"))
117                {
118                    InitNASeq(this_elem, MASK);
119                }
120                else if (Find(in_line, "TEXT"))
121                {
122                    InitNASeq(this_elem, TEXT);
123                }
124                else if (Find(in_line, "PROT"))
125                {
126                    InitNASeq(this_elem, PROTEIN);
127                }
128                else
129                    InitNASeq(this_elem, DNA);
130
131                strncpy_terminate(this_elem->short_name, fields[1], SIZE_SHORT_NAME);
132                AsciiTime(&(this_elem->t_stamp.origin), fields[n-1]);
133                this_elem->attr = DEFAULT_X_ATTR;
134
135                if (Find(in_line, "Circular"))
136                    this_elem->attr |= IS_CIRCULAR;
137
138                gencomments = NULL;
139                genclen = 0;
140            }
141            else if (Find(in_line, "DEFINITION"))
142                strncpy_terminate(this_elem->description, in_line+12, SIZE_DESCRIPTION);
143
144            else if (Find(in_line, "AUTHOR"))
145                strncpy_terminate(this_elem->authority, in_line+12, SIZE_AUTHORITY);
146
147            else if (Find(in_line, "  ORGANISM"))
148                strncpy_terminate(this_elem->seq_name, in_line+12, SIZE_SEQ_NAME);
149
150            else if (Find(in_line, "ACCESSION"))
151                strncpy_terminate(this_elem->id, in_line+12, SIZE_ID);
152
153            else if (Find(in_line, "ORIGIN"))
154            {
155                done = false;
156                len = 0;
157                for (; !done && fgets(in_line, GBUFSIZ, file) != 0;)
158                {
159                    if (in_line[0] != '/')
160                    {
161                        if (buflen == 0)
162                        {
163                            buflen = GBUFSIZ;
164                            buffer = Calloc(sizeof(char), buflen);
165                        }
166
167                        else if (len+strlen(in_line) >= buflen)
168                        {
169                            buflen += GBUFSIZ;
170                            buffer = Realloc((char*)buffer, sizeof(char)*buflen);
171                            for (size_t j=buflen-GBUFSIZ ; j<buflen; j++) buffer[j] = '\0';
172                        }
173                        // Search for the fist column of data (whitespace-number-whitespace)data
174                        if (start_col == -1)
175                        {
176                            for (start_col=0; in_line[start_col] == ' ' || in_line[start_col] == '\t'; start_col++) ;
177                            for (start_col++; strchr("1234567890", in_line[start_col]) != NULL; start_col++) ;
178                            for (start_col++; in_line[start_col] == ' ' || in_line[start_col] == '\t'; start_col++) ;
179                        }
180                        for (int j=start_col; (c = in_line[j]) != '\0'; j++)
181                        {
182                            if ((c != '\n') && ((j-start_col + 1) % 11 != 0)) {
183                                buffer[len++] = c;
184                            }
185                        }
186                    }
187                    else
188                    {
189                        AppendNA((NA_Base*)buffer, len, &(dataset.element[curelem]));
190                        for (size_t j=0; j<len; j++) buffer[j] = '\0';
191                        len = 0;
192                        done = true;
193                        dataset.element[curelem].comments = gencomments;
194                        dataset.element[curelem].comments_len= genclen - 1;
195                        dataset.element[curelem].comments_maxlen = genclen;
196
197                        gencomments = NULL;
198                        genclen = 0;
199                    }
200                }
201                /* Test if sequence should be converted by the translation table
202                 * If it looks like a protein...
203                 */
204                if (dataset.element[curelem].rmatrix && !IS_REALLY_AA) {
205                    IS_REALLY_AA = CheckType((char*)dataset.element[curelem]. sequence, dataset.element[curelem].seqlen);
206
207                    if (!IS_REALLY_AA)
208                        Ascii2NA((char*)dataset.element[curelem].sequence,
209                                 dataset.element[curelem].seqlen,
210                                 dataset.element[curelem].rmatrix);
211                    else {
212                        // Force the sequence to be AA
213                        dataset.element[curelem].elementtype = PROTEIN;
214                        dataset.element[curelem].rmatrix = NULL;
215                        dataset.element[curelem].tmatrix = NULL;
216                        dataset.element[curelem].col_lut = Default_PROColor_LKUP;
217                    }
218                }
219            }
220            else if (Find(in_line, "ZZZZZ"))
221            {
222                Cfree(gencomments);
223                genclen = 0;
224            }
225            else
226            {
227                if (gencomments == NULL)
228                {
229                    gencomments = strdup(in_line);
230                    genclen = strlen(gencomments)+1;
231                }
232                else
233                {
234                    genclen += strlen(in_line)+1;
235                    gencomments = Realloc((char*)gencomments, genclen * sizeof(char));
236                    strncat(gencomments, in_line, GBUFSIZ);
237                    strncat(gencomments, "\n", GBUFSIZ);
238                }
239            }
240        }
241        Cfree(buffer);
242        fclose(file);
243    }
244    for (size_t j=0; j<dataset.numelements; j++) {
245        dataset.maxlen = std::max(dataset.maxlen,
246                                  dataset.element[j].seqlen+dataset.element[j].offset);
247    }
248    return error;
249}
250
251int WriteGen(NA_Alignment& aln, char *filename, int method)
252{
253    int i;
254    size_t j;
255    int k;
256    NA_Sequence *this_elem;
257
258    FILE *file = fopen(filename, "w");
259    if (file == NULL)
260    {
261        Warning("Cannot open file for output");
262        return (1);
263    }
264
265    for (j=0; j<aln.numelements; j++)
266    {
267        this_elem = &(aln.element[j]);
268        if (method == ALL)
269        {
270            fprintf(file,
271                    "LOCUS       %10s%8d bp    %4s  %10s   %2d%5s%4d\n",
272                    this_elem->short_name, this_elem->seqlen+this_elem->offset,
273                    (this_elem->elementtype == DNA) ? "DNA" :
274                    (this_elem->elementtype == RNA) ? "RNA" :
275                    (this_elem->elementtype == MASK) ? "MASK" :
276                    (this_elem->elementtype == PROTEIN) ? "PROT" : "TEXT",
277                    this_elem->attr & IS_CIRCULAR ? "Circular" : "",
278                    this_elem->t_stamp.origin.dd,
279                    GDEmonth[this_elem->t_stamp.origin.mm-1],
280                    (this_elem->t_stamp.origin.yy>1900) ? this_elem->t_stamp.origin.yy :
281                    this_elem->t_stamp.origin.yy+1900);
282
283            if (this_elem->description[0])
284                fprintf(file, "DEFINITION  %s\n", this_elem->description);
285            if (this_elem->seq_name[0])
286                fprintf(file, "  ORGANISM  %s\n", this_elem->seq_name);
287            if (this_elem->id[0])
288                fprintf(file, " ACCESSION  %s\n", this_elem->id);
289            if (this_elem->authority[0])
290                fprintf(file, "   AUTHORS  %s\n", this_elem->authority);
291            if (this_elem->comments)
292                fprintf(file, "%s\n", this_elem->comments);
293            fprintf(file, "ORIGIN");
294            if (this_elem->tmatrix)
295            {
296                for (i=0, k=0; k<this_elem->seqlen+this_elem->offset; k++)
297                {
298                    if (i%60 == 0)
299                        fprintf(file, "\n%9d", i+1);
300                    if (i%10 == 0)
301                        fprintf(file, " ");
302                    fprintf(file, "%c", this_elem->tmatrix
303                            [getelem(this_elem, k)]);
304                    i++;
305                }
306            }
307            else
308            {
309                for (i=0, k=0; k<this_elem->seqlen+this_elem->offset; k++)
310                {
311                    if (i%60 == 0)
312                        fprintf(file, "\n%9d", i+1);
313                    if (i%10 == 0)
314                        fprintf(file, " ");
315                    fprintf(file, "%c", getelem(this_elem, k));
316                    i++;
317                }
318            }
319            fprintf(file, "\n//\n");
320        }
321    }
322    fclose(file);
323    return (0);
324}
325
326
327void SetTime(void *b)
328{
329    ARB_TIME *a=(ARB_TIME*)b;
330    struct tm *tim;
331    time_t clock;
332
333    clock = time(0);
334    tim = localtime(&clock);
335
336    a->yy = tim->tm_year;
337    a->mm = tim->tm_mon+1;
338    a->dd = tim->tm_mday;
339    a->hr = tim->tm_hour;
340    a->mn = tim->tm_min;
341    a->sc = tim->tm_sec;
342    return;
343}
344
Note: See TracBrowser for help on using the repository browser.