source: tags/ms_r18q1/RNA3D/RNA3D_StructureData.cxx

Last change on this file was 16766, checked in by westram, 4 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.5 KB
Line 
1#include "RNA3D_GlobalHeader.hxx"
2#include "RNA3D_Global.hxx"
3#include "RNA3D_OpenGLGraphics.hxx"
4#include "RNA3D_Graphics.hxx"
5#include "RNA3D_StructureData.hxx"
6
7#include <cerrno>
8#include <string>
9#include <iostream>
10#include <fstream>
11
12#include <arbdbt.h>
13#include <aw_msg.hxx>
14#include <aw_awars.hxx>
15#include <aw_root.hxx>
16#include <ed4_extern.hxx>
17#include <ed4_plugins.hxx>
18#include <BI_basepos.hxx>
19
20
21#define COLORLINK (ED4_G_SBACK_0 - RNA3D_GC_SBACK_0)  // to link to the colors defined in primary editor ed4_defs.hxx
22#define SAICOLORS (ED4_G_CBACK_0 - RNA3D_GC_CBACK_0)
23
24using namespace std;
25
26static Struct3Dinfo   *start3D                 = NULp;
27static Struct2Dinfo   *start2D                 = NULp;
28static Struct2Dplus3D *start2D3D               = NULp;
29static HelixNrInfo    *start                   = NULp;
30static CurrSpecies    *startSp                 = NULp;
31static bool            bOldSpeciesDataExists   = false;
32static Insertions     *startIns                = NULp;
33static bool            bOldInsertionDataExists = false;
34static bool            bStartPosStored         = false;
35static bool            bEndPosStored           = false;
36
37static char *find_data_file(const char *name) {
38    char *fname = GB_lib_file(false, "rna3d/", name);
39    if (!fname) throw string("file not found: ")+name;
40    return fname;
41}
42
43static void throw_IO_error(const char *filename) __ATTR__NORETURN;
44static void throw_IO_error(const char *filename) {
45    string error = string("IO-Error: ")+strerror(errno)+" ('"+filename+"')";
46    throw error;
47}
48
49GBDATA *Structure3D::gb_main = NULp;
50
51Structure3D::Structure3D(ED4_plugin_host& host_)
52    : strCen(new Vector3(0.0, 0.0, 0.0)),
53      iInterval(25),
54      iMapSAI(0),
55      iMapSearch(0),
56      iMapEnable(0),
57      iStartPos(0),
58      iEndPos(0),
59      iEColiStartPos(0),
60      iEColiEndPos(0),
61      iTotalSubs(0),
62      iTotalDels(0),
63      iTotalIns(0),
64      LSU_molID(3), // default molecule to load in case of 23S rRNA : 1C2W
65      HelixBase(500),
66      EColiRef(NULp),
67      Host(host_),
68      GRAPHICS(new OpenGLGraphics)
69{}
70
71Structure3D::~Structure3D() {
72    delete GRAPHICS;
73    delete strCen;
74}
75
76
77void Structure3D::StoreCoordinates(float x, float y, float z, char base, unsigned int pos) {
78    Struct3Dinfo *data, *temp;
79    data = new Struct3Dinfo;
80    data->x    = x;
81    data->y    = y;
82    data->z    = z;
83    data->base = base;
84    data->pos  = pos;
85    data->next = NULp;
86
87    if (!start3D)
88        start3D = data;
89    else {
90        temp = start3D;
91        // We know this is not NULp - list not empty!
92        while (temp->next) {
93            temp = temp->next;  // Move to next link in chain
94        }
95        temp->next = data;
96    }
97}
98
99// ------------------Selecting rRNA type--------------------------//
100
101int Structure3D::FindTypeOfRNA() { // @@@ should better return an enum type
102    int rnaType = 0;
103    GB_push_transaction(gb_main);  // opening a transaction
104
105    GBDATA *gbTemplate = GBT_find_SAI(gb_main, "ECOLI");
106    if (!gbTemplate) {
107        aw_message("SAI:ECOLI not found");
108    }
109    else {
110        char   *ali_name          = GBT_get_default_alignment(gb_main);
111        GBDATA *gbAlignment       = GB_entry(gbTemplate, ali_name);
112        GBDATA *gbTemplateSeqData = gbAlignment ? GB_entry(gbAlignment, "data") : NULp;
113
114        if (!gbTemplateSeqData) {
115          aw_message("Could not find species in the database!");
116        }
117        else {
118            const char *pTemplateSeqData  = GB_read_char_pntr(gbTemplateSeqData);
119            int iSeqLen = strlen(pTemplateSeqData);
120            int iBaseCount = 0;
121            for (int i = 0; i<iSeqLen; i++) {
122                if ((pTemplateSeqData[i] != '.') && (pTemplateSeqData[i] != '-')) {
123                    iBaseCount++;
124                }
125            }
126
127            if (iBaseCount > 2000)      rnaType = LSU_23S;
128            else if (iBaseCount > 300)  rnaType = SSU_16S;
129            else                        rnaType = LSU_5S;
130        }
131    }
132    GB_pop_transaction(gb_main);
133
134    return rnaType;
135}
136
137// ----------Delete old molecule data-------------------------------//
138// Struct2Dinfo, Struct2Dplus3D, Struct3Dinfo, HelixNrInfo, Insertions
139
140void Structure3D::DeleteOldMoleculeData() {
141    // Struct2Dinfo -> start2D
142    {
143        Struct2Dinfo *tmp, *data;
144        for (data = start2D; data; data = tmp) {
145            tmp = data->next;
146            delete data;
147        }
148        start2D = NULp;
149    }
150
151    // Struct2Dplus3D ->start2D3D
152    {
153        Struct2Dplus3D *tmp, *data;
154        for (data = start2D3D; data; data = tmp) {
155            tmp = data->next;
156            delete data;
157        }
158        start2D3D = NULp;
159    }
160
161    // Struct3Dinfo ->start3D
162    {
163        Struct3Dinfo *tmp, *data;
164        for (data = start3D; data; data = tmp) {
165            tmp = data->next;
166            delete data;
167        }
168        start3D = NULp;
169    }
170
171    // HelixNrInfo -> start
172    {
173        HelixNrInfo *tmp, *data;
174        for (data = start; data; data = tmp) {
175            tmp = data->next;
176            delete data;
177        }
178        start = NULp;
179    }
180
181    // Delete insertions data
182    DeleteOldInsertionData();
183
184    // Delete mapped species data
185    DeleteOldSpeciesData ();
186}
187
188// =========== Reading 3D Coordinates from PDB file ====================//
189
190static char globalComment[1000];
191
192void Structure3D::ReadCoOrdinateFile() {
193    static char *DataFile = NULp;
194    int          rnaType  = FindTypeOfRNA();
195
196    RNA3D->bDisplayComments = true; // displaying comments in main window
197
198    switch (rnaType) {
199    case LSU_23S:
200        switch (LSU_molID) {
201        case _1PNU:
202            DataFile = find_data_file("Ecoli_1PNU_23S_rRNA.pdb");
203            sprintf(globalComment, "The 3D molecule rendered from PDB entry : 1PNU at 8.7 Angstrom.");
204            break;
205        case _1VOR:
206            DataFile = find_data_file("Ecoli_1VOR_23S_rRNA.pdb");
207            sprintf(globalComment, "The 3D molecule is rendered from PDB entry [1VOR] with 11.5 Angstrom resolution.");
208            break;
209        case _1C2W:
210            DataFile = find_data_file("Ecoli_1C2W_23S_rRNA.pdb");
211            sprintf(globalComment, "The 3D molecule is rendered from PDB entry [1C2W] with 7.5 Angstrom resolution.");
212            break;
213        }
214        break;
215    case LSU_5S:
216        DataFile = find_data_file("Ecoli_1C2X_5S_rRNA.pdb");
217        sprintf(globalComment, "The 3D molecule is rendered from PDB entry [1C2X] with 7.5 Angstrom resolution.");
218        break;
219    case SSU_16S:
220        DataFile = find_data_file("Ecoli_1M5G_16S_rRNA.pdb");
221        sprintf(globalComment, "The 3D molecule is rendered from PDB entry [1M5G] with 11.5 Angstrom resolution.");
222        break;
223    }
224
225    ifstream readData;
226    readData.open(DataFile, ios::in);
227    if (!readData.is_open()) {
228        throw_IO_error(DataFile);
229    }
230
231    int          cntr                 = 0;
232    unsigned int last3Dpos            = 0;
233    static bool  bEColiStartPosStored = false;
234
235    while (!readData.eof()) {
236        char buf[256];
237        readData.getline(buf, 100);
238
239        string line(buf);
240        if ((line.find("ATOM") != string::npos) || (line.find("HETATM") != string::npos)) {
241            string atom = line.substr(77, 2);
242            if (atom.find("P") != string::npos) {
243                char     base = line[19];
244                unsigned pos  = atoi(line.substr(22, 4).c_str());
245
246                // special filter for 23S rRNA structure (IVOR/IPNU)
247                // IVOR/IPNU contains artifacts and are not mentioned in any of the
248                // remarks of PDB file
249                if (last3Dpos != pos && !(pos >= 3093)) {
250                    float X = atof(line.substr(31, 8).c_str());
251                    float Y = atof(line.substr(39, 8).c_str());
252                    float Z = atof(line.substr(47, 8).c_str());
253
254                    StoreCoordinates(X, Y, Z, base, pos);
255                    last3Dpos  = pos;
256                    strCen->x += X; strCen->y += Y; strCen->z += Z;
257                    cntr++;
258                }
259                if (!bEColiStartPosStored) {
260                    iEColiStartPos = pos;
261                    bEColiStartPosStored = true;
262                }
263                iEColiEndPos = pos;
264            }
265        }
266    }
267
268    cout<<"--------------------------------------------------"<<endl
269        <<globalComment<<endl
270        <<"Structure contains "<<iEColiEndPos-iEColiStartPos<<"( "<<iEColiStartPos<<" - "
271        <<iEColiEndPos<<" ) positions."<<endl
272        <<"---------------------------------------------------"<<endl;
273
274    strCen->x = strCen->x/cntr;
275    strCen->y = strCen->y/cntr;
276    strCen->z = strCen->z/cntr;
277
278    readData.close();
279    free(DataFile);
280}
281
282
283void Structure3D::Store2Dinfo(char *info, int pos, int helixNr) {
284    Struct2Dinfo *data = new Struct2Dinfo;
285
286    data->base    = info[0];
287    data->mask    = info[1];
288    data->code    = info[2];
289    data->pos     = pos;
290    data->helixNr = helixNr;
291
292    data->next = NULp;
293
294    if (!start2D) {
295        start2D = data;
296    }
297    else {
298        Struct2Dinfo *temp = start2D;
299        // We know this is not NULp - list not empty!
300        while (temp->next) {
301            temp = temp->next;  // Move to next link in chain
302        }
303        temp->next = data;
304    }
305}
306
307// =========== Reading Secondary Structure Data from Ecoli Secondary Structure Mask file ====================//
308
309void Structure3D::GetSecondaryStructureInfo() {
310    static char *DataFile = NULp;
311    int          rnaType  = FindTypeOfRNA();
312
313    switch (rnaType) {
314        case LSU_23S:
315            DataFile = find_data_file("SecondaryStructureModel_23SrRNA.data");
316            break;
317        case LSU_5S:
318            DataFile = find_data_file("SecondaryStructureModel_5SrRNA.data");
319            break;
320        case SSU_16S:
321            DataFile = find_data_file("SecondaryStructureModel_16SrRNA.data");
322            break;
323    }
324
325    char  buf[256];
326
327    int  pos         = 0;
328    int  helixNr     = 0;
329    int  lastHelixNr = 0;
330    char info[4];
331    info[3]          = '\0';
332    bool insideHelix = false;
333    bool skip        = false;
334
335    ifstream readData;
336    readData.open(DataFile, ios::in);
337    if (!readData.is_open()) {
338        throw_IO_error(DataFile);
339    }
340
341    while (!readData.eof()) {
342        readData.getline(buf, 100);
343        char *tmp;
344        tmp = strtok(buf, " ");
345        for (int i = 0; tmp; tmp = strtok(NULp, " "), i++) {
346            switch (i) {
347                case 0: pos = atoi(tmp);      break;
348                case 1: info[0] = tmp[0];     break;
349                case 2: info[1] = tmp[0];     break;
350                case 3: info[2] = tmp[0];     break;
351                case 4: helixNr = atoi(tmp); break;
352            }
353        }
354
355        bool is_SE = (info[2] == 'S') || (info[2] == 'E'); // 'S' = helix start 'E' = helix end
356        if (is_SE) {
357            if (helixNr>0) lastHelixNr = helixNr;
358
359            if (!insideHelix) insideHelix = true;
360            else {
361                Store2Dinfo(info, pos, lastHelixNr);
362                skip        = true;
363                insideHelix = false;
364            }
365        }
366        if (insideHelix) Store2Dinfo(info, pos, lastHelixNr);
367        else {
368            if (skip) skip = false;
369            else Store2Dinfo(info, pos, 0);
370        }
371    }
372    readData.close();
373    free(DataFile);
374}
375
376void Structure3D::Store2D3Dinfo(Struct2Dinfo *s2D, Struct3Dinfo *s3D) {
377    Struct2Dplus3D *data, *temp;
378    data = new Struct2Dplus3D;
379    data->base    = s2D->base;
380    data->mask    = s2D->mask;
381    data->code    = s2D->code;
382    data->pos     = s2D->pos;
383    data->helixNr = s2D->helixNr;
384    data->x       = s3D->x;
385    data->y       = s3D->y;
386    data->z       = s3D->z;
387    data->next    = NULp;
388
389    if (!start2D3D)
390        start2D3D = data;
391    else {
392        temp = start2D3D;
393        while (temp->next) {
394            temp = temp->next;
395        }
396        temp->next = data;
397    }
398}
399
400// =========== Combining Secondary Structure Data with 3D Coordinates =======================//
401
402void Structure3D::Combine2Dand3DstructureInfo() {
403    Struct3Dinfo *temp3D;
404    Struct2Dinfo *temp2D;
405    int cntr = 0;
406
407    cout<<"Missing Base Positions : "<<endl;
408
409    temp3D = start3D;
410    temp2D = start2D;
411    while (temp3D && temp2D) {
412        if (temp2D->pos == temp3D->pos) {
413            Store2D3Dinfo(temp2D, temp3D);
414        }
415        else {
416            while (temp2D->pos != temp3D->pos) {
417                cout<<temp2D->pos<<", "; // missing base positions
418                cntr++;
419                temp2D = temp2D->next;
420            }
421            Store2D3Dinfo(temp2D, temp3D);
422        }
423
424        temp3D = temp3D->next;
425        temp2D = temp2D->next;
426    }
427    cout<<endl<<"Total No. of bases missing = "<<cntr<<endl;
428
429    // printing comments in the main window
430    {
431        RNA3D->bDisplayComments = true;
432        char buf[256];
433        sprintf(buf, "Total No. of bases missing = %d. See the console messages for the actual missing base positions.", cntr);
434        strcat(globalComment, buf);
435    }
436}
437
438void Structure3D::PointsToQuads(float x, float y, float z) {
439
440    if (RNA3D->bPointSpritesSupported) {
441        glVertex3f(x, y, z);
442    }
443    else {
444        glBegin(GL_QUADS);
445        // Front Face
446        glTexCoord2f(0, 0); glVertex3f(x - 1, y + 1, z + 1);
447        glTexCoord2f(1, 0); glVertex3f(x + 1, y + 1, z + 1);
448        glTexCoord2f(1, 1); glVertex3f(x + 1, y - 1, z + 1);
449        glTexCoord2f(0, 1); glVertex3f(x - 1, y - 1, z + 1);
450
451        // Back Face
452        glTexCoord2f(0, 0); glVertex3f(x + 1, y + 1, z - 1);
453        glTexCoord2f(1, 0); glVertex3f(x - 1, y + 1, z - 1);
454        glTexCoord2f(1, 1); glVertex3f(x - 1, y - 1, z - 1);
455        glTexCoord2f(0, 1); glVertex3f(x + 1, y - 1, z - 1);
456
457        // Top Face
458        glTexCoord2f(0, 0); glVertex3f(x + 1, y + 1, z + 1);
459        glTexCoord2f(1, 0); glVertex3f(x - 1, y + 1, z + 1);
460        glTexCoord2f(1, 1); glVertex3f(x - 1, y + 1, z - 1);
461        glTexCoord2f(0, 1); glVertex3f(x + 1, y + 1, z - 1);
462
463        // Bottom Face
464        glTexCoord2f(0, 0); glVertex3f(x + 1, y - 1, z - 1);
465        glTexCoord2f(1, 0); glVertex3f(x - 1, y - 1, z - 1);
466        glTexCoord2f(1, 1); glVertex3f(x - 1, y - 1, z + 1);
467        glTexCoord2f(0, 1); glVertex3f(x + 1, y - 1, z + 1);
468
469        // Left Face
470        glTexCoord2f(0, 0); glVertex3f(x + 1, y + 1, z + 1);
471        glTexCoord2f(1, 0); glVertex3f(x + 1, y + 1, z - 1);
472        glTexCoord2f(1, 1); glVertex3f(x + 1, y - 1, z - 1);
473        glTexCoord2f(0, 1); glVertex3f(x + 1, y - 1, z + 1);
474
475        // Right Face
476        glTexCoord2f(0, 0); glVertex3f(x - 1, y + 1, z - 1);
477        glTexCoord2f(1, 0); glVertex3f(x - 1, y + 1, z + 1);
478        glTexCoord2f(1, 1); glVertex3f(x - 1, y - 1, z + 1);
479        glTexCoord2f(0, 1); glVertex3f(x - 1, y - 1, z - 1);
480
481        glEnd();
482    }
483}
484
485void Structure3D::PositionsToCoordinatesDispList(int listID, int *pos, int len) {
486    Struct2Dplus3D *t;
487    int tmpPos = 0;
488
489    glNewList(listID, GL_COMPILE);
490    {
491        if (RNA3D->bPointSpritesSupported) {
492            glBegin(GL_POINTS);
493        }
494        for (int i = 0; i < len; i++) {
495            tmpPos = pos[i];
496            t = start2D3D;
497            while (t) {
498                if (t->pos == tmpPos) {
499                    PointsToQuads(t->x, t->y, t->z);
500                    break;
501                }
502                t = t->next;
503            }
504        }
505        if (RNA3D->bPointSpritesSupported) {
506            glEnd();
507        }
508    }
509    glEndList();
510}
511
512void Structure3D::GenerateSecStructureNonHelixRegions() {
513    Struct2Dplus3D *t;
514    const int MAX_BASE = 1000;
515    int baseA[MAX_BASE], baseG[MAX_BASE], baseC[MAX_BASE], baseU[MAX_BASE];
516    int a, g, c, u; a=g=c=u=0;
517
518    {
519        t = start2D3D;
520        while (t) {
521            if (t->helixNr == 0) {
522                switch (t->base) {
523                    case 'A': baseA[a++] = t->pos; break;
524                    case 'G': baseG[g++] = t->pos; break;
525                    case 'C': baseC[c++] = t->pos; break;
526                    case 'U': baseU[u++] = t->pos; break;
527                }
528            }
529            t = t->next;
530        }
531    }
532
533    PositionsToCoordinatesDispList(NON_HELIX_A, baseA, a);
534    PositionsToCoordinatesDispList(NON_HELIX_G, baseG, g);
535    PositionsToCoordinatesDispList(NON_HELIX_C, baseC, c);
536    PositionsToCoordinatesDispList(NON_HELIX_U, baseU, u);
537}
538
539void Structure3D::GenerateSecStructureHelixRegions() {
540    Struct2Dplus3D *t;
541    const int MAX_BASE = 1000;
542    int baseA[MAX_BASE], baseG[MAX_BASE], baseC[MAX_BASE], baseU[MAX_BASE];
543    int a, g, c, u; a=g=c=u=0;
544
545    {
546        t = start2D3D;
547        while (t) {
548            if (t->helixNr > 0) {
549                if ((t->mask == '[') || (t->mask == ']') || (t->mask == '<') || (t->mask == '>')) {
550                    switch (t->base) {
551                        case 'A': baseA[a++] = t->pos; break;
552                        case 'G': baseG[g++] = t->pos; break;
553                        case 'C': baseC[c++] = t->pos; break;
554                        case 'U': baseU[u++] = t->pos; break;
555                    }
556                }
557            }
558            t = t->next;
559        }
560    }
561
562    PositionsToCoordinatesDispList(HELIX_A, baseA, a);
563    PositionsToCoordinatesDispList(HELIX_G, baseG, g);
564    PositionsToCoordinatesDispList(HELIX_C, baseC, c);
565    PositionsToCoordinatesDispList(HELIX_U, baseU, u);
566}
567
568void Structure3D::GenerateSecStructureUnpairedHelixRegions() {
569    Struct2Dplus3D *t;
570    const int MAX_BASE = 500;
571    int baseA[MAX_BASE], baseG[MAX_BASE], baseC[MAX_BASE], baseU[MAX_BASE];
572    int a, g, c, u; a=g=c=u=0;
573
574    {
575        t = start2D3D;
576        while (t) {
577            if (t->helixNr > 0) {
578                if (t->mask == '.') {
579                    switch (t->base) {
580                        case 'A': baseA[a++] = t->pos; break;
581                        case 'G': baseG[g++] = t->pos; break;
582                        case 'C': baseC[c++] = t->pos; break;
583                        case 'U': baseU[u++] = t->pos; break;
584                    }
585                }
586            }
587            t = t->next;
588        }
589    }
590
591    PositionsToCoordinatesDispList(UNPAIRED_HELIX_A, baseA, a);
592    PositionsToCoordinatesDispList(UNPAIRED_HELIX_G, baseG, g);
593    PositionsToCoordinatesDispList(UNPAIRED_HELIX_C, baseC, c);
594    PositionsToCoordinatesDispList(UNPAIRED_HELIX_U, baseU, u);
595}
596
597// ==============================================================================
598// Tertiary Interactions of 16S ribosomal RNA model of E.coli.
599// Reference : http://www.rna.icmb.utexas.edu/
600// Year of Last Update : 2001.
601// Pseudoknots and Triple Base pairs are extracted and displayed in
602// the 3D model.
603// ==============================================================================
604
605void Structure3D::GenerateTertiaryInteractionsDispLists() {
606    Struct2Dplus3D *t;
607    static char    *DataFile = NULp;
608    int             rnaType  = FindTypeOfRNA();
609
610    switch (rnaType) {
611        case LSU_23S:
612            DataFile = find_data_file("TertiaryInteractions_23SrRNA.data");
613            break;
614        case LSU_5S:
615            cout<<"There are no tertiary interactions observed in 5S rRNA model!"<<endl;
616            return;
617        case SSU_16S:
618            DataFile = find_data_file("TertiaryInteractions_16SrRNA.data");
619            break;
620    }
621
622    cout<<"Tertiary Interactions are fetched from comparative RNA website [http://www.rna.icmb.utexas.edu]."<<endl
623        <<"The same are located in the file \""<<DataFile<<"\"."<<endl;
624
625    char  buf[256];
626
627    ifstream readData;
628    readData.open(DataFile, ios::in);
629    if (!readData.is_open()) {
630        throw_IO_error(DataFile);
631    }
632
633    int K[100];
634    int R[100];
635    int k, r; k = r = 0;
636
637    while (!readData.eof()) {
638        readData.getline(buf, 100);
639        char *tmp;
640        tmp = strtok(buf, " ");
641        if (tmp) {
642            if (strcmp(tmp, "KNOT") == 0) {
643                tmp = strtok(NULp, ":");
644                while (tmp) {
645                    K[k++] = atoi(tmp);
646                    tmp = strtok(NULp, ":");
647                }
648            }
649            else if (strcmp(tmp, "TRIPLE") == 0) {
650                tmp = strtok(NULp, ":");
651                while (tmp) {
652                    R[r++] = atoi(tmp);
653                    tmp = strtok(NULp, ":");
654                }
655            }
656        }
657    }
658    readData.close();
659
660    glNewList(ECOLI_TERTIARY_INTRACTION_PSEUDOKNOTS, GL_COMPILE);
661    {
662        for (int i = 0; i < k;) {
663            glBegin(GL_LINES);
664            for (int j = 0; j < 2; j++) {
665                t = start2D3D;
666                while (t) {
667                    if (t->pos == K[i]) {
668                        glVertex3f(t->x, t->y, t->z);
669                        i++;
670                        break;
671                    }
672                    t = t->next;
673                }
674            }
675            glEnd();
676        }
677    }
678    glEndList();
679
680    glNewList(ECOLI_TERTIARY_INTRACTION_TRIPLE_BASES, GL_COMPILE);
681    {
682        for (int i = 0; i < r;) {
683            glBegin(GL_LINE_STRIP);
684            for (int j = 0; j < 3; j++) {
685                t = start2D3D;
686                while (t) {
687                    if (t->pos == R[i]) {
688                        glVertex3f(t->x, t->y, t->z);
689                        i++;
690                        break;
691                    }
692                    t = t->next;
693                }
694            }
695            glEnd();
696        }
697    }
698    glEndList();
699    free(DataFile);
700}
701
702// ==============================================================================
703
704void Structure3D::StoreHelixNrInfo(float x, float y, float z, int helixNr) {
705    HelixNrInfo *data, *temp;
706    data = new HelixNrInfo;
707    data->helixNr = helixNr;
708    data->x       = x;
709    data->y       = y;
710    data->z       = z;
711    data->next    = NULp;
712
713    if (!start)
714        start = data;
715    else {
716        temp = start;
717        while (temp->next) {
718            temp = temp->next;
719        }
720        temp->next = data;
721    }
722}
723
724void Structure3D::GenerateHelixDispLists(int HELIX_NR_ID, int HELIX_NR) {
725    Struct2Dinfo *temp2D;
726    Struct2Dplus3D *temp2D3D;
727
728    const int MAX = 500;
729
730    int thisStrandPos[MAX], otherStrandPos[MAX];
731    int i, j; i = j = 0;
732    {
733        temp2D = start2D;
734        while (temp2D) {
735            if (temp2D->helixNr == HELIX_NR) {
736                if ((temp2D->mask == '[') || (temp2D->mask == '<'))
737                    thisStrandPos[i++]  = temp2D->pos;
738                if ((temp2D->mask == ']') || (temp2D->mask == '>'))
739                    otherStrandPos[j++] = temp2D->pos;
740            }
741            temp2D = temp2D->next;
742        }
743    }
744
745    int tempPos = 0;
746    float x1, x2, y1, y2, z1, z2; x1=x2=y1=y2=z1=z2=0.0;
747
748    bool bThisStrand, bOtherStrand;
749    bThisStrand = bOtherStrand = false;
750
751    bool bMissingBasePair = false;
752
753    glNewList(HELIX_NR_ID, GL_COMPILE);
754    {
755        for (int k = 0, l = j-1; k < i && l >= 0; k++, l--) {
756            tempPos = thisStrandPos[k];
757            {
758                temp2D3D = start2D3D;
759                while (temp2D3D) {
760                    if (temp2D3D->pos == tempPos) {
761                        bThisStrand = true;
762                        x1=temp2D3D->x; y1=temp2D3D->y; z1=temp2D3D->z;
763                        break;
764                    }
765                    temp2D3D = temp2D3D->next;
766                }
767            }
768            tempPos = otherStrandPos[l];
769            {
770                temp2D3D = start2D3D;
771                while (temp2D3D) {
772                    if (temp2D3D->pos == tempPos) {
773                        bOtherStrand = true;
774                        x2=temp2D3D->x; y2=temp2D3D->y; z2=temp2D3D->z;
775                        break;
776                    }
777                    temp2D3D = temp2D3D->next;
778                }
779            }
780            // if bases present in both the strands then draw a bond
781            // and store the helix number information
782            if (bThisStrand && bOtherStrand) {
783                glVertex3f(x1, y1, z1);
784                glVertex3f(x2, y2, z2);
785
786                x1 = (x1+x2)/2; y1 = (y1+y2)/2; z1 = (z1+z2)/2;
787                StoreHelixNrInfo(x1, y1, z1, HELIX_NR);
788
789                bThisStrand = bOtherStrand = false;
790            }
791            else {
792                bMissingBasePair = true;
793                cout<<thisStrandPos[k]<<"-"<<tempPos<<"("<<HELIX_NR_ID<<"), ";
794            }
795        }
796        if (bMissingBasePair) cout<<endl;
797    }
798    glEndList();
799}
800
801void Structure3D::GenerateHelixNrDispList(int startHx, int endHx) {
802    HelixNrInfo *t;
803
804    glNewList(HELIX_NUMBERS, GL_COMPILE);
805    {
806        char POS[50];
807        for (int i = startHx; i <= endHx; i++) {
808            t = start;
809            while (t) {
810                if (t->helixNr == i) {
811                    sprintf(POS, "%d", t->helixNr);
812                    GRAPHICS->PrintString(t->x, t->y, t->z, POS, GLUT_BITMAP_8_BY_13);
813                }
814                t = t->next;
815            }
816        }
817    }
818    glEndList();
819
820    glNewList(HELIX_NUMBERS_POINTS, GL_COMPILE);
821    {
822        if (RNA3D->bPointSpritesSupported) {
823            glBegin(GL_POINTS);
824        }
825        for (int i = startHx; i <= endHx; i++) {
826            t = start;
827            while (t) {
828                if (t->helixNr == i) {
829                    PointsToQuads(t->x, t->y, t->z);
830                }
831                t = t->next;
832            }
833        }
834        if (RNA3D->bPointSpritesSupported) {
835            glEnd();
836        }
837    }
838    glEndList();
839}
840
841void Structure3D::GenerateDisplayLists() {
842
843    GenerateMoleculeSkeleton();
844    ComputeBasePositions();
845
846    int rnaType = FindTypeOfRNA();
847    switch (rnaType) {
848    case SSU_16S:
849        if (glIsList(HelixBase) != GL_TRUE) {
850            for (int i = 1; i <= 50; i++) {
851                GenerateHelixDispLists(HelixBase+i, i);
852            }
853        }
854        break;
855    case LSU_23S:
856        cout<<"=========================================================="<<endl
857            <<"Missing base pairs (bracket number indicates helix number) "<<endl;
858        if (glIsList(HelixBase) != GL_TRUE) {
859            for (int i = 1; i <= 101; i++) {
860                GenerateHelixDispLists(HelixBase+i, i);
861            }
862        }
863        cout<<"=========================================================="<<endl;
864        break;
865    case LSU_5S:
866        if (glIsList(HelixBase) != GL_TRUE) {
867            for (int i = 1; i <= 5; i++) {
868                GenerateHelixDispLists(HelixBase+i, i);
869            }
870        }
871        break;
872    }
873
874    GenerateSecStructureHelixRegions();
875    GenerateSecStructureNonHelixRegions();
876    GenerateSecStructureUnpairedHelixRegions();
877
878    GenerateTertiaryInteractionsDispLists();
879}
880
881void Structure3D::GenerateMoleculeSkeleton() {
882    Struct2Dplus3D *t;
883
884    glNewList(STRUCTURE_BACKBONE, GL_COMPILE);
885    {
886        glBegin(GL_LINE_STRIP);
887        t = start2D3D;
888        while (t) {
889            glVertex3f(t->x, t->y, t->z);
890            t = t->next;
891        }
892        glEnd();
893    }
894    glEndList();
895
896    glNewList(STRUCTURE_BACKBONE_CLR, GL_COMPILE);
897    {
898        glBegin(GL_LINE_STRIP);
899        t = start2D3D;
900        while (t) {
901            if (t->helixNr > 0) {
902                if ((t->mask == '[') || (t->mask == ']') || (t->mask == '<') || (t->mask == '>')) {
903                    GRAPHICS->SetColor(RNA3D_GC_BASES_HELIX);
904                    glVertex3f(t->x, t->y, t->z);
905                }
906                if (t->mask == '.') {
907                    GRAPHICS->SetColor(RNA3D_GC_BASES_UNPAIRED_HELIX);
908                    glVertex3f(t->x, t->y, t->z);
909                }
910            }
911            if (t->helixNr == 0) {
912                GRAPHICS->SetColor(RNA3D_GC_BASES_NON_HELIX);
913                glVertex3f(t->x, t->y, t->z);
914            }
915            t = t->next;
916        }
917        glEnd();
918    }
919    glEndList();
920}
921
922void Structure3D::GenerateCursorPositionDispList(long pos) {
923    Struct3Dinfo *temp;
924
925    glNewList(ECOLI_CURSOR_POSITION, GL_COMPILE);
926    {
927        glBegin(GL_POINTS);
928        temp = start3D;
929        while (temp) {
930            if (temp->pos == pos) {
931#ifdef DEBUG
932                cout<<"Cursor Position : "<<pos<<endl;
933#endif
934                glVertex3f(temp->x, temp->y, temp->z);
935                break;
936            }
937            temp = temp->next;
938        }
939        glEnd();
940    }
941    glEndList();
942}
943
944void Structure3D::ComputeBasePositions() {
945    Struct3Dinfo *temp;
946
947    char POS[50];
948    float spacer = 1.5;
949    int posSkip = iInterval;
950    if (posSkip <= 0) posSkip = 25; // default interval value
951
952    glNewList(STRUCTURE_POS, GL_COMPILE);
953    {
954        temp = start3D;
955        while (temp) {
956            if (temp->pos%posSkip == 0) {
957                sprintf(POS, "%d", temp->pos);
958                GRAPHICS->PrintString(temp->x-spacer, temp->y, temp->z-spacer, POS, GLUT_BITMAP_HELVETICA_10);
959            }
960            temp = temp->next;
961        }
962    }
963    glEndList();
964
965    glNewList(STRUCTURE_POS_ANCHOR, GL_COMPILE);
966    {
967        glLineWidth(0.5);
968        glBegin(GL_LINES);
969        temp = start3D;
970        while (temp) {
971            if (temp->pos%posSkip == 0) {
972                glVertex3f(temp->x, temp->y, temp->z);
973                glVertex3f(temp->x-spacer, temp->y, temp->z-spacer);
974            }
975            temp = temp->next;
976        }
977        glEnd();
978    }
979    glEndList();
980}
981
982void Structure3D::PrepareSecondaryStructureData() {
983//     const char
984//         outFile[]      = "/nfshome/yadhu/DataBase/Struct3D/model16SrRNA/test.data",
985//         EcoliFile[]    = "/nfshome/yadhu/DataBase/Struct3D/model16SrRNA/ECOLI_GAPS",
986//         HelixNrFile[]  = "/nfshome/yadhu/DataBase/Struct3D/model16SrRNA/HELIX_NR",
987//         HelixGapFile[] = "/nfshome/yadhu/DataBase/Struct3D/model16SrRNA/HELIX_GAP",
988//         ErrorMsg[] = "\n *** Error Opening File : ";
989
990    const char
991        outFile[]      = "/home/Yadhu/DataBase/Struct3D/model23SrRNA/test.data",
992        EcoliFile[]    = "/home/Yadhu/DataBase/Struct3D/model23SrRNA/Ecoli_23SrRNA_with_gaps.seq",
993        HelixNrFile[]  = "/home/Yadhu/DataBase/Struct3D/model23SrRNA/HelixNumbers.msk",
994        HelixGapFile[] = "/home/Yadhu/DataBase/Struct3D/model23SrRNA/SecStructureMask_with_gaps.msk",
995        ErrorMsg[] = "\n *** Error Opening File : ";
996
997//     const char
998//         outFile[]      = "/home/Yadhu/DataBase/Struct3D/model5SrRNA/test.data",
999//         EcoliFile[]    = "/home/Yadhu/DataBase/Struct3D/model5SrRNA/Ecoli_5SrRNA_with_gaps.seq",
1000//         HelixNrFile[]  = "/home/Yadhu/DataBase/Struct3D/model5SrRNA/HelixNumbers.msk",
1001//         HelixGapFile[] = "/home/Yadhu/DataBase/Struct3D/model5SrRNA/SecStructureMask_with_gaps.msk",
1002//         ErrorMsg[] = "\n *** Error Opening File : ";
1003
1004    int fileLen = 0;
1005    char *helixNrBuf, *helixGapBuf, *ecoliBuf;
1006
1007    ofstream out;
1008    out.open(outFile, ios::out);
1009    if (!out.is_open())  cerr<<ErrorMsg<<outFile<<endl;
1010
1011    ifstream inFile;
1012    {
1013        {
1014            inFile.open(EcoliFile, ios::binary);
1015            if (!inFile.is_open())  cerr<<ErrorMsg<<HelixGapFile<<endl;
1016
1017            inFile.seekg (0, ios::end);  // get length of file
1018            fileLen = inFile.tellg();
1019            inFile.seekg (0, ios::beg);
1020
1021            ecoliBuf = new char[fileLen];    // allocate memory:
1022
1023            inFile.read (ecoliBuf, fileLen);    // read data as a block:
1024            inFile.close();
1025        }
1026
1027        {
1028            inFile.open(HelixNrFile, ios::binary);
1029            if (!inFile.is_open())  cerr<<ErrorMsg<<HelixGapFile<<endl;
1030
1031            inFile.seekg (0, ios::end);  // get length of file
1032            fileLen = inFile.tellg();
1033            inFile.seekg (0, ios::beg);
1034
1035            helixNrBuf = new char[fileLen];    // allocate memory:
1036
1037            inFile.read (helixNrBuf, fileLen);    // read data as a block:
1038            inFile.close();
1039        }
1040
1041        {
1042            inFile.open(HelixGapFile, ios::binary);
1043            if (!inFile.is_open())  cerr<<ErrorMsg<<HelixNrFile<<endl;
1044
1045            inFile.seekg (0, ios::end);  // get length of file
1046            fileLen = inFile.tellg();
1047            inFile.seekg (0, ios::beg);
1048
1049            helixGapBuf = new char[fileLen];    // allocate memory:
1050
1051            inFile.read (helixGapBuf, fileLen);    // read data as a block:
1052            inFile.close();
1053        }
1054
1055        char helixNr[4]; helixNr[3] = '\0';
1056        int skip, gaps, k;
1057        skip = gaps = k = 0;
1058
1059        for (unsigned int i = 0; i < strlen(ecoliBuf); i++) {
1060            if (ecoliBuf[i] == '\n')    skip++;
1061            else {
1062                if ((ecoliBuf[i] == '.') || (ecoliBuf[i] == '-'))  gaps++;
1063                else {
1064                    int pos = (i - (skip + gaps)) + 1;
1065                    out<<pos<<"  "<<ecoliBuf[i]<<"  "<<helixGapBuf[i]<<"  ";
1066                    switch (helixGapBuf[i]) {
1067                        case '.': out<<"N  "; break;
1068                        case '[': out<<"S  "; break;
1069                        case ']': out<<"E  "; break;
1070                        case '<': out<<"H  "; break;
1071                        case '>': out<<"H  "; break;
1072                    }
1073                    for (unsigned int j = i; helixNrBuf[j] != '.'; j++) {
1074                        helixNr[k++] = helixNrBuf[j];
1075                    }
1076                    for (int l = k; l < 4; l++)  helixNr[l] = '\0';
1077                    k = 0;
1078                    if ((helixGapBuf[i] == '[') || (helixGapBuf[i] == ']')) {
1079                        out<<helixNr;
1080                    }
1081                    out<<endl;
1082                }
1083            }
1084        }
1085        delete [] helixNrBuf;
1086        delete [] ecoliBuf;
1087        delete [] helixGapBuf;
1088        out.close();
1089    }
1090}
1091
1092void Structure3D::StoreCurrSpeciesDifference(char base, int pos) {
1093    Struct2Dplus3D *st;
1094    st = start2D3D;
1095    while (st) {
1096        if (st->pos == pos) {
1097            CurrSpecies *data, *temp;
1098            data = new CurrSpecies;
1099            data->base = base;
1100            data->pos  = pos;
1101            data->x = st->x;
1102            data->y = st->y;
1103            data->z = st->z;
1104            data->next = NULp;
1105
1106            if (!startSp) {
1107                startSp = data;
1108                bOldSpeciesDataExists = true;
1109            }
1110            else {
1111                temp = startSp;
1112                while (temp->next) {
1113                    temp = temp->next;
1114                }
1115                temp->next = data;
1116            }
1117            break;
1118        }
1119        st = st->next;
1120    }
1121}
1122
1123void Structure3D::DeleteOldSpeciesData() {
1124    CurrSpecies *tmp, *data;
1125
1126    for (data = startSp; data; data = tmp) {
1127        tmp = data->next;
1128        delete data;
1129    }
1130    startSp = NULp;
1131
1132    bOldSpeciesDataExists = false;
1133}
1134
1135void Structure3D::GenerateBaseDifferencePositionDisplayList() {
1136    CurrSpecies *t;
1137    Struct3Dinfo *temp;
1138
1139    char POS[50];
1140    float spacer = 1.5;
1141
1142    glNewList(MAP_SPECIES_BASE_DIFFERENCE_POS, GL_COMPILE);
1143    {
1144        for (t = startSp; t; t = t->next) {
1145            for (temp = start3D; temp; temp = temp->next) {
1146                if (temp->pos == t->pos) {
1147                    sprintf(POS, "%d", temp->pos);
1148                    GRAPHICS->PrintString(temp->x-spacer, temp->y, temp->z-spacer, POS, GLUT_BITMAP_8_BY_13);
1149                }
1150            }
1151        }
1152    }
1153    glEndList();
1154
1155    glNewList(MAP_SPECIES_BASE_DIFFERENCE_POS_ANCHOR, GL_COMPILE);
1156    {
1157        glLineWidth(1.0);
1158        glBegin(GL_LINES);
1159
1160        for (t = startSp; t; t = t->next) {
1161            for (temp = start3D; temp; temp = temp->next) {
1162                if (temp->pos == t->pos) {
1163                    glVertex3f(temp->x, temp->y, temp->z);
1164                    glVertex3f(temp->x-spacer, temp->y, temp->z-spacer);
1165                }
1166            }
1167        }
1168        glEnd();
1169    }
1170    glEndList();
1171}
1172
1173void Structure3D::BuildDisplayList(int listID, int *pos, int len) {
1174    CurrSpecies *t;
1175    int tmpPos = 0;
1176
1177    glNewList(listID, GL_COMPILE);
1178    {
1179        if (RNA3D->bPointSpritesSupported) {
1180            glBegin(GL_POINTS);
1181        }
1182        for (int i = 0; i < len; i++) {
1183            tmpPos = pos[i];
1184            t = startSp;
1185            while (t) {
1186                if (t->pos == tmpPos) {
1187                    PointsToQuads(t->x, t->y, t->z);
1188                    break;
1189                }
1190                t = t->next;
1191            }
1192        }
1193        if (RNA3D->bPointSpritesSupported) {
1194            glEnd();
1195        }
1196    }
1197    glEndList();
1198}
1199
1200void Structure3D::GenerateBaseDifferenceDisplayList() {
1201    CurrSpecies *t;
1202
1203    glNewList(MAP_SPECIES_BASE_DIFFERENCE, GL_COMPILE);
1204    {
1205        if (RNA3D->bPointSpritesSupported) {
1206            glBegin(GL_POINTS);
1207        }
1208        t = startSp;
1209        while (t) {
1210            PointsToQuads(t->x, t->y, t->z);
1211            t = t->next;
1212        }
1213        if (RNA3D->bPointSpritesSupported) {
1214            glEnd();
1215        }
1216    }
1217    glEndList();
1218
1219    const int MAX_BASE = 400;
1220    int baseA[MAX_BASE], baseG[MAX_BASE], baseC[MAX_BASE],
1221        baseU[MAX_BASE], deletion[MAX_BASE], miss[MAX_BASE];
1222    int a, g, c, u, d, m; a=g=c=u=d=m=0;
1223
1224    {
1225        t = startSp;
1226        while (t) {
1227            switch (t->base) {
1228                case 'A': baseA[a++]     = t->pos; break;
1229                case 'G': baseG[g++]     = t->pos; break;
1230                case 'C': baseC[c++]     = t->pos; break;
1231                case 'U': baseU[u++]     = t->pos; break;
1232                case '-': deletion[d++] = t->pos; break;
1233                case '.': miss[m++] = t->pos; break;
1234            }
1235            t = t->next;
1236        }
1237        BuildDisplayList(MAP_SPECIES_BASE_A, baseA, a);
1238        BuildDisplayList(MAP_SPECIES_BASE_U, baseU, u);
1239        BuildDisplayList(MAP_SPECIES_BASE_G, baseG, g);
1240        BuildDisplayList(MAP_SPECIES_BASE_C, baseC, c);
1241        BuildDisplayList(MAP_SPECIES_DELETION, deletion, d);
1242        BuildDisplayList(MAP_SPECIES_MISSING, miss, m);
1243        GenerateBaseDifferencePositionDisplayList();
1244
1245        iTotalSubs = a+g+c+u;  // Summing up the substitutions/mutations
1246        iTotalDels = d;        // Storing total number of deletions
1247#ifdef DEBUG
1248        cout<<"Substitutions = "<<iTotalSubs<<endl;
1249        cout<<"Deletions     = "<<iTotalDels<<endl;
1250#endif
1251    }
1252}
1253
1254void Structure3D::StoreInsertions(char base, int pos) {
1255    Insertions *data, *temp;
1256    data = new Insertions;
1257    data->base = base;
1258    data->pos  = pos;
1259    data->next = NULp;
1260
1261    if (!startIns) {
1262        startIns = data;
1263        bOldInsertionDataExists = true;
1264    }
1265    else {
1266        temp = startIns;
1267        while (temp->next) {
1268            temp = temp->next;
1269        }
1270        temp->next = data;
1271    }
1272}
1273
1274void Structure3D::DeleteOldInsertionData() {
1275    Insertions *tmp, *data;
1276
1277    for (data = startIns; data; data = tmp) {
1278        tmp = data->next;
1279        delete data;
1280    }
1281    startIns = NULp;
1282
1283    bOldInsertionDataExists = false;
1284}
1285
1286void Structure3D::GenerateInsertionDisplayList() {
1287    Insertions   *ins;
1288    Struct3Dinfo *str;
1289    char inserts[500];
1290    int i, cntr;
1291    float spacer = 2.0;
1292    iTotalIns = 0;
1293
1294    glNewList(MAP_SPECIES_INSERTION_BASES, GL_COMPILE);
1295    {
1296        for (str = start3D; str; str = str->next) {
1297            i = cntr = 0;
1298            for (ins = startIns; ins; ins = ins->next) {
1299                if (str->pos == ins->pos) {
1300                    inserts[i++] = ins->base;
1301                    cntr++;
1302                }
1303            }
1304            if (cntr>0) {
1305                inserts[i] = '\0';
1306                char buffer[strlen(inserts) + 10];
1307                sprintf(buffer, "%d:%s", cntr, inserts);
1308                GRAPHICS->PrintString(str->x, str->y+spacer, str->z,
1309                                      buffer, GLUT_BITMAP_8_BY_13);
1310
1311                iTotalIns += cntr; // Summing up the insertions
1312            }
1313        }
1314    }
1315    glEndList();
1316
1317#ifdef DEBUG
1318    cout<<"Insertions = "<<iTotalIns<<endl;
1319#endif
1320
1321    glNewList(MAP_SPECIES_INSERTION_BASES_ANCHOR, GL_COMPILE);
1322    {
1323        glLineWidth(1.0);
1324        glBegin(GL_LINES);
1325
1326        for (str = start3D; str; str = str->next) {
1327            for (ins = startIns; ins; ins = ins->next) {
1328                if (str->pos == ins->pos) {
1329                    glVertex3f(str->x, str->y, str->z);
1330                    glVertex3f(str->x, str->y+spacer, str->z);
1331                }
1332            }
1333        }
1334        glEnd();
1335    }
1336    glEndList();
1337
1338    glNewList(MAP_SPECIES_INSERTION_POINTS, GL_COMPILE);
1339    {
1340        if (RNA3D->bPointSpritesSupported) {
1341            glBegin(GL_POINTS);
1342        }
1343        for (str = start3D; str; str = str->next) {
1344            for (ins = startIns; ins; ins = ins->next) {
1345                if (str->pos == ins->pos) {
1346                    PointsToQuads(str->x, str->y, str->z);
1347                    break;
1348                }
1349            }
1350        }
1351        if (RNA3D->bPointSpritesSupported) {
1352            glEnd();
1353        }
1354    }
1355    glEndList();
1356}
1357
1358void Structure3D::MapCurrentSpeciesToEcoliTemplate(AW_root *awr) {
1359
1360    GB_push_transaction(gb_main);
1361
1362    GBDATA *gbTemplate = GBT_find_SAI(gb_main, "ECOLI");
1363
1364    if (!gbTemplate) {
1365        aw_message("SAI:ECOLI not found");
1366    }
1367    else {
1368        char *pSpeciesName = awr->awar(AWAR_SPECIES_NAME)->read_string();
1369        if (pSpeciesName) {
1370            Host.announce_current_species(pSpeciesName);
1371            GBDATA *gbSpecies = GBT_find_species(gb_main, pSpeciesName);
1372            if (gbSpecies) {
1373                char   *ali_name          = GBT_get_default_alignment(gb_main);
1374                GBDATA *gbAlignment       = GB_entry(gbTemplate, ali_name);
1375                GBDATA *gbTemplateSeqData = gbAlignment ? GB_entry(gbAlignment, "data") : NULp;
1376
1377                if (!gbTemplateSeqData) {
1378                    aw_message(GBS_global_string("Mapping impossible, since SAI:ECOLI has no data in alignment '%s'", ali_name));
1379                }
1380                else {
1381                    const char *pTemplateSeqData  = GB_read_char_pntr(gbTemplateSeqData);
1382
1383                    if (!RNA3D->bEColiRefInitialized) {
1384                        EColiRef = new BI_ecoli_ref;
1385                        EColiRef->init(gb_main);
1386                        RNA3D->bEColiRefInitialized = true;
1387                    }
1388
1389                    char buf[100];
1390                    char *pSpFullName = GB_read_string(GB_entry(gbSpecies, "full_name"));
1391                    sprintf(buf, "%s : %s", pSpeciesName, pSpFullName);
1392                    awr->awar(AWAR_3D_SELECTED_SPECIES)->write_string(buf);
1393
1394                    GBDATA *gbSeqData    = GBT_find_sequence(gbSpecies, ali_name);
1395                    const char *pSeqData = GB_read_char_pntr(gbSeqData);
1396
1397                    if (pSeqData && pTemplateSeqData) {
1398                        int iSeqLen = strlen(pTemplateSeqData);
1399
1400                        if (bOldSpeciesDataExists) {
1401                            DeleteOldSpeciesData();
1402                        }
1403
1404                        for (int i = 0, iSeqCount = 0; i<iSeqLen; i++) {
1405                            if ((pTemplateSeqData[i] != '.') && (pTemplateSeqData[i] != '-')) {
1406                                if (!bStartPosStored) {
1407                                    iStartPos = i;
1408                                    bStartPosStored = true;
1409                                }
1410                                if (pTemplateSeqData[i] != pSeqData[i]) {
1411                                    StoreCurrSpeciesDifference(pSeqData[i], iSeqCount);
1412                                }
1413                                iSeqCount++;
1414                            }
1415                        }
1416
1417                        for (int i = iSeqLen; i>0; i--) {
1418                            if ((pTemplateSeqData[i] != '.') && (pTemplateSeqData[i] != '-')) {
1419                                if (!bEndPosStored) {
1420                                    iEndPos = i;
1421                                    bEndPosStored = true;
1422                                    break;
1423                                }
1424                            }
1425                        }
1426
1427                        if (bOldInsertionDataExists) {
1428                            DeleteOldInsertionData();
1429                        }
1430
1431                        for (int i = iStartPos, iSeqCount = 0; i < iEndPos; i++) {
1432                            if ((pTemplateSeqData[i] != '.') && (pTemplateSeqData[i] != '-')) {
1433                                // Store EColi base positions : Insertion point!
1434                                iSeqCount++;
1435                            }
1436
1437                            if ((pTemplateSeqData[i] == '-') &&
1438                               (pSeqData[i]         != '-') &&
1439                               (pSeqData[i]         != '.'))
1440                            {
1441                                StoreInsertions(pSeqData[i], iSeqCount);
1442                            }
1443                        }
1444                    }
1445                    free(pSpFullName);
1446                }
1447                free(ali_name);
1448            }
1449            GenerateBaseDifferenceDisplayList();
1450            GenerateInsertionDisplayList();
1451        }
1452        free(pSpeciesName);
1453    }
1454    GB_pop_transaction(gb_main);
1455}
1456
1457static bool ValidSearchColor(int iColor, int mode) {
1458    bool isValid = false;
1459
1460    switch (mode) {
1461        case SAI:    isValid = (iColor >= RNA3D_GC_CBACK_0) && (iColor < RNA3D_GC_MAX); break;
1462        case SEARCH: isValid = (iColor >= RNA3D_GC_SBACK_0) && (iColor < RNA3D_GC_MAX); break;
1463    }
1464    return isValid;
1465}
1466
1467void Structure3D::MapSearchStringsToEcoliTemplate(AW_root * /* awr */) {
1468    const char *pSearchColResults = NULp;
1469
1470    if (iMapSearch) pSearchColResults = Host.get_search_background(iStartPos, iEndPos);
1471
1472    if (pSearchColResults) {
1473        int iColor = 0;
1474
1475        glNewList(MAP_SEARCH_STRINGS_TO_STRUCTURE, GL_COMPILE);
1476        {
1477            if (RNA3D->bPointSpritesSupported) glBegin(GL_POINTS);
1478            for (int i = iStartPos; i < iEndPos; i++) {
1479                if (RNA3D->bEColiRefInitialized) {
1480                    long absPos   = (long) i;
1481                    long EColiPos = EColiRef->abs_2_rel(absPos);
1482
1483                    for (Struct3Dinfo *t = start3D; t; t = t->next) {
1484                        if ((t->pos == EColiPos) && (pSearchColResults[i] >= 0)) {
1485                            iColor = pSearchColResults[i] - COLORLINK;
1486
1487                            if (ValidSearchColor(iColor, SEARCH)) {
1488                                RNA3D->cGraphics->SetColor(iColor);
1489                                PointsToQuads(t->x, t->y, t->z);
1490                            }
1491                            break;
1492                        }
1493                    }
1494                }
1495            }
1496            if (RNA3D->bPointSpritesSupported) glEnd();
1497        }
1498        glEndList();
1499
1500        glNewList(MAP_SEARCH_STRINGS_BACKBONE, GL_COMPILE);
1501        {
1502            int iLastClr = 0; int iLastPos = 0; Vector3 vLastPt;
1503            glBegin(GL_LINES);
1504            for (int i = iStartPos; i < iEndPos; i++) {
1505                if (RNA3D->bEColiRefInitialized) {
1506                    long absPos   = (long) i;
1507                    long EColiPos = EColiRef->abs_2_rel(absPos);
1508
1509                    for (Struct3Dinfo *t = start3D; t; t = t->next) {
1510                        if ((t->pos == EColiPos) && (pSearchColResults[i] >= 0)) {
1511                            iColor = pSearchColResults[i] - COLORLINK;
1512
1513                            if (ValidSearchColor(iColor, SEARCH)) {
1514                                if ((iLastClr == iColor) && (iLastPos == EColiPos-1)) {
1515                                    RNA3D->cGraphics->SetColor(iColor);
1516                                    glVertex3f(vLastPt.x, vLastPt.y, vLastPt.z);
1517                                    glVertex3f(t->x, t->y, t->z);
1518                                }
1519                                iLastPos = EColiPos;
1520                                iLastClr = iColor;
1521                                vLastPt.x = t->x; vLastPt.y = t->y; vLastPt.z = t->z;
1522                            }
1523                            break;
1524                        }
1525                    }
1526                }
1527            }
1528            glEnd();
1529        }
1530        glEndList();
1531
1532        RNA3D->bMapSearchStringsDispListCreated = true;
1533    }
1534    else cout<<"BuildColorString did not get the colors : SAI cannot be Visualized!"<<endl;
1535}
1536
1537void Structure3D::MapSaiToEcoliTemplate() {
1538    const char *pSearchColResults = NULp;
1539
1540    if (iMapSAI && Host.SAIs_visualized()) {
1541        pSearchColResults = Host.get_SAI_background(iStartPos, iEndPos); // returns 0 if sth went wrong
1542    }
1543
1544    if (pSearchColResults) {
1545        int iColor = 0;
1546
1547        glNewList(MAP_SAI_TO_STRUCTURE, GL_COMPILE);
1548        {
1549            if (RNA3D->bPointSpritesSupported) glBegin(GL_POINTS);
1550
1551            for (int i = iStartPos; i < iEndPos; i++) {
1552                if (RNA3D->bEColiRefInitialized) {
1553                    long absPos   = (long) i;
1554                    long EColiPos = EColiRef->abs_2_rel(absPos);
1555
1556                    for (Struct3Dinfo *t = start3D; t; t = t->next) {
1557                        if ((t->pos == EColiPos) && (pSearchColResults[i] >= 0)) {
1558                            iColor = pSearchColResults[i-1] - SAICOLORS;
1559
1560                            if (ValidSearchColor(iColor, SAI)) {
1561                                RNA3D->cGraphics->SetColor(iColor);
1562                                PointsToQuads(t->x, t->y, t->z);
1563                            }
1564                            break;
1565                        }
1566                    }
1567                }
1568            }
1569
1570            if (RNA3D->bPointSpritesSupported) glEnd();
1571        }
1572        glEndList();
1573
1574        RNA3D->bMapSaiDispListCreated = true;
1575    }
1576    else cout<<"ED4_getSaiColorString did not get the colors : SAI cannot be Visualized!"<<endl;
1577}
Note: See TracBrowser for help on using the repository browser.