clustermain.cpp

著名的基因芯片表达数据聚类软件。esein编写

    }

    for (i=0;i<2*Rows-1;i++)
    {
        GeneID[i] = "";
    }

    delete ArrayID;
    delete GeneID;

    delete ColumnOrder;
    delete RowOrder;

    StatusBar1->SimpleText = "Done Clustering";
}
//---------------------------------------------------------------------------

/* Average Linkage Hierarchical Clustering Code */

void    TMainForm::Cluster(double **ClusterData, bool **ClusterMask,
        int ClusterRows, int ClusterColumns,
        TStringList *TreeFile, double *Order, double *RowWeight, double *ColumnWeight,
        AnsiString *ID, int DistFunction,
        bool CalculateWeights, double WeightCutoff, double WeightPower,
        TStringList *ClusterOrder, TStatusBar *StatusBar,
        bool ReturnTNode, TNode *TopNode)
{
        int i,j,k,l,m;


        unsigned short *minval;
        int *minpair;
        int *elements;
        AnsiString OutString;
        unsigned short **Dist;
        int **NodeElement;
        double *NodeDistance;
        double *CalcWeight;
        int element1, element2;

        Dist = new unsigned short*[2*ClusterRows-1];
        minval = new unsigned short[2*ClusterRows-1];
        minpair = new int[2*ClusterRows-1];
        elements = new int[2*ClusterRows-1];
        NodeElement = new int*[2*ClusterRows-1];
        NodeDistance = new double[2*ClusterRows-1];
        CalcWeight = new double[2*ClusterRows-1];

        bool SaveDistances = false;

        double WeightDist;

        for (i=0;i<ClusterRows;i++)
        {
            Dist[i] = new unsigned short[2*ClusterRows-1];
        }

        for (i=0;i<2*ClusterRows-1;i++)
        {
            minval[i] = 32769;
            elements[i] = 1;
            NodeElement[i] = new int[2];
            CalcWeight[i] = 1;
            NodeDistance[i] = 0.0;
        }


        bool Centered = false;

        if ((DistFunction == 1) || (DistFunction == 3) )
        {
            Centered = true;
        }

        bool Absolute = false;

        if ((DistFunction == 3) || (DistFunction == 4) )
        {
            Absolute = true;
        }

        bool CalcInferredDistances = false;
        TStringList **MemberList;
        unsigned short **TreeDist;

        if (CalcInferredDistances == true)
        {
            MemberList = new TStringList*[2*ClusterRows-1];
            TreeDist = new unsigned short*[ClusterRows];
            for (i=0;i<2*ClusterRows-1;i++)
            {
                MemberList[i] = new TStringList();
                if (i<ClusterRows)
                {
                    MemberList[i]->Add(AnsiString(i));
                    TreeDist[i] = new unsigned short[ClusterRows];
                }
            }
        }

        int TotalDistance = ClusterRows * (ClusterRows - 1) / 2;

        /* First step is to compute the distance matrix
        This is the slowest step */

        for (i=0;i<ClusterRows;i++)
        {
            Dist[i][i] = 0.0;
            for (j=0;j<i;j++)
            {

                Dist[i][j] = Distance(DistFunction,ClusterData,ClusterMask,ColumnWeight,
                    i,j,ClusterColumns);
                Dist[j][i] = Dist[i][j];

                if (CalculateWeights == true)
                {

                    WeightDist = ( (double) (32768 - Dist[i][j])) /32768.0;
                    WeightDist = max(0.0, (WeightDist - WeightCutoff)/(1.0 - WeightCutoff));
                    WeightDist = pow(WeightDist,WeightPower);

                    CalcWeight[i] += WeightDist;
                    CalcWeight[j] += WeightDist;
                }

                if (Dist[i][j] < minval[i])
                {
                    minval[i] = Dist[i][j];
                    minpair[i] = j;
                }

                if (Dist[i][j] < minval[j])
                {
                    minval[j] = Dist[i][j];
                    minpair[j] = i;
                }
            }
            StatusBar->SimpleText = "Calculating Distances " + AnsiString(i*(i+1)/2) + " of "
                + AnsiString(TotalDistance);
            Application->ProcessMessages();
        }

        TStringList *MinList = new TStringList();

        for (i=0;i<ClusterRows;i++)
        {
            MinList->Add(AnsiString(minval[i]));
        }

        MinList->SaveToFile("minlist.txt");

        delete MinList;

        bool *Active;
        Active = new bool[2*ClusterRows-1];

        for (i=0;i<ClusterRows;i++)
        {
            Active[i] = true;
        }

        for (i=ClusterRows;i<2*ClusterRows-1;i++)
        {
            Active[i] = false;
        }

        unsigned short minminval;

        int min1, min2;

        int nodeindex;

        //int node;

        /* Now we join nodes */

        /* If we are going to return TNode, need to set up TNodes */

        TNode **Nodes;
        if (ReturnTNode)
        {
           Nodes = new TNode*[2*ClusterRows-1];

           for (i=0;i<2*ClusterRows-1;i++)
           {
               Nodes[i] = new TNode();

               if (i < ClusterRows)
               {
                   Nodes[i]->IsNode = false;
                   Nodes[i]->ID = AnsiString(i);
               }
               else
               {
                   Nodes[i]->IsNode = true;
               }
           }
        }

        for (nodeindex=ClusterRows; nodeindex<2*ClusterRows-1;nodeindex++)

        {
            ID[nodeindex] = "NODE" + AnsiString(nodeindex-ClusterRows+1) + "X";

            minminval = 32769;

            for (i=0;i<2*ClusterRows-1;i++)
            {
                if (Active[i] == true)
                {
                    if (minval[i] < minminval)
                    {
                        minminval = minval[i];
                        min1 = i;
                        min2 = minpair[i];
                    }
                }
            }

            /* When we get here min1 and min2 are the elements that are to be joined */




        Active[min1] = false;
        Active[min2] = false;

        if (Order[min1] > Order[min2])
        {
            NodeElement[nodeindex][0] = min1;
            NodeElement[nodeindex][1] = min2;
        }
        else
        {
            NodeElement[nodeindex][1] = min1;
            NodeElement[nodeindex][0] = min2;
        }

        NodeDistance[nodeindex] = ((double) minminval)/16384.0;
        NodeDistance[nodeindex] = max(NodeDistance[nodeindex],NodeDistance[min1]);
        NodeDistance[nodeindex] = max(NodeDistance[nodeindex],NodeDistance[min2]);

        if (ReturnTNode)
        {
               Nodes[nodeindex]->Child1 = Nodes[NodeElement[nodeindex][0]];
               Nodes[nodeindex]->Child2 = Nodes[NodeElement[nodeindex][1]];
               Nodes[nodeindex]->Length = 1.0 - NodeDistance[nodeindex];
        }

        if (CalcInferredDistances == true)
        {
            int elem1,elem2,ielem1,ielem2;

            for (elem1=0;elem1<MemberList[min1]->Count;elem1++)
            {
                for (elem2=0;elem2<MemberList[min2]->Count;elem2++)
                {
                    ielem1 = (MemberList[min1]->Strings[elem1]).ToInt();
                    ielem2 = (MemberList[min2]->Strings[elem2]).ToInt();
                    TreeDist[ielem1][ielem2] = (NodeDistance[nodeindex] * 16384);
                    TreeDist[ielem2][ielem1] = (NodeDistance[nodeindex] * 16384);
                }
            }

            MemberList[nodeindex]->AddStrings(MemberList[min1]);
            MemberList[nodeindex]->AddStrings(MemberList[min2]);

            delete MemberList[min1];
            delete MemberList[min2];
        }

        OutString = ID[nodeindex] + AnsiString("\t") + ID[NodeElement[nodeindex][0]]
                + AnsiString("\t") + ID[NodeElement[nodeindex][1]] + AnsiString("\t")
                + AnsiString(1.0 - NodeDistance[nodeindex]);

        TreeFile->Add(OutString);

        StatusBar->SimpleText = "Generating Node " + AnsiString(nodeindex-ClusterRows+1) + " of " + AnsiString(ClusterRows-1);
        Application->ProcessMessages();

        elements[nodeindex] = elements[min1] + elements[min2];

        for (k=0;k<ClusterColumns;k++)
        {
            if ( (ClusterMask[min1][k] == true) && (ClusterMask[min2][k] == true) )
            {
                ClusterData[nodeindex][k] =
                    (elements[min1] * ClusterData[min1][k] +
                        elements[min2] * ClusterData[min2][k]) /
                    (elements[min1] + elements[min2]);
                ClusterMask[nodeindex][k] = true;
            }
            else if (ClusterMask[min1][k] == true)
            {
                ClusterData[nodeindex][k] = ClusterData[min1][k];
                ClusterMask[nodeindex][k] = true;
            }
            else if (ClusterMask[min2][k] == true)
            {
                ClusterData[nodeindex][k] = ClusterData[min2][k];
                ClusterMask[nodeindex][k] = true;
            }
            else
            {
                ClusterMask[nodeindex][k] = false;
            }
        }

        Order[nodeindex] = (elements[min1]*Order[min1] + elements[min2]*Order[min2])/(elements[min1] + elements[min2]);
        Dist[nodeindex] = new unsigned short[2*ClusterRows-1];

        if (SaveDistances == false)
        {
            delete Dist[min1];
            delete Dist[min2];
        }

        for (l=0;l<2*ClusterRows-1;l++)
        {
            if (Active[l] == true)
            {

                Dist[nodeindex][l] =
                    Distance(DistFunction,ClusterData,ClusterMask,ColumnWeight,
                    nodeindex,l,ClusterColumns);
                Dist[l][nodeindex] = Dist[nodeindex][l];

                if (Dist[nodeindex][l] < minval[nodeindex])
                {
                    minval[nodeindex] = Dist[nodeindex][l];
                    minpair[nodeindex] = l;
                }

                //Rescan for min if minpair was one of the fused elements

                if ( (minpair[l] == min1) || (minpair[l] == min2) )
                {
                    minval[l] = 32769;

                    for (m=0;m<2*ClusterRows-1;m++)
                    {
                        if ( (Active[m] == true) && (l != m) )
                        {
                            if (Dist[l][m] < minval[l])
                            {
                                minval[l] = Dist[l][m];
                                minpair[l] = m;
                            }
                        }
                    }
                }

                //otherwise check new distance to see if it is new min

                else if (Dist[nodeindex][l] < minval[l])
                {
                    minval[l] = Dist[nodeindex][l];
                    minpair[l] = nodeindex;
                }
            }
        }

        Active[nodeindex] = true;
        }

        if (SaveDistances == false)
        {
            delete Dist[2*ClusterRows-2];
        }

        TList *NewList = new TList();
        TList *OutList = new TList();

        int **Elements;
        Elements = new int*[2*ClusterRows-1];
        for (i=0;i<2*ClusterRows-1;i++)
        {
            Elements[i] = new int;
            *Elements[i] = i;
        }

        OutList->Clear();

        OutList->Add((void *)Elements[2*ClusterRows-2]);

        int Position;

        bool Replaced;

        Replaced = true;

        while (Replaced == true)
        {
            Replaced = false;

            for (Position = 0; Position < OutList->Count; Position ++)
            {
                nodeindex =  *(int *) OutList->Items[Position];
                if (nodeindex >= 0)
                {
                    if (nodeindex < ClusterRows)
                    {
                        NewList->Add((void *) Elements[nodeindex]);
                    }
                    else
                    {
                        element1 =  NodeElement[nodeindex][0];
                        element2 =  NodeElement[nodeindex][1];

                        NewList->Add((void *) Elements[element1]);
                        NewList->Add((void *) Elements[element2]);

                        Replaced = true;
                    }
                }
            }

            OutList->Clear();
            for (i=0;i<NewList->Count;i++)
            {
                OutList->Add(NewList->Items[i]);
            }
            NewList->Clear();
        }

        int index;

        ClusterOrder->Clear();

        if (SaveDistances == true)
        {
            if (DistFileDialogBox->Execute())
            {

                TFileStream *DistFile = new TFileStream(DistFileDialogBox->FileName,fmCreate);

                int ii, ij;

                for (i=0;i<OutList->Count;i++)
                {
                    ii = *(int *)OutList->Items[i];
                    for (j=0;j<OutList->Count;j++)
                    {
                        {
                            ij = *(int *)OutList->Items[j];
                            DistFile->Write(&Dist[ii][ij],2);
                        }
                    }
                }
                delete DistFile;

                for (i=0;i<2*ClusterRows-1;i++)
                {
                    delete Dist[i];
                }
            }
        }

        if (CalcInferredDistances == true)
        {
            if (DistFileDialogBox->Execute())
            {

                TFileStream *DistFile = new TFileStream(DistFileDialogBox->FileName,fmCreate);

                int ii, ij;

                for (i=0;i<OutList->Count;i++)
                {
                    ii = *(int *)O