cav_alndisplay.cpp

来自「ncbi源码」· C++ 代码 · 共 1,340 行 · 第 1/4 页

int AlignmentDisplay::DumpFASTA(int firstCol, int lastCol, int nColumns,
    bool doLowercase, CNcbiOstream& os) const
{
    if (firstCol < 0 || lastCol >= GetWidth() || firstCol > lastCol || nColumns < 1) {
        ERR_POST(Error << "AlignmentDisplay::DumpFASTA() - nonsensical display region parameters");
        return CAV_ERROR_BAD_PARAMS;
    }

    // output each alignment row
    for (int row=0; row<GetNRows(); ++row) {

        // create title line
        os << '>';
        const Sequence *seq = indexAlnLocToSeqLocRows[row]->sequence;
        bool prevID = false;
        if (seq->gi != Sequence::NOT_SET) {
            os << "gi|" << seq->gi;
            prevID = true;
        }
        if (seq->pdbID.size() > 0) {
            if (prevID) os << '|';
            if (seq->pdbID == "query" || seq->pdbID == "consensus") {
                os << "lcl|" << seq->pdbID;
            } else {
                os << "pdb|" << seq->pdbID;
                if (seq->pdbChain != ' ')
                    os << '|' << (char) seq->pdbChain << " Chain "
                       << (char) seq->pdbChain << ',';
            }
            prevID = true;
        }
        if (seq->accession.size() > 0) {
            if (prevID) os << '|';
            os << seq->accession;
            prevID = true;
        }
        if (seq->description.size() > 0)
            os << ' ' << seq->description;
        os << '\n';

        // split alignment up into "paragraphs", each with nColumns
        int paragraphStart, nParags = 0, i;
        for (paragraphStart=0; (firstCol+paragraphStart)<=lastCol; paragraphStart+=nColumns, ++nParags) {
            for (i=0; i<nColumns && (firstCol+paragraphStart+i)<=lastCol; ++i) {
                char ch = textRows[row]->GetCharAt(firstCol+paragraphStart+i);
                if (!doLowercase) ch = toupper(ch);
                os << ch;
            }
            os << '\n';
        }
    }

    return CAV_SUCCESS;
}

const string& AlignmentDisplay::GetColumnColor(int alnLoc, double conservationThreshhold) const
{
    // standard probabilities (calculated by BLAST using BLOSUM62 - see conservation_colorer.cpp in Cn3D++)
    typedef map < char , double > Char2Double;
    static Char2Double StandardProbabilities;
    if (StandardProbabilities.size() == 0) {
        StandardProbabilities['A'] = 0.07805;
        StandardProbabilities['C'] = 0.01925;
        StandardProbabilities['D'] = 0.05364;
        StandardProbabilities['E'] = 0.06295;
        StandardProbabilities['F'] = 0.03856;
        StandardProbabilities['G'] = 0.07377;
        StandardProbabilities['H'] = 0.02199;
        StandardProbabilities['I'] = 0.05142;
        StandardProbabilities['K'] = 0.05744;
        StandardProbabilities['L'] = 0.09019;
        StandardProbabilities['M'] = 0.02243;
        StandardProbabilities['N'] = 0.04487;
        StandardProbabilities['P'] = 0.05203;
        StandardProbabilities['Q'] = 0.04264;
        StandardProbabilities['R'] = 0.05129;
        StandardProbabilities['S'] = 0.07120;
        StandardProbabilities['T'] = 0.05841;
        StandardProbabilities['V'] = 0.06441;
        StandardProbabilities['W'] = 0.01330;
        StandardProbabilities['X'] = 0;
        StandardProbabilities['Y'] = 0.03216;
    }

    // if this column isn't completely aligned, use plain color
    int row;
    for (row=0; row<GetNRows(); ++row)
        if (!IsAligned(textRows[row]->GetCharAt(alnLoc))) return plainColor;

    // create profile (residue frequencies) for this column
    Char2Double profile;
    Char2Double::iterator p, pe;
    for (row=0; row<GetNRows(); ++row) {
        char ch = toupper(textRows[row]->GetCharAt(alnLoc));
        switch (ch) {
            case 'A': case 'R': case 'N': case 'D': case 'C':
            case 'Q': case 'E': case 'G': case 'H': case 'I':
            case 'L': case 'K': case 'M': case 'F': case 'P':
            case 'S': case 'T': case 'W': case 'Y': case 'V':
                break;
            default:
                ch = 'X'; // make all but natural aa's just 'X'
        }
        if ((p=profile.find(ch)) != profile.end())
            p->second += 1.0/GetNRows();
        else
            profile[ch] = 1.0/GetNRows();
    }

    // check for identity...
    if (conservationThreshhold == SHOW_IDENTITY) {
        if (profile.size() == 1 && profile.begin()->first != 'X')
            return conservedColor;
        else
            return blockColor;
    }

    // ... or calculate information content for this column (calculated in bits -> logs of base 2)
    double information = 0.0;
    pe = profile.end();
    for (p=profile.begin(); p!=pe; ++p) {
        static const double ln2 = log(2.0), threshhold = 0.0001;
        double expFreq = StandardProbabilities[p->first];
        if (expFreq > threshhold) {
            float freqRatio = p->second / expFreq;
            if (freqRatio > threshhold)
                information += p->second * log(freqRatio) / ln2;
        }
    }

    // if conserved, use conservation color
    if (information > conservationThreshhold) return conservedColor;

    // if is unconserved, use block color to show post-IBM block locations
    return blockColor;
}


///// Row class methods /////

void TextRow::InsertGaps(int nGaps, int beforePos)
{
    if (beforePos < 0 || beforePos > Length()) {
        ERR_POST(Error << "TextRow::InsertGaps() - beforePos out of range");
        return;
    }

    chars.insert(beforePos, nGaps, '-');
}

void TextRow::DeleteGaps(int nGaps, int startPos)
{
    if (startPos < 0 || startPos+nGaps-1 > Length()) {
        ERR_POST(Error << "TextRow::DeleteGaps() - startPos out of range");
        return;
    }

    // check to make sure they're all gaps
    for (int i=0; i<nGaps; ++i) {
        if (!IsGap(chars[startPos + i])) {
            ERR_POST(Error << "TextRow::DeleteGaps() - trying to delete non-gap");
            return;
        }
    }

    chars.erase(startPos, nGaps);
}

// find out how many gaps (up to maxGaps) are present from alnLoc to the next aligned
// residue to the right; set startPos to the first gap, nGaps to # gaps starting at startPos
bool TextRow::IsSqueezable(int alnLoc, int *nGaps, int *startPos, int maxGaps) const
{
    if (alnLoc < 0 || alnLoc >= chars.size()) {
        ERR_POST(Error << "TextRow::IsSqueezable() - alnLoc out of range");
        return false;
    }

    // skip unaligned residues
    while (alnLoc < chars.size() && IsUnaligned(chars[alnLoc])) ++alnLoc;
    if (alnLoc == chars.size() || IsAligned(chars[alnLoc])) return false;

    // count gaps
    *startPos = alnLoc;
    for (*nGaps=1, ++alnLoc; alnLoc < chars.size() && IsGap(chars[alnLoc]); (*nGaps)++, ++alnLoc)
        if (*nGaps == maxGaps) break;
    return true;
}


IndexAlnLocToSeqLocRow::IndexAlnLocToSeqLocRow(const Sequence *seq, int length) :
    sequence(seq)
{
    if (length > 0) seqLocs.resize(length, -1);
}

void IndexAlnLocToSeqLocRow::InsertGaps(int nGaps, int beforePos)
{
    if (nGaps <= 0 || seqLocs.size() == 0) return;

    if (beforePos < 0 || beforePos > Length()) {
        ERR_POST(Error << "IndexAlnLocToSeqLocRow::InsertGaps() - beforePos out of range");
        return;
    }

    IntVec::iterator s = seqLocs.begin();
    for (int i=0; i<beforePos; ++i) ++s;
    seqLocs.insert(s, nGaps, -1);
}

void IndexAlnLocToSeqLocRow::ReIndex(const TextRow& textRow)
{
    seqLocs.resize(textRow.Length());
    int seqLoc = 0;
    for (int i=0; i<textRow.Length(); ++i) {
        if (IsGap(textRow.GetCharAt(i)))
            seqLocs[i] = -1;
        else
            seqLocs[i] = seqLoc++;
    }
    if (seqLoc != sequence->sequenceString.size())
        ERR_POST(Error << "IndexAlnLocToSeqLocRow::ReIndex() - wrong sequence length");
}

END_NCBI_SCOPE


/*
* ---------------------------------------------------------------------------
* $Log: cav_alndisplay.cpp,v $
* Revision 1000.3  2004/06/01 19:41:15  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.5
*
* Revision 1.5  2004/05/21 21:42:51  gorelenk
* Added PCH ncbi_pch.hpp
*
* Revision 1.4  2004/03/15 18:51:27  thiessen
* prefer prefix vs. postfix ++/-- operators
*
* Revision 1.3  2003/11/15 13:16:04  thiessen
* fix uid link urls
*
* Revision 1.2  2003/06/02 16:06:41  dicuccio
* Rearranged src/objects/ subtree.  This includes the following shifts:
*     - src/objects/asn2asn --> arc/app/asn2asn
*     - src/objects/testmedline --> src/objects/ncbimime/test
*     - src/objects/objmgr --> src/objmgr
*     - src/objects/util --> src/objmgr/util
*     - src/objects/alnmgr --> src/objtools/alnmgr
*     - src/objects/flat --> src/objtools/flat
*     - src/objects/validator --> src/objtools/validator
*     - src/objects/cddalignview --> src/objtools/cddalignview
* In addition, libseq now includes six of the objects/seq... libs, and libmmdb
* replaces the three libmmdb? libs.
*
* Revision 1.1  2003/03/19 19:04:12  thiessen
* move again
*
* Revision 1.2  2003/03/19 13:25:09  thiessen
* fix log ambiguity
*
* Revision 1.1  2003/03/19 05:33:43  thiessen
* move to src/app/cddalignview
*
* Revision 1.23  2003/02/03 17:52:03  thiessen
* move CVS Log to end of file
*
* Revision 1.22  2003/01/21 18:01:07  thiessen
* add condensed alignment display
*
* Revision 1.21  2002/12/09 13:31:04  thiessen
* use new query.fcgi for genpept links
*
* Revision 1.20  2002/11/08 19:38:11  thiessen
* add option for lowercase unaligned in FASTA
*
* Revision 1.19  2002/04/25 13:14:55  thiessen
* fix range test
*
* Revision 1.18  2002/02/12 13:08:20  thiessen
* annot description optional
*
* Revision 1.17  2002/02/12 12:54:10  thiessen
* feature legend at bottom; annot only where aligned
*
* Revision 1.16  2002/02/08 19:53:17  thiessen
* add annotation to text/HTML displays
*
* Revision 1.15  2001/03/02 01:19:24  thiessen
* add FASTA output
*
* Revision 1.14  2001/02/16 19:18:47  thiessen
* change color scheme again
*
* Revision 1.13  2001/02/15 19:23:44  thiessen
* add identity coloring
*
* Revision 1.12  2001/02/15 18:08:15  thiessen
* change color scheme
*
* Revision 1.11  2001/02/14 16:06:09  thiessen
* add block and conservation coloring to HTML display
*
* Revision 1.10  2001/02/14 03:16:27  thiessen
* fix seqloc markers and right-justification of left tails
*
* Revision 1.9  2001/02/14 00:06:49  thiessen
* filter out consensus
*
* Revision 1.8  2001/01/29 23:55:09  thiessen
* add AlignmentDisplay verification
*
* Revision 1.7  2001/01/29 18:13:33  thiessen
* split into C-callable library + main
*
* Revision 1.6  2001/01/25 00:51:20  thiessen
* add command-line args; can read asn data from stdin
*
* Revision 1.5  2001/01/23 23:17:48  thiessen
* fix uid link problems
*
* Revision 1.4  2001/01/23 20:42:00  thiessen
* add description
*
* Revision 1.3  2001/01/23 17:34:12  thiessen
* add HTML output
*
* Revision 1.2  2001/01/22 15:55:11  thiessen
* correctly set up ncbi namespacing
*
* Revision 1.1  2001/01/22 13:15:23  thiessen
* initial checkin
*
*/