su_alignment_set.cpp

ncbi源码

    multiple->GetUngappedAlignedBlocks(&blocks);

    // create Seq-aligns; if there's only one row (the master), then cheat and create an alignment
    // of the master with itself, because asn data doesn't take well to single-row "alignment"
    if (multiple->NRows() > 1) {
        for (unsigned int row=1; row<multiple->NRows(); ++row, ++sa) {
            sa->Reset(CreatePairwiseSeqAlignFromMultipleRow(multiple, blocks,
                (rowOrder ? (*rowOrder)[row] : row)));
        }
    } else
        sa->Reset(CreatePairwiseSeqAlignFromMultipleRow(multiple, blocks, 0));

    auto_ptr<AlignmentSet> newAlignmentSet;
    try {
        newAlignmentSet.reset(new AlignmentSet(newSeqAnnots, sequenceSet));
    } catch (CException& e) {
        ERROR_MESSAGE(
            "AlignmentSet::CreateFromMultiple() - failed to create AlignmentSet from new asn object: "
            << e.GetMsg());
        return NULL;
    }

    *newAsnAlignmentData = newSeqAnnots;
    return newAlignmentSet.release();
}


///// MasterSlaveAlignment methods /////

MasterSlaveAlignment::MasterSlaveAlignment(const ncbi::objects::CSeq_align& seqAlign,
    const Sequence *masterSequence, const SequenceSet& sequenceSet) :
        m_master(masterSequence), m_slave(NULL)
{
    // resize alignment and block vector
    m_masterToSlave.resize(m_master->Length(), eUnaligned);
    m_blockStructure.resize(m_master->Length(), eUnaligned);

    // find slave sequence for this alignment, and order (master or slave first)
    const CSeq_id& frontSeqId = seqAlign.GetSegs().IsDendiag() ?
        seqAlign.GetSegs().GetDendiag().front()->GetIds().front().GetObject() :
        seqAlign.GetSegs().GetDenseg().GetIds().front().GetObject();
    const CSeq_id& backSeqId = seqAlign.GetSegs().IsDendiag() ?
        seqAlign.GetSegs().GetDendiag().front()->GetIds().back().GetObject() :
        seqAlign.GetSegs().GetDenseg().GetIds().back().GetObject();

    bool masterFirst = true;
    SequenceSet::SequenceList::const_iterator s, se = sequenceSet.m_sequences.end();
    for (s=sequenceSet.m_sequences.begin(); s!=se; ++s) {
        if (m_master->MatchesSeqId(frontSeqId) &&
            (*s)->MatchesSeqId(backSeqId)) {
            break;
        } else if ((*s)->MatchesSeqId(frontSeqId) &&
                   m_master->MatchesSeqId(backSeqId)) {
            masterFirst = false;
            break;
        }
    }
    if (s == se) {
        ERROR_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - couldn't find matching sequences; "
            << "both " << frontSeqId.AsFastaString() << " and "
            << backSeqId.AsFastaString() << " must be in the sequence list for this file!");
        THROW_MESSAGE("couldn't find matching sequences");
    } else {
        m_slave = *s;
    }

    int masterRes, slaveRes, blockNum = 0;
	unsigned int i;

    // unpack dendiag alignment
    if (seqAlign.GetSegs().IsDendiag()) {

        CSeq_align::C_Segs::TDendiag::const_iterator d , de = seqAlign.GetSegs().GetDendiag().end();
        for (d=seqAlign.GetSegs().GetDendiag().begin(); d!=de; ++d, ++blockNum) {
            const CDense_diag& block = d->GetObject();

            if (block.GetDim() != 2 || block.GetIds().size() != 2 || block.GetStarts().size() != 2)
                THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - incorrect dendiag block dimensions");

            // make sure identities of master and slave sequences match in each block
            if ((masterFirst &&
                    (!m_master->MatchesSeqId(block.GetIds().front().GetObject()) ||
                     !m_slave->MatchesSeqId(block.GetIds().back().GetObject()))) ||
                (!masterFirst &&
                    (!m_master->MatchesSeqId(block.GetIds().back().GetObject()) ||
                     !m_slave->MatchesSeqId(block.GetIds().front().GetObject()))))
                THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - mismatched Seq-id in dendiag block");

            // finally, actually unpack the data into the alignment vector
            for (i=0; i<block.GetLen(); ++i) {
                if (masterFirst) {
                    masterRes = block.GetStarts().front() + i;
                    slaveRes = block.GetStarts().back() + i;
                } else {
                    masterRes = block.GetStarts().back() + i;
                    slaveRes = block.GetStarts().front() + i;
                }
                if ((unsigned) masterRes >= m_master->Length() || (unsigned) slaveRes >= m_slave->Length())
                    THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - seqloc in dendiag block > length of sequence!");
                m_masterToSlave[masterRes] = slaveRes;
                m_blockStructure[masterRes] = blockNum;
            }
        }
    }

    // unpack denseg alignment
    else if (seqAlign.GetSegs().IsDenseg()) {

        const CDense_seg& block = seqAlign.GetSegs().GetDenseg();

        if (!block.IsSetDim() && block.GetDim() != 2 ||
                block.GetIds().size() != 2 ||
                block.GetStarts().size() != ((unsigned int) 2 * block.GetNumseg()) ||
                block.GetLens().size() != ((unsigned int) block.GetNumseg()))
            THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - incorrect denseg block dimension");

        // make sure identities of master and slave sequences match in each block
        if ((masterFirst &&
                (!m_master->MatchesSeqId(block.GetIds().front().GetObject()) ||
                 !m_slave->MatchesSeqId(block.GetIds().back().GetObject()))) ||
            (!masterFirst &&
                (!m_master->MatchesSeqId(block.GetIds().back().GetObject()) ||
                 !m_slave->MatchesSeqId(block.GetIds().front().GetObject()))))
            THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - mismatched Seq-id in denseg block");

        // finally, actually unpack the data into the alignment vector
        CDense_seg::TStarts::const_iterator starts = block.GetStarts().begin();
        CDense_seg::TLens::const_iterator lens, le = block.GetLens().end();
        for (lens=block.GetLens().begin(); lens!=le; ++lens) {
            if (masterFirst) {
                masterRes = *(starts++);
                slaveRes = *(starts++);
            } else {
                slaveRes = *(starts++);
                masterRes = *(starts++);
            }
            if (masterRes != -1 && slaveRes != -1) { // skip gaps
                if ((masterRes + *lens - 1) >= m_master->Length() ||
                        (slaveRes + *lens - 1) >= m_slave->Length())
                    THROW_MESSAGE("MasterSlaveAlignment::MasterSlaveAlignment() - "
                        "seqloc in denseg block > length of sequence!");
                for (i=0; i<*lens; ++i) {
                    m_masterToSlave[masterRes + i] = slaveRes + i;
                    m_blockStructure[masterRes + i] = blockNum;
                }
                ++blockNum; // a "block" of a denseg is an aligned (non-gap) segment
            }
        }
    }
}

END_SCOPE(struct_util)


/*
* ---------------------------------------------------------------------------
* $Log: su_alignment_set.cpp,v $
* Revision 1000.0  2004/06/01 18:14:18  gouriano
* PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.6
*
* Revision 1.6  2004/05/26 14:45:11  gorelenk
* UNALIGNED->eUnaligned
*
* Revision 1.5  2004/05/26 01:58:05  ucko
* Move #include <corelib/ncbi_limits.hpp> to su_alignment_set.hpp.
*
* Revision 1.4  2004/05/25 21:23:03  ucko
* Remove definition of MasterSlaveAlignment::UNALIGNED (now part of an enum)
*
* Revision 1.3  2004/05/25 16:12:30  thiessen
* fix GCC warnings
*
* Revision 1.2  2004/05/25 15:52:17  thiessen
* add BlockMultipleAlignment, IBM algorithm
*
* Revision 1.1  2004/05/24 23:04:05  thiessen
* initial checkin
*
*/