insertSeqSA.cpp

/*
 * inserts sequences into the SA
 * returns number of SA indexes inserted
 */
#include "insertSeqSA.h"
#include "ErrorWarning.h"
#include "SuffixArrayFuns.h"
#include "SequenceFuns.h"
#include "serviceFuns.cpp"
#include "streamFuns.h"
#include "binarySearch2.h"
#include "funCompareUintAndSuffixes.h"
#include "funCompareUintAndSuffixesMemcmp.h"
#include <cmath>
#include "genomeSAindex.h"
#include "sortSuffixesBucket.h"

uint insertSeqSA(PackedArray & SA, PackedArray & SA1, PackedArray & SAi, char * G, char * G1, uint64 nG, uint64 nG1, uint64 nG2, Parameters & P, Genome &mapGen)
{//insert new sequences into the SA

    uint GstrandBit1 = (uint) floor(log(nG+nG1)/log(2))+1;
    if (GstrandBit1<32) GstrandBit1=32; //TODO: use simple access function for SA
    if ( GstrandBit1+1 != SA.wordLength)
    {//sequence is too long - GstrandBit changed
        ostringstream errOut;
        errOut << "EXITING because of FATAL ERROR: cannot insert sequence on the fly because of strand GstrandBit problem\n";
        errOut << "SOLUTION: please contact STAR author at https://groups.google.com/forum/#!forum/rna-star\n";
        exitWithError(errOut.str(),std::cerr, P.inOut->logMain, EXIT_CODE_GENOME_FILES, P);
    };

    uint N2bit= 1LLU << (SA.wordLength-1);
    uint strandMask=~N2bit;
    for (uint64 isa=0;isa<SA.length; isa++)
    {
        uint64 ind1=SA[isa];
        if ( (ind1 & N2bit)>0 )
        {//- strand
            if ( (ind1 & strandMask)>=nG2 )
            {//the first nG bases
                ind1+=nG1; //reverse complementary indices are all shifted by the length of the sequence
                SA.writePacked(isa,ind1);
            };
        } else
        {//+ strand
            if ( ind1>=nG )
            {//the last nG2 bases
                ind1+=nG1; //reverse complementary indices are all shifted by the length of the sequence
                SA.writePacked(isa,ind1);
            };
        };
    };

    char** seq1=new char*[2];

    #define GENOME_endFillL 16
    char* seqq=new char [4*nG1+3*GENOME_endFillL];//ends shouldbe filled with 5 to mark boundaries

    seq1[0]=seqq+GENOME_endFillL;//TODO: avoid defining an extra array, use reverse search
    seq1[1]=seqq+2*GENOME_endFillL+2*nG1;

    memset(seqq,GENOME_spacingChar,GENOME_endFillL);
    memset(seqq+2*nG1+GENOME_endFillL,GENOME_spacingChar,GENOME_endFillL);
    memset(seqq+4*nG1+2*GENOME_endFillL,GENOME_spacingChar,GENOME_endFillL);

    memcpy(seq1[0], G1, nG1);
    for (uint ii=0; ii<nG1; ii++)
    {//reverse complement sequence
        seq1[0][2*nG1-1-ii]=seq1[0][ii]<4 ? 3-seq1[0][ii] : seq1[0][ii];
    };
    complementSeqNumbers(seq1[0], seq1[1], 2*nG1);//complement

    uint64* indArray=new uint64[nG1*2*2+2];// for each base, 1st number - insertion place in SA, 2nd number - index, *2 for reverse compl


    #pragma omp parallel num_threads(P.runThreadN)
    #pragma omp for schedule (dynamic,1000)
    for (uint ii=0; ii<2*nG1; ii++) {//find insertion points for each of the sequences

        if (seq1[0][ii]>3)
        {//no index for suffices starting with N
            indArray[ii*2]=-1;
        } else
        {
            indArray[ii*2] =  suffixArraySearch1(mapGen, seq1, ii, 10000, nG, (ii<nG1 ? true:false), 0, SA.length-1, 0) ;
            indArray[ii*2+1] = ii;
        };
    };

    uint64 nInd=0;//true number of new indices
    for (uint ii=0; ii<2*nG1; ii++) {//remove entries that cannot be inserted, this cannot be done in the parallel cycle above
        if (indArray[ii*2]!= (uint) -1) {
            indArray[nInd*2]=indArray[ii*2];
            indArray[nInd*2+1]=indArray[ii*2+1];
            ++nInd;
        };
    };

    time_t rawtime;
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished SA search, number of new SA indices = "<<nInd<<endl;

    /*//old-debug
    uint64* indArray1=new uint64[nG1*2*2+2];
    memcpy((void*) indArray1, (void*) indArray, 8*(nG1*2*2+2));
    g_funCompareUintAndSuffixes_G=seq1[0];
    qsort((void*) indArray1, nInd, 2*sizeof(uint64), funCompareUintAndSuffixes);
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished qsort - old " <<endl;
    */
    
    g_funCompareUintAndSuffixesMemcmp_G=seq1[0];
    g_funCompareUintAndSuffixesMemcmp_L=mapGen.pGe.gSuffixLengthMax/sizeof(uint64_t);
    qsort((void*) indArray, nInd, 2*sizeof(uint64_t), funCompareUintAndSuffixesMemcmp);  
    
//     qsort((void*) indArray, nInd, 2*sizeof(uint64), funCompareUint2);
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished qsort" <<endl;


    
    /*//new sorting, 2-step: qsort for indArray, bucket sort for suffixes
    qsort((void*) indArray, nInd, 2*sizeof(uint64), funCompareUint2);
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished qsort"<<nInd<<endl;
    
    sortSuffixesBucket(seq1[0], (void*) indArray, nInd, 2*sizeof(uint64));
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished ordering suffixes"<<nInd<<endl;
    */
    
    /* //debug 
    for (int ii=0;ii<2*nInd;ii++)
    {
        if (indArray[ii]!=indArray1[ii]) 
        {
            cout << ii <<"   "<< indArray[ii]  <<"   "<< indArray1[ii] <<endl;
        };
    };
    */
    
    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished sorting SA indices"<<endl;

    indArray[2*nInd]=-999; //mark the last junction
    indArray[2*nInd+1]=-999; //mark the last junction

    SA1.defineBits(SA.wordLength,SA.length+nInd);

    /*testing
    PackedArray SAo;
    SAo.defineBits(mapGen.GstrandBit+1,mapGen.nSA+nInd);
    SAo.allocateArray();
    ifstream oldSAin("./DirTrue/SA");
    oldSAin.read(SAo.charArray,SAo.lengthByte);
    oldSAin.close();
    */

    uint isa1=0, isa2=0;
    for (uint isa=0;isa<SA.length;isa++) {
        while (isa==indArray[isa1*2]) {//insert new index before the existing index
            uint ind1=indArray[isa1*2+1];
            if (ind1<nG1) {
                ind1+=nG;
            } else {//reverse strand
                ind1=(ind1-nG1+nG2) | N2bit;
            };
            SA1.writePacked(isa2,ind1);
            /*testing
            if (SA1[isa2]!=SAo[isa2]) {
               cout <<isa2 <<" "<< SA1[isa2]<<" "<<SAo[isa2]<<endl;
               //sleep(100);
            };
            */
            ++isa2; ++isa1;

        };

        SA1.writePacked(isa2,SA[isa]); //TODO make sure that the first sj index is not before the first array index
            /*testing
            if (SA1[isa2]!=SAo[isa2]) {
               cout <<isa2 <<" "<< SA1[isa2]<<" "<<SAo[isa2]<<endl;
               //sleep(100);
            };
            */
        ++isa2;
    };
    for (;isa1<nInd;isa1++)
    {//insert the last indices
        uint ind1=indArray[isa1*2+1];
        if (ind1<nG1)
        {
            ind1+=nG;
        } else
        {//reverse strand
            ind1=(ind1-nG1+nG2) | N2bit;
        };
        SA1.writePacked(isa2,ind1);
        ++isa2;
    };

    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished inserting SA indices" <<endl;

//     //SAi insertions
//     for (uint iL=0; iL < P.mapGen.gSAindexNbases; iL++) {
//         uint iSeq=0;
//         uint ind0=mapGen.genomeSAindexStart[iL]-1;//last index that was present in the old genome
//         for (uint ii=mapGen.genomeSAindexStart[iL];ii<mapGen.genomeSAindexStart[iL+1]; ii++) {//scan through the longest index
//             if (ii==798466)
//                 cout <<ii;
//
//             uint iSA1=SAi[ii];
//             uint iSA2=iSA1 & mapGen.SAiMarkNmask & mapGen.SAiMarkAbsentMask;
//
//             if ( iSeq<nInd && (iSA1 &  mapGen.SAiMarkAbsentMaskC)>0 )
//             {//index missing from the old genome
//                 uint iSeq1=iSeq;
//                 int64 ind1=funCalcSAi(seq1[0]+indArray[2*iSeq+1],iL);
//                 while (ind1 < (int64)(ii-mapGen.genomeSAindexStart[iL]) && indArray[2*iSeq]<iSA2) {
//                     ++iSeq;
//                     ind1=funCalcSAi(seq1[0]+indArray[2*iSeq+1],iL);
//                 };
//                 if (ind1 == (int64)(ii-mapGen.genomeSAindexStart[iL]) ) {
//                     SAi.writePacked(ii,indArray[2*iSeq]+iSeq+1);
//                     for (uint ii0=ind0+1; ii0<ii; ii0++) {//fill all the absent indices with this value
//                         SAi.writePacked(ii0,(indArray[2*iSeq]+iSeq+1) | mapGen.SAiMarkAbsentMaskC);
//                     };
//                     ++iSeq;
//                     ind0=ii;
//                 } else {
//                     iSeq=iSeq1;
//                 };
//             } else
//             {//index was present in the old genome
//                 while (iSeq<nInd && indArray[2*iSeq]+1<iSA2) {//for this index insert "smaller" junctions
//                     ++iSeq;
//                 };
//
//                 while (iSeq<nInd && indArray[2*iSeq]+1==iSA2) {//special case, the index falls right behind SAi
//                     if (funCalcSAi(seq1[0]+indArray[2*iSeq+1],iL) >= (int64) (ii-mapGen.genomeSAindexStart[iL]) ) {//this belongs to the next index
//                         break;
//                     };
//                     ++iSeq;
//                 };
//
//                 SAi.writePacked(ii,iSA1+iSeq);
//
//                 for (uint ii0=ind0+1; ii0<ii; ii0++) {//fill all the absent indices with this value
//                     SAi.writePacked(ii0,(iSA2+iSeq) | mapGen.SAiMarkAbsentMaskC);
//                 };
//                 ind0=ii;
//             };
//         };
//
//     };
// //     time ( &rawtime );    cout << timeMonthDayTime(rawtime) << "SAi first" <<endl;
//
//     for (uint isj=0;isj<nInd;isj++) {
//         int64 ind1=0;
//         for (uint iL=0; iL < P.mapGen.gSAindexNbases; iL++) {
//             uint g=(uint) seq1[0][indArray[2*isj+1]+iL];
//             ind1 <<= 2;
//             if (g>3) {//this iSA contains N, need to mark the previous
//                 for (uint iL1=iL; iL1 < P.mapGen.gSAindexNbases; iL1++) {
//                     ind1+=3;
//                     int64 ind2=mapGen.genomeSAindexStart[iL1]+ind1;
//                     for (; ind2>=0; ind2--) {//find previous index that is not absent
//                         if ( (SAi[ind2] & mapGen.SAiMarkAbsentMaskC)==0 ) {
//                             break;
//                         };
//                     };
//                     SAi.writePacked(ind2,SAi[ind2] | mapGen.SAiMarkNmaskC);
//                     ind1 <<= 2;
//                 };
//                 break;
//             } else {
//                 ind1 += g;
//             };
//         };
//     };
//     time ( &rawtime );
//     P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished SAi" <<endl;
//
//     /* testing
//     PackedArray SAio=SAi;
//     SAio.allocateArray();
//     ifstream oldSAiin("./DirTrue/SAindex");
//     oldSAiin.read(SAio.charArray,8*(P.mapGen.gSAindexNbases+2));//skip first bytes
//     oldSAiin.read(SAio.charArray,SAio.lengthByte);
//     oldSAiin.close();
//
//     for (uint iL=0; iL < P.mapGen.gSAindexNbases; iL++) {
//         for (uint ii=mapGen.genomeSAindexStart[iL];ii<mapGen.genomeSAindexStart[iL+1]; ii++) {//scan through the longets index
//                 if ( SAio[ii]!=SAi[ii] ) {
//                     cout <<iL<<" "<<ii<<" "<<SAio[ii]<<" "<<SAi[ii]<<endl;
//                 };
//         };
//     };
//     */

    //change parameters, most parameters are already re-defined in sjdbPrepare.cpp
    SA.defineBits(mapGen.GstrandBit+1,SA.length+nInd);//same as SA2
    SA.pointArray(SA1.charArray);
    mapGen.nSA=SA.length;
    mapGen.nSAbyte=SA.lengthByte;

    //generate SAi
    genomeSAindex(G,SA,P,SAi,mapGen);

    time ( &rawtime );
    P.inOut->logMain  << timeMonthDayTime(rawtime) << "   Finished SAi" <<endl;


//     mapGen.sjGstart=mapGen.chrStart[mapGen.nChrReal];
//     memcpy(G+mapGen.chrStart[mapGen.nChrReal],seq1[0], nseq1[0]);


    return nInd;
};