block_multiple_alignment.cpp

ncbi源码

    AlignedBlockMap correspondingNewBlocks;

    for (b=blocks.begin(); b!=be; ++b) {
        if (!(*b)->IsAligned()) continue;
        for (nb=newAlignment->blocks.begin(); nb!=nbe; ++nb) {
            if (!(*nb)->IsAligned()) continue;

            newRange = (*nb)->GetRangeOfRow(0);
            thisRange = (*b)->GetRangeOfRow(0);
            if (newRange->from <= thisRange->from && newRange->to >= thisRange->to) {
                correspondingNewBlocks[dynamic_cast<UngappedAlignedBlock*>(*b)] =
                    dynamic_cast<const UngappedAlignedBlock*>(*nb);
                break;
            }
        }
        if (nb == nbe) return false;    // no corresponding block found
    }

    // add slave sequences from new alignment; also copy scores/status
    SequenceList *modSequences = const_cast<SequenceList*>(sequences);
    int i, nNewRows = newAlignment->sequences->size() - 1;
    modSequences->resize(modSequences->size() + nNewRows);
    rowDoubles.resize(rowDoubles.size() + nNewRows);
    rowStrings.resize(rowStrings.size() + nNewRows);
    for (i=0; i<nNewRows; ++i) {
        (*modSequences)[modSequences->size() + i - nNewRows] = (*(newAlignment->sequences))[i + 1];
        SetRowDouble(NRows() + i - nNewRows, newAlignment->GetRowDouble(i + 1));
        SetRowStatusLine(NRows() + i - nNewRows, newAlignment->GetRowStatusLine(i + 1));
    }

    // now that we know blocks are compatible, add new rows at end of this alignment, containing
    // all rows (except master) from new alignment; only that part of the aligned blocks from
    // the new alignment that intersect with the aligned blocks from this object are added, so
    // that this object's block structure is unchanged

    // resize all aligned blocks
    for (b=blocks.begin(); b!=be; ++b) (*b)->AddRows(nNewRows);

    // set ranges of aligned blocks, and add them to the list
    AlignedBlockMap::const_iterator ab, abe = correspondingNewBlocks.end();
    const Block::Range *thisMaster, *newMaster;
    for (ab=correspondingNewBlocks.begin(); ab!=abe; ++ab) {
        thisMaster = ab->first->GetRangeOfRow(0);
        newMaster = ab->second->GetRangeOfRow(0);
        for (i=0; i<nNewRows; ++i) {
            newRange = ab->second->GetRangeOfRow(i + 1);
            ab->first->SetRangeOfRow(NRows() + i - nNewRows,
                newRange->from + thisMaster->from - newMaster->from,
                newRange->to + thisMaster->to - newMaster->to);
        }
    }

    // delete then recreate the unaligned blocks, in case a merge requires the
    // creation of a new unaligned block
    for (b=blocks.begin(); b!=be; ) {
        if (!(*b)->IsAligned()) {
            BlockList::iterator bb(b);
            ++bb;
            delete *b;
            blocks.erase(b);
            b = bb;
        } else
            ++b;
    }
    InitCache();

    // update this alignment, but leave row scores/status alone
    if (!AddUnalignedBlocks() || !UpdateBlockMapAndColors(false)) {
        ERRORMSG("BlockMultipleAlignment::MergeAlignment() - internal update after merge failed");
        return false;
    }
    return true;
}

int BlockMultipleAlignment::ShowGeometryViolations(bool showGV)
{
    geometryViolations.clear();

    if (!showGV || !GetMaster()->molecule || GetMaster()->molecule->parentSet->isAlphaOnly) {
        showGeometryViolations = false;
        return 0;
    }

    Threader::GeometryViolationsForRow violations;
    int nViolations = alignmentManager->threader->GetGeometryViolations(this, &violations);
    if (violations.size() != NRows()) {
        ERRORMSG("BlockMultipleAlignment::ShowGeometryViolations() - wrong size list");
        showGeometryViolations = false;
        return 0;
    }

    geometryViolations.resize(NRows());
    for (int row=0; row<NRows(); ++row) {
        geometryViolations[row].resize(GetSequenceOfRow(row)->Length(), false);
        Threader::IntervalList::const_iterator i, ie = violations[row].end();
        for (i=violations[row].begin(); i!=ie; ++i)
            for (int l=i->first; l<=i->second; ++l)
                geometryViolations[row][l] = true;
    }

    showGeometryViolations = true;
    return nViolations;
}

CSeq_align * CreatePairwiseSeqAlignFromMultipleRow(const BlockMultipleAlignment *multiple,
    const BlockMultipleAlignment::UngappedAlignedBlockList& blocks, int slaveRow)
{
    if (!multiple || slaveRow < 0 || slaveRow >= multiple->NRows()) {
        ERRORMSG("CreatePairwiseSeqAlignFromMultipleRow() - bad parameters");
        return NULL;
    }

    CSeq_align *seqAlign = new CSeq_align();
    seqAlign->SetType(CSeq_align::eType_partial);
    seqAlign->SetDim(2);

    CSeq_align::C_Segs::TDendiag& denDiags = seqAlign->SetSegs().SetDendiag();
    denDiags.resize((blocks.size() > 0) ? blocks.size() : 1);

    CSeq_align::C_Segs::TDendiag::iterator d, de = denDiags.end();
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b = blocks.begin();
    const Block::Range *range;
    for (d=denDiags.begin(); d!=de; ++d, ++b) {

        CDense_diag *denDiag = new CDense_diag();
        d->Reset(denDiag);
        denDiag->SetDim(2);
        denDiag->SetIds().resize(2);

        // master row
        denDiag->SetIds().front().Reset(multiple->GetSequenceOfRow(0)->CreateSeqId());
        if (blocks.size() > 0) {
            range = (*b)->GetRangeOfRow(0);
            denDiag->SetStarts().push_back(range->from);
        } else
            denDiag->SetStarts().push_back(0);

        // slave row
        denDiag->SetIds().back().Reset(multiple->GetSequenceOfRow(slaveRow)->CreateSeqId());
        if (blocks.size() > 0) {
            range = (*b)->GetRangeOfRow(slaveRow);
            denDiag->SetStarts().push_back(range->from);
        } else
            denDiag->SetStarts().push_back(0);

        // block width
        denDiag->SetLen((blocks.size() > 0) ? (*b)->width : 0);
    }

    return seqAlign;
}

bool BlockMultipleAlignment::MarkBlock(int column)
{
    if (column >= 0 && column < blockMap.size() && blockMap[column].block->IsAligned()) {
        TRACEMSG("marked block for realignment");
        markBlocks[blockMap[column].block] = true;
        return true;
    }
    return false;
}

bool BlockMultipleAlignment::ClearMarks(void)
{
    if (markBlocks.size() == 0) return false;
    markBlocks.clear();
    return true;
}

bool BlockMultipleAlignment::HighlightAlignedColumnsOfMasterRange(int from, int to) const
{
    const Sequence *master = GetMaster();

    // must do one column at a time, rather than range, in case there are inserts wrt the master
    bool anyError = false;
    for (int i=from; i<=to; ++i) {

        // sanity check
        if (i < 0 || i >= master->Length() || !IsAligned(0, i)) {
            WARNINGMSG("Can't highlight alignment at master residue " << (i+1));
            anyError = true;
            // highlight unaligned residues, but master only
            if (i >= 0 && i < master->Length())
                GlobalMessenger()->AddHighlights(GetSequenceOfRow(0), i, i);
            continue;
        }

        // get block and offset
        const Block *block = GetBlock(0, i);
        int blockOffset = i - block->GetRangeOfRow(0)->from;

        // highlight aligned residue in each row
        for (int row=0; row<NRows(); ++row) {
            int slaveIndex = block->GetRangeOfRow(row)->from + blockOffset;
            GlobalMessenger()->AddHighlights(GetSequenceOfRow(row), slaveIndex, slaveIndex);
        }
    }

    return !anyError;
}

int BlockMultipleAlignment::NAlignedBlocks(void) const
{
    int n = 0;
    BlockList::const_iterator b, be = blocks.end();
    for (b=blocks.begin(); b!=be; ++b) if ((*b)->IsAligned()) ++n;
    return n;
}

int BlockMultipleAlignment::GetAlignmentIndex(int row, int seqIndex, eUnalignedJustification justification)
{
    if (row < 0 || row >= NRows() || seqIndex < 0 || seqIndex >= GetSequenceOfRow(row)->Length()) {
        ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - coordinate out of range");
        return -1;
    }

    int alignmentIndex, blockColumn;
    const Block *block = NULL;
    const Block::Range *range;

    for (alignmentIndex=0; alignmentIndex<totalWidth; ++alignmentIndex) {

        // check each block to see if index is in range
        if (block != blockMap[alignmentIndex].block) {
            block = blockMap[alignmentIndex].block;

            range = block->GetRangeOfRow(row);
            if (seqIndex >= range->from && seqIndex <= range->to) {

                // override requested justification for end blocks
                if (block == blocks.back()) // also true if there's a single aligned block
                    justification = eLeft;
                else if (block == blocks.front())
                    justification = eRight;

                // search block columns to find index (inefficient, but avoids having to write
                // inverse functions of Block::GetIndexAt()
                for (blockColumn=0; blockColumn<block->width; ++blockColumn) {
                    if (seqIndex == block->GetIndexAt(blockColumn, row, justification))
                        return alignmentIndex + blockColumn;
                }
                ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - can't find index in block");
                return -1;
            }
        }
    }

    // should never get here...
    ERRORMSG("BlockMultipleAlignment::GetAlignmentIndex() - confused");
    return -1;
}

// creates a SeqAlign from a BlockMultipleAlignment
SeqAlignPtr BlockMultipleAlignment::CreateCSeqAlign(void) const
{
    // one SeqAlign (chained into a linked list) for each slave row
    SeqAlignPtr prevSap = NULL, firstSap = NULL;
    for (int row=1; row<NRows(); ++row) {

        SeqAlignPtr sap = SeqAlignNew();
        if (prevSap) prevSap->next = sap;
        prevSap = sap;
		if (!firstSap) firstSap = sap;

        sap->type = SAT_PARTIAL;
        sap->dim = 2;
        sap->segtype = SAS_DENDIAG;

        DenseDiagPtr prevDd = NULL;
        UngappedAlignedBlockList blocks;
        GetUngappedAlignedBlocks(&blocks);
        UngappedAlignedBlockList::const_iterator b, be = blocks.end();

        for (b=blocks.begin(); b!=be; ++b) {
            DenseDiagPtr dd = DenseDiagNew();
            if (prevDd) prevDd->next = dd;
            prevDd = dd;
            if (b == blocks.begin()) sap->segs = dd;

            dd->dim = 2;
            GetSequenceOfRow(0)->AddCSeqId(&(dd->id), false);      // master
            GetSequenceOfRow(row)->AddCSeqId(&(dd->id), false);    // slave
            dd->len = (*b)->width;

            dd->starts = (Int4Ptr) MemNew(2 * sizeof(Int4));
            const Block::Range *range = (*b)->GetRangeOfRow(0);
            dd->starts[0] = range->from;
            range = (*b)->GetRangeOfRow(row);
            dd->starts[1] = range->from;
        }
    }

	return firstSap;
}


///// UngappedAlignedBlock methods /////

char UngappedAlignedBlock::GetCharacterAt(int blockColumn, int row) const
{
    return (*(parentAlignment->sequences))[row]->sequenceString[GetIndexAt(blockColumn, row)];
}

Block * UngappedAlignedBlock::Clone(const BlockMultipleAlignment *newMultiple) const
{
    UngappedAlignedBlock *copy = new UngappedAlignedBlock(newMultiple);
    const Block::Range *range;
    for (int row=0; row<NSequences(); ++row) {
        range = GetRangeOfRow(row);
        copy->SetRangeOfRow(row, range->from, range->to);
    }
    copy->width = width;
    return copy;
}

void UngappedAlignedBlock::DeleteRow(int row)
{
    RangeList::iterator r = ranges.begin();
    for (int i=0; i<row; ++i) ++r;
    ranges.erase(r);
}

void UngappedAlignedBlock::DeleteRows(vector < bool >& removeRows, int nToRemove)
{
    VectorRemoveElements(ranges, removeRows, nToRemove);
}


///// UnalignedBlock methods /////

int UnalignedBlock::GetIndexAt(int blockColumn, int row,
        BlockMultipleAlignment::eUnalignedJustification justification) const
{
    const Block::Range *range = GetRangeOfRow(row);
    int seqIndex, rangeWidth, rangeMiddle, extraSpace;

    switch (justification) {
        case BlockMultipleAlignment::eLeft:
            seqIndex = range->from + blockColumn;
            break;
        case BlockMultipleAlignment::eRight:
            seqIndex = range->to - width + blockColumn + 1;
            break;
        case BlockMultipleAlignment::eCenter:
            rangeWidth = (range->to - range->from + 1);
            extraSpace = (width - rangeWidth) / 2;
            if (blockColumn < extraSpace || blockColumn >= extraSpace + rangeWidth)
                seqIndex = -1;
            else
                seqIndex = range->from + blockColumn - extraSpace;
            break;
        case BlockMultipleAlignment::eSplit:
            rangeWidth = (range->to - range->from + 1);