blast_psi_cxx.cpp

ncbi源码

        } else if (subject_offset == GAP_IN_ALIGNMENT) {

            // gap in subject, initialize appropriately
            for (TSeqPos i = 0; i < lengths[segmt_idx]; i++) {
                PsiDesc& pos_desc =
                    m_AlignmentData->desc_matrix[seq_index][query_offset++];
                pos_desc.letter = GAP;
                pos_desc.used = true;
                pos_desc.e_value = kDefaultEvalueForPosition;
            }

        } else {

            // Aligned segments without any gaps
            for (TSeqPos i = 0; i < lengths[segmt_idx]; i++) {
                PsiDesc& pos_desc =
                    m_AlignmentData->desc_matrix[seq_index][query_offset++];
                pos_desc.letter = seq.first.get()[subj_seq_idx++];
                pos_desc.used = true;
                pos_desc.e_value = e_value;
            }
        }

    }

}

// Should be called once per HSP
// @todo merge with blast_objmgr_tools next week
CScoreMatrixBuilder::TSeqPair 
CScoreMatrixBuilder::x_GetSubjectSequence(const CDense_seg& ds, 
                                          CScope& scope) 
{
    ASSERT(ds.GetDim() == 2);
    Uint1* retval = NULL;
    TSeqPos subjlen = 0;                    // length of the return value
    TSeqPos subj_start = kInvalidSeqPos;    // start of subject alignment
    bool subj_start_found = false;
    TSeqPos subj_index = 1;                 // index into starts vector

    const vector<TSignedSeqPos>& starts = ds.GetStarts();
    const vector<TSeqPos>& lengths = ds.GetLens();

    for (int i = 0; i < ds.GetNumseg(); i++) {

        if (starts[subj_index] != (TSignedSeqPos)GAP_IN_ALIGNMENT) {
            if ( !subj_start_found ) {
                subj_start = starts[subj_index];
                subj_start_found = true;
            }
            subjlen += lengths[i];
        }

        subj_index += ds.GetDim();
    }
    ASSERT(subj_start_found);

    CSeq_loc seqloc;
    seqloc.SetInt().SetFrom(subj_start);
    seqloc.SetInt().SetTo(subj_start+subjlen-1);
    seqloc.SetInt().SetStrand(eNa_strand_unknown);
    seqloc.SetInt().SetId(const_cast<CSeq_id&>(*ds.GetIds().back()));

    CBioseq_Handle bh = scope.GetBioseqHandle(seqloc);
    CSeqVector sv = bh.GetSequenceView(seqloc,
                                       CBioseq_Handle::eViewConstructed,
                                       CBioseq_Handle::eCoding_Ncbi);
    sv.SetCoding(CSeq_data::e_Ncbistdaa);
    retval = s_ExtractSequenceFromSeqVector(sv);

    return make_pair(AutoPtr<Uint1, ArrayDeleter<Uint1> >(retval), subjlen);
}

//////////////////////////////////////////////////////////////////////////////
// Static function definitions

// Returns the evalue from this score object
static double
s_GetEvalue(const CScore& score)
{
    string score_type = score.GetId().GetStr();
    if (score.GetValue().IsReal() && 
       (score_type == "e_value" || score_type == "sum_e")) {
        return score.GetValue().GetReal();
    }
    return numeric_limits<double>::max();
}

static Uint1* 
s_ExtractSequenceFromSeqVector(CSeqVector& sv) 
{
    Uint1* retval = NULL;

    try {
        retval = new Uint1[sv.size()];
    } catch (const std::bad_alloc&) {
        ostringstream os;
        os << sv.size();
        string msg = string("Could not allocate ") + os.str() + " bytes";
        NCBI_THROW(CBlastException, eOutOfMemory, msg);
    }
    for (TSeqPos i = 0; i < sv.size(); i++) {
        retval[i] = sv[i];
    }

    return retval;
}

static double 
s_GetLowestEvalue(const CDense_seg::TScores& scores)
{
    double retval = numeric_limits<double>::max();
    double tmp;

    ITERATE(CDense_seg::TScores, i, scores) {
        if ( (tmp = s_GetEvalue(**i)) < retval) {
            retval = tmp;
        }
    }
    return retval;
}

//////////////////////////////////////////////////////////////////////////////
// Debugging code

// Pretty print sequence
static void 
s_DBG_printSequence(const Uint1* seq, TSeqPos len, ostream& out,
                    bool show_markers, TSeqPos chars_per_line)
{
    TSeqPos nlines = len/chars_per_line;

    for (TSeqPos line = 0; line < nlines + 1; line++) {

        // print chars_per_line residues/bases
        for (TSeqPos i = (chars_per_line*line); 
             i < chars_per_line*(line+1) && (i < len); i++) {
            out << AMINOACID_TO_NCBISTDAA[seq[i]];
        }
        out << endl;

        if ( !show_markers )
            continue;

        // print the residue/base markers
        for (TSeqPos i = (chars_per_line*line); 
             i < chars_per_line*(line+1) && (i < len); i++) {
            if (i == 0 || ((i%10) == 0)) {
                out << i;
                ostringstream os;
                os << i;
                TSeqPos marker_length = os.str().size();
                i += (marker_length-1);
            } else {
                out << " ";
            }
        }
        out << endl;
    }
}

std::string
PsiAlignmentData2String(const PsiAlignmentData* alignment)
{
    if ( !alignment ) {
        return "";
    }

    stringstream ss;

    for (TSeqPos i = 0; i < alignment->dimensions->num_seqs+1; i++) {
        ss << "Seq " << setw(3) << i;
        if ( !alignment->use_sequences[i] ) {
            ss << " NOT USED";
        }
        ss << endl;
        for (TSeqPos j = 0; j < alignment->dimensions->query_sz; j++) {
            if ( !alignment->use_sequences[i] ) {
                continue;
            }

            ss << setw(3) << j << " {" 
               << AMINOACID_TO_NCBISTDAA[alignment->desc_matrix[i][j].letter]
               << ",";
            if (alignment->desc_matrix[i][j].used)
                ss << "used";
            else
                ss << "unused";
            ss << "," << alignment->desc_matrix[i][j].e_value
               << "," << alignment->desc_matrix[i][j].extents.left
               << "," << alignment->desc_matrix[i][j].extents.right
               << "}" << endl;
        }
        ss << "*****************************************************" << endl;
    }
    ss << endl;

    // Print the number of matching sequences per column
    ss << "Matching sequences (alignment->match_seqs)" << endl;
    for (TSeqPos i = 0; i < alignment->dimensions->query_sz; i++) { 
        ss << alignment->match_seqs[i] << " "; 
    }
    ss << endl;
    ss << "*****************************************************" << endl;

    // Print the number of distinct residues per position
    ss << "Residue counts in the matrix " << endl;
    for (TSeqPos i = 0; i < alignment->dimensions->query_sz; i++) {
        for (TSeqPos j = 0; j < PSI_ALPHABET_SIZE; j++) {
            ss << alignment->res_counts[i][j] << " ";
        }
    }

    return ss.str();
}

ostream&
operator<<(ostream& out, const CScoreMatrixBuilder& smb)
{
    TSeqPos query_sz = smb.m_AlignmentData->dimensions->query_sz;

    // Print query sequence in IUPAC
    out << "QUERY:\n";
    s_DBG_printSequence(smb.m_AlignmentData->query, query_sz, out);

    // Print description of alignment and associated data
    out << PsiAlignmentData2String(smb.m_AlignmentData.get());
    out << endl;

    return out;
}

// End debugging code
//////////////////////////////////////////////////////////////////////////////

END_SCOPE(blast)
END_NCBI_SCOPE

/*
* ===========================================================================
*
* $Log: blast_psi_cxx.cpp,v $
* Revision 1000.0  2004/06/01 18:13:18  gouriano
* PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.3
*
* Revision 1.3  2004/05/28 17:42:02  camacho
* Fix compiler warning
*
* Revision 1.2  2004/05/28 17:15:43  camacho
* Fix NCBI {BEGIN,END}_SCOPE macro usage, remove warning
*
* Revision 1.1  2004/05/28 16:41:39  camacho
* Initial revision
*
*
* ===========================================================================
*/