conservation_colorer.cpp

来自「ncbi源码」· C++ 代码 · 共 641 行 · 第 1/2 页

    typedef map < const UngappedAlignedBlock * , FloatVector > BlockRowScores;
    BlockRowScores blockFitScores, blockZFitScores, blockRowFitScores;
    float minBlockFit, maxBlockFit, minBlockZFit, maxBlockZFit, minBlockRowFit, maxBlockRowFit;

    BlockMap::const_iterator b, be = blocks.end();
    for (b=blocks.begin(); b!=be; ++b) {
        blockFitScores[b->first].resize(nRows, 0.0f);

        for (int blockColumn=0; blockColumn<b->first->width; ++blockColumn) {
            profileColumn = b->second[blockColumn];
            ColumnProfile& profile = alignmentProfile[profileColumn];
            pe = profile.end();

            // fit scores
            for (p=profile.begin(); p!=pe; ++p) {
                int& fit = fitScores[profileColumn][p->first];
                fit = GetPSSMScore(alignment->GetPSSM(),
                    p->first, b->first->GetRangeOfRow(0)->from + blockColumn);
                if (blockColumn == 0 && b == blocks.begin() && p == profile.begin()) {
                    minFit = maxFit = fit;
                } else {
                    if (fit < minFit) minFit = fit;
                    else if (fit > maxFit) maxFit = fit;
                }
            }

            // add up residue fit scores to get block fit scores
            for (row=0; row<nRows; ++row) {
                char ch = ScreenResidueCharacter(b->first->GetCharacterAt(blockColumn, row));
                blockFitScores[b->first][row] += fitScores[profileColumn][ch];
            }
        }

        // find average/min/max block fit
        float average = 0.0f;
        for (row=0; row<nRows; ++row) {
            float& score = blockFitScores[b->first][row];
            score /= b->first->width;   // average fit score across the block for this row
            average += score;
            if (row == 0 && b == blocks.begin()) {
                minBlockFit = maxBlockFit = score;
            } else {
                if (score < minBlockFit) minBlockFit = score;
                else if (score > maxBlockFit) maxBlockFit = score;
            }
        }
        average /= nRows;

        // calculate block Z scores from block fit scores
        if (nRows >= 2) {
            // calculate standard deviation of block fit score over all rows of this block
            float stdDev = 0.0f;
            for (row=0; row<nRows; ++row)
                stdDev += (blockFitScores[b->first][row] - average) *
                          (blockFitScores[b->first][row] - average);
            stdDev /= nRows - 1;
            stdDev = sqrt(stdDev);
            if (stdDev > 1e-10) {
                // calculate Z scores for each row
                blockZFitScores[b->first].resize(nRows);
                for (row=0; row<nRows; ++row)
                    blockZFitScores[b->first][row] = (blockFitScores[b->first][row] - average) / stdDev;
            }
        }
    }

    // calculate row fit scores based on Z-scores for each block across a given row
    if (blocks.size() >= 2) {
        for (b=blocks.begin(); b!=be; ++b)
            blockRowFitScores[b->first].resize(nRows, kMin_Float);

        // calculate row average, standard deviation, and Z-scores
        for (row=0; row<nRows; ++row) {
            float average = 0.0f;
            for (b=blocks.begin(); b!=be; ++b)
                average += blockFitScores[b->first][row];
            average /= blocks.size();
            float stdDev = 0.0f;
            for (b=blocks.begin(); b!=be; ++b)
                stdDev += (blockFitScores[b->first][row] - average) *
                          (blockFitScores[b->first][row] - average);
            stdDev /= blocks.size() - 1;
            stdDev = sqrt(stdDev);
            if (stdDev > 1e-10) {
                for (b=blocks.begin(); b!=be; ++b)
                    blockRowFitScores[b->first][row] = (blockFitScores[b->first][row] - average) / stdDev;
            }
        }
    }

    // now assign colors
    double scale;
    fitColors.resize(nRows * nColumns);
    for (profileColumn=0; profileColumn<nColumns; ++profileColumn) {
        // fit
        CharIntMap::const_iterator c, ce = fitScores[profileColumn].end();
        for (c=fitScores[profileColumn].begin(); c!=ce; ++c) {
            if (maxFit == minFit)
                scale = 1.0;
            else
                scale = 1.0 * (c->second - minFit) / (maxFit - minFit);
            fitColors[profileColumn][c->first] = GlobalColors()->Get(Colors::eConservationMap, scale);
        }
    }

    // block fit
    blockFitColors.clear();
    for (b=blocks.begin(); b!=be; ++b) {
        blockFitColors[b->first].resize(nRows);
        for (row=0; row<nRows; ++row) {
            if (maxBlockFit == minBlockFit)
                scale = 1.0;
            else
                scale = 1.0 * (blockFitScores[b->first][row] - minBlockFit) / (maxBlockFit - minBlockFit);
            blockFitColors[b->first][row] = GlobalColors()->Get(Colors::eConservationMap, scale);
        }
    }

    // block Z fit
    blockZFitColors.clear();
    for (b=blocks.begin(); b!=be; ++b) {
        blockZFitColors[b->first].resize(nRows, Vector(0,0,0));
        if (blockZFitScores.find(b->first) != blockZFitScores.end()) {  // if this column has scores
            for (row=0; row<nRows; ++row) {                             // normalize colors per column
                float zScore = blockZFitScores[b->first][row];
                if (row == 0) {
                    minBlockZFit = maxBlockZFit = zScore;
                } else {
                    if (zScore < minBlockZFit) minBlockZFit = zScore;
                    else if (zScore > maxBlockZFit) maxBlockZFit = zScore;
                }
            }
            for (row=0; row<nRows; ++row) {
                if (maxBlockZFit == minBlockZFit)
                    scale = 1.0;
                else
                    scale = 1.0 * (blockZFitScores[b->first][row] - minBlockZFit) /
                                        (maxBlockZFit - minBlockZFit);
                blockZFitColors[b->first][row] = GlobalColors()->Get(Colors::eConservationMap, scale);
            }
        }
    }

    // block row fit
    blockRowFitColors.clear();
    for (b=blocks.begin(); b!=be; ++b)
        blockRowFitColors[b->first].resize(nRows, Vector(0,0,0));
    if (blocks.size() >= 2) {
        for (row=0; row<nRows; ++row) {
            if (blockRowFitScores.begin()->second[row] != kMin_Float) { // if this row has fit scores
                for (b=blocks.begin(); b!=be; ++b) {                    // normalize colors per row
                    float zScore = blockRowFitScores[b->first][row];
                    if (b == blocks.begin()) {
                        minBlockRowFit = maxBlockRowFit = zScore;
                    } else {
                        if (zScore < minBlockRowFit) minBlockRowFit = zScore;
                        else if (zScore > maxBlockRowFit) maxBlockRowFit = zScore;
                    }
                }
                for (b=blocks.begin(); b!=be; ++b) {
                    if (maxBlockRowFit == minBlockRowFit)
                        scale = 1.0;
                    else
                        scale = 1.0 * (blockRowFitScores[b->first][row] - minBlockRowFit) /
                                            (maxBlockRowFit - minBlockRowFit);
                    blockRowFitColors[b->first][row] = GlobalColors()->Get(Colors::eConservationMap, scale);
                }
            }
        }
    }

    fitColorsCurrent = true;
}

void ConservationColorer::GetProfileIndexAndResidue(
    const UngappedAlignedBlock *block, int blockColumn, int row,
    int *profileIndex, char *residue) const
{
    BlockMap::const_iterator b = blocks.find(block);
    *profileIndex = b->second[blockColumn];
    *residue = ScreenResidueCharacter(b->first->GetCharacterAt(blockColumn, row));
}

void ConservationColorer::Clear(void)
{
    nColumns = 0;
    blocks.clear();
    FreeColors();
}

void ConservationColorer::FreeColors(void)
{
    alignmentProfile.clear();
    identities.clear();
    varietyColors.clear();
    weightedVarietyColors.clear();
    informationContentColors.clear();
    fitColors.clear();
    blockFitColors.clear();
    blockZFitColors.clear();
    blockRowFitColors.clear();
    basicColorsCurrent = fitColorsCurrent = false;
}

END_SCOPE(Cn3D)


/*
* ---------------------------------------------------------------------------
* $Log: conservation_colorer.cpp,v $
* Revision 1000.3  2004/06/01 18:28:29  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.35
*
* Revision 1.35  2004/05/21 21:41:39  gorelenk
* Added PCH ncbi_pch.hpp
*
* Revision 1.34  2004/03/15 18:23:01  thiessen
* prefer prefix vs. postfix ++/-- operators
*
* Revision 1.33  2004/02/19 17:04:53  thiessen
* remove cn3d/ from include paths; add pragma to disable annoying msvc warning
*
* Revision 1.32  2003/12/11 19:35:54  thiessen
* use built-in blosum62 matrix
*
* Revision 1.31  2003/11/04 18:09:17  thiessen
* rearrange headers for OSX build
*
* Revision 1.30  2003/02/13 18:31:39  thiessen
* separate basic from fit color calculation
*
* Revision 1.29  2003/02/12 15:36:14  thiessen
* clean up fit scoring code
*
* Revision 1.28  2003/02/06 16:39:53  thiessen
* add block row fit coloring
*
* Revision 1.27  2003/02/03 19:20:03  thiessen
* format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
*
* Revision 1.26  2003/01/31 17:18:58  thiessen
* many small additions and changes...
*
* Revision 1.25  2003/01/30 14:00:23  thiessen
* add Block Z Fit coloring
*
* Revision 1.24  2003/01/28 21:07:56  thiessen
* add block fit coloring algorithm; tweak row dragging; fix style bug
*
* Revision 1.23  2002/07/23 16:40:39  thiessen
* remove static decl
*
* Revision 1.22  2002/07/23 15:46:50  thiessen
* print out more BLAST info; tweak save file name
*
* Revision 1.21  2002/04/27 16:32:12  thiessen
* fix small leaks/bugs found by BoundsChecker
*
* Revision 1.20  2002/03/01 15:47:10  thiessen
* fix bug in fit coloring for non-standard residues
*
* Revision 1.19  2002/02/13 19:45:29  thiessen
* Fit coloring by info content contribution
*
* Revision 1.18  2001/09/27 15:37:58  thiessen
* decouple sequence import and BLAST
*
* Revision 1.17  2001/08/24 00:41:35  thiessen
* tweak conservation colors and opengl font handling
*
* Revision 1.16  2001/08/13 22:30:59  thiessen
* add structure window mouse drag/zoom; add highlight option to render settings
*
* Revision 1.15  2001/06/04 14:58:00  thiessen
* add proximity sort; highlight sequence on browser launch
*
* Revision 1.14  2001/06/02 17:22:45  thiessen
* fixes for GCC
*
* Revision 1.13  2001/05/31 18:47:08  thiessen
* add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
*
* Revision 1.12  2001/05/25 01:38:16  thiessen
* minor fixes for compiling on SGI
*
* Revision 1.11  2001/05/22 22:37:13  thiessen
* check data registry path
*
* Revision 1.10  2001/05/22 19:09:30  thiessen
* many minor fixes to compile/run on Solaris/GTK
*
* Revision 1.9  2001/03/22 00:33:16  thiessen
* initial threading working (PSSM only); free color storage in undo stack
*
* Revision 1.8  2001/02/22 00:30:06  thiessen
* make directories global ; allow only one Sequence per StructureObject
*
* Revision 1.7  2001/02/15 00:03:54  thiessen
* print out total information content
*
* Revision 1.6  2001/02/13 20:33:49  thiessen
* add information content coloring
*
* Revision 1.5  2000/11/30 15:49:38  thiessen
* add show/hide rows; unpack sec. struc. and domain features
*
* Revision 1.4  2000/10/16 20:03:07  thiessen
* working block creation
*
* Revision 1.3  2000/10/16 14:25:48  thiessen
* working alignment fit coloring
*
* Revision 1.2  2000/10/04 17:41:30  thiessen
* change highlight color (cell background) handling
*
* Revision 1.1  2000/09/20 22:24:48  thiessen
* working conservation coloring; split and center unaligned justification
*
*/