//////////////////////////////////////////////////////////////////////////////////

//

// STAMP version 1

//

// Written By: Shaun Mahony

//

// Tree.cpp

//

// Started: 5th Dec 2005

//

// Copyright 2007-2015 Shaun Mahony

//

// This file is part of STAMP.

//

// STAMP is free software; you can redistribute it and/or modify

// it under the terms of the GNU General Public License as published by

// the Free Software Foundation; either version 2 of the License, or

// (at your option) any later version.

//

// STAMP is distributed in the hope that it will be useful,

// but WITHOUT ANY WARRANTY; without even the implied warranty of

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

// GNU General Public License for more details.

//

// You should have received a copy of the GNU General Public License

// along with STAMP; if not, write to the Free Software

// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

//

////////////////////////////////////////////////////////////////////////////////////

#include "Tree.h"

//Destructor

TreeNode::~TreeNode()

{

Child* tmpChild;

if(profile!=NULL){delete profile;}

tmpChild = progeny;

while(progeny!=NULL){

progeny = tmpChild->next;

delete tmpChild;

tmpChild=progeny;

}

if(alignment!=NULL)

delete alignment;

}

//kill the child nodes

void Tree::KillChildren(TreeNode* n){

Child* tmpChild;

tmpChild = n->progeny;

while(n->progeny!=NULL){

n->progeny = tmpChild->next;

delete tmpChild;

tmpChild=n->progeny;

}n->progeny=NULL;

}

//Tree node Random init

void TreeNode::RandomInit()

{

double Ac, Gc, Cc, Tc;

int INIT_SIZE=8;

//Random initialisation

profile = new Motif(INIT_SIZE);

for(int k=0; k<INIT_SIZE; k++)

{

Ac=(double)rand()/RAND_MAX;

Cc=(double)rand()/RAND_MAX;

Gc=(double)rand()/RAND_MAX;

Tc=(double)rand()/RAND_MAX;

double ttl = Ac+Cc+Gc+Tc;

Ac=Ac/ttl;

Cc=Cc/ttl;

Gc=Gc/ttl;

Tc=Tc/ttl;

profile->f[k][0]=Ac;

profile->f[k][1]=Cc;

profile->f[k][2]=Gc;

profile->f[k][3]=Tc;

}

}

//UPGMA's tree building method

void UPGMA::BuildTree(PlatformSupport* p, bool treeTest)

{

Plat = p;

treeTesting=treeTest;

int numMotifs = Plat->GetMatCount();

Motif** motifSet = Plat->inputMotifs;

AlignRec** pairwiseAlign = Plat->pairwiseAlign;

int i, j, k, z;

TreeNode* tmpNode;

TreeNode** nodes = new TreeNode*[numMotifs];

bool* active = new bool[numMotifs];

double minCH = 10000000;

nodesMinCH = numMotifs;

// double* error = new double[numMotifs];////////////////Necessary only for Davies-Bouldin

for(i=0; i<numMotifs; i++){

nodes[i] = new TreeNode();

nodes[i]->profile = new Motif(motifSet[i]->GetLen());

motifSet[i]->CopyMotif(nodes[i]->profile);

nodes[i]->leaf=true;

nodes[i]->members=1;

nodes[i]->leafID=i;

active[i]=true;

PostorderListChildren(nodes[i], nodes[i]);

}double generalMean=0, count=0;

double** pairwiseDist = new double* [numMotifs];

for(i=0; i<numMotifs; i++){

pairwiseDist[i]=new double [numMotifs];

for(j=0; j<numMotifs; j++){

if(i!=j)

{ pairwiseDist[i][j]=pairwiseAlign[i][j].dist;

generalMean+=(pairwiseAlign[i][j].dist*pairwiseAlign[i][j].dist);

count++;

}

}

}generalMean=generalMean/count;

//This loop iterates until all nodes are on the tree

for(z=0; z<numMotifs-1; z++){

double minDist=10000;

int minNodeA, minNodeB;

//Step 1: find the two nodes for which dij is minimal

for(i=0; i<numMotifs; i++){

if(active[i]){

for(j=0; j<numMotifs; j++){

if(i!=j && active[j] && pairwiseDist[i][j]<minDist){

minDist = pairwiseDist[i][j];

minNodeA=i; minNodeB=j;

}

}

}

}

//Step 2: put minNodes A&B into new Node

tmpNode = new TreeNode();

tmpNode->left = nodes[minNodeA];

tmpNode->right = nodes[minNodeB];

nodes[minNodeB]->parent=tmpNode;

tmpNode->members = nodes[minNodeA]->members + nodes[minNodeB]->members;

tmpNode->height = minDist/2.0;

nodes[minNodeA]->edge=tmpNode->height-nodes[minNodeA]->height;

nodes[minNodeB]->edge=tmpNode->height-nodes[minNodeB]->height;

tmpNode->leaf=false;

tmpNode->nodeID = z;

PostorderListChildren(tmpNode, tmpNode);

if(treeTesting)

IRAlignment(tmpNode);

else

BuildFBP(tmpNode, pairwiseAlign,z);

char tmpname[50];

sprintf(tmpname, "Internal_%d", z);

strcpy(tmpNode->profile->name, tmpname);

//Step 3: Add tmpNode into minNodeA's spot, inactivate minNodeB's spot, update dists for new node

nodes[minNodeA] = tmpNode;

active[minNodeB]=false;

for(j=0; j<numMotifs; j++){

if(minNodeA!=j && active[j])

{ pairwiseDist[minNodeA][j] = (pairwiseDist[minNodeA][j]*(double)nodes[minNodeA]->members + pairwiseDist[minNodeB][j]*(double)nodes[minNodeB]->members)/((double)nodes[minNodeA]->members + (double)nodes[minNodeB]->members);

pairwiseDist[j][minNodeA] = (pairwiseDist[j][minNodeA]*(double)nodes[minNodeA]->members + pairwiseDist[j][minNodeB]*(double)nodes[minNodeB]->members)/((double)nodes[minNodeA]->members + (double)nodes[minNodeB]->members);

}

}

//test the exclusivity of each node

if(treeTesting){

Child* tmpA; Child* tmpB;

double nodeCount=0;

double totalHomo=0, currHomo=0;

for(j=0; j<numMotifs; j++){

if(active[j]){nodeCount++;}

}

nodeCount=0;

//Homogeneity of the node contents

for(j=0; j<numMotifs; j++){

currHomo=0;

if(active[j]){

nodeCount++;

if(nodes[j]->members==1){

currHomo++;

}else{

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){

for(tmpB=nodes[j]->progeny; tmpB!=NULL; tmpB=tmpB->next)

if(strcmp(tmpA->m->famName, tmpB->m->famName)==0)

currHomo++;

}

}

currHomo=currHomo/((double)nodes[j]->members*(double)nodes[j]->members);

totalHomo+=currHomo;

}

}

totalHomo=totalHomo/nodeCount;

///////////Calinski & Harabasz///////////////////////////

//Internal similarity of the nodes & between clusters

double internalSOSE=0, betweenSOSE=0, Ak=0;

double currE=0, currCount, tmpAvgDist;double tmpE;

for(j=0; j<numMotifs; j++){

if(active[j]){

if(nodes[j]->members==1)

{ currE=0;currCount=1;count=1;

}else{currCount=0;count=0;currE=0;

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){count++;

for(tmpB=nodes[j]->progeny; tmpB->mID!=tmpA->mID; tmpB=tmpB->next)

if(tmpA->mID!=tmpB->mID)

{

currCount++; currE+=pairwiseAlign[tmpA->mID][tmpB->mID].dist*pairwiseAlign[tmpB->mID][tmpA->mID].dist;

}

}

}

internalSOSE+=currE;

Ak+=(count-1)*(generalMean-(currE/currCount));

}

}Ak = Ak/((double)numMotifs-nodeCount);

for(int w=0; w<numMotifs; w++){

for(int v=0; v<w; v++){

if(v!=w && active[w] && active[v]){

tmpAvgDist = (pairwiseDist[v][w]+pairwiseDist[w][v])/2;

betweenSOSE+=(tmpAvgDist*tmpAvgDist);

}

}

}

internalSOSE = log(internalSOSE);

betweenSOSE = log(betweenSOSE);

if(nodeCount>1){

//This works!

//printf("%.0lf\t%lf\t%.10lf\t%lf\n", nodeCount, betweenSOSE, internalSOSE, (((betweenSOSE+nodeCount)/(nodeCount-1))/((internalSOSE+nodeCount)/(numMotifs-nodeCount))));

//printf("%.0lf\t%lf\t%.10lf\t%lf\n", nodeCount, betweenSOSE, internalSOSE, ((betweenSOSE/(nodeCount-1))+(nodeCount/numMotifs))/((internalSOSE/(numMotifs-nodeCount))+(nodeCount/numMotifs)));

//printf("%.0lf\t%lf\t%lf\t%lf\t%lf\n", nodeCount, totalHomo, (betweenSOSE/(nodeCount-1)), internalSOSE/(numMotifs-nodeCount), (((betweenSOSE)/(nodeCount-1))/((internalSOSE)/(numMotifs-nodeCount))));

double currCH = (((betweenSOSE)/(nodeCount-1))/((internalSOSE)/(numMotifs-nodeCount)));

if(!silence){

printf("\t%.0lf\t%lf\n", nodeCount, currCH);

}

if(currCH < minCH){minCH=currCH; nodesMinCH = nodeCount;}

//HC p10

if((generalMean-Ak)>0){

double HC = (generalMean+(((double)numMotifs-nodeCount)/(nodeCount-1))*Ak)/(generalMean-Ak);

}

double numerate=(generalMean+(((double)numMotifs-nodeCount)/(nodeCount-1))*Ak);

double denominate=generalMean-Ak;

}

//Level-specific code!

if(nodeCount==-1){

FILE* tmpF = fopen("FBPmotifs.transfac", "w");

for(j=0; j<numMotifs; j++){

if(active[j]){

//nodes[j]->profile->PrintMotif(tmpF);

PrintNode(nodes[j], tmpF);

}

}fclose(tmpF);

//pseudo LOOCV

printf("Non-matches\n");

count=0; double correct=0;

for(j=0; j<numMotifs; j++){

if(active[j]){

if(nodes[j]->members==1){

}else{

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){

count++;

printf("Testing %s_%s\n", tmpA->m->famName, tmpA->m->GetName());

//PPAAlignment(nodes[j], nodes[j], tmpA->mID);

IRAlignment(nodes[j], tmpA->mID);

double LOOMin = 100000000, LOOMax=0; int LOOMinIndex = 0;

double currPred=0;

for(int q=0; q<numMotifs; q++){//Find closest node

currPred=0;

if(active[q]){

int t_i1, t_i2, t_al; bool t_f1, t_f2;

double currS = Aman->AlignMotifs2D(tmpA->m, nodes[q]->profile, t_i1, t_i2, t_al, t_f1, t_f2);

double currP = Plat->Score2PVal(tmpA->m->GetLen(), nodes[q]->profile->GetLen(), currS);

if(currP>LOOMax){

LOOMax=currP; LOOMinIndex=q;

}

}

}

if(LOOMinIndex==j){

correct++; nodes[j]->profile->PrintMotifConsensus(); printf("Correct... score = %.4e\n", 1-LOOMax);

}else{

printf("%s_%s\t%s\n", tmpA->m->famName, tmpA->m->GetName(), nodes[LOOMinIndex]->profile->GetName());

printf("\t%s\t", tmpA->m->GetName());

tmpA->m->PrintMotifConsensus();

printf("\n\t%s\t", nodes[j]->profile->GetName());

nodes[j]->profile->PrintMotifConsensus();

printf("\n\t%s\t", nodes[LOOMinIndex]->profile->GetName());

nodes[LOOMinIndex]->profile->PrintMotifConsensus();

printf("\n");

printf("Incorrect... score = %.4e\n", 1-LOOMax);

}

}

//PPAAlignment(nodes[j], nodes[j]);

IRAlignment(nodes[j]);

}

}

}

printf("\nPerformance: %lf (%.0lf from %.0lf)\n", correct/count, correct, count);

}

//

/*

//////////////////Gap Statistic///////////////////////////////

//Internal similarity of the nodes & between clusters

double internalE=0, betweenE=0, Ak=0;

double currE=0, currCount, tmpAvgDist;double tmpE;

for(j=0; j<numMotifs; j++){

if(active[j]){

if(nodes[j]->members==1)

{ currE=0;currCount=1;count=1;

}else{currCount=0;count=0;currE=0;

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){count++;

for(tmpB=nodes[j]->progeny; tmpB->mID!=tmpA->mID; tmpB=tmpB->next)

if(tmpA->mID!=tmpB->mID)

{ //tmpE=pairwiseAlign[tmpA->mID][tmpB->mID].p_value;

currCount++; currE+=pairwiseAlign[tmpA->mID][tmpB->mID].dist*pairwiseAlign[tmpB->mID][tmpA->mID].dist;}

}

}

internalE+=currE;//(currE/(currCount));

}

}printf("Gap\t%.0lf\t%lf\t%lf\n", nodeCount, internalE, log(internalE));

//Davies-Bouldin

//first compute e's

double currE=0, currCount, tmpAvgDist;double tmpE;

for(j=0; j<numMotifs; j++){

if(active[j]){

if(nodes[j]->members==1)

{ currE=0;currCount=1;count=1;

}else{currCount=0;count=0;currE=0;

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){count++;

for(tmpB=nodes[j]->progeny; tmpB!=NULL; tmpB=tmpB->next)

if(tmpA->mID!=tmpB->mID)

{ currCount++; currE+=pairwiseAlign[tmpA->mID][tmpB->mID].dist;

}

}

}

error[j]=(currE/currCount);

}

}

//Next find the maximal R_jk's

double DB=0, Rk=0, currR; int MaxRk=0;

for(k=0; k<numMotifs; k++){

if(active[k]){

Rk=0;

for(j=0; j<numMotifs; j++){

if(active[j] && k!=j){

currR = (error[j]+error[k])/((pairwiseDist[j][k]+pairwiseDist[k][j])/2);

if(currR>Rk){

Rk=currR;

MaxRk = j;

}

}

}

DB+=Rk;

}

}DB = DB/nodeCount;

printf("%.0lf\t%lf\n", nodeCount, DB);

*/

}

}

//Tree building finished, the last node (i.e. tmpNode), should be the root

root = tmpNode;

root->nodeID=numMotifs-1;

root->edge = tmpNode->height;

delete [] active;

}

//UPGMA's tree building method (LOOCV)

void UPGMA::LOOCVBuildTree(PlatformSupport* p, bool treeTest)

{

Plat = p;

treeTesting=treeTest;

int numMotifs = Plat->GetMatCount();

Motif** motifSet = Plat->inputMotifs;

AlignRec** pairwiseAlign = Plat->pairwiseAlign;

int i, j, k, z, x, q;

//Matrices here for keeping average performance

double* homogeneities = new double[numMotifs-1];

double* predictivity = new double[numMotifs-1];

double* HCList = new double[numMotifs-1]; double HCMin=1000;double HC, HCMinPred, HCMinIndex;

double totalHCMin=0, totalHCMinPred=0;

for(q=0; q<numMotifs-1; q++){

homogeneities[q]=0;

predictivity[q]=0;

HCList[q]=0;

}

for(x=0; x<numMotifs; x++){HCMin=1000;HCMinIndex=-1;

TreeNode* tmpNode;

TreeNode** nodes = new TreeNode*[numMotifs];

bool* active = new bool[numMotifs];

for(i=0; i<numMotifs; i++){

nodes[i] = new TreeNode();

nodes[i]->profile = new Motif(motifSet[i]->GetLen());

motifSet[i]->CopyMotif(nodes[i]->profile);

nodes[i]->leaf=true;

nodes[i]->members=1;

nodes[i]->leafID=i;

active[i]=true;

}

double** pairwiseDist = new double* [numMotifs];

for(i=0; i<numMotifs; i++){

pairwiseDist[i]=new double [numMotifs];

for(j=0; j<numMotifs; j++){

if(i!=j)

{ pairwiseDist[i][j]=pairwiseAlign[i][j].dist;

}

}

}

//This loop iterates until all nodes are on the tree

for(z=0; z<numMotifs-2; z++){

double minDist=10000;

int minNodeA, minNodeB;

//Step 1: find the two nodes for which dij is minimal

for(i=0; i<numMotifs; i++){

if(i!=x){

if(active[i]){

for(j=0; j<numMotifs; j++){

if(j!=x){

if(i!=j && active[j] && pairwiseDist[i][j]<minDist){

minDist = pairwiseDist[i][j];

minNodeA=i; minNodeB=j;

}}

}

}}

}

//Step 2: put minNodes A&B into new Node

tmpNode = new TreeNode();

tmpNode->left = nodes[minNodeA];

tmpNode->right = nodes[minNodeB];

nodes[minNodeB]->parent=tmpNode;

tmpNode->members = nodes[minNodeA]->members + nodes[minNodeB]->members;

tmpNode->height = minDist/2.0;

nodes[minNodeA]->edge=tmpNode->height-nodes[minNodeA]->height;

nodes[minNodeB]->edge=tmpNode->height-nodes[minNodeB]->height;

tmpNode->leaf=false;

tmpNode->nodeID = z;

PostorderListChildren(tmpNode, tmpNode);

//BuildFBP(tmpNode, pairwiseAlign,z);

PPAAlignment(tmpNode, tmpNode);

//IRAlignment(tmpNode);

//Step 3: Add tmpNode into minNodeA's spot, inactivate minNodeB's spot, update dists for new node

nodes[minNodeA] = tmpNode;

active[minNodeB]=false;

for(j=0; j<numMotifs; j++){

if(minNodeA!=j && active[j])

{ pairwiseDist[minNodeA][j] = (pairwiseDist[minNodeA][j]*(double)nodes[minNodeA]->members + pairwiseDist[minNodeB][j]*(double)nodes[minNodeB]->members)/((double)nodes[minNodeA]->members + (double)nodes[minNodeB]->members);

pairwiseDist[j][minNodeA] = (pairwiseDist[j][minNodeA]*(double)nodes[minNodeA]->members + pairwiseDist[j][minNodeB]*(double)nodes[minNodeB]->members)/((double)nodes[minNodeA]->members + (double)nodes[minNodeB]->members);

}

}

//test the exclusivity of each node

if(treeTesting){

Child* tmpA; Child* tmpB;

double nodeCount=0;

double totalHomo=0, currHomo=0;//Homogeneity of the node contents

for(j=0; j<numMotifs; j++){

if(j!=x){

currHomo=0;

if(active[j]){

nodeCount++;

if(nodes[j]->members==1)

currHomo++;

else{

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next)

for(tmpB=nodes[j]->progeny; tmpB!=NULL; tmpB=tmpB->next)

if(strcmp(tmpA->m->famName, tmpB->m->famName)==0)

currHomo++;

}

currHomo=currHomo/((double)nodes[j]->members*(double)nodes[j]->members);

totalHomo+=currHomo;

}}

}

totalHomo=totalHomo/nodeCount;

homogeneities[(int)nodeCount]+=totalHomo;

//printf("%.0lf\t%lf\n", nodeCount, totalHomo);

///////////Calinski & Harabasz///////////////////////////

double internalSOSE=0, betweenSOSE=0;

double currE=0, currCount, tmpAvgDist;double tmpE;

for(j=0; j<numMotifs; j++){

if(j!=x && active[j]){

if(nodes[j]->members==1)

{ currE=0;currCount=1;

}else{currCount=0;currE=0;

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next){

for(tmpB=nodes[j]->progeny; tmpB->mID!=tmpA->mID; tmpB=tmpB->next)

if(tmpA->mID!=tmpB->mID)

{ currCount++; currE+=pairwiseAlign[tmpA->mID][tmpB->mID].dist*pairwiseAlign[tmpB->mID][tmpA->mID].dist;

}

}

}internalSOSE+=currE;

}

}

for(int w=0; w<numMotifs; w++){

for(int v=0; v<w; v++){

if(v!=x && w!=x && v!=w && active[w] && active[v]){

tmpAvgDist = (pairwiseDist[v][w]+pairwiseDist[w][v])/2;

betweenSOSE+=(tmpAvgDist*tmpAvgDist);

}

}

}

internalSOSE = log(internalSOSE);

betweenSOSE = log(betweenSOSE);

if(nodeCount>1 && nodeCount<(numMotifs-1)){

HC = (((betweenSOSE)/(nodeCount-1))/((internalSOSE)/((numMotifs-1)-nodeCount)));

HCList[(int)nodeCount]+=HC;

if(HC<HCMin){

HCMin=HC;

HCMinIndex=nodeCount;

}

//printf("%.0lf\t%lf\t%lf\t%lf\t%lf\n", nodeCount, totalHomo, (betweenSOSE/(nodeCount-1)), internalSOSE/(numMotifs-nodeCount), (((betweenSOSE)/(nodeCount-1))/((internalSOSE)/(numMotifs-nodeCount))));

}

//now test the predictive capacity of the current tree against the left-out motif

double LOOMin = 100000000, LOOMax=0;; int LOOMinIndex = 0;

double currPred=0;

for(j=0; j<numMotifs; j++){//Find closest node

if(j!=x){

currPred=0;

if(active[j]){

//This is one way of doing it...

/*if(pairwiseDist[x][j]<LOOMin){

LOOMin=pairwiseDist[x][j];LOOMinIndex=j;

}if(pairwiseDist[j][x]<LOOMin){

LOOMin=pairwiseDist[j][x];LOOMinIndex=j;

}*/

//This is another; align versus the FBPs

int t_i1, t_i2, t_al; bool t_f1, t_f2;

double currS = Aman->AlignMotifs2D(motifSet[x], nodes[j]->profile, t_i1, t_i2, t_al, t_f1, t_f2);

double currP = Plat->Score2PVal(motifSet[x]->GetLen(), nodes[j]->profile->GetLen(), currS);

if(currP>LOOMax){

LOOMax=currP; LOOMinIndex=j;

}

}

}

}

if(nodes[LOOMinIndex]->members==1 && strcmp(motifSet[LOOMinIndex]->famName, motifSet[x]->famName)==0){

currPred++;

}else{

for(tmpA=nodes[LOOMinIndex]->progeny; tmpA!=NULL; tmpA=tmpA->next)

if(strcmp(tmpA->m->famName, motifSet[x]->famName)==0)

currPred++;

}

currPred = currPred/((double)nodes[LOOMinIndex]->members);

predictivity[(int)nodeCount]+=currPred;

//Catch the predictivity for the predicted number of clusters!

if(nodeCount==HCMinIndex){

HCMinPred=currPred;

}

}

}

//Tree building finished, the last node (i.e. tmpNode), should be the root

root = tmpNode;

root->nodeID=numMotifs-1;

root->edge = tmpNode->height;

delete [] active;

PostorderDeleteTree(root);

printf("%.0lf\t%lf\n", HCMinIndex, HCMinPred);

totalHCMin+=HCMinIndex; totalHCMinPred+=HCMinPred;

}

printf("Homogeneities\n");

for(q=0; q<numMotifs-1; q++){

printf("%d\t%lf\t%lf\n", q, homogeneities[q]/((double)numMotifs), predictivity[q]/((double)numMotifs));

}printf("Average Clusters Predicted: %lf\nAverage Performance: %lf\n", totalHCMin/(double)numMotifs, totalHCMinPred/(double)numMotifs);

delete [] homogeneities;

delete [] predictivity;

delete [] HCList;

BuildTree(p, false);

}

// *********************************************************************************** //

//Neighbour-joinings's tree building method

void Neighbourjoin::BuildTree(PlatformSupport* p, bool treeTest)

{

Plat = p;

treeTesting=treeTest;

int numMotifs = Plat->GetMatCount();

Motif** motifSet = Plat->inputMotifs;

AlignRec** pairwiseAlign = Plat->pairwiseAlign;

int i, j, k, z;

int minNodeA, minNodeB;

TreeNode* tmpNode;

TreeNode** nodes = new TreeNode*[numMotifs];

bool* active = new bool[numMotifs];

double numActive = (double)numMotifs;

double* r = new double[numMotifs];

//double r_i, r_j, min_r_i, min_r_j;

for(i=0; i<numMotifs; i++){

nodes[i] = new TreeNode();

nodes[i]->profile = new Motif(motifSet[i]->GetLen());

motifSet[i]->CopyMotif(nodes[i]->profile);

nodes[i]->leaf=true;

nodes[i]->members=1;

active[i]=true;

nodes[i]->leafID=i;

r[i]=0;

}

double** pairwiseDist = new double* [numMotifs];

for(i=0; i<numMotifs; i++){

pairwiseDist[i]=new double [numMotifs];

for(j=0; j<numMotifs; j++){

if(i!=j)

{ //pairwiseDist[i][j]=(1-pairwiseAlign[i][j].p_value);

pairwiseDist[i][j]=pairwiseAlign[i][j].dist;

r[i]+=pairwiseDist[i][j];

}

}

r[i]=r[i]/(numActive-2.0);

}

//Test for additivity!

/* int a,b,c,d;

for(d=0; d<numMotifs; d++)

for(c=0; c<d; c++)

for(b=0; b<c; b++)

for(a=0; a<b; a++){

printf("A:%d B:%d C:%d D:%d\tAB+CD:%lf\tAC+BD:%lf\tAD+BC:%lf\n", a, b, c, d, pairwiseDist[a][b]+pairwiseDist[c][d], pairwiseDist[a][c]+pairwiseDist[b][d], pairwiseDist[a][d]+pairwiseDist[b][c]);

}

*/

//This loop iterates until all nodes are on the tree

for(z=0; z<numMotifs-2; z++){

double minDist=1000000;

//Step 1: find the two nodes for which Dij is minimal

for(i=0; i<numMotifs; i++){

if(active[i]){

for(j=0; j<i; j++){

if(i!=j && active[j]){//&& (pairwiseDist[i][j]-(r[i]+r[j]))<minDist){

if(pairwiseDist[i][j]-(r[i]+r[j])<minDist){

minDist = (pairwiseDist[i][j]-(r[i]+r[j]));

minNodeA=i; minNodeB=j;

}

}

}

}

}

// printf("Got Here: %d, %d, %lf, %.0lf, %d, %d\n",minNodeA, minNodeB, pairwiseDist[minNodeA][minNodeB]-(r[minNodeA]+r[minNodeB]), numActive, numMotifs, z);

//Step 2: put minNodes A&B into new Node

tmpNode = new TreeNode();

tmpNode->left = nodes[minNodeA];

nodes[minNodeA]->parent=tmpNode;

nodes[minNodeA]->edge=(pairwiseDist[minNodeA][minNodeB]/2)+((r[minNodeA]-r[minNodeB])/2);

tmpNode->right = nodes[minNodeB];

nodes[minNodeB]->parent=tmpNode;

nodes[minNodeB]->edge=pairwiseDist[minNodeA][minNodeB]-nodes[minNodeA]->edge;

tmpNode->members = nodes[minNodeA]->members + nodes[minNodeB]->members;

// tmpNode->height = minDist/2.0;

tmpNode->leaf=false;

tmpNode->nodeID=z;

PostorderListChildren(tmpNode, tmpNode);

BuildFBP(tmpNode, pairwiseAlign,z);

//Step 3: Add tmpNode into minNodeA's spot, inactivate minNodeB's spot, update dists for new node

nodes[minNodeA] = tmpNode;

active[minNodeB]=false; numActive--;

double AvsB = pairwiseDist[minNodeA][minNodeB];

double BvsA = pairwiseDist[minNodeB][minNodeA];

for(j=0; j<numMotifs; j++){

if(minNodeA!=j && active[j])

{ pairwiseDist[minNodeA][j] = (pairwiseDist[minNodeA][j]+pairwiseDist[minNodeB][j]-AvsB)/2;

pairwiseDist[j][minNodeA] = (pairwiseDist[j][minNodeA]+pairwiseDist[j][minNodeB]-BvsA)/2;

}

}

for(i=0; i<numMotifs; i++){

r[i]=0;

for(j=0; j<numMotifs; j++){

if(i!=j && active[i] && active[j])

r[i]+=pairwiseDist[i][j];

}r[i]=r[i]/(numActive-2);

}

//test the exclusivity of each node

if(treeTesting){

Child* tmpA; Child* tmpB;

double nodeCount=0;

double totalHomo=0, currHomo=0;//Homogeneity of the node contents

for(j=0; j<numMotifs; j++){

currHomo=0;

if(active[j]){

nodeCount++;

if(nodes[j]->members==1)

currHomo++;

else{

for(tmpA=nodes[j]->progeny; tmpA!=NULL; tmpA=tmpA->next)

for(tmpB=nodes[j]->progeny; tmpB!=NULL; tmpB=tmpB->next)

if(strcmp(tmpA->m->famName, tmpB->m->famName)==0)

currHomo++;

}

currHomo=currHomo/((double)nodes[j]->members*(double)nodes[j]->members);

totalHomo+=currHomo;

}

}

totalHomo=totalHomo/nodeCount;

printf("%.0lf\t%lf\n", nodeCount, totalHomo);

}

}

bool oneFound=false;

for(i=0; i<numMotifs; i++){

if(active[i] && !oneFound)

{ minNodeA=i; oneFound=true;}

else if(active[i])

minNodeB=i;

}

tmpNode = new TreeNode();

tmpNode->leaf=false;

tmpNode->nodeID=numMotifs-1;

tmpNode->left = nodes[minNodeA];

nodes[minNodeA]->parent=tmpNode;

nodes[minNodeA]->edge=pairwiseDist[minNodeA][minNodeB]/2;

tmpNode->right = nodes[minNodeB];

nodes[minNodeB]->parent=tmpNode;

nodes[minNodeB]->edge=pairwiseDist[minNodeA][minNodeB]/2;

tmpNode->members = nodes[minNodeA]->members + nodes[minNodeB]->members;

PostorderListChildren(tmpNode, tmpNode);

BuildFBP(tmpNode, pairwiseAlign,numMotifs-1);

//Tree building finished, the last node (i.e. tmpNode), should be the root

root = tmpNode;

root->edge = pairwiseDist[minNodeA][minNodeB];

delete [] active;

delete [] r;

}

// *********************************************************************************** //

//Private method: Print out the tree (now with full names and internal node labels)

void Tree::PostorderPrintTree(TreeNode* n, FILE* out, FILE* orderMat)

{

double dist;

if(n->left != NULL){ fprintf(out, "("); PostorderPrintTree(n->left, out, orderMat); fprintf(out, ",");}

if(n->right != NULL){

PostorderPrintTree(n->right, out, orderMat);

dist=n->edge;

// Internal nodes now get labels (Internal_N)

fprintf(out, ")%s:%lf", n->profile ? n->profile->name : "", fabs(dist));

}if(n->leaf)

{

if(orderMat!=NULL)

n->profile->PrintMotif(orderMat);

// Full name — no more 20-char Phylip truncation

fprintf(out, "%s:%lf", n->profile->name, fabs(n->edge));

}else{

if(orderMat!=NULL)

n->profile->PrintMotif(orderMat);

}

}

//Private method: Print enhanced tree (labeled Newick) — same logic, used for webmode stdout

void Tree::PostorderPrintEnhancedTree(TreeNode* n, FILE* out)

{

double dist;

if(n->left != NULL){ fprintf(out, "("); PostorderPrintEnhancedTree(n->left, out); fprintf(out, ",");}

if(n->right != NULL){

PostorderPrintEnhancedTree(n->right, out);

dist=n->edge;

fprintf(out, ")%s:%lf", n->profile ? n->profile->name : "", fabs(dist));

}if(n->leaf)

{

fprintf(out, "%s:%lf", n->profile->name, fabs(n->edge));

}

}

//Private method: Postorder print internal profiles

void Tree::PostorderPrintInternalProfiles(TreeNode* n, FILE* out)

{

if(n->left != NULL) PostorderPrintInternalProfiles(n->left, out);

if(n->right != NULL) PostorderPrintInternalProfiles(n->right, out);

if(!n->leaf && n->profile != NULL){

fprintf(out, ">>INTERNAL_NODE\t%s\t%d\n", n->profile->name, n->nodeID);

fprintf(out, "DE\t%s\n", n->profile->name);

for(int i = 0; i < n->profile->len; i++){

fprintf(out, "%d\t", i);

for(int j = 0; j < B; j++)

fprintf(out, "%.4lf\t", n->profile->f[i][j]);

fprintf(out, "%c\n", n->profile->ColConsensus(i));

}

fprintf(out, "XX\n");

}

}

//Print the tree

void Tree::PrintTree(char* outFile)

{

char outName[STR_LEN];

sprintf(outName, "%s.tree", outFile);

FILE* out = fopen(outName, "w");

if(out==NULL){printf("error; can't open output file %s\n", outName); exit(1);}

//FILE* orderMat = fopen("orderMat.txt", "w+");

FILE* orderMat=NULL;

PostorderPrintTree(root, out, orderMat);

fprintf(out, ";");

if(orderMat!=NULL)

fclose(orderMat);

fclose(out);

}

//Print internal profiles to stdout in delimited section (webmode)

void Tree::PrintInternalProfiles(FILE* out)

{

if(root == NULL) return;

fprintf(out, ">>STAMP_INTERNAL_PROFILES_START\n");

PostorderPrintInternalProfiles(root, out);

fprintf(out, ">>STAMP_INTERNAL_PROFILES_END\n");

}

//Print labeled Newick tree to stdout in delimited section (webmode)

void Tree::PrintEnhancedTree(FILE* out)

{

if(root == NULL) return;

fprintf(out, ">>STAMP_LABELED_TREE_START\n");

PostorderPrintEnhancedTree(root, out);

fprintf(out, ";\n");

fprintf(out, ">>STAMP_LABELED_TREE_END\n");

}

//Write internal profiles to file (default mode)

void Tree::WriteInternalProfiles(char* outPrefix)

{

if(root == NULL) return;

char outFName[STR_LEN];

sprintf(outFName, "%s_internal_profiles.txt", outPrefix);

FILE* out = fopen(outFName, "w");

if(out == NULL){

fprintf(stderr, "Error: can't open output file %s\n", outFName);

return;

}

PostorderPrintInternalProfiles(root, out);

fclose(out);

}

//Tree deletion

void Tree::PostorderDeleteTree(TreeNode* n)

{

if(n->left != NULL){ PostorderDeleteTree(n->left);}

if(n->right != NULL){PostorderDeleteTree(n->right);}

delete n;

}

//Populate the list of progeny

void Tree::PostorderListChildren(TreeNode* n, TreeNode* start)

{

if(n->left != NULL){PostorderListChildren(n->left, start);}

if(n->right != NULL){PostorderListChildren(n->right, start);}

if(n->leaf)

{

Child* tmp=new Child();

tmp->next = start->progeny;

tmp->m=n->profile;

tmp->mID=n->leafID;

start->progeny = tmp;

}

}

//Print the tree's details

void Tree::PrintNodes(char* outFile)

{

char outName[STR_LEN];

sprintf(outName, "%s_order.transfac", outFile);

FILE* out = fopen(outName, "w");

if(out==NULL){printf("error; can't open output file\n"); exit(1);}

PostorderPrintNodes(root, out);

fclose(out);

}

//Print the matrices and children names for each node

void Tree::PostorderPrintNodes(TreeNode* n, FILE* out)

{

Child* tmp;

int i, j;

// for(tmp=n->progeny; tmp!=NULL; tmp=tmp->next){

// fprintf(out, "%s\t", tmp->m->name);//printf("%s\n", tmp->m->name);

// }

if(n->leaf)

{

fprintf(out, "\nDE %s\n", n->profile->name);

for(i=0; i<n->profile->len; i++)

{

fprintf(out, "%d\t", i);

for(j=0; j<B; j++)

fprintf(out, "%.4lf\t", n->profile->f[i][j]);

fprintf(out, "%c\n", n->profile->ColConsensus(i));

}

fprintf(out, "XX\n\n");/**/

//fprintf(out, "\n\n");

}

if(n->left != NULL){PostorderPrintNodes(n->left, out);}

if(n->right != NULL){PostorderPrintNodes(n->right, out);}

}

//Print the matrices and children names for a node

void Tree::PrintNode(TreeNode* n, FILE* out)

{

Child* tmp;

int i, j;

fprintf(out, "XX\tCluster FBP\n");

fprintf(out, "DE %s\n", n->profile->name);

for(i=0; i<n->profile->len; i++)

{

fprintf(out, "%d\t", i);

for(j=0; j<B; j++)

fprintf(out, "%.4lf\t", n->profile->f[i][j]);

fprintf(out, "%c\n", n->profile->ColConsensus(i));

}

fprintf(out, "XX\n");/**/

fprintf(out, "XX\tCluster_Members:\t");

for(tmp=n->progeny; tmp!=NULL; tmp=tmp->next){

fprintf(out, "%s\t", tmp->m->name);//printf("%s\n", tmp->m->name);

}

fprintf(out, "\n\n");

}

//Build FBPs for a node based on its list of children (GAPS NOT SUPPORTED YET!)

void Tree::BuildFBP(TreeNode* n, AlignRec** pairwiseAlign, int nameID)

{

Child* curr; Child* currB;

int i, j, k, m, a, x, y, currCount;

int maxLen3 = 3*maxLen;

int aStart = maxLen;

int anchor; double currSim; double maxSim=0; double count=0;

double** alignmentMat = new double*[maxLen3];

double* alignCnt = new double [maxLen3];

for(k=0; k<maxLen3; k++)

{ alignmentMat[k] = new double[B];

alignCnt[k]=0;

for(m=0; m<B; m++)

{ alignmentMat[k][m]=0;

}

}

for(curr=n->progeny; curr!=NULL; curr=curr->next){count++;};

double* weights = new double[(int)count];

//first find the anchor (maxSim to all others)

i=0;

for(curr=n->progeny; curr!=NULL; curr=curr->next){

currSim=0;

for(currB=n->progeny; currB!=NULL; currB=currB->next){

if(curr->mID != currB->mID){

currSim+=pairwiseAlign[curr->mID][currB->mID].p_value;

}

}

currSim = currSim/(count-1);

weights[i]=currSim;

if(currSim>maxSim){

maxSim=currSim;

anchor = curr->mID;

}

i++;

}

//Add in the hits according to weight

curr = n->progeny; currCount=0;

while(curr!=NULL){

if(curr->mID != anchor)

{

Motif* rev = new Motif(curr->m->len);

curr->m->RevCompMotif(rev);

// ------->

j=pairwiseAlign[anchor][curr->mID].i2;

for(i=aStart+pairwiseAlign[anchor][curr->mID].i1; (i<maxLen3 && j<curr->m->len); i++)

{

for(k=0; k<B; k++)

{ if(pairwiseAlign[anchor][curr->mID].forward1)

alignmentMat[i][k]+=(curr->m->f[j][k] * weights[currCount]);

else

alignmentMat[i][k]+=(rev->f[j][k] * weights[currCount]);

}

alignCnt[i]++;

j++;

}

// <-------

j=pairwiseAlign[anchor][curr->mID].i2-1;

for(i=aStart+(pairwiseAlign[anchor][curr->mID].i1-1); (i>=0 && j>=0); i--)

{

for(k=0; k<B; k++)

{ if(pairwiseAlign[anchor][curr->mID].forward1)

alignmentMat[i][k]+=(curr->m->f[j][k] * weights[currCount]);

else

alignmentMat[i][k]+=(rev->f[j][k] * weights[currCount]);

}