root/trunk/CONVERTALN/gcg.cxx

Revision 8607, 5.6 KB (checked in by westram, 5 weeks ago)

merge from e4fix [8135] [8136] [8137] [8138] [8139] [8140] [8141] [8142] [8143] [8144] [8222]
this revives the reverted patches [8129] [8130] [8131] [8132]

  • fixes
    • some free/delete mismatches
    • wrong definition of ORF objects (Level was no bit value)
    • amino consensus (failed for columns only containing 'C')
  • rename
    • AA_sequence_term -> orf_term
    • ED4_sequence_terminal_basic -> ED4_abstract_sequence_terminal
  • cleaned up hierarchy dumps
  • tweaked is_terminal()/to_terminal()
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#include "genbank.h"
2#include "embl.h"
3#include "macke.h"
4
5static void gcg_doc_out(const char *line, Writer& writer) {
6    // Output non-sequence data(document) of gcg format.
7    int indi, len;
8    int previous_is_dot;
9
10    ca_assert(writer.ok());
11
12    for (indi = 0, len = str0len(line), previous_is_dot = 0; indi < len; indi++) {
13        if (previous_is_dot) {
14            if (line[indi] == '.')
15                writer.out(' ');
16            else
17                previous_is_dot = 0;
18        }
19        writer.out(line[indi]);
20        if (line[indi] == '.')
21            previous_is_dot = 1;
22    }
23}
24
25static int gcg_checksum(const char *Str, int numofstr) {
26    // Calculate gcg_checksum for GCG format.
27    int cksum = 0;
28    int count = 0;
29    for (int indi = 0; indi < numofstr; indi++) {
30        if (!is_gapchar(Str[indi])) {
31            count++;
32            cksum = ((cksum + count * toupper(Str[indi])) % 10000);
33            if (count == 57) count = 0;
34        }
35    }
36    return cksum;
37}
38
39static void gcg_out_origin(const Seq& seq, Writer& write) {
40    // Output sequence data in gcg format.
41    int         indi, indj, indk;
42    const char *sequence = seq.get_seq();
43
44    for (indi = 0, indj = 0, indk = 1; indi < seq.get_len(); indi++) {
45        if (!is_gapchar(sequence[indi])) {
46            if ((indk % 50) == 1) write.outf("%8d  ", indk);
47            write.out(sequence[indi]);
48            indj++;
49            if (indj == 10) {
50                write.out(' ');
51                indj = 0;
52            }
53            if ((indk % 50) == 0) write.out("\n\n");
54            indk++;
55        }
56    }
57    if ((indk % 50) != 1) write.out(" \n");
58}
59
60static void gcg_seq_out(const Seq& seq, Writer& write, const char *key) {
61    // Output sequence data in gcg format.
62    write.outf("\n%s  Length: %d  %s  Type: N  Check: %d  ..\n\n",
63            key,
64            seq.get_len()-seq.count_gaps(),
65            gcg_date(today_date()),
66            gcg_checksum(seq.get_seq(), seq.get_len()));
67    gcg_out_origin(seq, write);
68}
69
70class GcgWriter;
71
72class GcgCommentWriter : public Writer {
73    GcgWriter& gcg_writer;
74
75    char linebuf[LINESIZE];
76    int  used;
77public:
78    GcgCommentWriter(GcgWriter& gcg_writer_)
79        : gcg_writer(gcg_writer_),
80          used(0)
81    {}
82    ~GcgCommentWriter() {
83        ca_assert(used == 0); // trailing \n has not been written
84    }
85
86    bool ok() const { return true; }
87    void throw_write_error() { ca_assert(0); }
88    void out(char ch);
89    const char *name() const { return "comment-writer"; }
90};
91
92class GcgWriter : public FileWriter { // derived from a Noncopyable
93    char *species_name;
94    bool  seq_written;             // if true, any further sequences are ignored
95
96    GcgCommentWriter writer;
97
98public:
99    GcgWriter(const char *outname)
100        : FileWriter(outname),
101          species_name(NULL),
102          seq_written(false),
103          writer(*this)
104    {}
105    ~GcgWriter() {
106        free(species_name);
107        FileWriter::seq_done(seq_written);
108    }
109
110    void set_species_name(const char *next_name) {
111        if (!seq_written) species_name = nulldup(next_name);
112        else warningf(111, "Species '%s' dropped (GCG allows only 1 sequence per file)", next_name);
113    }
114
115    void add_comment(const char *comment) {
116        if (!seq_written) gcg_doc_out(comment, *this);
117    }
118
119    Writer& comment_writer() {
120        ca_assert(!seq_written);
121        return writer;
122    }
123
124    void write_seq_data(const Seq& seq) {
125        if (!seq_written) {
126            ca_assert(species_name); // you have to call set_species_name() before!
127            gcg_seq_out(seq, *this, species_name);
128            seq_written = true;
129        }
130    }
131};
132
133void GcgCommentWriter::out(char ch) {
134    linebuf[used++] = ch;
135    ca_assert(used<LINESIZE);
136    if (ch == '\n') {
137        linebuf[used] = 0;
138        gcg_writer.add_comment(linebuf);
139        used = 0;
140    }
141}
142
143static void macke_to_gcg(const char *inf, const char *outf) {
144    MackeReader reader(inf);
145    GcgWriter   out(outf);
146
147    Seq seq;
148    if (!reader.read_one_entry(seq)) return;
149
150    Macke& macke = dynamic_cast<Macke&>(reader.get_data());
151    out.set_species_name(macke.get_id());
152    macke_seq_info_out(macke, out);
153    out.write_seq_data(seq);
154
155    reader.ignore_rest_of_file();
156}
157
158static void genbank_to_gcg(const char *inf, const char *outf) {
159    FormatReaderPtr reader = FormatReader::create(FormattedFile(inf, GENBANK));
160    GcgWriter       write(outf);
161
162    GenBank gbk;
163    Seq     seq;
164
165    {
166        GenbankReader& greader = dynamic_cast<GenbankReader&>(*reader);
167        if (!GenbankParser(gbk, seq, greader).parse_entry()) return;
168    }
169
170    genbank_out_header(gbk, seq, write.comment_writer());
171    genbank_out_base_count(seq, write.comment_writer());
172    write.out("ORIGIN\n");
173    write.set_species_name(gbk.get_id());
174    write.write_seq_data(seq);
175
176    reader->ignore_rest_of_file();
177}
178
179static void embl_to_gcg(const char *inf, const char *outf) {
180    EmblSwissprotReader reader(inf);
181    GcgWriter           write(outf);
182
183    Embl embl;
184    Seq  seq;
185
186    if (!EmblParser(embl, seq, reader).parse_entry()) return;
187
188    embl_out_header(embl, seq, write);
189    write.set_species_name(embl.get_id());
190    write.write_seq_data(seq);
191
192    reader.ignore_rest_of_file();
193}
194
195void to_gcg(const FormattedFile& in, const char *outf) {
196    // Convert from whatever to GCG format
197    // @@@ use InputFormat ?
198
199    switch (in.type()) {
200        case MACKE:     macke_to_gcg(in.name(), outf); break;
201        case GENBANK:   genbank_to_gcg(in.name(), outf); break;
202        case EMBL:
203        case SWISSPROT: embl_to_gcg(in.name(), outf); break;
204        default:
205            throw_conversion_not_supported(in.type(), GCG);
206            break;
207    }
208}
209
Note: See TracBrowser for help on using the browser.