← Back to branch summary

~ubuntu-branches/ubuntu/trusty/clustalx/trusty

« back to all changes in this revision

Viewing changes to clustalW/tree/Tree.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy, Steffen Moeller, Charles Plessy
  • Date: 2009-10-21 13:25:44 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091021132544-r4hbcnjxp354wxh0
Tags: 2.0.12-1
https://launchpad.net/bugs/423648
https://launchpad.net/bugs/393769
http://bugs.debian.org/550893
* New upstream release (LP: #423648, #393769):
  - Uses Qt instead of lesstif.
  - Includes new code for UPGMA guide trees.
  - Includes iterative alignment facility.

[ Steffen Moeller ]
* New upstream release.
* Updated watch file (Closes: #550893).
* Removed LICENSE from debian/clustalx.docs
* rename to clustalx seems no longer required in debian/rules
* moved clustalx.1 into debian folder (eases working with svn-buildpackage)
* added German translation to desktop file

[ Charles Plessy ]
* Updated my email address.
* debian/copyright made machine-readable.
* Added various informations in debian/upstream-metadata.yaml.
* Switched to Debhelper 7.
  (debian/rules, debian/control, debian/patches, debian/compat)
* Removed useless Debhelper file debian/clustalx.dirs.
* Updated package description.
* Hardcoded the localisation of accessory files in /usr/share/clustalx.
  (debian/patches/hardcode-accessory-file-locations.patch)
* Documented in debian/README.source that the documentation for quilt
  is in /usr/share/doc/quilt.
* Added upstream changelog downloaded from upstream website
  (debian/rules, debian/CHANGELOG.upstream).
* Incremented Standards-Version to reflect conformance with Policy 3.8.3
  (debian/control, no other changes needed).
* Updated homepage in debian/clustalw.1.

Show diffs side-by-side

added added

removed removed

Lines of Context:
clustalW/tree/Tree.cpp
 
1
/**
 
2
 * Author: Mark Larkin
 
3
 * 
 
4
 * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.  
 
5
 */
 
6
#ifdef HAVE_CONFIG_H
 
7
    #include "config.h"
 
8
#endif
 
9
#include <sstream>
 
10
#include "Tree.h"
 
11
 
 
12
namespace clustalw
 
13
{
 
14
 
 
15
/**
 
16
 * 
 
17
 * @param firstSeq 
 
18
 * @param lastSeq 
 
19
 * @param sweight 
 
20
 */
 
21
 
 
22
void Tree::calcSeqWeights(int firstSeq, int lastSeq, vector<int>* sweight)
 
23
{
 
24
    if((int)sweight->size() < lastSeq - 1)
 
25
    {
 
26
        sweight->resize(lastSeq - 1);
 
27
    }
 
28
    
 
29
    int i, _nSeqs;
 
30
    int temp, sum;
 
31
    int *weight;
 
32
    /*
 
33
     * If there are more than three sequences....
 
34
     */
 
35
    _nSeqs = lastSeq - firstSeq;
 
36
    if ((_nSeqs >= 2) && (clustalw::userParameters->getDistanceTree() == true) && 
 
37
        (clustalw::userParameters->getNoWeights() == false))
 
38
    {
 
39
        /*
 
40
         * Calculate sequence weights based on Phylip tree.
 
41
         */
 
42
 
 
43
        weight = new int[lastSeq + 1];
 
44
        
 
45
        for (i = firstSeq; i < lastSeq; i++)
 
46
        {
 
47
            weight[i] = calcWeight(i);
 
48
        }
 
49
 
 
50
        /*
 
51
         * Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
 
52
         */
 
53
 
 
54
        sum = 0;
 
55
        for (i = firstSeq; i < lastSeq; i++)
 
56
        {
 
57
            sum += weight[i];
 
58
        }
 
59
 
 
60
        if (sum == 0)
 
61
        {
 
62
            for (i = firstSeq; i < lastSeq; i++)
 
63
            {
 
64
                weight[i] = 1;
 
65
            }
 
66
            sum = i;
 
67
        }
 
68
 
 
69
        for (i = firstSeq; i < lastSeq; i++)
 
70
        {
 
71
            (*sweight)[i] = (weight[i] * INT_SCALE_FACTOR) / sum;
 
72
            if ((*sweight)[i] < 1)
 
73
            {
 
74
                (*sweight)[i] = 1;
 
75
            }
 
76
        }
 
77
        
 
78
        delete []weight;
 
79
        weight = NULL;
 
80
    }
 
81
 
 
82
    else
 
83
    {
 
84
        /*
 
85
         * Otherwise, use identity weights.
 
86
         */
 
87
        temp = INT_SCALE_FACTOR / _nSeqs;
 
88
        // AW 2009-05-28: goes wrong if we have more than
 
89
        // INT_SCALE_FACTOR seqs. if so, set to 1, just as above
 
90
        if (temp < 1)
 
91
            temp = 1;
 
92
 
 
93
        for (i = firstSeq; i < lastSeq; i++)
 
94
        {
 
95
            (*sweight)[i] = temp;
 
96
        }
 
97
    }
 
98
 
 
99
}
 
100
 
 
101
/**
 
102
 * 
 
103
 * @param alignPtr 
 
104
 * @param treeFileName 
 
105
 * @param firstSeq 
 
106
 * @param lastSeq 
 
107
 * @return 
 
108
 */
 
109
int Tree::readTree(clustalw::Alignment* alignPtr, const string& treeFileName, int firstSeq, int lastSeq)
 
110
{
 
111
 
 
112
    if(alignPtr == NULL || firstSeq < 0 || lastSeq < 1)
 
113
    {
 
114
        return 0;
 
115
    }
 
116
    char c;
 
117
    string name1 = "", name2 = "";
 
118
    int i, j, k;
 
119
    bool found;
 
120
 
 
121
    numSeq = 0;
 
122
    nnodes = 0;
 
123
    ntotal = 0;
 
124
    rootedTree = true;
 
125
 
 
126
    // Need to check what happens if I try to open a file that doesnt exist!
 
127
    if(treeFileName.empty())
 
128
    {
 
129
        return 0; // Cannot open, empty string!
 
130
    }
 
131
    else
 
132
    {
 
133
        file.open(treeFileName.c_str(), ifstream::in);
 
134
        if(!file.is_open())
 
135
        {
 
136
            clustalw::utilityObject->error("Cannot open output file [%s]\n", treeFileName.c_str());
 
137
        }
 
138
    }
 
139
 
 
140
    skipSpace(&file);
 
141
    charFromFile = file.get();
 
142
    
 
143
    if (charFromFile != '(')
 
144
    {
 
145
        clustalw::utilityObject->error("Wrong format in tree file %s\n", treeFileName.c_str());
 
146
        return 0;
 
147
    }
 
148
    file.seekg(0, ios::beg); // Put get pointer back to the begining!
 
149
 
 
150
    clustalw::userParameters->setDistanceTree(true);
 
151
 
 
152
    
 
153
    // Allocate memory for tree
 
154
    
 
155
    // Allocate memory.
 
156
    nptr = new TreeNode*[3 * (lastSeq - firstSeq + 1)];
 
157
    ptrs = new TreeNode*[3 * (lastSeq - firstSeq + 1)];
 
158
    lptr = new TreeNode*[lastSeq - firstSeq + 1];
 
159
    olptr = new TreeNode*[lastSeq + 1];
 
160
    
 
161
    seqTree = avail();
 
162
    setInfo(seqTree, NULL, 0, string(""), 0.0);
 
163
 
 
164
    createTree(seqTree, NULL, &file);
 
165
    file.close();
 
166
 
 
167
    if (numSeq != lastSeq - firstSeq)
 
168
    {
 
169
        stringstream ss;
 
170
        ss << "tree not compatible with alignment\n(" << lastSeq - firstSeq 
 
171
        << " sequences in alignment and "<< numSeq <<" in tree\n";
 
172
        string errorMes;
 
173
        ss >> errorMes;
 
174
        clustalw::utilityObject->error(errorMes.c_str());
 
175
        return 0;
 
176
    }
 
177
   
 
178
    // If the tree is unrooted, reroot the tree - ie. minimise the difference
 
179
    // between the mean root->leaf distances for the left and right branches of
 
180
    // the tree.     
 
181
 
 
182
    if (clustalw::userParameters->getDistanceTree() == false)
 
183
    {
 
184
        if (rootedTree == false)
 
185
        {
 
186
            clustalw::utilityObject->error("input tree is unrooted and has no distances.\nCannot align sequences");
 
187
            return 0;
 
188
        }
 
189
    }
 
190
 
 
191
    if (rootedTree == false)
 
192
    {
 
193
        root = reRoot(seqTree, lastSeq - firstSeq + 1);
 
194
    }
 
195
    else
 
196
    {
 
197
        root = seqTree;
 
198
    }
 
199
   
 
200
    // calculate the 'order' of each node.
 
201
    int nameLength;
 
202
    string nameSeq;
 
203
    orderNodes();
 
204
 
 
205
    if (numSeq >= 2)
 
206
    {
 
207
        // If there are more than three sequences....
 
208
        // assign the sequence nodes (in the same order as in the alignment file)
 
209
        for (i = firstSeq; i < lastSeq; i++)
 
210
        {
 
211
            nameLength = alignPtr->getName(i + 1).length();
 
212
            nameSeq = alignPtr->getName(i + 1);
 
213
            name1 = "";
 
214
            name2 = "";
 
215
            if (nameLength > clustalw::MAXNAMES)
 
216
            {
 
217
                stringstream ss;
 
218
                ss << "name " << nameSeq << " is too long for PHYLIP tree format (max " 
 
219
                   << clustalw::MAXNAMES << " chars)";
 
220
                string msg; 
 
221
                ss >> msg; 
 
222
                clustalw::utilityObject->warning(msg.c_str());// Mark change 17-5-07
 
223
            }
 
224
 
 
225
            for (k = 0; k < nameLength && k < clustalw::MAXNAMES; k++)
 
226
            {
 
227
                c = nameSeq[k];
 
228
                if ((c > 0x40) && (c < 0x5b)) 
 
229
                {
 
230
                    c = c | 0x20;
 
231
                }
 
232
                if (c == ' ')
 
233
                {
 
234
                    c = '_';
 
235
                }
 
236
                name2 += c;
 
237
            }
 
238
            found = false;
 
239
 
 
240
            for (j = 0; j < numSeq; j++)
 
241
            {
 
242
                name1 = "";
 
243
                for (k = 0; k < (int)lptr[j]->name.length() && k < clustalw::MAXNAMES; k++)
 
244
                {
 
245
                    c = lptr[j]->name[k];
 
246
                    if ((c > 0x40) && (c < 0x5b))
 
247
                    {
 
248
                        c = c | 0x20;
 
249
                    }
 
250
 
 
251
                    name1 += c;
 
252
                }
 
253
 
 
254
                if (name1.compare(name2) == 0)
 
255
                {
 
256
                    olptr[i] = lptr[j];
 
257
                    found = true;
 
258
                }
 
259
            }
 
260
 
 
261
            if (found == false)
 
262
            {
 
263
                utilityObject->error("tree not compatible with alignment:\n %s not found\n", name2.c_str());
 
264
                return 0;
 
265
            }
 
266
        }
 
267
 
 
268
    }
 
269
    return (1);
 
270
}
 
271
 
 
272
/**
 
273
 * 
 
274
 * @param firstSeq 
 
275
 * @param lastSeq 
 
276
 * @return 
 
277
 */
 
278
auto_ptr<AlignmentSteps> Tree::createSets(int firstSeq, int lastSeq)
 
279
{
 
280
    auto_ptr<AlignmentSteps> progAlignSteps;
 
281
    progAlignSteps.reset(new AlignmentSteps);
 
282
    
 
283
    int i, j, _nSeqs;
 
284
 
 
285
    numSets = 0;
 
286
    _nSeqs = lastSeq - firstSeq;
 
287
    if (_nSeqs >= 2)
 
288
    {
 
289
        // If there are more than three sequences....
 
290
        groups = new int[_nSeqs + 1];
 
291
        groupSeqs(root, groups, _nSeqs, progAlignSteps.get());
 
292
 
 
293
        delete []groups;
 
294
        groups  = NULL;
 
295
 
 
296
    }
 
297
 
 
298
    else
 
299
    {
 
300
        groups = new int[_nSeqs + 1];
 
301
        for (i = 0; i < _nSeqs - 1; i++)
 
302
        {
 
303
            for (j = 0; j < _nSeqs; j++)
 
304
                if (j <= i)
 
305
                {
 
306
                    groups[j] = 1;
 
307
                }
 
308
                else if (j == i + 1)
 
309
                {
 
310
                    groups[j] = 2;
 
311
                }
 
312
                else
 
313
                {
 
314
                    groups[j] = 0;
 
315
                }
 
316
 
 
317
            progAlignSteps->saveSet(_nSeqs, groups);
 
318
        }
 
319
        delete []groups;
 
320
        groups  = NULL;
 
321
    }    
 
322
    
 
323
    return progAlignSteps;
 
324
}
 
325
 
 
326
/**
 
327
 * calcSimilarities changes the distMat.
 
328
 * @param alignPtr 
 
329
 * @param distMat 
 
330
 * @return 
 
331
 */
 
332
int Tree::calcSimilarities(clustalw::Alignment* alignPtr, clustalw::DistMatrix* distMat)
 
333
{
 
334
    int depth = 0, i, j, k, n;
 
335
    bool found;
 
336
    int nerrs, seq1[MAXERRS], seq2[MAXERRS];
 
337
    TreeNode *p,  **pathToRoot;
 
338
    float dist;
 
339
    float *distToNode, badDist[MAXERRS];
 
340
    double **dmat;
 
341
    ostringstream err1;
 
342
    char reply;
 
343
    int nSeqs = alignPtr->getNumSeqs();
 
344
    
 
345
    pathToRoot = new TreeNode*[nSeqs];
 
346
    distToNode = new float[nSeqs];
 
347
    dmat = new double*[nSeqs];
 
348
 
 
349
    for (i = 0; i < nSeqs; i++)
 
350
    {
 
351
        dmat[i] = new double[nSeqs];
 
352
    }
 
353
 
 
354
    if (nSeqs >= 2)
 
355
    {
 
356
        /*
 
357
         * for each leaf, determine all nodes between the leaf and the root;
 
358
         */
 
359
        for (i = 0; i < nSeqs; i++)
 
360
        {
 
361
            depth = 0;dist = 0.0;
 
362
            p = olptr[i];
 
363
            while (p != NULL)
 
364
            {
 
365
                pathToRoot[depth] = p;
 
366
                dist += p->dist;
 
367
                distToNode[depth] = dist;
 
368
                p = p->parent;
 
369
                depth++;
 
370
            }
 
371
 
 
372
            /*
 
373
             * for each pair....
 
374
             */
 
375
            for (j = 0; j < i; j++)
 
376
            {
 
377
                p = olptr[j];
 
378
                dist = 0.0;
 
379
                /*
 
380
                 * find the common ancestor.
 
381
                 */
 
382
                found = false;
 
383
                n = 0;
 
384
                while ((found == false) && (p->parent != NULL))
 
385
                {
 
386
                    for (k = 0; k < depth; k++)
 
387
                    if (p->parent == pathToRoot[k])
 
388
                    {
 
389
                        found = true;
 
390
                        n = k;
 
391
                    }
 
392
 
 
393
                    dist += p->dist;
 
394
                    p = p->parent;
 
395
                }
 
396
 
 
397
                dmat[i][j] = dist + distToNode[n - 1];
 
398
            }
 
399
        }
 
400
 
 
401
        nerrs = 0;
 
402
        for (i = 0; i < nSeqs; i++)
 
403
        {
 
404
            dmat[i][i] = 0.0;
 
405
            for (j = 0; j < i; j++)
 
406
            {
 
407
                if (dmat[i][j] < 0.01)
 
408
                {
 
409
                    dmat[i][j] = 0.01;
 
410
                }
 
411
                if (dmat[i][j] > 1.0)
 
412
                {
 
413
                    if (dmat[i][j] > 1.1 && nerrs < MAXERRS)
 
414
                    {
 
415
                        seq1[nerrs] = i;
 
416
                        seq2[nerrs] = j;
 
417
                        badDist[nerrs] = dmat[i][j];
 
418
                        nerrs++;
 
419
                    }
 
420
                    dmat[i][j] = 1.0;
 
421
                }
 
422
            }
 
423
        }
 
424
        if (nerrs > 0)
 
425
        {
 
426
            string errMess = "The following sequences are too divergent to be aligned:\n";
 
427
            
 
428
            for (i = 0; i < nerrs && i < 5; i++)
 
429
            {
 
430
                err1 << "           " << alignPtr->getName(seq1[i] + 1) << " and " 
 
431
                     << alignPtr->getName(seq2[i] + 1) << " (distance " 
 
432
                     << setprecision(3) << badDist[i] << ")\n"; 
 
433
            }
 
434
            errMess += err1.str();
 
435
            errMess += "(All distances should be between 0.0 and 1.0)\n";
 
436
            errMess += "This may not be fatal but you have been warned!\n";
 
437
            errMess += "SUGGESTION: Remove one or more problem sequences and try again";
 
438
            
 
439
            if (clustalw::userParameters->getInteractive())
 
440
            {
 
441
                reply = clustalw::utilityObject->promptForYesNo(errMess.c_str(), "Continue ");
 
442
            }
 
443
            else
 
444
            {
 
445
                reply = 'y';
 
446
            }
 
447
            if ((reply != 'y') && (reply != 'Y'))
 
448
            {
 
449
                return 0;
 
450
            }
 
451
        }
 
452
    }
 
453
    else
 
454
    {
 
455
        for (i = 0; i < nSeqs; i++)
 
456
        {
 
457
            for (j = 0; j < i; j++)
 
458
            {
 
459
                dmat[i][j] = (*distMat)(i + 1, j + 1);
 
460
            }
 
461
        }
 
462
    }
 
463
    delete []pathToRoot;
 
464
    delete []distToNode;
 
465
    
 
466
    double value;
 
467
    for (i = 0; i < nSeqs; i++)
 
468
    {
 
469
        distMat->SetAt(i + 1, i + 1, 0.0);
 
470
        for (j = 0; j < i; j++)
 
471
        {
 
472
            value = 100.0 - (dmat[i][j]) * 100.0;
 
473
            distMat->SetAt(i + 1, j + 1, value);
 
474
            distMat->SetAt(j + 1, i + 1, value);
 
475
        }
 
476
    }
 
477
 
 
478
    for (i = 0; i < nSeqs; i++)
 
479
    {
 
480
        delete [] dmat[i];
 
481
    }
 
482
 
 
483
    delete []dmat;
 
484
 
 
485
    return 1;
 
486
    
 
487
}
 
488
 
 
489
/** *************************************************************************
 
490
 * Private functions!!!!!!!!!!!!!!!                                         *
 
491
 ****************************************************************************/
 
492
 
 
493
/**
 
494
 * 
 
495
 * @param ptree 
 
496
 * @param parent 
 
497
 * @param file 
 
498
 */
 
499
void Tree::createTree(clustalw::TreeNode* ptree, clustalw::TreeNode* parent, ifstream* file)
 
500
{
 
501
    TreeNode* p;
 
502
 
 
503
    int i, type;
 
504
    float dist;
 
505
    string name;
 
506
 
 
507
    
 
508
    // is this a node or a leaf ?
 
509
    skipSpace(file);
 
510
    charFromFile = file->get();
 
511
    if (charFromFile == '(')
 
512
    {
 
513
        // this must be a node....
 
514
        type = NODE;
 
515
        name = "";
 
516
        ptrs[ntotal] = nptr[nnodes] = ptree;
 
517
        nnodes++;
 
518
        ntotal++;
 
519
 
 
520
        createNode(ptree, parent);
 
521
 
 
522
        p = ptree->left;
 
523
        createTree(p, ptree, file);
 
524
 
 
525
        if (charFromFile == ',')
 
526
        {
 
527
            p = ptree->right;
 
528
            createTree(p, ptree, file);
 
529
            if (charFromFile == ',')
 
530
            {
 
531
                ptree = insertNode(ptree);
 
532
                ptrs[ntotal] = nptr[nnodes] = ptree;
 
533
                nnodes++;
 
534
                ntotal++;
 
535
                p = ptree->right;
 
536
                createTree(p, ptree, file);
 
537
                rootedTree = false;
 
538
            }
 
539
        }
 
540
 
 
541
        skipSpace(file);
 
542
        charFromFile = file->get();
 
543
    }
 
544
    
 
545
    // ...otherwise, this is a leaf
 
546
    
 
547
    else
 
548
    {
 
549
        type = LEAF;
 
550
        ptrs[ntotal++] = lptr[numSeq++] = ptree;        
 
551
        // get the sequence name
 
552
        name = "";
 
553
        name += charFromFile;
 
554
        charFromFile = file->get();
 
555
        
 
556
        i = 1;
 
557
        while ((charFromFile != ':') && (charFromFile != ',') && (charFromFile != ')'))
 
558
        {
 
559
            if (i < MAXNAMES)
 
560
            {
 
561
                name += charFromFile;
 
562
                i++;
 
563
            }
 
564
            charFromFile = file->get();
 
565
        }
 
566
 
 
567
        if (charFromFile != ':')
 
568
        {
 
569
            clustalw::userParameters->setDistanceTree(false);
 
570
            dist = 0.0;
 
571
        }
 
572
    }
 
573
    
 
574
    // get the distance information
 
575
    
 
576
    dist = 0.0;
 
577
    if (charFromFile == ':')
 
578
    {
 
579
        skipSpace(file);
 
580
        (*file) >> dist;
 
581
        skipSpace(file);
 
582
        charFromFile = file->get();
 
583
    }
 
584
    setInfo(ptree, parent, type, name, dist);
 
585
}
 
586
 
 
587
/**
 
588
 * 
 
589
 * @param pptr 
 
590
 * @param parent 
 
591
 */
 
592
void Tree::createNode(TreeNode* pptr, TreeNode* parent)
 
593
{
 
594
    TreeNode* t;
 
595
    pptr->parent = parent;
 
596
    t = avail();
 
597
    pptr->left = t;
 
598
    t = avail();
 
599
    pptr->right = t;
 
600
}
 
601
 
 
602
/**
 
603
 * 
 
604
 * @param pptr 
 
605
 * @return 
 
606
 */
 
607
TreeNode* Tree::insertNode(TreeNode* pptr)
 
608
{
 
609
 
 
610
    TreeNode* newnode;
 
611
 
 
612
    newnode = avail();
 
613
    createNode(newnode, pptr->parent);
 
614
 
 
615
    newnode->left = pptr;
 
616
    pptr->parent = newnode;
 
617
 
 
618
    setInfo(newnode, pptr->parent, NODE, "", 0.0);
 
619
 
 
620
    return newnode;
 
621
}
 
622
 
 
623
/**
 
624
 * 
 
625
 * @param p 
 
626
 */
 
627
void Tree::clearTree(clustalw::TreeNode* p)
 
628
{
 
629
    clearTreeNodes(p);
 
630
    delete [] nptr;
 
631
    nptr  = NULL;
 
632
    delete [] ptrs;
 
633
    ptrs  = NULL;
 
634
    delete [] lptr;
 
635
    lptr  = NULL;
 
636
    delete [] olptr;
 
637
    olptr  = NULL;
 
638
}
 
639
 
 
640
/**
 
641
 * 
 
642
 * @param p 
 
643
 */
 
644
void Tree::clearTreeNodes(clustalw::TreeNode* p)
 
645
{
 
646
    if (p == NULL)
 
647
    {
 
648
        p = root;
 
649
    }
 
650
    if (p->left != NULL)
 
651
    {
 
652
        clearTreeNodes(p->left);
 
653
    }
 
654
    if (p->right != NULL)
 
655
    {
 
656
        clearTreeNodes(p->right);
 
657
    }
 
658
    p->left = NULL;
 
659
    p->right = NULL;
 
660
    
 
661
    delete p;
 
662
    p  = NULL;
 
663
}
 
664
 
 
665
 
 
666
 
 
667
 
 
668
/**
 
669
 * 
 
670
 * @param 
 
671
 * @param 
 
672
 * @return 
 
673
 */
 
674
void Tree::debugPrintAllNodes(int nseqs)
 
675
{
 
676
  clustalw::TreeNode *p;
 
677
  int i;
 
678
   float diff, maxDist;
 
679
 
 
680
    cerr << "\nDEBUG: reportAllNodes\n";
 
681
    for (i = 0; i < ntotal; i++) {
 
682
        p = ptrs[i];
 
683
        //        ios::sync_with_stdio();
 
684
 
 
685
        // same design as TreeNode
 
686
        if (p->parent == NULL)
 
687
            diff = calcRootMean(p, &maxDist);
 
688
        else
 
689
            diff = calcMean(p, &maxDist, nseqs);
 
690
        fprintf(stdout,
 
691
                "i=%d p=%p: parent=%p left=%p right=%p dist=%f diff=%f\n",
 
692
                i, (void*)p, (void*)p->parent, (void*)p->left, (void*)p->right,
 
693
                p->dist, diff);
 
694
    }
 
695
}
 
696
 
 
697
 
 
698
 
 
699
 
 
700
 
 
701
/**
 
702
 * 
 
703
 * @param ptree 
 
704
 * @param nseqs 
 
705
 * @return 
 
706
 */
 
707
clustalw::TreeNode* Tree::reRoot(clustalw::TreeNode* ptree, int nseqs)
 
708
{
 
709
    clustalw::TreeNode *p, *rootNode, *rootPtr;
 
710
    float diff, minDiff = 0.0, minDepth = 1.0, maxDist;
 
711
    int i;
 
712
    bool first = true;
 
713
 
 
714
    // find the difference between the means of leaf->node
 
715
    // distances on the left and on the right of each node
 
716
    rootPtr = ptree;
 
717
    for (i = 0; i < ntotal; i++)
 
718
    {
 
719
        p = ptrs[i];
 
720
        if (p->parent == NULL)
 
721
        {
 
722
            /* AW Bug 94: p->parent must be chosen as rootNode
 
723
               (non-optimized executables (-O0) never do), otherwise
 
724
               insertRoot fails.
 
725
               Is p->parent == NULL valid at all?
 
726
               Simplest thing for now is to continue here. Tree code
 
727
               needs serious dismantling anyway. See debugPrintAllNodes
 
728
            */
 
729
 
 
730
            continue;
 
731
            //diff = calcRootMean(p, &maxDist);
 
732
        }
 
733
        else
 
734
        {
 
735
            diff = calcMean(p, &maxDist, nseqs);
 
736
        }
 
737
 
 
738
        if ((diff == 0) || ((diff > 0) && (diff < 2 *p->dist)))
 
739
        {
 
740
            if ((maxDist < minDepth) || (first == true))
 
741
            {
 
742
                first = false;
 
743
                rootPtr = p;
 
744
                minDepth = maxDist;
 
745
                minDiff = diff;
 
746
            }
 
747
        }
 
748
 
 
749
    }
 
750
 
 
751
    
 
752
    // insert a new node as the ancestor of the node which produces the shallowest
 
753
    // tree.
 
754
    /* AW Bug 94: could also be prevented here */
 
755
    if (rootPtr == ptree)
 
756
    {
 
757
        minDiff = rootPtr->left->dist + rootPtr->right->dist;
 
758
        rootPtr = rootPtr->right;
 
759
    }
 
760
    rootNode = insertRoot(rootPtr, minDiff);
 
761
 
 
762
    diff = calcRootMean(rootNode, &maxDist);
 
763
 
 
764
    return rootNode;
 
765
}
 
766
 
 
767
/**
 
768
 * 
 
769
 * @param p 
 
770
 * @param diff 
 
771
 * @return 
 
772
 */
 
773
clustalw::TreeNode* Tree::insertRoot(clustalw::TreeNode* p, float diff)
 
774
{
 
775
    clustalw::TreeNode *newp, *prev, *q, *t;
 
776
    float dist, prevDist, td;
 
777
 
 
778
    newp = avail();
 
779
 
 
780
    if (p->parent==NULL) {
 
781
        // AW bug 94: question remains if access here should be handled differently
 
782
        cerr << "\n\n*** INTERNAL ERROR: Tree::insertRoot: TreeNode p->parent is NULL\n";
 
783
        cerr << "To help us fix this bug, please send sequence file and used options to clustalw@ucd.ie\n";
 
784
        exit(1);
 
785
    }
 
786
 
 
787
    t = p->parent;
 
788
    prevDist = t->dist;
 
789
 
 
790
    p->parent = newp;
 
791
 
 
792
    dist = p->dist;
 
793
 
 
794
    p->dist = diff / 2;
 
795
 
 
796
    if (p->dist < 0.0)
 
797
    {
 
798
        p->dist = 0.0;
 
799
    }
 
800
 
 
801
    if (p->dist > dist)
 
802
    {
 
803
        p->dist = dist;
 
804
    }
 
805
 
 
806
    t->dist = dist - p->dist;
 
807
 
 
808
    newp->left = t;
 
809
    newp->right = p;
 
810
    newp->parent = NULL;
 
811
    newp->dist = 0.0;
 
812
    newp->leaf = NODE;
 
813
 
 
814
    if (t->left == p)
 
815
    {
 
816
        t->left = t->parent;
 
817
    }
 
818
    else
 
819
    {
 
820
        t->right = t->parent;
 
821
    }
 
822
 
 
823
    prev = t;
 
824
    q = t->parent;
 
825
 
 
826
    t->parent = newp;
 
827
 
 
828
    while (q != NULL)
 
829
    {
 
830
        if (q->left == prev)
 
831
        {
 
832
            q->left = q->parent;
 
833
            q->parent = prev;
 
834
            td = q->dist;
 
835
            q->dist = prevDist;
 
836
            prevDist = td;
 
837
            prev = q;
 
838
            q = q->left;
 
839
        }
 
840
        else
 
841
        {
 
842
            q->right = q->parent;
 
843
            q->parent = prev;
 
844
            td = q->dist;
 
845
            q->dist = prevDist;
 
846
            prevDist = td;
 
847
            prev = q;
 
848
            q = q->right;
 
849
        }
 
850
    }
 
851
 
 
852
    /*
 
853
     * remove the old root node
 
854
     */
 
855
    q = prev;
 
856
    if (q->left == NULL)
 
857
    {
 
858
        dist = q->dist;
 
859
        q = q->right;
 
860
        q->dist += dist;
 
861
        q->parent = prev->parent;
 
862
 
 
863
        if (prev->parent->left == prev)
 
864
        {
 
865
            prev->parent->left = q;
 
866
        }
 
867
        else
 
868
        {
 
869
            prev->parent->right = q;
 
870
        }
 
871
 
 
872
        prev->right = NULL;
 
873
    }
 
874
    else
 
875
    {
 
876
        dist = q->dist;
 
877
        q = q->left;
 
878
        q->dist += dist;
 
879
        q->parent = prev->parent;
 
880
 
 
881
        if (prev->parent->left == prev)
 
882
        {
 
883
            prev->parent->left = q;
 
884
        }
 
885
        else
 
886
        {
 
887
            prev->parent->right = q;
 
888
        }
 
889
 
 
890
        prev->left = NULL;
 
891
    }
 
892
 
 
893
    return (newp);
 
894
}
 
895
 
 
896
/**
 
897
 * 
 
898
 * @param root 
 
899
 * @param maxDist 
 
900
 * @return 
 
901
 */
 
902
float Tree::calcRootMean(clustalw::TreeNode* root, float *maxDist)
 
903
{
 
904
    float dist, leftSum = 0.0, rightSum = 0.0, leftMean, rightMean, diff;
 
905
    clustalw::TreeNode* p;
 
906
    int i;
 
907
    int numLeft, numRight;
 
908
    int direction;
 
909
    
 
910
    // for each leaf, determine whether the leaf is left or right of the root.
 
911
 
 
912
    dist = (*maxDist) = 0;
 
913
    numLeft = numRight = 0;
 
914
    for (i = 0; i < numSeq; i++)
 
915
    {
 
916
        p = lptr[i];
 
917
        dist = 0.0;
 
918
        while (p->parent != root)
 
919
        {
 
920
            dist += p->dist;
 
921
            p = p->parent;
 
922
        }
 
923
 
 
924
        if (p == root->left)
 
925
        {
 
926
            direction = LEFT;
 
927
        }
 
928
        else
 
929
        {
 
930
            direction = RIGHT;
 
931
        }
 
932
 
 
933
        dist += p->dist;
 
934
 
 
935
        if (direction == LEFT)
 
936
        {
 
937
            leftSum += dist;
 
938
            numLeft++;
 
939
        }
 
940
        else
 
941
        {
 
942
            rightSum += dist;
 
943
            numRight++;
 
944
        }
 
945
        if (dist > (*maxDist))
 
946
        {
 
947
            *maxDist = dist;
 
948
        }
 
949
    }
 
950
 
 
951
    leftMean = leftSum / numLeft;
 
952
    rightMean = rightSum / numRight;
 
953
 
 
954
    diff = leftMean - rightMean;
 
955
    return diff;
 
956
}
 
957
 
 
958
/**
 
959
 * 
 
960
 * @param nptr 
 
961
 * @param maxDist 
 
962
 * @param nSeqs 
 
963
 * @return 
 
964
 */
 
965
float Tree::calcMean(clustalw::TreeNode* nptr, float *maxDist, int nSeqs)
 
966
{
 
967
    float dist, leftSum = 0.0, rightSum = 0.0, leftMean, rightMean, diff;
 
968
    clustalw::TreeNode* p;  
 
969
    clustalw::TreeNode** pathToRoot;
 
970
    float *distToNode;
 
971
    int depth = 0, i, j, n = 0;
 
972
    int numLeft, numRight;
 
973
    int direction;
 
974
    bool found;
 
975
 
 
976
    pathToRoot = new clustalw::TreeNode*[nSeqs];
 
977
    distToNode = new float[nSeqs];
 
978
    // determine all nodes between the selected node and the root;
 
979
  
 
980
    depth = 0;
 
981
    (*maxDist) = dist = 0.0;
 
982
    numLeft = numRight = 0;
 
983
    p = nptr;
 
984
    while (p != NULL)
 
985
    {
 
986
        pathToRoot[depth] = p;
 
987
        dist += p->dist;
 
988
        distToNode[depth] = dist;
 
989
        p = p->parent;
 
990
        depth++;
 
991
    }
 
992
 
 
993
    // for each leaf, determine whether the leaf is left or right of the node.
 
994
    // (RIGHT = descendant, LEFT = not descendant)
 
995
 
 
996
    for (i = 0; i < numSeq; i++)
 
997
    {
 
998
        p = lptr[i];
 
999
        if (p == nptr)
 
1000
        {
 
1001
            direction = RIGHT;
 
1002
            dist = 0.0;
 
1003
        }
 
1004
        else
 
1005
        {
 
1006
            direction = LEFT;
 
1007
            dist = 0.0;
 
1008
            
 
1009
            // find the common ancestor. 
 
1010
            
 
1011
            found = false;
 
1012
            n = 0;
 
1013
            while ((found == false) && (p->parent != NULL))
 
1014
            {
 
1015
                for (j = 0; j < depth; j++)
 
1016
                if (p->parent == pathToRoot[j])
 
1017
                {
 
1018
                    found = true;
 
1019
                    n = j;
 
1020
                }
 
1021
 
 
1022
                dist += p->dist;
 
1023
                p = p->parent;
 
1024
            }
 
1025
            if (p == nptr)
 
1026
            {
 
1027
                direction = RIGHT;
 
1028
            }
 
1029
        }
 
1030
 
 
1031
        if (direction == LEFT)
 
1032
        {
 
1033
            leftSum += dist;
 
1034
            leftSum += distToNode[n - 1];
 
1035
            numLeft++;
 
1036
        }
 
1037
        else
 
1038
        {
 
1039
            rightSum += dist;
 
1040
            numRight++;
 
1041
        }
 
1042
 
 
1043
        if (dist > (*maxDist))
 
1044
        {
 
1045
            *maxDist = dist;
 
1046
        }
 
1047
    }
 
1048
 
 
1049
    delete [] distToNode;
 
1050
    distToNode = NULL;
 
1051
    delete [] pathToRoot;
 
1052
    pathToRoot = NULL;
 
1053
    
 
1054
    leftMean = leftSum / numLeft;
 
1055
    rightMean = rightSum / numRight;
 
1056
 
 
1057
    diff = leftMean - rightMean;
 
1058
    return (diff);
 
1059
}
 
1060
 
 
1061
/**
 
1062
 * 
 
1063
 */
 
1064
void Tree::orderNodes()
 
1065
{
 
1066
    int i;
 
1067
    clustalw::TreeNode* p;
 
1068
 
 
1069
    for (i = 0; i < numSeq; i++)
 
1070
    {
 
1071
        p = lptr[i];
 
1072
        while (p != NULL)
 
1073
        {
 
1074
            p->order++;
 
1075
            p = p->parent;
 
1076
        }
 
1077
    }
 
1078
}
 
1079
 
 
1080
/**
 
1081
 * 
 
1082
 * @param leaf 
 
1083
 * @return 
 
1084
 */
 
1085
int Tree::calcWeight(int leaf)
 
1086
{
 
1087
    clustalw::TreeNode* p;
 
1088
    float weight = 0.0;
 
1089
 
 
1090
    p = olptr[leaf];
 
1091
    while (p->parent != NULL)
 
1092
    {
 
1093
        weight += p->dist / p->order;
 
1094
        p = p->parent;
 
1095
    }
 
1096
 
 
1097
    weight *= 100.0;
 
1098
 
 
1099
    return ((int)weight);
 
1100
}
 
1101
 
 
1102
/**
 
1103
 * skipSpace is used to skip all the spaces at the begining of a file. The next read
 
1104
 * will give a character other than a space.
 
1105
 * @param file 
 
1106
 */
 
1107
void Tree::skipSpace(ifstream* file)
 
1108
{
 
1109
    char c;
 
1110
 
 
1111
    do
 
1112
    {
 
1113
        c = file->get();
 
1114
    }
 
1115
    while (isspace(c));
 
1116
 
 
1117
    file->putback(c);
 
1118
}
 
1119
 
 
1120
/**
 
1121
 * 
 
1122
 * @param p 
 
1123
 * @param nextGroups 
 
1124
 * @param nSeqs 
 
1125
 * @param stepsPtr 
 
1126
 */
 
1127
void Tree::groupSeqs(clustalw::TreeNode* p, int *nextGroups, int nSeqs, AlignmentSteps* stepsPtr)
 
1128
{
 
1129
    int i;
 
1130
    int *tmpGroups;
 
1131
 
 
1132
    tmpGroups = new int[nSeqs + 1];
 
1133
    for (i = 0; i < nSeqs; i++)
 
1134
    {
 
1135
        tmpGroups[i] = 0;
 
1136
    }
 
1137
 
 
1138
    if (p->left != NULL)
 
1139
    {
 
1140
        if (p->left->leaf == NODE)
 
1141
        {
 
1142
            groupSeqs(p->left, nextGroups, nSeqs, stepsPtr);
 
1143
 
 
1144
            for (i = 0; i < nSeqs; i++)
 
1145
                if (nextGroups[i] != 0)
 
1146
                {
 
1147
                    tmpGroups[i] = 1;
 
1148
                }
 
1149
        }
 
1150
        else
 
1151
        {
 
1152
            markGroup1(p->left, tmpGroups, nSeqs);
 
1153
        }
 
1154
 
 
1155
    }
 
1156
 
 
1157
    if (p->right != NULL)
 
1158
    {
 
1159
        if (p->right->leaf == NODE)
 
1160
        {
 
1161
            groupSeqs(p->right, nextGroups, nSeqs, stepsPtr);
 
1162
            for (i = 0; i < nSeqs; i++)
 
1163
                if (nextGroups[i] != 0)
 
1164
                {
 
1165
                    tmpGroups[i] = 2;
 
1166
                }
 
1167
        }
 
1168
        else
 
1169
        {
 
1170
            markGroup2(p->right, tmpGroups, nSeqs);
 
1171
        }
 
1172
        stepsPtr->saveSet(nSeqs, tmpGroups);
 
1173
    }
 
1174
 
 
1175
    for (i = 0; i < nSeqs; i++)
 
1176
    {
 
1177
        nextGroups[i] = tmpGroups[i];
 
1178
    }
 
1179
 
 
1180
    delete [] tmpGroups;
 
1181
    tmpGroups  = NULL;
 
1182
}
 
1183
 
 
1184
/**
 
1185
 * 
 
1186
 * @param p 
 
1187
 * @param groups 
 
1188
 * @param n 
 
1189
 */
 
1190
void Tree::markGroup1(clustalw::TreeNode* p, int *groups, int n)
 
1191
{
 
1192
    int i;
 
1193
 
 
1194
    for (i = 0; i < n; i++)
 
1195
    {
 
1196
        if (olptr[i] == p)
 
1197
        {
 
1198
            groups[i] = 1;
 
1199
        }
 
1200
        else
 
1201
        {
 
1202
            groups[i] = 0;
 
1203
        }
 
1204
    }
 
1205
}
 
1206
 
 
1207
/**
 
1208
 * 
 
1209
 * @param p 
 
1210
 * @param groups 
 
1211
 * @param n 
 
1212
 */
 
1213
void Tree::markGroup2(clustalw::TreeNode* p, int *groups, int n)
 
1214
{
 
1215
    int i;
 
1216
 
 
1217
    for (i = 0; i < n; i++)
 
1218
    {
 
1219
        if (olptr[i] == p)
 
1220
        {
 
1221
            groups[i] = 2;
 
1222
        }
 
1223
        else if (groups[i] != 0)
 
1224
        {
 
1225
            groups[i] = 1;
 
1226
        }
 
1227
    }
 
1228
}
 
1229
 
 
1230
/**
 
1231
 * 
 
1232
 * @return 
 
1233
 */
 
1234
clustalw::TreeNode* Tree::avail()
 
1235
{
 
1236
    clustalw::TreeNode* p;
 
1237
    p = new TreeNode; 
 
1238
    p->left = NULL;
 
1239
    p->right = NULL;
 
1240
    p->parent = NULL;
 
1241
    p->dist = 0.0;
 
1242
    p->leaf = 0;
 
1243
    p->order = 0;
 
1244
    return (p);
 
1245
}
 
1246
 
 
1247
/**
 
1248
 * 
 
1249
 * @param p 
 
1250
 * @param parent 
 
1251
 * @param pleaf 
 
1252
 * @param pname 
 
1253
 * @param pdist 
 
1254
 */
 
1255
void Tree::setInfo(TreeNode* p, TreeNode* parent, int pleaf, string pname, float
 
1256
                     pdist)
 
1257
{
 
1258
    p->parent = parent;
 
1259
    p->leaf = pleaf;
 
1260
    p->dist = pdist;
 
1261
    p->order = 0;
 
1262
    p->name = pname;
 
1263
    if (p->leaf == true)
 
1264
    {
 
1265
        p->left = NULL;
 
1266
        p->right = NULL;
 
1267
    }
 
1268
}
 
1269
          
 
1270
}