alnvwr.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: alnvwr.cpp,v $
 * PRODUCTION Revision 1000.2  2004/06/01 19:40:58  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.22
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: alnvwr.cpp,v 1000.2 2004/06/01 19:40:58 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================
*
* Author:  Kamen Todorov, NCBI
*
* File Description:
*   Various alignment viewers. Demonstration of CAlnMap/CAlnVec usage.
*
* ===========================================================================
*/
#include <ncbi_pch.hpp>
#include <corelib/ncbiapp.hpp>
#include <corelib/ncbiargs.hpp>
#include <corelib/ncbienv.hpp>

#include <serial/iterator.hpp>
#include <serial/objistr.hpp>
#include <serial/objostr.hpp>
#include <serial/serial.hpp>

#include <objects/seq/Bioseq.hpp>
#include <objects/seqloc/Seq_id.hpp>
#include <objects/seqloc/Textseq_id.hpp>
#include <objects/seqset/Seq_entry.hpp>

#include <objects/seqalign/Seq_align.hpp>
#include <objects/seqalign/Seq_align_set.hpp>
#include <objects/seqalign/Std_seg.hpp>
#include <objects/seq/Seq_annot.hpp>
#include <objects/submit/Seq_submit.hpp>

#include <objtools/data_loaders/genbank/gbloader.hpp>
#include <objmgr/object_manager.hpp>
#include <objmgr/scope.hpp>
#include <objmgr/seq_vector.hpp>

#include <objtools/alnmgr/alnmix.hpp>

USING_SCOPE(ncbi);
USING_SCOPE(objects);

void LogTime(const string& s)
{

    static time_t prev_t;
    time_t        t = time(0);

    if (prev_t==0) {
        prev_t=t;
    }
    
    cout << s << " " << (int)(t-prev_t) << endl;
}

class CAlnMgrTestApp : public CNcbiApplication
{
    virtual void     Init(void);
    virtual int      Run(void);
    void             LoadDenseg(void);
    void             View1();
    void             View2(int screen_width);
    void             View3(int screen_width);
    void             View4(int screen_width);
    void             View5();
    void             View6();
    void             View7();
    void             View8(int aln_pos);
    void             View9(int row0, int row1);
    void             GetSeqPosFromAlnPosDemo();
private:
    CRef<CObjectManager> m_ObjMgr;
    CRef<CScope>  m_Scope;
    CRef<CAlnVec> m_AV;
};

void CAlnMgrTestApp::Init(void)
{
    // Create command-line argument descriptions class
    auto_ptr<CArgDescriptions> arg_desc(new CArgDescriptions);

    // Specify USAGE context
    arg_desc->SetUsageContext
        (GetArguments().GetProgramBasename(),
         "Alignment manager demo program");

    // Describe the expected command-line arguments
    arg_desc->AddDefaultKey
        ("in", "InputFile",
         "Name of file to read the Dense-seg from (standard input by default)",
         CArgDescriptions::eInputFile, "-", CArgDescriptions::fPreOpen);

    arg_desc->AddOptionalKey
        ("se_in", "SeqEntryInputFile",
         "An optional Seq-entry file to load a local top level seq entry from.",
         CArgDescriptions::eInputFile, CArgDescriptions::fPreOpen);

    arg_desc->AddOptionalKey
        ("log", "LogFile",
         "Name of log file to write to",
         CArgDescriptions::eOutputFile, CArgDescriptions::fPreOpen);

    arg_desc->AddOptionalKey
        ("a", "AnchorRow",
         "Anchor row (zero based)",
         CArgDescriptions::eInteger);

    arg_desc->AddKey
        ("v", "",
         "View format:\n"
         "1. CSV table\n"
         "2. Popset style using chunks\n"
         "3. Popset style\n"
         "4. Popset style speed optimized\n"
         "5. Print segments\n"
         "6. Print chunks\n"
         "7. Alternative ways to get sequence\n"
         "8. Demonstrate obtaining column vector in two alternative ways.\n"
         "   (Use numeric param n to choose alignment position)\n"
         "9. Print relative residue index mapping for two rows.\n"
         "   (Use row0 and row1 params to choose the rows)\n",
         CArgDescriptions::eInteger);

    arg_desc->AddDefaultKey
        ("w", "ScreenWidth",
         "Screen width for some of the viewers",
         CArgDescriptions::eInteger, "60");

    arg_desc->AddDefaultKey
        ("n", "Number",
         "Generic Numeric Parameter, used by some viewers",
         CArgDescriptions::eInteger, "0");

    arg_desc->AddDefaultKey
        ("row0", "Row0",
         "Generic Row Parameter, used by some viewers",
         CArgDescriptions::eInteger, "0");

    arg_desc->AddDefaultKey
        ("row1", "Row1",
         "Generic Row Parameter, used by some viewers",
         CArgDescriptions::eInteger, "1");

    arg_desc->AddDefaultKey
        ("cf", "GetChunkFlags",
         "Flags for GetChunks (CAlnMap::TGetChunkFlags)",
         CArgDescriptions::eInteger, "0");

    // Setup arg.descriptions for this application
    SetupArgDescriptions(arg_desc.release());
}


void CAlnMgrTestApp::LoadDenseg(void)
{
    CArgs args = GetArgs();

    CNcbiIstream& is = args["in"].AsInputFile();
    bool done = false;
    
    string asn_type;
    {{
        auto_ptr<CObjectIStream> in
            (CObjectIStream::Open(eSerial_AsnText, is));

        asn_type = in->ReadFileHeader();
        in->Close();
        is.seekg(0);
    }}

    auto_ptr<CObjectIStream> in
        (CObjectIStream::Open(eSerial_AsnText, is));
    
    //create scope
    {{
        m_ObjMgr = new CObjectManager;
        
        m_ObjMgr->RegisterDataLoader
            (*new CGBDataLoader("ID", NULL, 2),
             CObjectManager::eDefault);

        m_Scope = new CScope(*m_ObjMgr);
        m_Scope->AddDefaults();
    }}

    if (asn_type == "Dense-seg") {
        CRef<CDense_seg> ds(new CDense_seg);
        *in >> *ds;
        m_AV = new CAlnVec(*ds, *m_Scope);
    } else if (asn_type == "Seq-submit") {
        CRef<CSeq_submit> ss(new CSeq_submit);
        *in >> *ss;
        CTypesIterator i;
        CType<CDense_seg>::AddTo(i);
        CType<CSeq_entry>::AddTo(i);
        int tse_cnt = 0;
        for (i = Begin(*ss); i; ++i) {
            if (CType<CDense_seg>::Match(i)) {
                m_AV = new CAlnVec(*(CType<CDense_seg>::Get(i)), *m_Scope);
            } else if (CType<CSeq_entry>::Match(i)) {
                if ( !(tse_cnt++) ) {
                    m_Scope->AddTopLevelSeqEntry
                        (*(CType<CSeq_entry>::Get(i)));
                }
            }
        }
    } else {
        cerr << "Cannot read: " << asn_type;
        return;
    }

    if ( args["se_in"] ) {
        CNcbiIstream& se_is = args["se_in"].AsInputFile();
        bool done = false;
    
        string se_asn_type;
        {{
            auto_ptr<CObjectIStream> se_in
                (CObjectIStream::Open(eSerial_AsnText, se_is));
            
            se_asn_type = se_in->ReadFileHeader();
            se_in->Close();
            se_is.seekg(0);
        }}
        
        auto_ptr<CObjectIStream> se_in
            (CObjectIStream::Open(eSerial_AsnText, se_is));
        
        if (se_asn_type == "Seq-entry") {
            CRef<CSeq_entry> se (new CSeq_entry);
            *se_in >> *se;
            m_Scope->AddTopLevelSeqEntry(*se);
        } else {
            cerr << "se_in only accepts a Seq-entry asn text file.";
            return;
        }
    }
}

void CAlnMgrTestApp::View1()
{
    cout << ",";
    for (int seg=0; seg<m_AV->GetNumSegs(); seg++) {
        cout << "," << m_AV->GetLen(seg) << ",";
    }
    cout << endl;
    for (int row=0; row<m_AV->GetNumRows(); row++) {
        cout << row << ",";
        for (int seg=0; seg<m_AV->GetNumSegs(); seg++) {
            cout << m_AV->GetStart(row, seg) << "," 
                 << m_AV->GetStop(row, seg) << ",";
        }
        cout << endl;
    }
}

void CAlnMgrTestApp::View2(int screen_width)
{
    int aln_pos = 0;
    CAlnMap::TSignedRange rng;

    do {
        // create range
        rng.Set(aln_pos, aln_pos + screen_width - 1);

        string aln_seq_str;
        aln_seq_str.reserve(screen_width + 1);
        // for each sequence
        for (CAlnMap::TNumrow row = 0; row < m_AV->GetNumRows(); row++) {
            cout << m_AV->GetSeqId(row).AsFastaString()
                 << "\t" 
                 << m_AV->GetSeqPosFromAlnPos(row, rng.GetFrom(),
                                              CAlnMap::eLeft)
                 << "\t"
                 << m_AV->GetAlnSeqString(aln_seq_str, row, rng)
                 << "\t"
                 << m_AV->GetSeqPosFromAlnPos(row, rng.GetTo(),
                                              CAlnMap::eLeft)
                 << endl;
        }
        cout << endl;
        aln_pos += screen_width;
    } while (aln_pos < m_AV->GetAlnStop());
}


void CAlnMgrTestApp::View3(int screen_width)
{
    TSeqPos aln_len = m_AV->GetAlnStop() + 1;
    const CAlnMap::TNumrow nrows = m_AV->GetNumRows();
    const CAlnMap::TNumseg nsegs = m_AV->GetNumSegs();
    const CDense_seg::TStarts& starts = m_AV->GetDenseg().GetStarts();
    const CDense_seg::TLens& lens = m_AV->GetDenseg().GetLens();

    vector<string> buffer(nrows);
    for (CAlnMap::TNumrow row = 0; row < nrows; row++) {

        // allocate space for the row
        buffer[row].reserve(aln_len + 1);
        string buff;

        int seg, pos, left_seg = -1, right_seg = -1;
        TSignedSeqPos start;
        TSeqPos len;

        // determine the ending right seg
        for (seg = nsegs - 1, pos = seg * nrows + row;
             seg >= 0; --seg, pos -= nrows) {
            if (starts[pos] >= 0) {
                right_seg = seg;
                break;
            }
        }

        for (seg = 0, pos = row;  seg < nsegs; ++seg, pos += nrows) {
            len = lens[seg];
            if ((start = starts[pos]) >= 0) {

                left_seg = seg; // ending left seg is at most here

                m_AV->GetSeqString(buff, row, start, start + len - 1);
                buffer[row] += buff;
            } else {
                // add appropriate number of gap/end chars
                char* ch_buff = new char[len+1];
                char fill_ch;
                if (left_seg < 0  ||  seg > right_seg  &&  right_seg > 0) {
                    fill_ch = m_AV->GetEndChar();
                } else {
                    fill_ch = m_AV->GetGapChar(row);
                }
                memset(ch_buff, fill_ch, len);
                ch_buff[len] = 0;
                buffer[row] += ch_buff;
                delete[] ch_buff;
            }
        }
    }

    TSeqPos pos = 0;
    do {
        for (CAlnMap::TNumrow row = 0; row < nrows; row++) {
            cout << m_AV->GetSeqId(row).AsFastaString()
                 << "\t"
                 << m_AV->GetSeqPosFromAlnPos(row, pos, CAlnMap::eLeft)
                 << "\t"
                 << buffer[row].substr(pos, screen_width)
                 << "\t"
                 << m_AV->GetSeqPosFromAlnPos(row, pos + screen_width - 1,
                                              CAlnMap::eLeft)
                 << endl;
        }
        cout << endl;