seqport_util.cpp

ncbi源码

    m_Ncbi4naRev = InitNcbi4naRev();


    // Initialize map tables

    m_Ncbi2naIupacna = InitMaps(eSeq_code_type_ncbi2na,
                                eSeq_code_type_iupacna);

    m_Ncbi2naNcbi4na = InitMaps(eSeq_code_type_ncbi2na,
                                eSeq_code_type_ncbi4na);

    m_Ncbi4naIupacna = InitMaps(eSeq_code_type_ncbi4na,
                                eSeq_code_type_iupacna);

    m_IupacnaNcbi2na = InitMaps(eSeq_code_type_iupacna,
                                eSeq_code_type_ncbi2na);

    m_IupacnaNcbi4na = InitMaps(eSeq_code_type_iupacna,
                                eSeq_code_type_ncbi4na);

    m_Ncbi4naNcbi2na = InitMaps(eSeq_code_type_ncbi4na,
                                eSeq_code_type_ncbi2na);

    m_IupacaaNcbieaa = InitMaps(eSeq_code_type_iupacaa,
                                eSeq_code_type_ncbieaa);

    m_NcbieaaIupacaa = InitMaps(eSeq_code_type_ncbieaa,
                                eSeq_code_type_iupacaa);

    m_IupacaaNcbistdaa = InitMaps(eSeq_code_type_iupacaa,
                                  eSeq_code_type_ncbistdaa);

    m_NcbieaaNcbistdaa = InitMaps(eSeq_code_type_ncbieaa,
                                  eSeq_code_type_ncbistdaa);

    m_NcbistdaaNcbieaa = InitMaps(eSeq_code_type_ncbistdaa,
                                  eSeq_code_type_ncbieaa);

    m_NcbistdaaIupacaa = InitMaps(eSeq_code_type_ncbistdaa,
                                  eSeq_code_type_iupacaa);

    // Initialize fast conversion tables
    //m_FastNcbi2naIupacna = InitFastNcbi2naIupacna();
    //m_FastNcbi2naNcbi4na = InitFastNcbi2naNcbi4na();
    //m_FastNcbi4naIupacna = InitFastNcbi4naIupacna();
    m_FastIupacnaNcbi2na = InitFastIupacnaNcbi2na();
    m_FastIupacnaNcbi4na = InitFastIupacnaNcbi4na();
    m_FastNcbi4naNcbi2na = InitFastNcbi4naNcbi2na();
    
    // Initialize tables for conversion of index to codes or names
    InitIndexCodeName();

    // Initialize m_Masks used for random ambiguity resolution
    m_Masks = CSeqportUtil_implementation::InitMasks();

    // Initialize m_DetectAmbigNcbi4naNcbi2na used for ambiguity
    // detection and reporting
    m_DetectAmbigNcbi4naNcbi2na = InitAmbigNcbi4naNcbi2na();

    // Initialize m_DetectAmbigIupacnaNcbi2na used for ambiguity detection
    // and reporting
    m_DetectAmbigIupacnaNcbi2na = InitAmbigIupacnaNcbi2na();

}

// Destructor. All memory allocated on the
// free store is wrapped in smart pointers.
// Therefore, the destructor does not need
// to deallocate memory.
CSeqportUtil_implementation::~CSeqportUtil_implementation()
{
    return;
}


/////////////////////////////////////////////////////////////////////////////
//  PRIVATE::
//


// Helper function to initialize m_SeqCodeSet from sm_StrAsnData
CRef<CSeq_code_set> CSeqportUtil_implementation::Init()
{
    // Compose a long-long string
    string str;
    for (size_t i = 0;  sm_StrAsnData[i];  i++) {
        str += sm_StrAsnData[i];
    }

    // Create an in memory stream on sm_StrAsnData
    CNcbiIstrstream is(str.c_str(), str.length());

    auto_ptr<CObjectIStream>
        asn_codes_in(CObjectIStream::Open(eSerial_AsnText, is));

    // Create a CSeq_code_set
    CRef<CSeq_code_set> ptr_seq_code_set(new CSeq_code_set());

    // Initialize the newly created CSeq_code_set
    *asn_codes_in >> *ptr_seq_code_set;

    // Return a newly created CSeq_code_set
    return ptr_seq_code_set;
}


// Function to initialize code tables
CRef<CSeqportUtil_implementation::CCode_table>
CSeqportUtil_implementation::InitCodes(ESeq_code_type code_type)
{
    // Get list of code tables
    const list<CRef<CSeq_code_table> >& code_list = m_SeqCodeSet->GetCodes();

    // Get table for code_type
    list<CRef<CSeq_code_table> >::const_iterator i_ct;
    for(i_ct = code_list.begin(); i_ct != code_list.end(); ++i_ct)
        if((*i_ct)->GetCode() == code_type)
            break;


    if(i_ct == code_list.end())
        throw runtime_error("Requested code table not found");

    // Get table data
    const list<CRef<CSeq_code_table::C_E> >& table_data = (*i_ct)->GetTable();
    SIZE_TYPE size = table_data.size();
    int start_at = (*i_ct)->GetStart_at();
    CRef<CCode_table> codeTable(new CCode_table(size, start_at));

    // Initialize codeTable to 255
    for(int i=0; i<256; i++)
        codeTable->m_Table[i] = '\xff';

    // Copy table data to codeTable
    int nIdx = start_at;
    list<CRef<CSeq_code_table::C_E> >::const_iterator i_td;
    for(i_td = table_data.begin(); i_td != table_data.end(); ++i_td) {
        codeTable->m_Table[nIdx] =  *((*i_td)->GetSymbol().c_str());
        if(codeTable->m_Table[nIdx] == '\x00')
            codeTable->m_Table[nIdx++] = '\xff';
        else
            nIdx++;
    }

    // Return codeTable
    return codeTable;
}


// Function to initialize iupacna complement table
CRef<CSeqportUtil_implementation::CCode_comp>
CSeqportUtil_implementation::InitIupacnaComplement()
{

    // Get list of code tables
    const list<CRef<CSeq_code_table> >& code_list = m_SeqCodeSet->GetCodes();

    // Get table for code_type iupacna
    list<CRef<CSeq_code_table> >::const_iterator i_ct;
    for(i_ct = code_list.begin(); i_ct != code_list.end(); ++i_ct)
        if((*i_ct)->GetCode() == eSeq_code_type_iupacna)
            break;


    if(i_ct == code_list.end())
        throw runtime_error("Code table for Iupacna not found");

    // Check that complements are set
    if(!(*i_ct)->IsSetComps())
        throw runtime_error("Complement data is not set for iupacna table");

    // Get complement data, start at and size of complement data
    const list<int>& comp_data = (*i_ct)->GetComps();
    int start_at = (*i_ct)->GetStart_at();

    // Allocate memory for complement data
    CRef<CCode_comp> compTable(new CCode_comp(256, start_at));

    // Initialize compTable to 255 for illegal codes
    for(unsigned int i = 0; i<256; i++)
        compTable->m_Table[i] = (char) 255;

    // Loop trhough the complement data and set compTable
    list<int>::const_iterator i_comp;
    unsigned int nIdx = start_at;
    for(i_comp = comp_data.begin(); i_comp != comp_data.end(); ++i_comp)
        compTable->m_Table[nIdx++] = (*i_comp);

    // Return the complement data
    return compTable;

}


// Function to initialize ncbi2na complement table
CRef<CSeqportUtil_implementation::CCode_comp>
CSeqportUtil_implementation::InitNcbi2naComplement()
{

    // Get list of code tables
    const list<CRef<CSeq_code_table> >& code_list = m_SeqCodeSet->GetCodes();

    // Get table for code_type ncbi2na
    list<CRef<CSeq_code_table> >::const_iterator i_ct;
    for(i_ct = code_list.begin(); i_ct != code_list.end(); ++i_ct)
        if((*i_ct)->GetCode() == eSeq_code_type_ncbi2na)
            break;

    if(i_ct == code_list.end())
        throw runtime_error("Code table for Iupacna not found");

    // Check that complements are set
    if(!(*i_ct)->IsSetComps())
        throw runtime_error("Complement data is not set for ncbi2na table");

    // Get complement data, start at and size of complement data
    const list<int>& comp_data = (*i_ct)->GetComps();
    int start_at = (*i_ct)->GetStart_at();

    // Allocate memory for complement data
    CRef<CCode_comp> compTable(new CCode_comp(256, start_at));

    // Put complement data in an array
    char compArray[4];
    int nIdx = start_at;
    list<int>::const_iterator i_comp;
    for(i_comp = comp_data.begin(); i_comp != comp_data.end(); ++i_comp)
        compArray[nIdx++] = (*i_comp);

    // Set compTable
    for(unsigned int i = 0; i < 4; i++)
        for(unsigned int j = 0; j < 4; j++)
            for(unsigned int k = 0; k < 4; k++)
                for(unsigned int l = 0; l < 4; l++)
                    {
                        nIdx = i<<6 | j<<4 | k<<2 | l;
                        char c1 = compArray[i] << 6;
                        char c2 = compArray[j] << 4;
                        char c3 = compArray[k] << 2;
                        char c4 = compArray[l];
                        compTable->m_Table[nIdx] = c1 | c2 | c3 | c4;
                    }

    // Return complement data
    return compTable;

}


// Function to initialize ncbi4na complement table
CRef<CSeqportUtil_implementation::CCode_comp>
CSeqportUtil_implementation::InitNcbi4naComplement()
{

    // Get list of code tables
    const list<CRef<CSeq_code_table> >& code_list = m_SeqCodeSet->GetCodes();

    // Get table for code_type ncbi2na
    list<CRef<CSeq_code_table> >::const_iterator i_ct;
    for(i_ct = code_list.begin(); i_ct != code_list.end(); ++i_ct)
        if((*i_ct)->GetCode() == eSeq_code_type_ncbi4na)
            break;

    if(i_ct == code_list.end())
        throw runtime_error("Code table for Iupacna not found");

    // Check that complements are set
    if(!(*i_ct)->IsSetComps())
        throw runtime_error("Complement data is not set for iupacna table");

    // Get complement data, start at and size of complement data
    const list<int>& comp_data = (*i_ct)->GetComps();
    int start_at = (*i_ct)->GetStart_at();

    // Allocate memory for complement data
    CRef<CCode_comp> compTable(new CCode_comp(256, start_at));


    // Put complement data in an array
    char compArray[16];
    int nIdx = start_at;
    list<int>::const_iterator i_comp;
    for(i_comp = comp_data.begin(); i_comp != comp_data.end(); ++i_comp)
        compArray[nIdx++] = (*i_comp);

    // Set compTable
    for(unsigned int i = 0; i<16; i++)
        for(unsigned int j = 0; j < 16; j++)
            {
                nIdx = i<<4 | j;
                char c1 = compArray[i] << 4;
                char c2 = compArray[j];
                compTable->m_Table[nIdx] = c1 | c2;
            }

    // Return complement data
    return compTable;

}


// Function to initialize m_Ncbi2naRev
CRef<CSeqportUtil_implementation::CCode_rev> CSeqportUtil_implementation::InitNcbi2naRev()
{

    // Allocate memory for reverse table
    CRef<CCode_rev> revTable(new CCode_rev(256, 0));

    // Initialize table used to reverse a byte.
    for(unsigned int i = 0; i < 4; i++)
        for(unsigned int j = 0; j < 4; j++)
            for(unsigned int k = 0; k < 4; k++)
                for(unsigned int l = 0; l < 4; l++)
                    revTable->m_Table[64*i + 16*j + 4*k + l] =
                        64*l + 16*k + 4*j +i;

    // Return the reverse table
    return revTable;
}


// Function to initialize m_Ncbi4naRev
CRef<CSeqportUtil_implementation::CCode_rev> CSeqportUtil_implementation::InitNcbi4naRev()
{

    // Allocate memory for reverse table
    CRef<CCode_rev> revTable(new CCode_rev(256, 0));

    // Initialize table used to reverse a byte.
    for(unsigned int i = 0; i < 16; i++)
        for(unsigned int j = 0; j < 16; j++)
            revTable->m_Table[16*i + j] = 16*j + i;

    // Return the reverse table
    return revTable;
}



// Function to initialize map tables
CRef<CSeqportUtil_implementation::CMap_table> 
CSeqportUtil_implementation::InitMaps
(ESeq_code_type from_type,
 ESeq_code_type to_type)
{

    // Get list of map tables
    const list< CRef< CSeq_map_table > >& map_list = m_SeqCodeSet->GetMaps();

    // Get requested map table
    list<CRef<CSeq_map_table> >::const_iterator i_mt;
    for(i_mt = map_list.begin(); i_mt != map_list.end(); ++i_mt)
        if((*i_mt)->GetFrom() == from_type && (*i_mt)->GetTo() == to_type)
            break;

    if(i_mt == map_list.end())
        throw runtime_error("Requested map table not found");

    // Get the map table
    const list<int>& table_data = (*i_mt)->GetTable();

    // Create a map table reference
    SIZE_TYPE size = table_data.size();
    int start_at = (*i_mt)->GetStart_at();
    CRef<CMap_table> mapTable(new CMap_table(size,start_at));

    // Copy the table data to mapTable
    int nIdx = start_at;
    list<int>::const_iterator i_td;
    for(i_td = table_data.begin(); i_td != table_data.end(); ++i_td)
        {
            mapTable->m_Table[nIdx++] = *i_td;
        }

    return mapTable;
}


// Functions to initialize fast conversion tables
// Function to initialize FastNcib2naIupacna
/*
CRef<CSeqportUtil_implementation::CFast_table4> CSeqportUtil_implementation::InitFastNcbi2naIupacna()
{

    CRef<CFast_table4> fastTable(new CFast_table4(256,0));
    unsigned char i,j,k,l;
    for(i = 0; i < 4; i++)
        for(j = 0; j < 4; j++)
            for(k = 0; k < 4; k++)
                for(l = 0; l < 4; l++)
                    {
                        unsigned char aByte = (i<<6) | (j<<4) | (k<<2) | l;
                        char chi = m_Ncbi2naIupacna->m_Table[i];
                        char chj = m_Ncbi2naIupacna->m_Table[j];
                        char chk = m_Ncbi2naIupacna->m_Table[k];
                        char chl = m_Ncbi2naIupacna->m_Table[l];

                        // Note high order bit pair corresponds to low order
                        // byte etc., on Unix machines.