seqport_util.cpp

ncbi源码

                        char *pt = 
                            reinterpret_cast<char*>(&fastTable->m_Table[aByte]);
                        *(pt++) = chi;
                        *(pt++) = chj;
                        *(pt++) = chk;
                        *(pt) = chl;                       
                     }
    return fastTable;
}
*/

// Function to initialize FastNcib2naNcbi4na
CRef<CSeqportUtil_implementation::CFast_table2> CSeqportUtil_implementation::InitFastNcbi2naNcbi4na()
{

    CRef<CFast_table2> fastTable(new CFast_table2(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;
                    unsigned char chi = m_Ncbi2naNcbi4na->m_Table[i];
                    unsigned char chj = m_Ncbi2naNcbi4na->m_Table[j];
                    unsigned char chk = m_Ncbi2naNcbi4na->m_Table[k];
                    unsigned char chl = m_Ncbi2naNcbi4na->m_Table[l];
                    char *pt = 

                        reinterpret_cast<char*>(&fastTable->m_Table[aByte]);
                    *(pt++) = (chi << 4) | chj;
                    *pt = (chk << 4) | chl;
                }
    return fastTable;
}


// Function to initialize FastNcib4naIupacna
CRef<CSeqportUtil_implementation::CFast_table2> CSeqportUtil_implementation::InitFastNcbi4naIupacna()
{

    CRef<CFast_table2> fastTable(new CFast_table2(256,0));
    unsigned char i,j;
    for(i = 0; i < 16; i++)
        for(j = 0; j < 16; j++) {
            unsigned char aByte = (i<<4) | j;
            unsigned char chi = m_Ncbi4naIupacna->m_Table[i];
            unsigned char chj = m_Ncbi4naIupacna->m_Table[j];

            // Note high order nible corresponds to low order byte
            // etc., on Unix machines.
            char *pt = reinterpret_cast<char*>(&fastTable->m_Table[aByte]);
            *(pt++) = chi;
            *pt = chj;
        }
    return fastTable;
}


// Function to initialize m_FastIupacnancbi2na
CRef<CSeqportUtil_implementation::CFast_4_1> CSeqportUtil_implementation::InitFastIupacnaNcbi2na()
{

    int start_at = m_IupacnaNcbi2na->m_StartAt;
    int size = m_IupacnaNcbi2na->m_Size;
    CRef<CFast_4_1> fastTable(new CFast_4_1(4,0,256,0));
    for(int ch = 0; ch < 256; ch++) {
        if((ch >= start_at) && (ch < (start_at + size)))
            {
                unsigned char uch = m_IupacnaNcbi2na->m_Table[ch];
                uch &= '\x03';
                for(unsigned int pos = 0; pos < 4; pos++)
                    fastTable->m_Table[pos][ch] = uch << (6-2*pos);
            }
        else
            for(unsigned int pos = 0; pos < 4; pos++)
                fastTable->m_Table[pos][ch] = '\x00';
    }
    return fastTable;
}


// Function to initialize m_FastIupacnancbi4na
CRef<CSeqportUtil_implementation::CFast_2_1> CSeqportUtil_implementation::InitFastIupacnaNcbi4na()
{

    int start_at = m_IupacnaNcbi4na->m_StartAt;
    int size = m_IupacnaNcbi4na->m_Size;
    CRef<CFast_2_1> fastTable(new CFast_2_1(2,0,256,0));
    for(int ch = 0; ch < 256; ch++) {
        if((ch >= start_at) && (ch < (start_at + size)))
            {
                unsigned char uch = m_IupacnaNcbi4na->m_Table[ch];
                for(unsigned int pos = 0; pos < 2; pos++)
                    fastTable->m_Table[pos][ch] = uch << (4-4*pos);
            }
        else
            {
                fastTable->m_Table[0][ch] = 0xF0;
                fastTable->m_Table[1][ch] = 0x0F;
            }
    }
    return fastTable;
}


// Function to initialize m_FastNcbi4naNcbi2na
CRef<CSeqportUtil_implementation::CFast_2_1> CSeqportUtil_implementation::InitFastNcbi4naNcbi2na()
{

    int start_at = m_Ncbi4naNcbi2na->m_StartAt;
    int size = m_Ncbi4naNcbi2na->m_Size;
    CRef<CFast_2_1> fastTable(new CFast_2_1(2,0,256,0));
    for(int n1 = 0; n1 < 16; n1++)
        for(int n2 = 0; n2 < 16; n2++) {
            int nIdx = 16*n1 + n2;
            unsigned char u1, u2;
            if((n1 >= start_at) && (n1 < start_at + size))
                u1 = m_Ncbi4naNcbi2na->m_Table[n1] & 3;
            else
                u1 = '\x00';
            if((n2 >= start_at) && (n2 < start_at + size))
                u2 = m_Ncbi4naNcbi2na->m_Table[n2] & 3;
            else
                u2 = '\x00';
            fastTable->m_Table[0][nIdx] = (u1<<6) | (u2<<4);
            fastTable->m_Table[1][nIdx] = (u1<<2) | u2;
        }

    return fastTable;
}


// Function to initialize m_IndexString and m_StringIndex
void CSeqportUtil_implementation::InitIndexCodeName()
{
    typedef list<CRef<CSeq_code_table> >      Ttables;
    typedef list<CRef<CSeq_code_table::C_E> > Tcodes;
    
    m_IndexString[kName].resize(kNumCodes);
    m_IndexString[kSymbol].resize(kNumCodes);
    m_IndexComplement.resize(kNumCodes);
    m_StringIndex.resize(kNumCodes);
    m_StartAt.resize(kNumCodes);

    bool found[kNumCodes];
    for (unsigned int ii = 0; ii < kNumCodes; ii++) {
        found[ii] = false;
    }
    ITERATE (Ttables, it, m_SeqCodeSet->GetCodes()) {
        const ESeq_code_type& code = (*it)->GetCode();
        if (!found[code-1]) {
            found[code-1] = true;
            m_StartAt[code-1] = (*it)->IsSetStart_at() ?
                (*it)->GetStart_at() : 0;
            TIndex i = m_StartAt[code-1];
            ITERATE(Tcodes, is, (*it)->GetTable()) {                
                m_IndexString[kSymbol][code-1].push_back((*is)->GetSymbol());
                m_IndexString[kName][code-1].push_back((*is)->GetName());
                m_StringIndex[code-1].insert
                    (make_pair((*is)->GetSymbol(), i++));
            }
            if ( (*it)->IsSetComps() ) {
                ITERATE (list<int>, ic, (*it)->GetComps()) {
                    m_IndexComplement[code-1].push_back(*ic);
                }
            }
        }
    }
    
     
}


// Function to initialize m_Masks
CRef<CSeqportUtil_implementation::SMasksArray> CSeqportUtil_implementation::InitMasks()
{

    unsigned int i, j, uCnt;
    unsigned char cVal, cRslt;
    CRef<SMasksArray> aMask(new SMasksArray);

    // Initialize possible masks for converting ambiguous
    // ncbi4na bytes to unambiguous bytes
    static const unsigned char mask[16] = {
        0x11, 0x12, 0x14, 0x18,
        0x21, 0x22, 0x24, 0x28,
        0x41, 0x42, 0x44, 0x48,
        0x81, 0x82, 0x84, 0x88
    };

    static const unsigned char maskUpper[4] = { 0x10, 0x20, 0x40, 0x80 };
    static const unsigned char maskLower[4] = { 0x01, 0x02, 0x04, 0x08 };

    // Loop through possible ncbi4na bytes and
    // build masks that convert it to unambiguous na
    for(i = 0; i < 256; i++) {
        cVal = i;
        uCnt = 0;

        // Case where both upper and lower nible > 0
        if(((cVal & '\x0f') != 0) && ((cVal & '\xf0') != 0))
            for(j = 0; j < 16; j++) {
                cRslt = cVal & mask[j];
                if(cRslt == mask[j])
                    aMask->m_Table[i].cMask[uCnt++] = mask[j];
            }

        // Case where upper nible = 0 and lower nible > 0
        else if((cVal & '\x0f') != 0)
            for(j = 0; j < 4; j++)
                {
                    cRslt = cVal & maskLower[j];
                    if(cRslt == maskLower[j])
                        aMask->m_Table[i].cMask[uCnt++] = maskLower[j];
                }


        // Case where lower nible = 0 and upper nible > 0
        else if((cVal & '\xf0') != 0)
            for(j = 0; j < 4; j++)
                {
                    cRslt = cVal & maskUpper[j];
                    if(cRslt == maskUpper[j])
                        aMask->m_Table[i].cMask[uCnt++] = maskUpper[j];
                }

        // Both upper and lower nibles = 0
        else
            aMask->m_Table[i].cMask[uCnt++] = '\x00';

        // Number of distict masks for ncbi4na byte i
        aMask->m_Table[i].nMasks = uCnt;

        // Fill out the remainder of cMask array with copies
        // of first uCnt masks
        for(j = uCnt; j < 16 && uCnt > 0; j++)
            aMask->m_Table[i].cMask[j] = aMask->m_Table[i].cMask[j % uCnt];

    }

    return aMask;
}


// Function to initialize m_DetectAmbigNcbi4naNcbi2na used for
// ambiguity detection
CRef<CSeqportUtil_implementation::CAmbig_detect> CSeqportUtil_implementation::InitAmbigNcbi4naNcbi2na()
{

    unsigned char low, high, ambig;

    // Create am new CAmbig_detect object
    CRef<CAmbig_detect> ambig_detect(new CAmbig_detect(256,0));

    // Loop through low and high order nibles and assign
    // values as follows: 0 - no ambiguity, 1 - low order nible ambigiguous
    // 2 - high order ambiguous, 3 -- both high and low ambiguous.

    // Loop for low order nible
    for(low = 0; low < 16; low++) {
        // Determine if low order nible is ambiguous
        if((low == 1) || (low ==2) || (low == 4) || (low == 8))
            ambig = 0;  // Not ambiguous
        else
            ambig = 1;  // Ambiguous

        // Loop for high order nible
        for(high = 0; high < 16; high++) {

            // Determine if high order nible is ambiguous
            if((high != 1) && (high != 2) && (high != 4) && (high != 8))
                ambig += 2;  // Ambiguous

            // Set ambiguity value
            ambig_detect->m_Table[16*high + low] = ambig;

            // Reset ambig
            ambig &= '\xfd';  // Set second bit to 0
        }
    }

    return ambig_detect;
}


// Function to initialize m_DetectAmbigIupacnaNcbi2na used for ambiguity
// detection
CRef<CSeqportUtil_implementation::CAmbig_detect> CSeqportUtil_implementation::InitAmbigIupacnaNcbi2na()
{

    // Create am new CAmbig_detect object
    CRef<CAmbig_detect> ambig_detect(new CAmbig_detect(256,0));

    // 0 implies no ambiguity. 1 implies ambiguity
    // Initialize to 0
    for(unsigned int i = 0; i<256; i++)
        ambig_detect->m_Table[i] = 0;

    // Set iupacna characters that are ambiguous when converted
    // to ncib2na
    ambig_detect->m_Table[66] = 1; // B
    ambig_detect->m_Table[68] = 1; // D
    ambig_detect->m_Table[72] = 1; // H
    ambig_detect->m_Table[75] = 1; // K
    ambig_detect->m_Table[77] = 1; // M
    ambig_detect->m_Table[78] = 1; // N
    ambig_detect->m_Table[82] = 1; // R
    ambig_detect->m_Table[83] = 1; // S
    ambig_detect->m_Table[86] = 1; // V
    ambig_detect->m_Table[87] = 1; // W
    ambig_detect->m_Table[89] = 1; // Y

    return ambig_detect;
}

/*
struct SSeqDataToSeqUtil
{
    CSeq_data::E_Choice  seq_data_coding;
    CSeqConvert::TCoding seq_convert_coding;
};


static SSeqDataToSeqUtil s_SeqDataToSeqUtilMap[] = {
    { CSeq_data::e_Iupacna,   CSeqUtil::e_Iupacna },
    { CSeq_data::e_Iupacaa,   CSeqUtil::e_Iupacna },
    { CSeq_data::e_Ncbi2na,   CSeqUtil::e_Ncbi2na },
    { CSeq_data::e_Ncbi4na,   CSeqUtil::e_Ncbi4na },
    { CSeq_data::e_Ncbi8na,   CSeqUtil::e_Ncbi8na },
    { CSeq_data::e_Ncbi8aa,   CSeqUtil::e_Ncbi8aa },
    { CSeq_data::e_Ncbieaa,   CSeqUtil::e_Ncbieaa },
    { CSeq_data::e_Ncbistdaa, CSeqUtil::e_Ncbistdaa }
};
*/

static CSeqUtil::TCoding s_SeqDataToSeqUtil[] = {
    CSeqUtil::e_not_set,
    CSeqUtil::e_Iupacna,
    CSeqUtil::e_Iupacaa,
    CSeqUtil::e_Ncbi2na,
    CSeqUtil::e_Ncbi4na,
    CSeqUtil::e_Ncbi8na,
    CSeqUtil::e_not_set,
    CSeqUtil::e_Ncbi8aa,
    CSeqUtil::e_Ncbieaa,
    CSeqUtil::e_not_set,
    CSeqUtil::e_Ncbistdaa
};


// Convert from one coding scheme to another. The following
// 12 conversions are supported: ncbi2na<=>ncbi4na;
// ncbi2na<=>iupacna; ncbi4na<=>iupacna; ncbieaa<=>ncbistdaa;
// ncbieaa<=>iupacaa; ncbistdaa<=>iupacaa. Convert is
// really just a dispatch function--it calls the appropriate
// priviate conversion function.
TSeqPos CSeqportUtil_implementation::x_ConvertAmbig
(const CSeq_data&      in_seq,
 CSeq_data*            out_seq,
 CSeq_data::E_Choice   to_code,
 TSeqPos               uBeginIdx,
 TSeqPos               uLength,
 CRandom::TValue       seed)
    const
{
    CSeq_data::E_Choice from_code = in_seq.Which();

    if(to_code == CSeq_data::e_not_set || from_code == CSeq_data::e_not_set)
        throw std::runtime_error("to_code or from_code not set");

    if ( to_code != CSeq_data::e_Ncbi2na ) {
        throw std::runtime_error("to_code is not Ncbi2na");
    }

    switch (from_code) {
    case CSeq_data::e_Iupacna:
        return MapIupacnaToNcbi2na(in_seq, out_seq,
            uBeginIdx, uLength, true, seed);
    case CSeq_data::e_Ncbi4na:
        return MapNcbi4naToNcbi2na(in_seq, out_seq,
            uBeginIdx, uLength, true, seed);
    default:
        throw runtime_error("Requested conversion not implemented");
    }
}

// Convert from one coding scheme to another. The following
// 12 conversions are supported: ncbi2na<=>ncbi4na;
// ncbi2na<=>iupacna; ncbi4na<=>iupacna; ncbieaa<=>ncbistdaa;
// ncbieaa<=>iupacaa; ncbistdaa<=>iupacaa. Convert is
// really just a dispatch function--it calls the appropriate
// priviate conversion function.
TSeqPos CSeqportUtil_implementation::Convert
(const CSeq_data&      in_seq,
 CSeq_data*