modified: makefile

[GalaxyCodeBases.git] / BGI / SOAPdenovo2 / standardPregraph / prlReadFillGap.c

/*
 * prlReadFillGap.c
 *
 * Copyright (c) 2008-2012 BGI-Shenzhen <soap at genomics dot org dot cn>.
 *
 * This file is part of SOAPdenovo.
 *
 * SOAPdenovo is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * SOAPdenovo is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with SOAPdenovo.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include "stdinc.h"
#include "newhash.h"
#include "kmerhash.h"
#include "extfunc.h"
#include "extvab.h"
#include "zlib.h"

#define RDBLOCKSIZE 50
#define CTGappend 50

static Kmer MAXKMER;
static int Ncounter;
static int allGaps;

// for multi threads
static int * counters;
static pthread_mutex_t mutex;
static int scafBufSize = 100;
static boolean * flagBuf;
static unsigned char * thrdNoBuf;
static STACK ** ctgStackBuffer;
static int scafCounter;
static int scafInBuf;
static char * scaffBuffer;

static void MarkCtgOccu ( unsigned int ctg );

/*
static void printRead(int len,char *seq)
{
    int j;
    fprintf(stderr,">read\n");
    for(j=0;j<len;j++)
        fprintf(stderr,"%c",int2base((int)getCharInTightString(seq,j)));
    fprintf(stderr,"\n");
}
*/
static void attach1read2contig ( unsigned int ctgID, int len, int pos, long long starter )
{
        unsigned int ctg = index_array[ctgID];  //new index in contig array

        if ( isLargerThanTwin ( ctg ) )
        {
                ctg = getTwinCtg ( ctg ); // put all reads in one contig of a twin
                pos = contig_array[ctg].length + overlaplen - pos - len;
        }

        if ( !contig_array[ctg].closeReads )
        {
                contig_array[ctg].closeReads = ( STACK * ) createStack ( RDBLOCKSIZE, sizeof ( READNEARBY ) );
        }

        READNEARBY * rd = ( READNEARBY * ) stackPush ( contig_array[ctg].closeReads );
        rd->len = len;
        rd->dis = pos;
        rd->seqStarter = starter;
}

static void convertIndex ()
{
        int * length_array = ( int * ) ckalloc ( ( num_ctg + 1 ) * sizeof ( int ) );
        unsigned int i;

        for ( i = 1; i <= num_ctg; i++ )
        {
                length_array[i] = 0;
        }

        for ( i = 1; i <= num_ctg; i++ )
        {
                if ( index_array[i] > 0 )
                {
                        length_array[index_array[i]] = i;
                }
        }

        for ( i = 1; i <= num_ctg; i++ )
        {
                index_array[i] = length_array[i];
        }           //contig i with new index: index_array[i]

        free ( ( void * ) length_array );
}

static long long getRead1by1_gz ( gzFile * fp, DARRAY * readSeqInGap )
{
        long long readCounter = 0;

        if ( !fp )
        {
                return readCounter;
        }

        int len, ctgID, pos;
        long long starter;
        char * pt;
        char * freadBuf = ( char * ) ckalloc ( ( maxReadLen / 4 + 1 ) * sizeof ( char ) );

        while ( gzread ( fp, &len, sizeof ( int ) ) == 4 )
        {
                if ( gzread ( fp, &ctgID, sizeof ( int ) ) != 4 )
                {
                        break;
                }

                if ( gzread ( fp, &pos, sizeof ( int ) ) != 4 )
                {
                        break;
                }

                if ( gzread ( fp, freadBuf, sizeof ( char ) * ( unsigned ) ( len / 4 + 1 ) ) != ( unsigned ) ( len / 4 + 1 ) )
                {
                        break;
                }

                //put seq to dynamic array
                starter = readSeqInGap->item_c;

                if ( !darrayPut ( readSeqInGap, starter + len / 4 ) ) // make sure there's room for this seq
                {
                        break;
                }

                pt = ( char * ) darrayPut ( readSeqInGap, starter );
                bcopy ( freadBuf, pt, len / 4 + 1 );
                attach1read2contig ( ctgID, len, pos, starter );
                readCounter++;
        }

        free ( ( void * ) freadBuf );
        return readCounter;
}

static long long getRead1by1 ( FILE * fp, DARRAY * readSeqInGap )
{
        long long readCounter = 0;

        if ( !fp )
        {
                return readCounter;
        }

        int len, ctgID, pos;
        long long starter;
        char * pt;
        char * freadBuf = ( char * ) ckalloc ( ( maxReadLen / 4 + 1 ) * sizeof ( char ) );

        while ( fread ( &len, sizeof ( int ), 1, fp ) == 1 )
        {
                if ( fread ( &ctgID, sizeof ( int ), 1, fp ) != 1 )
                {
                        break;
                }

                if ( fread ( &pos, sizeof ( int ), 1, fp ) != 1 )
                {
                        break;
                }

                if ( fread ( freadBuf, sizeof ( char ), len / 4 + 1, fp ) != ( unsigned ) ( len / 4 + 1 ) )
                {
                        break;
                }

                //put seq to dynamic array
                starter = readSeqInGap->item_c;

                if ( !darrayPut ( readSeqInGap, starter + len / 4 ) ) // make sure there's room for this seq
                {
                        break;
                }

                pt = ( char * ) darrayPut ( readSeqInGap, starter );
                bcopy ( freadBuf, pt, len / 4 + 1 );
                attach1read2contig ( ctgID, len, pos, starter );
                readCounter++;
        }

        free ( ( void * ) freadBuf );
        return readCounter;
}

// Darray *readSeqInGap
static boolean loadReads4gap ( char * graphfile )
{
        FILE * fp1, *fp2;
        gzFile * fp;
        char name[1024];
        long long readCounter;

        if ( COMPATIBLE_MODE == 1 )
        {
                sprintf ( name, "%s.readInGap", graphfile );
                fp1 = fopen ( name, "rb" );
        }
        else
        {
                sprintf ( name, "%s.readInGap.gz", graphfile );
                fp = gzopen ( name, "rb" );
        }

        sprintf ( name, "%s.longReadInGap", graphfile );
        fp2 = fopen ( name, "rb" );

        if ( COMPATIBLE_MODE == 1 && !fp1 && !fp2 )
        {
                fprintf ( stderr, "Can't open %s.readInGap and %s.longReadInGap!\n", graphfile, graphfile );
                return 0;
        }
        else if ( COMPATIBLE_MODE == 0 && !fp && !fp2 )
        {
                fprintf ( stderr, "Can't open %s.readInGap.gz and %s.longReadInGap!\n", graphfile, graphfile );
                return 0;
        }

        if ( !orig2new )
        {
                convertIndex ();
                orig2new = 1;
        }

        readSeqInGap = ( DARRAY * ) createDarray ( 1000000, sizeof ( char ) );

        if ( COMPATIBLE_MODE == 1 && fp1 )
        {
                readCounter = getRead1by1 ( fp1, readSeqInGap );
                fprintf ( stderr, "Loaded %lld reads from %s.readInGap.\n", readCounter, graphfile );
                fclose ( fp1 );
        }
        else if ( COMPATIBLE_MODE == 0 && fp )
        {
                readCounter = getRead1by1_gz ( fp, readSeqInGap );
                fprintf ( stderr, "Loaded %lld reads from %s.readInGap.\n", readCounter, graphfile );
                gzclose ( fp );
        }

        if ( fp2 )
        {
                readCounter = getRead1by1 ( fp2, readSeqInGap );
                fprintf ( stderr, "Loaded %lld reads from %s.LongReadInGap.\n", readCounter, graphfile );
                fclose ( fp2 );
        }

        return 1;
}

static void debugging1 ()
{
        unsigned int i;

        if ( orig2new )
        {
                unsigned int * length_array = ( unsigned int * ) ckalloc ( ( num_ctg + 1 ) * sizeof ( unsigned int ) );

                //use length_array to change info in index_array
                for ( i = 1; i <= num_ctg; i++ )
                {
                        length_array[i] = 0;
                }

                for ( i = 1; i <= num_ctg; i++ )
                {
                        if ( index_array[i] > 0 )
                        {
                                length_array[index_array[i]] = i;
                        }
                }

                for ( i = 1; i <= num_ctg; i++ )
                {
                        index_array[i] = length_array[i];
                }       //contig i with original index: index_array[i]

                orig2new = 0;
        }

        READNEARBY * rd;
        int j;
        char * pt;

        for ( i = 1; i <= num_ctg; i++ )
        {
                if ( !contig_array[i].closeReads )
                {
                        continue;
                }

                if ( index_array[i] != 735 )
                {
                        continue;
                }

                fprintf ( stderr, "Contig %d, len %d: \n", index_array[i], contig_array[i].length );
                stackBackup ( contig_array[i].closeReads );

                while ( ( rd = ( READNEARBY * ) stackPop ( contig_array[i].closeReads ) ) != NULL )
                {
                        fprintf ( stderr, "%d\t%d\t%lld\t", rd->dis, rd->len, rd->seqStarter );
                        pt = ( char * ) darrayGet ( readSeqInGap, rd->seqStarter );

                        for ( j = 0; j < rd->len; j++ )
                        {
                                fprintf ( stderr, "%c", int2base ( ( int ) getCharInTightString ( pt, j ) ) );
                        }

                        fprintf ( stderr, "\n" );
                }

                stackRecover ( contig_array[i].closeReads );
        }
}

static void initiateCtgInScaf ( CTGinSCAF * actg )
{
        actg->cutTail = 0;
        actg->cutHead = overlaplen;
        actg->gapSeqLen = 0;
}

static int procGap ( char * line, STACK * ctgsStack )
{
        char * tp;
        int length, i, seg;
        unsigned int ctg;
        CTGinSCAF * ctgPt;
        tp = strtok ( line, " " );
        tp = strtok ( NULL, " " );  //length
        length = atoi ( tp );
        tp = strtok ( NULL, " " );  //seg
        seg = atoi ( tp );

        if ( !seg )
        {
                return length;
        }

        for ( i = 0; i < seg; i++ )
        {
                tp = strtok ( NULL, " " );
                ctg = atoi ( tp );
                MarkCtgOccu ( ctg );
                ctgPt = ( CTGinSCAF * ) stackPush ( ctgsStack );
                initiateCtgInScaf ( ctgPt );
                ctgPt->ctgID = ctg;
                ctgPt->start = 0;
                ctgPt->end = 0;
                ctgPt->scaftig_start = 0;
                ctgPt->mask = 1;
        }

        return length;
}

static void debugging2 ( int index, STACK * ctgsStack )
{
        CTGinSCAF * actg;
        stackBackup ( ctgsStack );
        fprintf ( stderr, ">scaffold%d\t%d 0.0\n", index, ctgsStack->item_c );

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                fprintf ( stderr, "%d\t%d\t%d\t%d\n", actg->ctgID, actg->start, actg->end, actg->scaftig_start );
        }

        stackRecover ( ctgsStack );
}

static int cmp_reads ( const void * a, const void * b )
{
        READNEARBY * A, *B;
        A = ( READNEARBY * ) a;
        B = ( READNEARBY * ) b;

        if ( A->dis > B->dis )
        {
                return 1;
        }
        else if ( A->dis == B->dis )
        {
                return 0;
        }
        else
        {
                return -1;
        }
}

static void cutRdArray ( READNEARBY * rdArray, int gapStart, int gapEnd, int * count, int arrayLen, READNEARBY * cutArray )
{
        int i;
        int num = 0;

        for ( i = 0; i < arrayLen; i++ )
        {
                if ( rdArray[i].dis > gapEnd )
                {
                        break;
                }

                if ( ( rdArray[i].dis + rdArray[i].len ) >= gapStart )
                {
                        cutArray[num].dis = rdArray[i].dis;
                        cutArray[num].len = rdArray[i].len;
                        cutArray[num++].seqStarter = rdArray[i].seqStarter;
                }
        }

        *count = num;
}

static void outputTightStr2Visual ( FILE * foc2, int ctg,  int * num, int start, int length, int outputlen, int revS )
{
        int i, end, column = 0, n = *num;
        char * tightStr = contig_array[ctg].seq;

        if ( n == 1 ) { fprintf ( foc2, ">%d\n", index_array[ctg] ); }
        else { fprintf ( foc2, ">%d-%d\n", index_array[ctg], n - 1 ); }

        if ( !revS )
        {
                end = start + outputlen <= length ? start + outputlen : length;

                for ( i = start; i < end; i++ )
                {
                        fprintf ( foc2, "%c", int2base ( ( int ) getCharInTightString ( tightStr, i ) ) );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                //column = 0;
                                fprintf ( foc2, "\n" );
                        }
                }

                if ( column % 100 != 0 )
                {
                        fprintf ( foc2, "\n" );
                }
        }
        else
        {
                end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

                for ( i = end; i <= length - 1 - start; i++ )
                {
                        fprintf ( foc2, "%c", int2base ( ( int ) getCharInTightString ( tightStr, i ) ) );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                fprintf ( foc2, "\n" );
                                //column = 0;
                        }
                }

                if ( column % 100 != 0 )
                {
                        fprintf ( foc2, "\n" );
                }
        }
}

void outputTightStr ( FILE * fp, char * tightStr, int start, int length, int outputlen, int revS, int * col )
{
        int i;
        int end;
        int column = *col;

        if ( !revS )
        {
                end = start + outputlen <= length ? start + outputlen : length;

                for ( i = start; i < end; i++ )
                {
                        fprintf ( fp, "%c", int2base ( ( int ) getCharInTightString ( tightStr, i ) ) );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                //column = 0;
                                fprintf ( fp, "\n" );
                        }
                }
        }
        else
        {
                end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

                for ( i = length - 1 - start; i >= end; i-- )
                {
                        fprintf ( fp, "%c", int2compbase ( ( int ) getCharInTightString ( tightStr, i ) ) );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                fprintf ( fp, "\n" );
                                //column = 0;
                        }
                }
        }

        *col = column;
}

static void outputTightStr2 ( char * tightStr, char * scaff, int start, int length, int outputlen, int revS, int * col )
{
        int i;
        int end;
        int column = *col;
        char a;

        if ( !revS )
        {
                end = start + outputlen <= length ? start + outputlen : length;

                for ( i = start; i < end; i++ )
                {
                        a = int2base ( ( int ) getCharInTightString ( tightStr, i ) );
                        //  fprintf (fp, "%c", int2base ((int) getCharInTightString (tightStr, i)));
                        scaff[column++] = a;
                }
        }
        else
        {
                end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

                for ( i = length - 1 - start; i >= end; i-- )
                {
                        a = int2compbase ( ( int ) getCharInTightString ( tightStr, i ) );
                        //  fprintf (fp, "%c", int2compbase ((int) getCharInTightString (tightStr, i)));
                        scaff[column++] = a;
                }
        }

        *col = column;
}

static void outputTightStrLowerCase2Visual ( FILE * foc2, int gapNum, char * tightStr, int start, int length, int outputlen )
{
        int i, end, column = 0;
        end = start + outputlen <= length ? start + outputlen : length;
        fprintf ( foc2, ">%d-0\n", gapNum );

        for ( i = start; i < end; i++ )
        {
                fprintf ( foc2, "%c", "actg"[ ( int ) getCharInTightString ( tightStr, i )] );

                if ( ( ++column ) % 100 == 0 )
                {
                        //column = 0;
                        fprintf ( foc2, "\n" );
                }
        }

        if ( column % 100 != 0 )
        {
                fprintf ( foc2, "\n" );
        }
}

static void outputTightStrLowerCase ( FILE * fp, char * tightStr, int start, int length, int outputlen, int revS, int * col )
{
        int i;
        int end;
        int column = *col;

        if ( !revS )
        {
                end = start + outputlen <= length ? start + outputlen : length;

                for ( i = start; i < end; i++ )
                {
                        fprintf ( fp, "%c", "actg"[ ( int ) getCharInTightString ( tightStr, i )] );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                //column = 0;
                                fprintf ( fp, "\n" );
                        }
                }
        }
        else
        {
                end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

                for ( i = length - 1 - start; i >= end; i-- )
                {
                        fprintf ( fp, "%c", "tgac"[ ( int ) getCharInTightString ( tightStr, i )] );

                        if ( ( ++column ) % 100 == 0 )
                        {
                                fprintf ( fp, "\n" );
                                //column = 0;
                        }
                }
        }

        *col = column;
}

static void outputTightStrLowerCase2 ( char * tightStr, char * scaff, int start, int length, int outputlen, int revS, int * col )
{
        int i;
        int end;
        int column = *col;
        char a;

        if ( !revS )
        {
                end = start + outputlen <= length ? start + outputlen : length;

                for ( i = start; i < end; i++ )
                {
                        a = "actg"[ ( int ) getCharInTightString ( tightStr, i )];
                        //  fprintf (fp, "%c", "actg"[(int) getCharInTightString (tightStr, i)]);
                        scaff[column++] = a;
                }
        }
        else
        {
                end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

                for ( i = length - 1 - start; i >= end; i-- )
                {
                        a = "tgac"[ ( int ) getCharInTightString ( tightStr, i )];
                        //  fprintf (fp, "%c", "tgac"[(int) getCharInTightString (tightStr, i)]);
                        scaff[column++] = a;
                }
        }

        *col = column;
}

static void outputNs ( FILE * fp, int gapN, int * col )
{
        int i, column = *col;

        for ( i = 0; i < gapN; i++ )
        {
                fprintf ( fp, "N" );

                if ( ( ++column ) % 100 == 0 )
                {
                        //column = 0;
                        fprintf ( fp, "\n" );
                }
        }

        *col = column;
}

static void outputNs2 ( char * scaff, int gapN, int * col )
{
        int i, column = *col;

        for ( i = 0; i < gapN; i++ )
        {
                scaff[column] = 'N';
                column++;
        }

        *col = column;
}

static void outputGapInfo ( unsigned int ctg1, unsigned int ctg2 )
{
        unsigned int bal_ctg1 = getTwinCtg ( ctg1 );
        unsigned int bal_ctg2 = getTwinCtg ( ctg2 );

        if ( isLargerThanTwin ( ctg1 ) )
        {
                fprintf ( stderr, "%d\t", index_array[bal_ctg1] );
        }
        else
        {
                fprintf ( stderr, "%d\t", index_array[ctg1] );
        }

        if ( isLargerThanTwin ( ctg2 ) )
        {
                fprintf ( stderr, "%d\n", index_array[bal_ctg2] );
        }
        else
        {
                fprintf ( stderr, "%d\n", index_array[ctg2] );
        }
}

static void output1gap ( FILE * fo, int scafIndex, CTGinSCAF * prevCtg, CTGinSCAF * actg, DARRAY * gapSeqArray )
{
        unsigned int ctg1, bal_ctg1, length1;
        int start1, outputlen1;
        unsigned int ctg2, bal_ctg2, length2;
        int start2, outputlen2;
        char * pt;
        int column = 0;
        ctg1 = prevCtg->ctgID;
        bal_ctg1 = getTwinCtg ( ctg1 );
        start1 = prevCtg->cutHead;
        length1 = contig_array[ctg1].length + overlaplen;

        if ( length1 - prevCtg->cutTail - start1 > CTGappend )
        {
                outputlen1 = CTGappend;
                start1 = length1 - prevCtg->cutTail - outputlen1;
        }
        else
        {
                outputlen1 = length1 - prevCtg->cutTail - start1;
        }

        ctg2 = actg->ctgID;
        bal_ctg2 = getTwinCtg ( ctg2 );
        start2 = actg->cutHead;
        length2 = contig_array[ctg2].length + overlaplen;

        if ( length2 - actg->cutTail - start2 > CTGappend )
        {
                outputlen2 = CTGappend;
        }
        else
        {
                outputlen2 = length2 - actg->cutTail - start2;
        }

        if ( isLargerThanTwin ( ctg1 ) )
        {
                fprintf ( fo, ">S%d_C%d_L%d_G%d", scafIndex, index_array[bal_ctg1], outputlen1, prevCtg->gapSeqLen );
        }
        else
        {
                fprintf ( fo, ">S%d_C%d_L%d_G%d", scafIndex, index_array[ctg1], outputlen1, prevCtg->gapSeqLen );
        }

        if ( isLargerThanTwin ( ctg2 ) )
        {
                fprintf ( fo, "_C%d_L%d\t", index_array[bal_ctg2], outputlen2 );
        }
        else
        {
                fprintf ( fo, "_C%d_L%d\t", index_array[ctg2], outputlen2 );
        }

        fprintf ( fo, "%d\n", start2 );

        if ( contig_array[ctg1].seq )
        {
                outputTightStr ( fo, contig_array[ctg1].seq, start1, length1, outputlen1, 0, &column );
        }
        else if ( contig_array[bal_ctg1].seq )
        {
                outputTightStr ( fo, contig_array[bal_ctg1].seq, start1, length1, outputlen1, 1, &column );
        }

        pt = ( char * ) darrayPut ( gapSeqArray, prevCtg->gapSeqOffset );
        outputTightStrLowerCase ( fo, pt, 0, prevCtg->gapSeqLen, prevCtg->gapSeqLen, 0, &column );

        if ( contig_array[ctg2].seq )
        {
                outputTightStr ( fo, contig_array[ctg2].seq, start2, length2, outputlen2, 0, &column );
        }
        else if ( contig_array[bal_ctg2].seq )
        {
                outputTightStr ( fo, contig_array[bal_ctg2].seq, start2, length2, outputlen2, 1, &column );
        }

        if ( column % 100 != 0 )
        {
                fprintf ( fo, "\n" );
        }
}

static void outputGapSeq ( FILE * fo, int index, STACK * ctgsStack, DARRAY * gapSeqArray )
{
        CTGinSCAF * actg, *prevCtg = NULL;
        stackRecover ( ctgsStack );
        //  fprintf (fo, ">scaffold%d\n", index);

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                /*  if (prevCtg)
                    {
                        if (actg->scaftig_start)
                        {
                            fprintf (fo, "0\t%d\t%d\n", prevCtg->mask, actg->mask);
                        }
                        else
                        {
                            fprintf (fo, "1\t%d\t%d\n", prevCtg->mask, actg->mask);
                        }
                    }
                */
                if ( prevCtg && prevCtg->gapSeqLen > 0 )
                        { output1gap ( fo, index, prevCtg, actg, gapSeqArray ); }

                prevCtg = actg;
        }
}

static void outputScafSeq ( FILE * fo, FILE * foc, FILE * foc2, FILE * fo3, int index, STACK * ctgsStack, DARRAY * gapSeqArray )
{
        CTGinSCAF * actg, *prevCtg = NULL;
        unsigned int ctg, bal_ctg, ctg_out, length;
        int start, outputlen, gapN;
        char * pt;
        int column = 0;
        long long cvgSum = 0;
        int lenSum = 0;
        int i, t, lu_len = 0, lu_end = 0;
        unsigned int ctg_start_pos = 0;
        char strand;
        unsigned int * pos_start = ( unsigned int * ) ckalloc ( 1000000 * sizeof ( unsigned int ) );
        unsigned int * pos_end = ( unsigned int * ) ckalloc ( 1000000 * sizeof ( unsigned int ) );
        //  char index_contig[num_ctg][20];
        char ** index_contig = ( char ** ) ckalloc ( 1000000 * sizeof ( char * ) );

        for ( i = 0; i < 1000000; i++ )
                { index_contig[i] = ( char * ) ckalloc ( 20 * sizeof ( char ) ); }

        char * orien_array;
        orien_array = ( char * ) ckalloc ( 1000000 * sizeof ( char ) );
        //  scaffBuffer = (char *) ckalloc (300000000 * sizeof (char));
        stackRecover ( ctgsStack );

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                if ( ! ( contig_array[actg->ctgID].cvg > 0 ) )
                {
                        continue;
                }

                lenSum += contig_array[actg->ctgID].length;
                cvgSum += contig_array[actg->ctgID].length * contig_array[actg->ctgID].cvg;
        }

        if ( lenSum > 0 )
        {
                fprintf ( fo, ">scaffold%d %4.1f\n", index, ( double ) cvgSum / lenSum );
        }
        else
        {
                fprintf ( fo, ">scaffold%d 0.0\n", index );
        }

        fprintf ( foc, ">scaffold%d\n", index );
        stackRecover ( ctgsStack );

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                ctg = actg->ctgID;
                bal_ctg = getTwinCtg ( ctg );
                length = contig_array[ctg].length + overlaplen;

                if ( prevCtg && actg->scaftig_start )
                {
                        gapN = actg->start - prevCtg->start - contig_array[prevCtg->ctgID].length;
                        gapN = gapN > 0 ? gapN : 1;
                        outputNs ( fo,  gapN, &column );
                        ctg_start_pos += gapN;
                        //outputGapInfo(prevCtg->ctgID,ctg);
                        Ncounter++;
                }

                if ( !prevCtg )
                {
                        start = 0;
                }
                else
                {
                        start = actg->cutHead;
                }

                outputlen = length - start - actg->cutTail;

                if ( contig_array[ctg].seq )
                {
                        outputTightStr ( fo, contig_array[ctg].seq, start, length, outputlen, 0, &column );
                        lu_end = start + outputlen > length ? length : start + outputlen;
                        lu_len = lu_end - start;
                        strand = '+';
                        fprintf ( foc, "%d\t", index_array[ctg] );
                }
                else if ( contig_array[bal_ctg].seq )
                {
                        outputTightStr ( fo, contig_array[bal_ctg].seq, start, length, outputlen, 1, &column );
                        lu_end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;
                        lu_len = length - start - lu_end;
                        strand = '-';
                        fprintf ( foc, "%d\t", index_array[bal_ctg] );
                }

                fprintf ( foc, "%u\t%c\t%d\n", ctg_start_pos, strand, lu_len );
                ctg_start_pos += lu_len;

                if ( actg->gapSeqLen < 1 )
                {
                        prevCtg = actg;
                        continue;
                }

                pt = ( char * ) darrayPut ( gapSeqArray, actg->gapSeqOffset );
                outputTightStrLowerCase ( fo, pt, 0, actg->gapSeqLen, actg->gapSeqLen, 0, &column );
                ctg_start_pos = ctg_start_pos + actg->gapSeqLen;
                prevCtg = actg;
        }

        if ( column % 100 != 0 )
        {
                fprintf ( fo, "\n" );
        }

        if ( visual )
        {
                scaffBuffer = ( char * ) ckalloc ( ( ctg_start_pos + 5 ) * sizeof ( char ) );
                prevCtg = NULL;
                column = 0;
                ctg_start_pos = 0;
                lenSum = 0;
                stackRecover ( ctgsStack );

                while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
                {
                        ctg = actg->ctgID;
                        bal_ctg = getTwinCtg ( ctg );
                        length = contig_array[ctg].length + overlaplen;

                        if ( prevCtg && actg->scaftig_start )
                        {
                                gapN = actg->start - prevCtg->start - contig_array[prevCtg->ctgID].length;
                                gapN = gapN > 0 ? gapN : 1;
                                outputNs2 ( scaffBuffer, gapN, &column );
                                ctg_start_pos += gapN;
                                //  Ncounter++;
                        }

                        if ( !prevCtg )
                        {
                                start = 0;
                        }
                        else
                        {
                                start = actg->cutHead;
                        }

                        outputlen = length - start - actg->cutTail;

                        if ( contig_array[ctg].seq )
                        {
                                t = ++contig_index_array[index_array[ctg]];
                                outputTightStr2 ( contig_array[ctg].seq, scaffBuffer, start, length, outputlen, 0, &column );
                                lu_end = start + outputlen > length ? length : start + outputlen;
                                lu_len = lu_end - start;
                                strand = '+';
                                //  fprintf (foc, "%d\t", index_array[ctg]);
                                lenSum++;

                                if ( ctg_start_pos - start >= 0 )
                                {
                                        pos_start[lenSum] = ctg_start_pos - start;
                                        pos_end[lenSum] = ctg_start_pos + length - start;
                                        orien_array[lenSum] = '+';

                                        if ( t == 1 )  { sprintf ( index_contig[lenSum], "%u", index_array[ctg] ); }
                                        else  { sprintf ( index_contig[lenSum], "%u-%d", index_array[ctg], t - 1 ); }

                                        fprintf ( fo3, "{AFG\nacc:%s\nclr:0,%d\n}\n", ( index_contig[lenSum] ), length );
                                        outputTightStr2Visual ( foc2, ctg, & ( contig_index_array[index_array[ctg]] ), 0, length, length, 0 );
                                }
                                else
                                {
                                        pos_start[lenSum] = 0;
                                        pos_end[lenSum] = ctg_start_pos + length - start;
                                        orien_array[lenSum] = '+';

                                        if ( t == 1 )  { sprintf ( index_contig[lenSum], "%u", index_array[ctg] ); }
                                        else  { sprintf ( index_contig[lenSum], "%u-%d", index_array[ctg], t - 1 ); }

                                        fprintf ( fo3, "{AFG\nacc:%s\nclr:0,%d\n}\n", ( index_contig[lenSum] ), length );
                                        outputTightStr2Visual ( foc2, ctg, & ( contig_index_array[index_array[ctg]] ), start - ctg_start_pos, length, length - start + ctg_start_pos, 0 );
                                }
                        }
                        else if ( contig_array[bal_ctg].seq )
                        {
                                t = ++contig_index_array[index_array[bal_ctg]];
                                outputTightStr2 ( contig_array[bal_ctg].seq, scaffBuffer, start, length, outputlen, 1, &column );
                                lu_end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;
                                lu_len = length - start - lu_end;
                                strand = '-';
                                //  fprintf (foc, "%d\t", index_array[bal_ctg]);
                                lenSum++;

                                if ( ctg_start_pos - start >= 0 )
                                {
                                        pos_start[lenSum] = ctg_start_pos - start;
                                        pos_end[lenSum] = ctg_start_pos + length - start;
                                        orien_array[lenSum] = '-';

                                        if ( t == 1 )  { sprintf ( index_contig[lenSum], "%u", index_array[bal_ctg] ); }
                                        else  { sprintf ( index_contig[lenSum], "%u-%d", index_array[bal_ctg], t - 1 ); }

                                        fprintf ( fo3, "{AFG\nacc:%s\nclr:0,%d\n}\n", ( index_contig[lenSum] ), length );
                                        outputTightStr2Visual ( foc2, bal_ctg, & ( contig_index_array[index_array[bal_ctg]] ), 0, length, length, 1 );
                                }
                                else
                                {
                                        pos_start[lenSum] = ctg_start_pos;
                                        pos_end[lenSum] = ctg_start_pos + lu_len;
                                        orien_array[lenSum] = '-';

                                        if ( t == 1 )  { sprintf ( index_contig[lenSum], "%u", index_array[bal_ctg] ); }
                                        else  { sprintf ( index_contig[lenSum], "%u-%d", index_array[bal_ctg], t - 1 ); }

                                        fprintf ( fo3, "{AFG\nacc:%s\nclr:0,%d\n}\n", ( index_contig[lenSum] ), length );
                                        outputTightStr2Visual ( foc2, bal_ctg, & ( contig_index_array[index_array[bal_ctg]] ), start, length, outputlen, 1 );
                                }
                        }

                        //  fprintf (foc, "%u\t%c\t%d\n", ctg_start_pos, strand, lu_len);
                        ctg_start_pos += lu_len;

                        if ( actg->gapSeqLen < 1 )
                        {
                                prevCtg = actg;
                                continue;
                        }

                        pt = ( char * ) darrayPut ( gapSeqArray, actg->gapSeqOffset );
                        outputTightStrLowerCase2 ( pt, scaffBuffer, 0, actg->gapSeqLen, actg->gapSeqLen, 0, &column );
                        outputTightStrLowerCase2Visual ( foc2, gapNum, pt, 0, actg->gapSeqLen, actg->gapSeqLen );
                        fprintf ( fo3, "{AFG\nacc:%d-0\nclr:0,%d\n}\n", gapNum, actg->gapSeqLen );
                        lenSum++;
                        pos_start[lenSum] = ctg_start_pos;
                        pos_end[lenSum] = ctg_start_pos + actg->gapSeqLen;
                        orien_array[lenSum] = '+';
                        sprintf ( index_contig[lenSum], "%d-0", gapNum );
                        gapNum++;
                        ctg_start_pos = ctg_start_pos + actg->gapSeqLen;
                        prevCtg = actg;
                }

                scaffNum++;
                fprintf ( fo3, "{CCO\nacc:%d\npla:P\nlen:%u\ncns:\n", scaffNum, ctg_start_pos );

                for ( i = 0; i < ctg_start_pos; i++ )
                {
                        if ( i != 0 && i % 100 == 0 && i < ctg_start_pos - 1 )  { fprintf ( fo3, "\n" ); }

                        fprintf ( fo3, "%c", scaffBuffer[i] );
                }

                fprintf ( fo3, "\n.\nqlt:\n" );

                for ( i = 0; i < ctg_start_pos; i++ )
                {
                        if ( i != 0 && i % 100 == 0 && i < ctg_start_pos - 1 )  { fprintf ( fo3, "\n" ); }

                        fprintf ( fo3, "D" );
                }

                fprintf ( fo3, "\n.\nnpc:%d\n", lenSum );

                for ( i = 1; i <= lenSum; i++ )
                {
                        if ( orien_array[i] == '+' ) { fprintf ( fo3, "{MPS\ntyp:R\nmid:%s\nsrc:\n.\npos:%u,%u\ndln:0\ndel:\n}\n", ( index_contig[i] ), pos_start[i], pos_end[i] ); }

                        if ( orien_array[i] == '-' ) { fprintf ( fo3, "{MPS\ntyp:R\nmid:%s\nsrc:\n.\npos:%u,%u\ndln:0\ndel:\n}\n", ( index_contig[i] ), pos_end[i], pos_start[i] ); }
                }

                fprintf ( fo3, "}\n" );
                free ( ( void * ) scaffBuffer );
        }

        free ( ( void * ) pos_start );
        free ( ( void * ) pos_end );
        free ( ( void * ) orien_array );

        for ( i = 0; i < 1000000; i++ )
        {
                free ( ( void * ) index_contig[i] );
        }

        free ( ( void * ) index_contig );
}

static void fill1scaf ( int index, STACK * ctgsStack, int thrdID );
static void check1scaf ( int t, int thrdID )
{
        if ( flagBuf[t] )
        {
                return;
        }

        boolean late = 0;
        pthread_mutex_lock ( &mutex );

        if ( !flagBuf[t] )
        {
                flagBuf[t] = 1;
                thrdNoBuf[t] = thrdID;
        }
        else
        {
                late = 1;
        }

        pthread_mutex_unlock ( &mutex );

        if ( late )
        {
                return;
        }

        counters[thrdID]++;
        fill1scaf ( scafCounter + t + 1, ctgStackBuffer[t], thrdID );
}

static void fill1scaf ( int index, STACK * ctgsStack, int thrdID )
{
        CTGinSCAF * actg, *prevCtg = NULL;
        READNEARBY * rdArray, *rdArray4gap, *rd;
        int numRd = 0, count, maxGLen = 0;
        unsigned int ctg, bal_ctg;
        STACK * rdStack;

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                if ( prevCtg )
                {
                        maxGLen = maxGLen < ( actg->start - prevCtg->end ) ? ( actg->start - prevCtg->end ) : maxGLen;
                }

                ctg = actg->ctgID;
                bal_ctg = getTwinCtg ( ctg );

                if ( actg->mask )
                {
                        prevCtg = actg;
                        continue;
                }

                if ( contig_array[ctg].closeReads )
                {
                        numRd += contig_array[ctg].closeReads->item_c;
                }
                else if ( contig_array[bal_ctg].closeReads )
                {
                        numRd += contig_array[bal_ctg].closeReads->item_c;
                }

                prevCtg = actg;
        }

        if ( numRd < 1 )
        {
                return;
        }

        rdArray = ( READNEARBY * ) ckalloc ( numRd * sizeof ( READNEARBY ) );
        rdArray4gap = ( READNEARBY * ) ckalloc ( numRd * sizeof ( READNEARBY ) );
        //fprintf(stderr,"scaffold%d reads4gap %d\n",index,numRd);
        // collect reads appended to contigs in this scaffold
        int numRd2 = 0;
        stackRecover ( ctgsStack );

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                ctg = actg->ctgID;
                bal_ctg = getTwinCtg ( ctg );

                if ( actg->mask )
                {
                        continue;
                }

                if ( contig_array[ctg].closeReads )
                {
                        rdStack = contig_array[ctg].closeReads;
                }
                else if ( contig_array[bal_ctg].closeReads )
                {
                        rdStack = contig_array[bal_ctg].closeReads;
                }
                else
                {
                        continue;
                }

                stackBackup ( rdStack );

                while ( ( rd = ( READNEARBY * ) stackPop ( rdStack ) ) != NULL )
                {
                        rdArray[numRd2].len = rd->len;
                        rdArray[numRd2].seqStarter = rd->seqStarter;

                        if ( isSmallerThanTwin ( ctg ) )
                        {
                                rdArray[numRd2++].dis = actg->start - overlaplen + rd->dis;
                        }
                        else
                                { rdArray[numRd2++].dis = actg->start - overlaplen + contig_array[ctg].length - rd->len - rd->dis; }
                }

                stackRecover ( rdStack );
        }

        if ( numRd2 != numRd )
        {
                fprintf ( stderr, "##reads numbers doesn't match, %d vs %d when scaffold %d.\n", numRd, numRd2, index );
        }

        qsort ( rdArray, numRd, sizeof ( READNEARBY ), cmp_reads );
        //fill gap one by one
        int gapStart, gapEnd;
        int numIn = 0;
        boolean flag;
        int buffer_size = maxReadLen > 100 ? maxReadLen : 100;
        int maxGSLen = maxGLen + GLDiff < 10 ? 10 : maxGLen + GLDiff;
        //fprintf(stderr,"maxGlen %d, maxGSlen %d\n",maxGLen,maxGSLen);
        char * seqGap = ( char * ) ckalloc ( maxGSLen * sizeof ( char ) ); // temp array for gap sequence
        Kmer * kmerCtg1 = ( Kmer * ) ckalloc ( buffer_size * sizeof ( Kmer ) );
        Kmer * kmerCtg2 = ( Kmer * ) ckalloc ( buffer_size * sizeof ( Kmer ) );
        char * seqCtg1 = ( char * ) ckalloc ( buffer_size * sizeof ( char ) );
        char * seqCtg2 = ( char * ) ckalloc ( buffer_size * sizeof ( char ) );
        prevCtg = NULL;
        stackRecover ( ctgsStack );

        while ( ( actg = stackPop ( ctgsStack ) ) != NULL )
        {
                if ( !prevCtg || !actg->scaftig_start )
                {
                        prevCtg = actg;
                        continue;
                }

                gapStart = prevCtg->end - 100;
                gapEnd = actg->start - overlaplen + 100;
                cutRdArray ( rdArray, gapStart, gapEnd, &count, numRd, rdArray4gap );
                numIn += count;
                /*
                   if(!count){
                   prevCtg = actg;
                   continue;
                   }
                 */
                int overlap;

                for ( overlap = overlaplen; overlap > 14; overlap -= 2 )
                {
                        flag = localGraph ( rdArray4gap, count, prevCtg, actg, overlaplen, kmerCtg1, kmerCtg2, overlap, darrayBuf[thrdID], seqCtg1, seqCtg2, seqGap );

                        //free_kmerset(kmerSet);

                        if ( flag == 1 )
                        {
                                /*
                                   fprintf(stderr,"Between ctg %d and %d, Found with %d\n",prevCtg->ctgID
                                   ,actg->ctgID,overlap);
                                 */
                                break;
                        }
                }

                /*
                   if(count==0)
                   printf("Gap closed without reads\n");
                   if(!flag)
                   fprintf(stderr,"Between ctg %d and %d, NO routes found\n",prevCtg->ctgID,actg->ctgID);
                 */
                prevCtg = actg;
        }

        //fprintf(stderr,"____scaffold%d reads in gap %d\n",index,numIn);
        free ( ( void * ) seqGap );
        free ( ( void * ) kmerCtg1 );
        free ( ( void * ) kmerCtg2 );
        free ( ( void * ) seqCtg1 );
        free ( ( void * ) seqCtg2 );
        free ( ( void * ) rdArray );
        free ( ( void * ) rdArray4gap );
}

static void reverseStack ( STACK * dStack, STACK * sStack )
{
        CTGinSCAF * actg, *ctgPt;
        emptyStack ( dStack );

        while ( ( actg = ( CTGinSCAF * ) stackPop ( sStack ) ) != NULL )
        {
                ctgPt = ( CTGinSCAF * ) stackPush ( dStack );
                ctgPt->ctgID = actg->ctgID;
                ctgPt->start = actg->start;
                ctgPt->end = actg->end;
                ctgPt->scaftig_start = actg->scaftig_start;
                ctgPt->mask = actg->mask;
                ctgPt->cutHead = actg->cutHead;
                ctgPt->cutTail = actg->cutTail;
                ctgPt->gapSeqLen = actg->gapSeqLen;
                ctgPt->gapSeqOffset = actg->gapSeqOffset;
        }

        stackBackup ( dStack );
}

#ifdef MER127
static Kmer tightStr2Kmer ( char * tightStr, int start, int length, int revS )
{
        int i;
        Kmer word;
        word.high1 = word.low1 = word.high2 = word.low2 = 0;

        if ( !revS )
        {
                if ( start + overlaplen > length )
                {
                        fprintf ( stderr, "The tightStr2Kmer A: no enough bases for kmer.\n" );
                        return word;
                }

                for ( i = start; i < start + overlaplen; i++ )
                {
                        word = KmerLeftBitMoveBy2 ( word );
                        word.low2 |= getCharInTightString ( tightStr, i );
                }
        }
        else
        {
                if ( length - start - overlaplen < 0 )
                {
                        fprintf ( stderr, "The tightStr2Kmer B: no enough bases for kmer.\n" );
                        return word;
                }

                for ( i = length - 1 - start; i > length - 1 - start - overlaplen; i-- )
                {
                        word = KmerLeftBitMoveBy2 ( word );
                        word.low2 |= int_comp ( getCharInTightString ( tightStr, i ) );
                }
        }

        return word;
}

static Kmer maxKmer ()
{
        Kmer word;
        word.high1 = word.low1 = word.high2 = word.low2 = 0;
        int i;

        for ( i = 0; i < overlaplen; i++ )
        {
                word = KmerLeftBitMoveBy2 ( word );
                word.low2 |= 0x3;
        }

        return word;
}
#else
static Kmer tightStr2Kmer ( char * tightStr, int start, int length, int revS )
{
        int i;
        Kmer word;
        word.high = word.low = 0;

        if ( !revS )
        {
                if ( start + overlaplen > length )
                {
                        fprintf ( stderr, "The tightStr2Kmer A: no enough bases for kmer.\n" );
                        return word;
                }

                for ( i = start; i < start + overlaplen; i++ )
                {
                        word = KmerLeftBitMoveBy2 ( word );
                        word.low |= getCharInTightString ( tightStr, i );
                }
        }
        else
        {
                if ( length - start - overlaplen < 0 )
                {
                        fprintf ( stderr, "The tightStr2Kmer B: no enough bases for kmer.\n" );
                        return word;
                }

                for ( i = length - 1 - start; i > length - 1 - start - overlaplen; i-- )
                {
                        word = KmerLeftBitMoveBy2 ( word );
                        word.low |= int_comp ( getCharInTightString ( tightStr, i ) );
                }
        }

        return word;
}

static Kmer maxKmer ()
{
        Kmer word;
        word.high = word.low = 0;
        int i;

        for ( i = 0; i < overlaplen; i++ )
        {
                word = KmerLeftBitMoveBy2 ( word );
                word.low |= 0x3;
        }

        return word;
}
#endif
static int contigCatch ( unsigned int prev_ctg, unsigned int ctg )
{
        if ( contig_array[prev_ctg].length == 0 || contig_array[ctg].length == 0 )
        {
                return 0;
        }

        Kmer kmerAtEnd, kmerAtStart;
        Kmer MaxKmer;
        unsigned int bal_ctg1 = getTwinCtg ( prev_ctg );
        unsigned int bal_ctg2 = getTwinCtg ( ctg );
        int i, start;
        int len1 = contig_array[prev_ctg].length + overlaplen;
        int len2 = contig_array[ctg].length + overlaplen;
        start = contig_array[prev_ctg].length;

        if ( contig_array[prev_ctg].seq )
        {
                kmerAtEnd = tightStr2Kmer ( contig_array[prev_ctg].seq, start, len1, 0 );
        }
        else
        {
                kmerAtEnd = tightStr2Kmer ( contig_array[bal_ctg1].seq, start, len1, 1 );
        }

        start = 0;

        if ( contig_array[ctg].seq )
        {
                kmerAtStart = tightStr2Kmer ( contig_array[ctg].seq, start, len2, 0 );
        }
        else
        {
                kmerAtStart = tightStr2Kmer ( contig_array[bal_ctg2].seq, start, len2, 1 );
        }

        MaxKmer = MAXKMER;

        for ( i = 0; i < 10; i++ )
        {
                if ( KmerEqual ( kmerAtStart, kmerAtEnd ) )
                {
                        break;
                }

                MaxKmer = KmerRightBitMoveBy2 ( MaxKmer );
                kmerAtEnd = KmerAnd ( kmerAtEnd, MaxKmer );
                kmerAtStart = KmerRightBitMoveBy2 ( kmerAtStart );
        }

        if ( i < 10 )
        {
                return overlaplen - i;
        }
        else
        {
                return 0;
        }
}

static void initStackBuf ( STACK ** ctgStackBuffer, int scafBufSize )
{
        int i;

        for ( i = 0; i < scafBufSize; i++ )
        {
                flagBuf[i] = 1;
                ctgStackBuffer[i] = ( STACK * ) createStack ( 100, sizeof ( CTGinSCAF ) );
        }
}
static void freeStackBuf ( STACK ** ctgStackBuffer, int scafBufSize )
{
        int i;

        for ( i = 0; i < scafBufSize; i++ )
        {
                freeStack ( ctgStackBuffer[i] );
        }
}

static void threadRoutine ( void * para )
{
        PARAMETER * prm;
        int i;
        prm = ( PARAMETER * ) para;

        //printf("%dth thread with threadID %d, hash_table %p\n",id,prm.threadID,prm.hash_table);
        while ( 1 )
        {
                if ( * ( prm->selfSignal ) == 1 )
                {
                        emptyDarray ( darrayBuf[prm->threadID] );

                        for ( i = 0; i < scafInBuf; i++ )
                        {
                                check1scaf ( i, prm->threadID );
                        }

                        * ( prm->selfSignal ) = 0;
                }
                else if ( * ( prm->selfSignal ) == 2 )
                {
                        * ( prm->selfSignal ) = 0;
                        break;
                }

                usleep ( 1 );
        }
}

static void creatThrds ( pthread_t * threads, PARAMETER * paras )
{
        unsigned char i;
        int temp;

        for ( i = 0; i < thrd_num; i++ )
        {
                if ( ( temp = pthread_create ( &threads[i], NULL, ( void * ) threadRoutine, & ( paras[i] ) ) ) != 0 )
                {
                        fprintf ( stderr, "Create threads failed.\n" );
                        exit ( 1 );
                }
        }

        fprintf ( stderr, "%d thread(s) initialized.\n", thrd_num );
}

static void sendWorkSignal ( unsigned char SIG, unsigned char * thrdSignals )
{
        int t;

        for ( t = 0; t < thrd_num; t++ )
        {
                thrdSignals[t + 1] = SIG;
        }

        while ( 1 )
        {
                usleep ( 10 );

                for ( t = 0; t < thrd_num; t++ )
                        if ( thrdSignals[t + 1] )
                        {
                                break;
                        }

                if ( t == thrd_num )
                {
                        break;
                }
        }
}

static void thread_wait ( pthread_t * threads )
{
        int i;

        for ( i = 0; i < thrd_num; i++ )
                if ( threads[i] != 0 )
                {
                        pthread_join ( threads[i], NULL );
                }
}

static void outputSeqs ( FILE * fo, FILE * foc, FILE * foc2, FILE * fo2, FILE * fo3, int scafInBuf )
{
        int i, thrdID;

        for ( i = 0; i < scafInBuf; i++ )
        {
                thrdID = thrdNoBuf[i];
                outputScafSeq ( fo, foc, foc2, fo3, scafCounter + i + 1, ctgStackBuffer[i], darrayBuf[thrdID] );
                outputGapSeq ( fo2, scafCounter + i + 1, ctgStackBuffer[i], darrayBuf[thrdID] );
        }
}

static void MaskContig ( unsigned int ctg )
{
        contig_array[ctg].mask = 1;
        contig_array[getTwinCtg ( ctg )].mask = 1;
}

static void MarkCtgOccu ( unsigned int ctg )
{
        contig_array[ctg].flag = 1;
        contig_array[getTwinCtg ( ctg )].flag = 1;
}

static void output_ctg ( unsigned int ctg, FILE * fo )
{
        if ( contig_array[ctg].length < 1 )
        {
                return;
        }

        int len;
        unsigned int bal_ctg = getTwinCtg ( ctg );
        len = contig_array[ctg].length + overlaplen;
        int col = 0;

        if ( contig_array[ctg].seq )
        {
                fprintf ( fo, ">C%d %4.1f\n", ctg, ( double ) contig_array[ctg].cvg );
                outputTightStr ( fo, contig_array[ctg].seq, 0, len, len, 0, &col );
        }
        else if ( contig_array[bal_ctg].seq )
        {
                fprintf ( fo, ">C%d %4.1f\n", bal_ctg, ( double ) contig_array[ctg].cvg );
                outputTightStr ( fo, contig_array[bal_ctg].seq, 0, len, len, 0, &col );
        }

        contig_array[ctg].flag = 1;
        contig_array[bal_ctg].flag = 1;

        if ( len % 100 != 0 )
        {
                fprintf ( fo, "\n" );
        }
}

void prlReadsCloseGap ( char * graphfile )
{
        //thrd_num=1;
        if ( fillGap )
        {
                boolean flag;
                fprintf ( stderr, "\nStart to load reads for gap filling. %d length discrepancy is allowed.\n", GLDiff );
                fprintf ( stderr, "...\n" );
                flag = loadReads4gap ( graphfile );

                if ( !flag )
                {
                        return;
                }
        }

        if ( orig2new )
        {
                convertIndex ();
                orig2new = 0;
        }

        FILE * fp, *fo, *fo2, *fo3 = NULL, *foc, *foc2 = NULL;
        char line[1024];
        CTGinSCAF * actg;
        STACK * ctgStack, *aStack;
        int index = 0, offset = 0, counter, overallLen;
        int i, starter, prev_start, gapLen, catchable, scafnum;
        unsigned int ctg, prev_ctg = 0;
        boolean IsPrevGap;
        pthread_t threads[thrd_num];
        unsigned char thrdSignal[thrd_num + 1];
        PARAMETER paras[thrd_num];

        for ( ctg = 1; ctg <= num_ctg; ctg++ )
        {
                contig_array[ctg].flag = 0;
        }

        MAXKMER = maxKmer ();
        ctgStack = ( STACK * ) createStack ( 1000, sizeof ( CTGinSCAF ) );
        sprintf ( line, "%s.scaf_gap", graphfile );
        fp = ckopen ( line, "r" );
        sprintf ( line, "%s.scafSeq", graphfile );
        fo = ckopen ( line, "w" );
        sprintf ( line, "%s.contigPosInscaff", graphfile );
        foc = ckopen ( line, "w" );

        if ( visual )
        {
                sprintf ( line, "%s.contig4asm", graphfile );
                foc2 = ckopen ( line, "w" );
                sprintf ( line, "%s.asm", graphfile );
                fo3 = ckopen ( line, "w" );
        }

        sprintf ( line, "%s.gapSeq", graphfile );
        fo2 = ckopen ( line, "w" );
        pthread_mutex_init ( &mutex, NULL );
        flagBuf = ( boolean * ) ckalloc ( scafBufSize * sizeof ( boolean ) );;
        thrdNoBuf = ( unsigned char * ) ckalloc ( scafBufSize * sizeof ( unsigned char ) );;
        memset ( thrdNoBuf, 0, scafBufSize * sizeof ( char ) );
        ctgStackBuffer = ( STACK ** ) ckalloc ( scafBufSize * sizeof ( STACK * ) );
        initStackBuf ( ctgStackBuffer, scafBufSize );
        darrayBuf = ( DARRAY ** ) ckalloc ( thrd_num * sizeof ( DARRAY * ) );
        counters = ( int * ) ckalloc ( thrd_num * sizeof ( int ) );
        contig_index_array = ( int * ) ckalloc ( ( num_ctg + 1 ) * sizeof ( int ) );

        for ( i = 0; i <= num_ctg; i++ )
        {
                contig_index_array[i] = 0;
        }

        for ( i = 0; i < thrd_num; i++ )
        {
                counters[i] = 0;
                darrayBuf[i] = ( DARRAY * ) createDarray ( 100000, sizeof ( char ) );
                thrdSignal[i + 1] = 0;
                paras[i].threadID = i;
                paras[i].mainSignal = &thrdSignal[0];
                paras[i].selfSignal = &thrdSignal[i + 1];
        }

        if ( fillGap )
        {
                creatThrds ( threads, paras );
        }

        Ncounter = scafCounter = scafInBuf = allGaps = 0;

        while ( fgets ( line, sizeof ( line ), fp ) != NULL )
        {
                if ( line[0] == '>' )
                {
                        if ( index )
                        {
                                aStack = ctgStackBuffer[scafInBuf];
                                flagBuf[scafInBuf++] = 0;
                                reverseStack ( aStack, ctgStack );

                                if ( scafInBuf == scafBufSize )
                                {
                                        if ( fillGap )
                                        {
                                                sendWorkSignal ( 1, thrdSignal );
                                        }

                                        outputSeqs ( fo, foc, foc2, fo2, fo3, scafInBuf );
                                        scafCounter += scafInBuf;
                                        scafInBuf = 0;
                                }

                                if ( index % 1000 == 0 )
                                {
                                        fprintf ( stderr, "%d scaffolds processed.\n", index );
                                }
                        }

                        //read next scaff
                        emptyStack ( ctgStack );
                        IsPrevGap = offset = prev_ctg = 0;
                        sscanf ( line + 9, "%d %d %d", &index, &counter, &overallLen );
                        continue;
                }

                if ( line[0] == 'G' ) // gap appears
                {
                        if ( fillGap )
                        {
                                gapLen = procGap ( line, ctgStack );
                                IsPrevGap = 1;
                        }

                        continue;
                }

                if ( line[0] >= '0' && line[0] <= '9' ) // a contig line
                {
                        sscanf ( line, "%d %d", &ctg, &starter );
                        actg = ( CTGinSCAF * ) stackPush ( ctgStack );
                        actg->ctgID = ctg;

                        if ( contig_array[ctg].flag )
                        {
                                MaskContig ( ctg );
                        }
                        else
                        {
                                MarkCtgOccu ( ctg );
                        }

                        initiateCtgInScaf ( actg );

                        if ( !prev_ctg )
                        {
                                actg->cutHead = 0;
                        }
                        else if ( !IsPrevGap )
                        {
                                allGaps++;
                        }

                        if ( !IsPrevGap )
                        {
                                if ( prev_ctg && ( starter - prev_start - ( int ) contig_array[prev_ctg].length ) < ( ( int ) overlaplen * 4 ) )
                                {
                                        /*
                                           if(fillGap)
                                           catchable = contigCatch(prev_ctg,ctg);
                                           else
                                         */
                                        catchable = 0;

                                        if ( catchable ) // prev_ctg and ctg overlap **bp
                                        {
                                                allGaps--;
                                                /*
                                                   if(isLargerThanTwin(prev_ctg))
                                                   fprintf(stderr,"%d #######  by_overlap\n",getTwinCtg(prev_ctg));
                                                   else
                                                   fprintf(stderr,"%d ####### by_overlap\n",prev_ctg);
                                                 */
                                                actg->scaftig_start = 0;
                                                actg->cutHead = catchable;
                                                offset += - ( starter - prev_start - contig_array[prev_ctg].length ) + ( overlaplen - catchable );
                                        }
                                        else
                                        {
                                                actg->scaftig_start = 1;
                                        }
                                }
                                else
                                {
                                        actg->scaftig_start = 1;
                                }
                        }
                        else
                        {
                                offset += - ( starter - prev_start - contig_array[prev_ctg].length ) + gapLen;
                                actg->scaftig_start = 0;
                        }

                        actg->start = starter + offset;
                        actg->end = actg->start + contig_array[ctg].length - 1;
                        actg->mask = contig_array[ctg].mask;
                        IsPrevGap = 0;
                        prev_ctg = ctg;
                        prev_start = starter;
                }
        }

        if ( index )
        {
                aStack = ctgStackBuffer[scafInBuf];
                flagBuf[scafInBuf++] = 0;
                reverseStack ( aStack, ctgStack );

                if ( fillGap )
                {
                        sendWorkSignal ( 1, thrdSignal );
                }

                outputSeqs ( fo, foc, foc2, fo2, fo3, scafInBuf );
        }

        if ( visual )
        {
                scafnum = scaffNum;

                for ( i = 1; i <= scafnum; i++ )
                {
                        fprintf ( fo3, "{SCF\nacc:%d\nnoc:0\n{CTP\nct1:%d\nct2:%d\nmea:0\nstd:0\nori:N\n}\n}\n", ++scaffNum, i, i );
                }
        }

        if ( fillGap )
        {
                sendWorkSignal ( 2, thrdSignal );
                thread_wait ( threads );
        }

        for ( ctg = 1; ctg <= num_ctg; ctg++ )
        {
                if ( ( contig_array[ctg].length + overlaplen ) < 100 || contig_array[ctg].flag )
                {
                        continue;
                }

                output_ctg ( ctg, fo );
        }

        fprintf ( stderr, "\nDone with %d scaffolds, %d gaps finished, %d gaps overall.\n", index, allGaps - Ncounter, allGaps );
        index = 0;

        for ( i = 0; i < thrd_num; i++ )
        {
                freeDarray ( darrayBuf[i] );
                index += counters[i];
        }

        if ( fillGap )
        {
                fprintf ( stderr, "Threads processed %d scaffolds.\n", index );
        }

        free ( ( void * ) darrayBuf );

        if ( readSeqInGap )
        {
                freeDarray ( readSeqInGap );
        }

        fclose ( fp );
        fclose ( fo );
        fclose ( foc );
        fclose ( fo2 );

        if ( visual )
        {
                fclose ( foc2 );
                fclose ( fo3 );
        }

        freeStack ( ctgStack );
        freeStackBuf ( ctgStackBuffer, scafBufSize );
        free ( ( void * ) flagBuf );
        free ( ( void * ) thrdNoBuf );
        free ( ( void * ) ctgStackBuffer );
}