main.cpp

#include <stdio.h>
#include <getopt.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <stdarg.h>
#include <pthread.h>

#include <vector>

#include "KmerCount.hpp"
#include "SeqSet.hpp"
#include "AlignAlgo.hpp"

char usage[] = "./trust4 [OPTIONS]:\n"
		"Required:\n"
		"\t-f STRING: fasta file containing the receptor genome sequence\n"
		"\t[Read file]\n"
		"\t-u STRING: path to single-end read file\n"
		"\t\tor\n"
		"\t-1 STRING -2 STRING: path to paried-end read files\n"
		"\t\tor\n"
		"\t-b STRING: path to BAM alignment file\n"
		"Optional:\n"
		"\t-o STRING: prefix of the output file (default: trust)\n"
		"\t-t INT: number of threads (default: 1)\n"
		"\t-c STRING: the path to the kmer count file\n"
		"\t-k INT: the starting k-mer size for indexing contigs (default: 9)\n"
		"\t--minHitLen INT: the minimal hit length for a valid overlap (default: auto)\n"
		"\t--skipMateExtension: skip the step of extension assemblies with mate-pair information\n"
		///"\t--noV: do not assemble the full length V gene (default: not used)\n"
		"\t--trimLevel INT: 0: no trim; 1: trim low quality; 2: trim unmatched (default: 1)\n"
		"\t--barcode STRING: the path to the barcode file (default: not used)\n"
		"\t--UMI STRING: the path to the UMI file (default: not used)\n"
		"\t--cgeneEnd INT: skipping reads mapped to C gene coordinate greater than INT (default: 200)\n"
		"\t--keepNoBarcode: assemble the reads with missing barcodes. (default: ignore the reads)\n" 
		"\t--contigMinCov INT: ignore contigs that have bases covered by fewer than INT reads (default: 0)\n" ;

int nucToNum[26] = { 0, -1, 1, -1, -1, -1, 2, 
	-1, -1, -1, -1, -1, -1, 0,
	-1, -1, -1, -1, -1, 3,
	-1, -1, -1, -1, -1, -1 } ;

char numToNuc[26] = {'A', 'C', 'G', 'T'} ;

char buffer[10240] = "" ;

static const char *short_options = "f:u:1:2:b:o:c:t:k:" ;
static struct option long_options[] = {
			{ "debug-ns", required_argument, 0, 10000 },
			{ "trimLevel", required_argument, 0, 10001 },
			{ "barcode", required_argument, 0, 10002 },
			{ "keepNoBarcode", no_argument, 0, 10003 },
			{ "UMI", required_argument, 0, 10004},
			{ "skipMateExtension", no_argument, 0, 10005},
			{ "minHitLen", required_argument, 0, 10006},
			{ "contigMinCov", required_argument, 0, 10007},
			{ "cgeneEnd", required_argument, 0, 10008},
			{ (char *)0, 0, 0, 0} 
			} ;

struct _sortRead
{
	char *id ;
	char *read ;
	char *qual ;

	// global k-mer count
	int minCnt ;
	int medianCnt ;
	float avgCnt ;

	// barcode-wise k-mer count
	int barcodeMinCnt ;
	int barcodeMedianCnt ;
	float barcodeAvgCnt ;

	//double avgQual ;
	//int germlineMatchCnt ; 

	int len ;

	int strand ;

	int mateIdx ; // the index of its mate pair.
	int info ; // some random information, such as it's index in orginal array.

	int barcode ; // this is for cell id
	int umi ; // this is for mrna id.

	struct _overlap geneOverlap[4] ; // Rough annotation of the read

	_sortRead()
	{
		id = read = qual = NULL ;
		minCnt = medianCnt = barcodeMinCnt = barcodeMedianCnt = 0 ;
		avgCnt = barcodeAvgCnt = 0 ;
		len = strand = 0 ;
		mateIdx = info = -1 ;
		barcode = umi = 1 ;
	}

	bool operator<( const struct _sortRead &b ) const 
	{
		if ( minCnt != b.minCnt )
			return minCnt > b.minCnt ;
		else if ( medianCnt != b.medianCnt )
			return medianCnt > b.medianCnt ;
		else if ( avgCnt != b.avgCnt )
			return avgCnt > b.avgCnt ;
		else if ( len != b.len )
			return len > b.len ;
		//else if (avgQual != b.avgQual)
		//	return avgQual > b.avgQual ;
		//else if (germlineMatchCnt != b.germlineMatchCnt)
		//	return germlineMatchCnt > b.germlineMatchCnt ;
		else
		{
			int tmp = strcmp( read, b.read ) ;
			if ( tmp != 0 )
				return tmp < 0;
			else
				return strcmp( id, b.id ) < 0 ;
		}
	}
} ;

bool CompReadWithBarcode(const struct _sortRead &a, const struct _sortRead &b)
{
	if (a.barcode != -1 && a.barcode != b.barcode)
		return a.barcode < b.barcode ;
	else if ( a.barcode != -1 && b.barcode != -1 && a.barcodeMinCnt != b.barcodeMinCnt )
		return a.barcodeMinCnt > b.barcodeMinCnt ;
	else
		return a < b ;
}

struct _quickAnnotateReadsThreadArg
{
	int tid ;
	int threadCnt ; 

	SeqSet *refSet ;
	std::vector<struct _sortRead> *pSortedReads ;
	int readCnt ;
} ;

// The structure used for threads related to kmer count, including read trimming and barcode-level kmer count
struct _kmerCountThreadArg
{
	int tid ;
	int threadCnt ;

	int countKmer ;
	KmerCount *pKmerCount ;
	SeqSet *pSeqSet ;
	std::vector<struct _sortRead> *pRead1Batch ;
	std::vector<struct _sortRead> *pRead2Batch ;
	
	std::vector<struct _sortRead> *pReads ;
	int readCnt ;
	int maxReadLen ;

	int trimLevel ;
} ;

struct _assignReadsThreadArg
{
	int tid ;
	int threadCnt ;

	SeqSet *seqSet ;
	std::vector<struct _assignRead> *pAssembledReads ;
	int assembledReadCnt ;
} ;


bool CompSortReadById( const struct _Read &a, const struct _Read &b )
{
	return strcmp( a.id, b.id ) < 0 ;
}

bool IsLowComplexity( const char *seq )
{
	int cnt[5] = {0, 0, 0, 0, 0} ;
	int i ;
	for ( i = 0 ; seq[i] ; ++i )
	{
		if ( seq[i] == 'N' )
			++cnt[4] ;
		else
			++cnt[ nucToNum[ seq[i] - 'A' ] ] ;
	}

	if ( cnt[0] >= i / 2 || cnt[1] >= i / 2 || cnt[2] >= i / 2 || cnt[3] >= i / 2 || cnt[4] >= i / 10 )
		return true ;

	int lowCnt = 0 ; 
	for ( i = 0 ; i < 4 ; ++i )
		if ( cnt[i] <= 2 )
			++lowCnt ;
	if ( lowCnt >= 2 )
		return true ;
	return false ;
}

void PrintLog( const char *fmt, ... )
{
	va_list args ;
	va_start( args, fmt ) ;
	vsprintf( buffer, fmt, args ) ;

	time_t mytime = time(NULL) ;
	struct tm *localT = localtime( &mytime ) ;
	char stime[500] ;
	strftime( stime, sizeof( stime ), "%c", localT ) ;
	fprintf( stderr, "[%s] %s\n", stime, buffer ) ;
}

// Process one read (pair), and append the result to reads
// This includes readthrough, mate pair merging, and kmer count
// The r1, r2's content can be modified in the process, including release
// origR1, origR2 and the content from the original input 
void ProcessRead(struct _sortRead &r1, struct _sortRead &r2, int countKmer, KmerCount &kmerCount, SeqSet &seqSet, std::vector<struct _sortRead> &reads)
{
	int j, k ;
	int rWeight = 1 ;
	// Check chimeric. After reverse-complement, the mate should be after the current read.
	if (r2.read != NULL)
	{
		int flen = strlen(r2.read) ;
		int slen = strlen(r1.read) ;
		seqSet.ReverseComplementInPlace(r2.read, flen ) ;
		if (r2.qual != NULL)
		{
			// Reverse the mateR's quality score
			for (j = 0, k = flen - 1 ; j < k ; ++j, --k)
			{
				char tmp = r2.qual[j] ;
				r2.qual[j] = r2.qual[k] ;
				r2.qual[k] = tmp ;
			}
		}
		int minOverlap = ( flen + slen ) / 10 ;
		int minOverlap2 = ( flen + slen ) / 20 ;
		if ( minOverlap > 31 )
			minOverlap = 31 ;
		if ( minOverlap2 > 31 )
			minOverlap2 = 31 ;
		int offset = -1 ;
		int bestMatchCnt = -1 ;

		/*if ( AlignAlgo::IsMateOverlap( mateR.read, flen, nr.read, slen, minOverlap, offset, bestMatchCnt ) >= 0 )
			{
			printf( "Outie\n%s\n%s\n", nr.read, mateR.read ) ;
			}
			else if ( AlignAlgo::IsMateOverlap( nr.read, slen, mateR.read, flen, minOverlap, offset, bestMatchCnt ) >= 0 )
			{
			printf( "Innie\n%s\n%s\n", nr.read, mateR.read ) ;
			}
			else
			printf( "Uncertain\n" ) ;*/

		int overlapSize = AlignAlgo::IsMateOverlap( r2.read, flen, r1.read, slen, minOverlap, offset, bestMatchCnt, false ) ;
		if ( overlapSize >= 0 )
		{
			// Wrong order (read-through) happened, 
			// Only keep the overlapped portion.
			r1.read[ overlapSize ] = '\0' ;
			if ( r1.qual != NULL )
			{
				r1.qual[ overlapSize ] = '\0' ;

				for ( j = 0 ; j < overlapSize ; ++j )
				{
					if ( r2.qual[j + offset] > r1.qual[j] || r1.read[j] == 'N' )
					{
						r1.read[j] = r2.read[j + offset] ;
						r1.qual[j] = r2.qual[j + offset] ;
					}
				}
			}

			// Release the mate.
			free( r2.read ) ; free( r2.id ) ; 
			if (r2.qual)
				free(r2.qual) ;
			r2.read = NULL ;
		}
		else if (  ( overlapSize = 
					AlignAlgo::IsMateOverlap( r1.read, slen, r2.read, flen, minOverlap2, offset, bestMatchCnt ) ) >= 0 )
		{
			// Merge paired-end reads
			if ( bestMatchCnt >= 0.95 * overlapSize )
			{
				char *r = (char *)malloc( sizeof( char ) * ( slen + flen + 1 ) ) ;
				char *q = (char *)malloc( sizeof( char ) * ( slen + flen + 1 ) ) ;
				for ( j = 0 ; j < flen ; ++j )
				{
					r[ offset + j ] = r2.read[j] ;
					q[ offset + j ] = r2.qual[j] ;
				}
				int len = offset + j ;
				for ( j = 0 ; j < slen && j < len ; ++j )
				{
					if ( j < offset || r1.qual[j] >= q[j] - 14 || r[j] == 'N' )
					{
						r[j] = r1.read[j] ;
						q[j] = r1.qual[j] ;
					}
				}

				//if ( j > len ) 
				//	len = j ;
				r[len] = q[len] = '\0' ;

				if ( 0 )//len > 4 )
				{
					for ( j = 2 ; j <= len ; ++j )
					{
						r[j - 2] = r[j] ;
						q[j - 2] = q[j] ;
					}
					r[len - 4] = '\0' ;
					q[len - 4] = '\0' ;
					len -= 4 ;
				}

				free( r2.read ) ; free( r2.id ) ; 
				if (r2.qual)
						free( r2.qual ) ;
				r2.read = NULL ;
				free( r1.read ) ; free( r1.qual ) ;
				r1.read = r ;
				r1.qual = q ;
				++rWeight ;
			}
			else
			{
				// Discard one of the read
				bool useFirst = true ;
				if ( r1.qual != NULL )
				{
					double avgQualR = 0, avgQualMate = 0 ;
					for ( j = offset ; j < slen ; ++j )
						avgQualR += r1.qual[j] - 32 ;
					for ( j = flen - 1 ; j >= flen - overlapSize ; --j )
						avgQualMate += r2.qual[j] - 32 ;

					avgQualR /= overlapSize ; avgQualMate /= overlapSize ;
					if ( avgQualR + 10 < avgQualMate )
						useFirst = false ;
				}
				if ( useFirst )
				{
					free( r2.read ) ; free( r2.id ) ; 
					if (r2.qual)
						free( r2.qual ) ;
					r2.read = NULL ;
				}
				else
				{
					free( r1.read ) ; 
					if (r1.qual)
						free(r1.qual) ;
					free( r2.id ) ;

					r1.read = r2.read ;
					seqSet.ReverseComplementInPlace(r1.read, flen) ;
					r1.qual = r2.qual ;
					/* Not sure why, keeping it reversed seems better
					if (r1.qual != NULL)
					{
						for (j = 0, k = flen - 1 ; j < k ; ++j, --k)
						{
							char tmp = r1.qual[j] ;
							r1.qual[j] = r1.qual[k] ;
							r1.qual[k] = tmp ;
						}
					}*/
					r2.read = NULL ;
				}
			}
		}
		else // Nothing changes
		{
			seqSet.ReverseComplementInPlace(r2.read, flen) ;
			if (r2.qual != NULL)
			{
				for (j = 0, k = flen - 1 ; j < k ; ++j, --k)
				{
					char tmp = r2.qual[j] ;
					r2.qual[j] = r2.qual[k] ;
					r2.qual[k] = tmp ;
				}
			}
		}
	} // End of testing read-through or overlapped paired-end.

	if ( !IsLowComplexity(r1.read) )
	{
		reads.push_back( r1 ) ;
		if ( countKmer )
			kmerCount.AddCount( r1.read ) ;

		if ( rWeight == 2 )
		{
			struct _sortRead wr ;
			wr = r1 ;
			wr.read = strdup( r1.read ) ;
			if ( r1.qual != NULL )
				wr.qual = strdup( r1.qual ) ;
			else
				wr.qual = NULL ;
			int len = strlen( r1.id ) ;
			wr.id = ( char * )malloc( sizeof( char ) * ( len + 3 ) ) ;
			strcpy( wr.id, r1.id ) ;
			wr.id[len] = '.' ;
			wr.id[len + 1] = '1' ;
			wr.id[len + 2] = '\0' ;

			reads.push_back( wr ) ;
			if ( countKmer )
				kmerCount.AddCount( wr.read ) ;
		}
	}
	else
	{
		free( r1.read ) ; free( r1.id ) ; 
		if (r1.qual)
			free( r1.qual ) ;
	}

	if (r2.read != NULL )
	{
		if ( !IsLowComplexity(r2.read) )
		{
			reads.push_back(r2) ;
			if ( countKmer )
				kmerCount.AddCount(r2.read ) ;
		}
		else
		{
			free( r2.read ) ; free( r2.id ) ; 
			if (r2.qual)
				free( r2.qual ) ;
		}
	}
}

void *ProcessReads_Thread(void *pArg)
{
	struct _kmerCountThreadArg &arg = *((struct _kmerCountThreadArg*)pArg) ;
	std::vector< struct _sortRead> &reads = *arg.pReads ;
	std::vector<struct _sortRead> &read1Batch = *arg.pRead1Batch ;	
	std::vector<struct _sortRead> &read2Batch = *arg.pRead2Batch ;	
	KmerCount &kmerCount = *arg.pKmerCount ;
	SeqSet &seqSet = *arg.pSeqSet ;

	int readCnt = read1Batch.size() ;
	int i, start, end ;
	
	int block = readCnt / arg.threadCnt ;
	if (readCnt % arg.threadCnt)
		++block ;
	start = arg.tid * block ;
	end = start + block - 1 ;
	if (end >= readCnt)
		end = readCnt - 1 ;

	for (i = start ; i <= end ; ++i)
		ProcessRead(read1Batch[i], read2Batch[i], arg.countKmer, kmerCount, seqSet, reads) ;

	pthread_exit(NULL) ;
}

void ProcessReadsWithThreads(std::vector<struct _sortRead> &r1Batch, std::vector<struct _sortRead> &r2Batch, int countKmer, KmerCount &kmerCount, SeqSet &seqSet, std::vector<struct _sortRead> &reads, int threadCnt)
{
	int i ;
	pthread_t *threads = new pthread_t[ threadCnt ] ;
	struct _kmerCountThreadArg *args = new struct _kmerCountThreadArg[threadCnt] ;
	pthread_attr_t attr ;

	pthread_attr_init( &attr ) ;
	pthread_attr_setdetachstate( &attr, PTHREAD_CREATE_JOINABLE ) ;

	for ( i = 0 ; i < threadCnt ; ++i )
	{
		args[i].tid = i ;
		args[i].threadCnt = threadCnt ;
		args[i].pRead1Batch = &r1Batch ;
		args[i].pRead2Batch = &r2Batch ;
		args[i].pReads = new std::vector<struct _sortRead> ;
		args[i].countKmer = countKmer ;
		args[i].pKmerCount = &kmerCount ;
		args[i].pSeqSet = &seqSet ;
		pthread_create( &threads[i], &attr, ProcessReads_Thread, (void *)( args + i ) ) ;
	}

	for ( i = 0 ; i < threadCnt ; ++i )
		pthread_join( threads[i], NULL ) ;

	for (i = 0 ; i < threadCnt ; ++i)
	{
		reads.insert(reads.end(), args[i].pReads->begin(), args[i].pReads->end()) ;
		delete args[i].pReads ;
	}

	delete[] threads ;
	delete[] args ;
}

void *QuickAnnotateReads_Thread( void *pArg )
{
	struct _quickAnnotateReadsThreadArg &arg = *( (struct _quickAnnotateReadsThreadArg *)pArg ) ;
	int start, end ;
	int i, j ;
	std::vector< struct _sortRead> &sortedReads = *arg.pSortedReads ;
	start = arg.readCnt / arg.threadCnt * arg.tid ;
	end = start + arg.readCnt / arg.threadCnt ;
	if ( arg.tid == arg.threadCnt - 1 )
		end = arg.readCnt ;
	struct _overlap geneOverlap[4] ;
	for ( i = start ; i < end ; ++i )
	{
		if ( i == start || strcmp( sortedReads[i].read, sortedReads[i - 1].read ) )
			arg.refSet->AnnotateRead( sortedReads[i].read, 0, geneOverlap, NULL, NULL ) ;
		for ( j = 0 ; j < 4 ; ++j )
			sortedReads[i].geneOverlap[j] = geneOverlap[j] ;
	}
	pthread_exit( NULL ) ;
}

void *GetReadCountStatsAndTrim_Thread(void *pArg)
{
	int i ;
	struct _kmerCountThreadArg &arg = *((struct _kmerCountThreadArg*)pArg) ;
	std::vector< struct _sortRead> &reads = *arg.pReads ;
	int readCnt = arg.readCnt ;
	
	KmerCount &kmerCount = *arg.pKmerCount ;
	int *buffer = new int[arg.maxReadLen] ;
	for (i = arg.tid ; i < readCnt ; i += arg.threadCnt)
	{
		if ( arg.trimLevel == 0 )
			kmerCount.GetCountStatsAndTrim( reads[i].read, NULL, 
					reads[i].minCnt, reads[i].medianCnt, reads[i].avgCnt, buffer ) ;
		else
			kmerCount.GetCountStatsAndTrim( reads[i].read, reads[i].qual, 
					reads[i].minCnt, reads[i].medianCnt, reads[i].avgCnt, buffer ) ;

		if ( reads[i].qual != NULL )
		{
			free( reads[i].qual ) ;
			reads[i].qual = NULL ;
		}

		if ( reads[i].read[0] == '\0' )
		{
			free( reads[i].read ) ;
			free( reads[i].id ) ;
			reads[i].read = NULL ;
		}
	}
	delete[] buffer ;
	pthread_exit(NULL) ;
}

void *BarcodeKmerCount_Thread(void *pArg)
{
	int i, j, k ;
	struct _kmerCountThreadArg &arg = *((struct _kmerCountThreadArg*)pArg) ;
	std::vector< struct _sortRead> &reads = *arg.pReads ;
	int readCnt = arg.readCnt ;

	KmerCount barcodeKmerCount(21, 23) ;
	barcodeKmerCount.SetBuffer( arg.maxReadLen ) ;
	for (i = 0 ; i < readCnt ; )
	{
		for (j = i + 1 ; j < readCnt ; ++j)
		{
			if (reads[i].barcode != reads[j].barcode)
				break ;
		}
		
		if (reads[i].barcode % arg.threadCnt != arg.tid)
		{
			i = j ;
			continue ;
		}

		for (k = i ; k < j ; ++k)
			barcodeKmerCount.AddCount(reads[k].read) ;

		for (k = i ; k < j ; ++k)
			barcodeKmerCount.GetCountStatsAndTrim(reads[k].read, NULL,
					reads[k].barcodeMinCnt, reads[k].barcodeMedianCnt,
					reads[k].barcodeAvgCnt ) ;  

		barcodeKmerCount.Clear() ;
		i = j ;
	}
	barcodeKmerCount.Release() ;
	pthread_exit(NULL) ;
}

void *AssignReads_Thread( void *pArg )
{
	struct _assignReadsThreadArg &arg = *( (struct _assignReadsThreadArg *)pArg ) ;
	int start, end ;
	int i ;
	std::vector< struct _assignRead> &assembledReads = *arg.pAssembledReads ;
	start = arg.assembledReadCnt / arg.threadCnt * arg.tid ;
	end = start + arg.assembledReadCnt / arg.threadCnt ;
	if ( arg.tid == arg.threadCnt - 1 )
		end = arg.assembledReadCnt ;
	struct _overlap assign ;
	for ( i = start ; i < end ; ++i )
	{
		if ( i == start || strcmp( assembledReads[i].read, assembledReads[i - 1].read ) )
			arg.seqSet->AssignRead( assembledReads[i].read, assembledReads[i].overlap.strand, 
				assembledReads[i].barcode, assign ) ;
		assembledReads[i].overlap = assign ;	
	}
	pthread_exit( NULL ) ;
}

int main( int argc, char *argv[] )
{
	int i, j, k ;

	if ( argc <= 1 )
	{
		fprintf( stderr, "%s", usage ) ;
		return 0 ;
	}

	int c, option_index ;
	option_index = 0 ;
	int indexKmerLength = 9 ;
	int changeKmerLengthThreshold = 4096 ;
	SeqSet seqSet( indexKmerLength ) ; // Only hold the novel seq.
	SeqSet refSet( 9 ) ;
	KmerCount kmerCount( 21 ) ;
	char outputPrefix[1024] = "trust" ;

	ReadFiles reads ;
	ReadFiles mateReads ;
	ReadFiles barcodeFile, umiFile ;
	bool countMyself = true ; // Count using our procedure, or read from some kmer count tool's output.
	int maxReadLen = -1 ;
	int firstReadLen = -1 ;
	int trimLevel = 1 ;
	int minHitLen = -1 ;
	bool hasMate = false ;
	bool hasBarcode = false ;
	bool hasUmi = false ;
	int constantGeneEnd = 200 ;
	bool keepMissingBarcode = false ;
	int threadCnt = 1 ;
	bool skipMateExtension = false ;
	int contigMinCov = 0 ;

	while ( 1 )
	{
		c = getopt_long( argc, argv, short_options, long_options, &option_index ) ;
		
		if ( c == -1 )
			break ;

		if ( c == 'f' )
		{
			//seqSet.InputRefFa( optarg ) ;
			refSet.InputRefFa( optarg ) ;
		}
		else if ( c == 'u' )
		{
			reads.AddReadFile( optarg, false ) ;
		}
		else if ( c == '1' )
		{
			reads.AddReadFile( optarg, true ) ;
			hasMate = true ;
		}
		else if ( c == '2' )
		{
			mateReads.AddReadFile( optarg, true ) ;
			hasMate = true ;
		}
		else if ( c == 'o' )
		{
			strcpy( outputPrefix, optarg ) ;
		}
		else if ( c == 'c' )
		{
			kmerCount.AddCountFromFile( optarg ) ;
			countMyself = false ;
			PrintLog( "Read in the kmer count information from %s", optarg ) ;
		}
		else if ( c == 't' )
		{
			threadCnt = atoi( optarg ) ;
		}
		else if ( c == 'k' )
		{
			indexKmerLength = atoi( optarg ) ;
			seqSet.ChangeKmerLength( indexKmerLength ) ;
		}
		else if ( c == 10000 ) //debug-ns
		{
			seqSet.InputNovelFa( optarg ) ;
		}
		else if ( c == 10001) //trimLevel
		{
			trimLevel = atoi( optarg ) ;
		}
		else if ( c == 10002 ) // barcode 
		{
			barcodeFile.AddReadFile( optarg, false ) ;
			hasBarcode = true ;
		}
		else if ( c == 10003 ) // keepNoBarcode
		{
			keepMissingBarcode = true ;
		}
		else if ( c == 10004 ) // umi
		{
			umiFile.AddReadFile( optarg, false ) ;
			hasUmi = true ;
		}
		else if ( c == 10005 ) // skipMateExtension
		{
			skipMateExtension = true ;
		}
		else if ( c == 10006 ) // minHitLen
		{
			minHitLen = atoi( optarg ) ;
		}
		else if ( c == 10007 )
		{
			contigMinCov = atoi( optarg ) ;
		}
		else if ( c == 10008)
		{
			constantGeneEnd = atoi( optarg ) ;
		}
		else
		{
			fprintf( stderr, "%s", usage ) ;
			return EXIT_FAILURE ;
		}
	}
	/*char align[10000] ;
	char t[] = "ACGGGGTTTTTT" ;
	char p[] = "ACTTTTTTGGGG" ;
	AlignAlgo::GlobalAlignment( t, strlen( t ), p, strlen( p ), align ) ;
	AlignAlgo::VisualizeAlignment( t, strlen( t ), p, strlen( p ), align ) ; */
	
	if ( refSet.Size() == 0 )
	{
		fprintf( stderr, "Need to use -f to specify the receptor genome sequence.\n" ) ;
		return EXIT_FAILURE ;
	}
	PrintLog( "Start to assemble reads." ) ;

	if ( trimLevel > 1 ) 
	{
		// We need a more careful rough annotation.
		refSet.ChangeKmerLength( 7 ) ;
	} 
	refSet.SetHitLenRequired( 17 ) ;

	std::vector< struct _sortRead > sortedReads ;
	std::map<std::string, int> barcodeStrToInt ;
	std::map<std::string, int> umiStrToInt ;
	std::vector<std::string> barcodeIntToStr ;
	
	// Relate to barcode release
	SimpleVector<int> barcodeTotalReadCount ; // Number of reads for this barcode
	SimpleVector<int> barcodeReadCount ;
	std::map<int, int> finishedBarcodes ; 

	std::vector<struct _sortRead> readBuffer, mateReadBuffer ;
	i = 0 ;
	if (threadCnt > 1)
		kmerCount.SetPthreadLocks() ;
	while ( reads.Next() )
	{
		/*struct _overlap geneOverlap[4] ;
		refSet.AnnotateRead( reads.seq, 1, geneOverlap, buffer ) ;
	
		if ( geneOverlap[0].seqIdx + geneOverlap[1].seqIdx + geneOverlap[2].seqIdx + geneOverlap[3].seqIdx == -4 )
		{
			continue ;
		}*/

		int barcode = -1 ;
		int umi = -1 ;
		if ( hasBarcode )
		{
			barcodeFile.Next() ;
			
			if ( !strcmp( barcodeFile.seq, "missing_barcode" ) && !keepMissingBarcode )
			{
				if ( hasMate )
					mateReads.Next() ;
				if ( hasUmi )
					umiFile.Next() ;
				continue ;
			}

			std::string s( barcodeFile.seq ) ;
			if ( barcodeStrToInt.find( s ) != barcodeStrToInt.end() )
				barcode = barcodeStrToInt[s] ;
			else
			{
				barcode = barcodeIntToStr.size() ;
				barcodeStrToInt[s] = barcode ;
				barcodeIntToStr.push_back( s ) ;
			}

			if (contigMinCov > 0)
			{
				if (barcode >= barcodeTotalReadCount.Size())
					barcodeTotalReadCount.PushBack(1) ;
				else
					++barcodeTotalReadCount[barcode] ;
			}
		}

		if ( hasUmi )
		{
			umiFile.Next() ;
			std::string s(umiFile.seq) ;
			if ( umiStrToInt.find(s) != umiStrToInt.end())
				umi = umiStrToInt[s] ;
			else
			{
				umi = umiStrToInt.size() ;
				umiStrToInt[s] = umi ;
			}
		}

		struct _sortRead nr, mateR ;
		nr.read = strdup( reads.seq ) ;
		nr.id = strdup( reads.id ) ;
		if ( reads.qual != NULL )
			nr.qual = strdup( reads.qual ) ;
		else
			nr.qual = NULL ;
		nr.barcode = barcode ;
		nr.umi = umi ;

		++i ;

		if ( threadCnt == 1 && countMyself && i % 100000 == 0 )
			PrintLog( "Read in and count kmers for %d reads.", i ) ;
		else if ( threadCnt == 1 && !countMyself && i % 1000000 == 0 )
			PrintLog( "Read in %d reads.", i ) ;
		
		if ( firstReadLen == -1 )
		{
			firstReadLen = strlen( reads.seq ) ;
		}

		mateR.read = NULL ;
		if ( mateReads.Next() )
		{
			mateR.read = strdup( mateReads.seq ) ;
			mateR.id = strdup( mateReads.id ) ;
			mateR.barcode = barcode ;
			mateR.umi = umi ;
			if ( mateReads.qual != NULL )
			{
				mateR.qual = strdup( mateReads.qual ) ;
			}
			else
				mateR.qual = NULL ;
			
			++i ;

			if ( threadCnt == 1 && countMyself && i % 100000 == 0 )
				PrintLog( "Read in and count kmers for %d reads.", i ) ;
			else if ( threadCnt == 1 && !countMyself && i % 1000000 == 0 )
				PrintLog( "Read in %d reads.", i ) ;
		}
		else if ( hasMate ) 
		{
			fprintf( stderr, "The two mate-pair read files have different number of reads.\n" ) ;
			exit( 1 ) ;
		}

		if (threadCnt == 1)
			ProcessRead(nr, mateR, countMyself, kmerCount, seqSet, sortedReads) ;
		else
		{
			readBuffer.push_back(nr) ;
			mateReadBuffer.push_back(mateR) ;

			if ((int)readBuffer.size() >= 100000 * threadCnt)
			{
				ProcessReadsWithThreads(readBuffer, mateReadBuffer, countMyself, kmerCount, seqSet, sortedReads, threadCnt) ;
				readBuffer.clear() ;
				mateReadBuffer.clear() ;
			
				if (countMyself)
					PrintLog( "Read in and count kmers for %d reads.", i ) ;
				else 
					PrintLog( "Read in %d reads.", i ) ;
			}
		}
	}

	if (readBuffer.size() > 0)
	{
		ProcessReadsWithThreads(readBuffer, mateReadBuffer, countMyself, kmerCount, seqSet, sortedReads, threadCnt) ;
	}

	int readCnt = sortedReads.size() ;
	
	maxReadLen = 0 ;
	for (i = 0 ; i < readCnt ; ++i)
	{
		sortedReads[i].len = strlen( sortedReads[i].read ) ;
		if (sortedReads[i].len > maxReadLen)
			maxReadLen = sortedReads[i].len ;
		//printf("%s %s %s\n", sortedReads[i].id, sortedReads[i].read, sortedReads[i].qual) ;
	}

	// empty file
	if ( readCnt <= 0 )
	{
		sprintf( buffer, "%s_raw.out", outputPrefix ) ;
		FILE *fp = fopen( buffer, "w" ) ;
		fclose(fp) ;

		sprintf( buffer, "%s_assembled_reads.fa", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
		fclose(fp) ;

		sprintf( buffer, "%s_final.out", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
		fclose(fp) ;

		return 0 ;
	}
#ifdef DEBUG
	printf( "Finish read in the reads and kmer count.\n") ;
#endif

	// Remove the reads from the barcode with too few reads
	if (contigMinCov > 0)
	{
		for (i = 0 ; i < readCnt ; ++i)
		{
			if (sortedReads[i].barcode != -1 && barcodeTotalReadCount[sortedReads[i].barcode] < contigMinCov)
			{
				free( sortedReads[i].read ) ;
				free( sortedReads[i].id ) ;
				if (sortedReads[i].qual != NULL)
				{
					free(sortedReads[i].qual) ;
					sortedReads[i].qual = NULL ;
				}
				sortedReads[i].read = NULL ;
			}
		}

		k = 0 ;
		for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( sortedReads[i].read == NULL )
				continue ;
			sortedReads[k] = sortedReads[i] ;
			++k ;
		}
		readCnt = k ;
		sortedReads.resize( k ) ;
	}
	kmerCount.SetBuffer( maxReadLen ) ;
	if (threadCnt == 1)
	{
		for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( trimLevel == 0 )
				kmerCount.GetCountStatsAndTrim( sortedReads[i].read, NULL, 
						sortedReads[i].minCnt, sortedReads[i].medianCnt, sortedReads[i].avgCnt ) ;
			else
				kmerCount.GetCountStatsAndTrim( sortedReads[i].read, sortedReads[i].qual, 
						sortedReads[i].minCnt, sortedReads[i].medianCnt, sortedReads[i].avgCnt ) ;

			if ( sortedReads[i].qual != NULL )
			{
				/*char *qual = sortedReads[i].qual ;
					double avgQual = 0 ;
					for (j = 0 ; qual[j] ; ++j)
					avgQual += qual[j] ;
					sortedReads[i].avgQual = avgQual / j ;*/

				free( sortedReads[i].qual ) ;
				sortedReads[i].qual = NULL ;
			}
			//else
			//	sortedReads[i].avgQual = 0 ;


			if ( sortedReads[i].read[0] == '\0' )
			{
				free( sortedReads[i].read ) ;
				free( sortedReads[i].id ) ;
				sortedReads[i].read = NULL ;
			}
		}
	}
	else
	{
		pthread_t *threads = new pthread_t[ threadCnt ] ;
		struct _kmerCountThreadArg *args = new struct _kmerCountThreadArg[threadCnt] ;
		pthread_attr_t attr ;

		pthread_attr_init( &attr ) ;
		pthread_attr_setdetachstate( &attr, PTHREAD_CREATE_JOINABLE ) ;

		for ( i = 0 ; i < threadCnt ; ++i )
		{
			args[i].tid = i ;
			args[i].threadCnt = threadCnt ;
			args[i].pReads = &sortedReads ;
			args[i].pKmerCount = &kmerCount ;
			args[i].readCnt = readCnt ;
			args[i].maxReadLen = maxReadLen ;
			args[i].trimLevel = trimLevel ;
			pthread_create( &threads[i], &attr, GetReadCountStatsAndTrim_Thread, (void *)( args + i ) ) ;
		}

		for ( i = 0 ; i < threadCnt ; ++i )
			pthread_join( threads[i], NULL ) ;

		delete[] threads ;
		delete[] args ;
	}
	
	k = 0 ;
	for ( i = 0 ; i < readCnt ; ++i )
	{
		if ( sortedReads[i].read == NULL )
			continue ;
		sortedReads[k] = sortedReads[i] ;
		sortedReads[k].len = strlen( sortedReads[i].read ) ;
		++k ;
	}
	readCnt = k ;
	sortedReads.resize( k ) ;

	PrintLog( "Found %i reads.", k ) ;
#ifdef DEBUG
	printf( "Finish put in the read kmer count.\n" ) ;
#endif
	
	kmerCount.Release() ;

	// Remaining intializations
	for ( i = 0 ; i < readCnt ; ++i )
	{
		sortedReads[i].info = i ;
		sortedReads[i].mateIdx = -1 ;
	}
	for ( i = 0 ; i < readCnt - 1 ; ++i )
	{
		if ( !strcmp( sortedReads[i].id, sortedReads[i + 1].id ) )
		{
			sortedReads[i].mateIdx = i + 1 ;
			sortedReads[i + 1].mateIdx = i ;
			++i ;	
		}
	}
	
	std::sort( sortedReads.begin(), sortedReads.end() ) ;
	PrintLog( "Finish sorting the reads." ) ;

	// Quickly annoate the reads.
	if ( trimLevel > 1 )
		refSet.SetRadius(0) ;
	if ( threadCnt <= 1 )
	{
		struct _overlap geneOverlap[4] ;
		for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( i == 0 || strcmp( sortedReads[i].read, sortedReads[i - 1].read ) )
				refSet.AnnotateRead( sortedReads[i].read, 0, geneOverlap, NULL, NULL ) ;
			for ( j = 0 ; j < 4 ; ++j )
				sortedReads[i].geneOverlap[j] = geneOverlap[j] ;
		}
	}
	else
	{
		pthread_t *threads = new pthread_t[ threadCnt ] ;
		struct _quickAnnotateReadsThreadArg *args = new struct _quickAnnotateReadsThreadArg[threadCnt] ;
		pthread_attr_t attr ;
		
		pthread_attr_init( &attr ) ;
		pthread_attr_setdetachstate( &attr, PTHREAD_CREATE_JOINABLE ) ;
		
		for ( i = 0 ; i < threadCnt ; ++i )
		{
			args[i].tid = i ;
			args[i].threadCnt = threadCnt ;
			args[i].refSet = &refSet ;
			args[i].pSortedReads = &sortedReads ;
			args[i].readCnt = readCnt ;
			pthread_create( &threads[i], &attr, QuickAnnotateReads_Thread, (void *)( args + i ) ) ;
		}

		for ( i = 0 ; i < threadCnt ; ++i )
			pthread_join( threads[i], NULL ) ;
	

		delete[] threads ;
		delete[] args ;
	}
	PrintLog( "Finish rough annotations." ) ;
	
	if (hasBarcode)
	{
		// Sort barcode afterwards so annotation can reuse same sequence from different barcode
		std::sort( sortedReads.begin(), sortedReads.end(), CompReadWithBarcode ) ;

		PrintLog("Get barcode-wise kmer count.") ;
		if (threadCnt == 1)
		{
			KmerCount barcodeKmerCount(21, 23) ;
			barcodeKmerCount.SetBuffer(maxReadLen) ;
			for (i = 0 ; i < readCnt ; )
			{
				for (j = i + 1 ; j < readCnt ; ++j)
				{
					if (sortedReads[i].barcode != sortedReads[j].barcode)
						break ;
				}

				for (k = i ; k < j ; ++k)
					barcodeKmerCount.AddCount(sortedReads[k].read) ;

				for (k = i ; k < j ; ++k)
				{
					barcodeKmerCount.GetCountStatsAndTrim(sortedReads[k].read, NULL,
							sortedReads[k].barcodeMinCnt, sortedReads[k].barcodeMedianCnt,
							sortedReads[k].barcodeAvgCnt ) ;  
				}

				barcodeKmerCount.Clear() ;
				i = j ;
			}
		}    
		else
		{
			pthread_t *threads = new pthread_t[ threadCnt ] ;
			struct _kmerCountThreadArg *args = new struct _kmerCountThreadArg[threadCnt] ;
			pthread_attr_t attr ;

			pthread_attr_init( &attr ) ;
			pthread_attr_setdetachstate( &attr, PTHREAD_CREATE_JOINABLE ) ;

			for ( i = 0 ; i < threadCnt ; ++i )
			{
				args[i].tid = i ;
				args[i].threadCnt = threadCnt ;
				args[i].pReads = &sortedReads ;
				args[i].readCnt = readCnt ;
				args[i].maxReadLen = maxReadLen ;
				pthread_create( &threads[i], &attr, BarcodeKmerCount_Thread, (void *)( args + i ) ) ;
			}

			for ( i = 0 ; i < threadCnt ; ++i )
				pthread_join( threads[i], NULL ) ;

			delete[] threads ;
			delete[] args ;
		}
		PrintLog("Finish barcode-wise kmer count.") ;

		// Resort to follow barcode k-mer count
		for (i = 0 ; i < readCnt ; )
		{
			for (j = i + 1 ; j < readCnt ; ++j)
				if (sortedReads[i].barcode != sortedReads[j].barcode)
					break ;
			if (j - i > 1)		
				std::sort(sortedReads.begin() + i, sortedReads.begin() + j,
						CompReadWithBarcode) ;
			i = j ;
		}
		PrintLog( "Finish re-sorting the reads based on barcode." ) ;
	}
	/*for (i = 0 ; i < readCnt ; ++i)
	{
		int germlineMatchCnt = 0 ;
		struct _overlap *overlaps = sortedReads[i].geneOverlap ;
		for (j = 0 ; j < 4 ; ++j)
			if (overlaps[j].seqIdx != -1)
				germlineMatchCnt += overlaps[j].matchCnt ;
		sortedReads[i].germlineMatchCnt = germlineMatchCnt ;
	}
	TODO: maybe just sort the elements with k-mer supports = 1.
	std::sort( sortedReads.begin(), sortedReads.end() ) ;
	PrintLog( "Finish re-sorting the reads.") ;*/

	// Connect mate pair information
	SimpleVector<int> originToSortedIdx ;
	SimpleVector<bool> goodCandidate ; // Use mate pair alignment information to infer 
					   // whether this read could be a good candidate
	originToSortedIdx.ExpandTo( readCnt ) ;
	goodCandidate.ExpandTo( readCnt ) ;
	for ( i = 0 ; i < readCnt ; ++i )
		originToSortedIdx[ sortedReads[i].info ] = i ;
	for ( i = 0 ; i < readCnt ; ++i )
	{
		if ( sortedReads[i].mateIdx != -1 )
			sortedReads[i].mateIdx = originToSortedIdx[ sortedReads[i].mateIdx ] ;
		goodCandidate[i] = false ;
	}


	if ( trimLevel > 1 && !hasBarcode )
	{
		// use V gene assignment as barcode to speed up.
		for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( sortedReads[i].geneOverlap[0].seqIdx != -1 && sortedReads[i].geneOverlap[0].similarity > 0.95 )
			{
				sortedReads[i].barcode = sortedReads[i].geneOverlap[0].seqIdx ;
				if ( sortedReads[i].mateIdx != -1 )
					sortedReads[ sortedReads[i].mateIdx ].barcode = sortedReads[i].geneOverlap[0].seqIdx ;
			}
		}
		/*for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( sortedReads[i].barcode == -1 )
			{
				free( sortedReads[i].id ) ;
				free( sortedReads[i].read ) ;
				if ( sortedReads[i].qual != NULL )
					free( sortedReads[i].qual ) ;
				sortedReads[i].read = NULL ;
			}
		}*/

		/*for ( i = 0 ; i < readCnt ; ++i )
		{
			if ( sortedReads[i].geneOverlap[0].seqIdx == -1 && sortedReads[i].geneOverlap[2].seqIdx == -1 )
			{
				free( sortedReads[i].id ) ;
				free( sortedReads[i].read ) ;
				if ( sortedReads[i].qual != NULL )
					free( sortedReads[i].qual ) ;
				sortedReads[i].read = NULL ;
			}
		}*/
	}

	// Remove the redudant sequence before V gene.
	for ( i = 0 ; i < readCnt ; ++i )
	{
		if ( sortedReads[i].read == NULL )
			continue ;
		struct _overlap *geneOverlap = sortedReads[i].geneOverlap ;
		if ( geneOverlap[0].seqIdx == -1 )
			continue ;
		bool mayTrim = false ;
		if ( geneOverlap[0].seqStart < 31 && geneOverlap[0].similarity > 0.9 )
			mayTrim  = true ;

		if ( geneOverlap[0].similarity > 0.95 && 
			geneOverlap[0].seqStart <= refSet.GetSeqConsensusLen( geneOverlap[0].seqIdx ) / 3 )
			mayTrim = true ; 
		
		if ( trimLevel > 1 )
			mayTrim = true ;

		if ( !mayTrim )
			continue ;
		
		int trimBase = geneOverlap[0].readStart ;
		if ( trimLevel > 1 && refSet.GetSeqName(geneOverlap[0].seqIdx)[0] == 'T') // Some bad annotation
		{	
			//Since TCR may not require full-length assembly, just V gene assignment, we can
			//  ignore the other part. Comment that out now. We may need those for better 
			//  allele identification or identify some variants
			/*if ( geneOverlap[0].seqEnd + 25 < refSet.GetSeqConsensusLen( geneOverlap[0].seqIdx ) )
				trimBase = geneOverlap[0].readEnd + 2 ;
			else*/
			if ( geneOverlap[0].similarity < 0.97 )
			{
				trimBase = ( geneOverlap[0].readStart + geneOverlap[0].readEnd ) / 2 ;
			}
		}
		//printf("%s\n", sortedReads[i].read) ;	
		//printf("0: %lf %d %d %s: %d %d\n", geneOverlap[0].similarity, mayTrim, trimBase,  
		//	refSet.GetSeqName(geneOverlap[0].seqIdx), geneOverlap[0].readStart, geneOverlap[0].readEnd) ;
		if ( trimBase <= 0 )
			continue ;

		if ( geneOverlap[2].seqIdx != -1 
			&& geneOverlap[2].readStart < trimBase && trimLevel <= 1 )
			continue ;
		if ( geneOverlap[3].seqIdx != -1 
			&& geneOverlap[3].readStart < trimBase && trimLevel <= 1 )
			continue ;

		if ( sortedReads[i].len - trimBase < 31 )
		{
			free( sortedReads[i].id ) ;
			free( sortedReads[i].read ) ;
			if ( sortedReads[i].qual != NULL )
				free( sortedReads[i].qual ) ;
			sortedReads[i].read = NULL ;
			continue ;
		}
		
		if ( geneOverlap[0].strand >= 0 )
		{
			// Remove the first bases
			int len = sortedReads[i].len ;
			for ( j = trimBase ; j <= len ; ++j )
			{
				sortedReads[i].read[j - trimBase] = sortedReads[i].read[j] ;
				if ( sortedReads[i].qual != NULL )
					sortedReads[i].qual[j - trimBase] = sortedReads[i].qual[j] ;

			}
		}
		else
		{
			// Remove the last bases.
			int len = sortedReads[i].len ;
			sortedReads[i].read[len - trimBase] = '\0' ;
			if ( sortedReads[i].qual != NULL )
				sortedReads[i].qual[len - trimBase] = '\0' ;
		}

		for ( j = 0 ; j < 4 ; ++j )
		{
			if ( geneOverlap[j].seqIdx == -1 )
				continue ;
			geneOverlap[j].readStart -= trimBase ;
			geneOverlap[j].readEnd -= trimBase ;
			if ( geneOverlap[j].readStart < 0 )
				geneOverlap[j].readStart = 0 ;
			if ( geneOverlap[j].readEnd < 0 )
			{
				geneOverlap[j].readEnd = 0 ;
				geneOverlap[j].seqIdx = -1 ;
			}
		}
		sortedReads[i].len -= trimBase ;
	}
	
	// Remove redundant sequences after constant gene ; a
	// and also constant gene for non IGH genes for long reads.
	for ( i = 0 ; i < readCnt ; ++i )
	{
		int len = sortedReads[i].len ;
		int gidx = 3 ;
		struct _overlap *geneOverlap = sortedReads[i].geneOverlap ;
		if ( sortedReads[i].read == NULL )
			continue ;
		for ( gidx = 2 ; gidx <= 3 ; ++gidx)
			if ( geneOverlap[gidx].seqIdx != -1 )
				break ;
		if ( gidx > 3 )
			continue ;
		
		if ( gidx == 2 && refSet.GetSeqName(geneOverlap[gidx].seqIdx)[2] == 'H' )
		{
			gidx = 3 ;
			if ( geneOverlap[gidx].seqIdx == -1 )
				continue ;
		}

		bool mayTrim = false ;
		if ( gidx == 3 && geneOverlap[3].seqStart < 9 && geneOverlap[3].similarity > 0.95 )
				//&& geneOverlap[3].readEnd + indexKmerLength >= sortedReads[i].len 
				//&& refSet.GetSeqName( geneOverlap[3].seqIdx )[2] != 'H' )
			mayTrim  = true ;
		if ( trimLevel > 1 )
			mayTrim = true ;
		if ( !mayTrim )
			continue ;
		
		int trimBase = len - geneOverlap[gidx].readEnd - 1 ;
		if ( trimLevel > 1 && refSet.GetSeqName(geneOverlap[gidx].seqIdx)[0] == 'T') // Some bad annotation
		{
			/*if ( geneOverlap[gidx].seqStart >= 25 )
				trimBase = len - ( geneOverlap[gidx].readStart - 1 ) ;
			else*/ 
			if ( geneOverlap[gidx].similarity < 0.97 )
			{
				trimBase = len - ( ( geneOverlap[gidx].readStart + geneOverlap[gidx].readEnd ) / 2 ) - 1 ;
			}
		}
		//printf("%s %d %s\n", sortedReads[i].id, gidx, sortedReads[i].read) ;
		//printf("%d: %lf %d %d %s: %d %d\n", gidx, geneOverlap[gidx].similarity, mayTrim, trimBase,  
		//	refSet.GetSeqName(geneOverlap[gidx].seqIdx), geneOverlap[gidx].readStart, geneOverlap[gidx].readEnd) ;
		/*if ( firstReadLen > 200 && gidx == 3 && trimLevel < 2 && refSet.GetSeqName( geneOverlap[3].seqIdx )[2] != 'H' )
		{
			trimBase = ( len - geneOverlap[3].readStart ) + geneOverlap[3].seqStart ;
		}*/
		
		if ( trimBase <= 0 )
			continue ;

		if ( gidx == 3 && geneOverlap[2].seqIdx != -1 
				&& geneOverlap[2].readStart + trimBase >= sortedReads[i].len && trimLevel <= 1 )
			continue ;
		if ( geneOverlap[0].seqIdx != -1 
				&& geneOverlap[0].readStart + trimBase >= sortedReads[i].len && trimLevel <= 1 )
			continue ;

		if ( sortedReads[i].len - trimBase < 31 )
		{
			free( sortedReads[i].id ) ;
			free( sortedReads[i].read ) ;
			if ( sortedReads[i].qual != NULL )
				free( sortedReads[i].qual ) ;
			sortedReads[i].read = NULL ;
			continue ;
		}

		if ( geneOverlap[gidx].strand < 0 )
		{
			// Remove the first bases
			for ( j = trimBase ; j <= len ; ++j )
			{
				sortedReads[i].read[j - trimBase] = sortedReads[i].read[j] ;
				if ( sortedReads[i].qual != NULL )
					sortedReads[i].qual[j - trimBase] = sortedReads[i].qual[j] ;

			}

			geneOverlap[3].seqIdx = -1 ;
		}
		else
		{
			// Remove the last bases.
			sortedReads[i].read[len - trimBase] = '\0' ;
			if ( sortedReads[i].qual != NULL )
				sortedReads[i].qual[len - trimBase] = '\0' ;
			geneOverlap[3].seqIdx = -1 ;
		}
		for ( j = 0 ; j < 4 ; ++j )
		{
			if ( geneOverlap[j].seqIdx == -1 )
				continue ;
			if ( geneOverlap[j].readStart + trimBase >= len )
			{
				geneOverlap[j].readStart = len - 1 ;
				geneOverlap[j].seqIdx = -1 ;
			}
			if ( geneOverlap[j].readEnd + trimBase >= len )
				geneOverlap[j].readEnd = len - 1 ;
		}
		//printf("%s\n", sortedReads[i].read);
		sortedReads[i].len -= trimBase ;
	}

	// Remove the reads that are too short if this is a long read data set.
	if ( firstReadLen > 200 )
	{
		for ( i = 0 ; i < readCnt ; ++i )
			if ( sortedReads[i].read != NULL && sortedReads[i].len < firstReadLen / 3 )
			{
				free( sortedReads[i].id ) ;
				free( sortedReads[i].read ) ;
				if ( sortedReads[i].qual != NULL )
					free( sortedReads[i].qual ) ;
				sortedReads[i].read = NULL ;
				continue ;
			}

		seqSet.SetIsLongSeqSet( true ) ;
	}
	
	// Remove the sequences whose J gene is not associate with C gene.
	/*for ( i = 0 ; i < readCnt ; ++i )
	  {
	  struct _overlap *geneOverlap = sortedReads[i].geneOverlap ;
	  if ( sortedReads[i].read == NULL )
	  continue ;

	  if ( geneOverlap[0].seqIdx == -1 && geneOverlap[2].seqIdx != -1 && geneOverlap[3].seqIdx == -1 )
	  {
	  if ( geneOverlap[2].similarity > 0.9 && 
	  sortedReads[i].len - geneOverlap[2].readEnd > 50 ) // With 50 base but could not identify C gene
	  {
	  free( sortedReads[i].id ) ;
	  free( sortedReads[i].read ) ;
	  if ( sortedReads[i].qual != NULL )
	  free( sortedReads[i].qual ) ;
	  sortedReads[i].read = NULL ;
	  }
	  }
	  }*/

	k = 0 ;
	for ( i = 0 ; i < readCnt ; ++i )
	{
		if ( sortedReads[i].read != NULL )
		{
			originToSortedIdx[i] = k ;
			sortedReads[k] = sortedReads[i] ;
			++k ;
		}
		else
			originToSortedIdx[i] = -1 ;
	}

	// Don't forget to update the mateidx.
	for ( i = 0 ; i < k ; ++i )
	{
		if ( sortedReads[i].mateIdx != -1 )
			sortedReads[i].mateIdx = originToSortedIdx[ sortedReads[i].mateIdx ] ;
	}

	originToSortedIdx.Release() ;
	sortedReads.resize( k ) ;
	readCnt = k ;

	std::vector<int> rescueReadIdx ;
	std::vector<int> assembledReadIdx ;

	int assembledReadCnt = 0 ;
	int prevAddRet = -1 ;
	
	/*if ( seqSet.GetSeqCnt() > 0 )
	{
		for ( i = 0 ; i < readCnt ; ++i )
			assembledReadIdx.push_back( i ) ;
		readCnt = 0 ;
	}*/
	
	if ( firstReadLen / 2 < 31 )
	{
		int l = firstReadLen / 2 ;
		if ( l < 21 )
			l = 21 ;
		seqSet.SetHitLenRequired( l ) ;
	}

	if ( hasBarcode )
	{
		if ( 1 )
			seqSet.SetHitLenRequired( 13 ) ;
		else
			seqSet.SetHitLenRequired( 17 ) ;

		if (keepMissingBarcode == false)
		{
			seqSet.SetConsiderBarcodeInIndexHash(true) ;
		}
	}
	
	if (minHitLen != -1)
	{
		seqSet.SetHitLenRequired( minHitLen ) ;
	}

	if ( firstReadLen > 200 || trimLevel > 1 )
		changeKmerLengthThreshold /= 2 ;

	if (hasBarcode)
	{
		barcodeTotalReadCount.ExpandTo(barcodeIntToStr.size()) ;
		barcodeReadCount.ExpandTo(barcodeIntToStr.size()) ;
		barcodeTotalReadCount.SetZero(0, barcodeIntToStr.size()) ;
		barcodeReadCount.SetZero(0, barcodeIntToStr.size()) ;
		for (i = 0 ; i < readCnt ; ++i)
		{
			if (sortedReads[i].barcode != -1)
				++barcodeTotalReadCount[ sortedReads[i].barcode ] ;
		}
	}
	
	for ( i = 0 ; i < readCnt ; ++i )
	{
		/*++assembledReadCnt ;
		assembledReadIdx.push_back( i ) ;
		continue ;*/
		static struct _overlap geneOverlap[4] ;

#ifdef DEBUG
		printf( "%s %s %d %lf\n", sortedReads[i].id, sortedReads[i].read, sortedReads[i].minCnt, sortedReads[i].avgCnt ) ;
		fflush( stdout ) ;
#endif
		int addRet = -1 ;
		
		if ( i == 0 || strcmp( sortedReads[i].read, sortedReads[i - 1].read ) 
			|| sortedReads[i].barcode != sortedReads[i - 1].barcode )
		{
			//printf( "new stuff\n" ) ;
			//buffer[0] = '\0' ;
			//refSet.AnnotateRead( sortedReads[i].read, 0, geneOverlap, NULL, NULL ) ;
			for ( j = 0 ; j < 4 ; ++j )
				geneOverlap[j] = sortedReads[i].geneOverlap[j] ;
			
			// If the order of V,D,J,C is wrong from this read, then we ignore this.
			//   probably from read through in cyclic fragment.
			// We do not need to worry about the strand here, since if the strand is -1, the coordinate
			//   is relative to the reversed-complementary read.
			bool filter = false ;
			int strand = 0 ;
			/*for ( j = 0 ; j < 4 ; ++j )
			{
				if ( geneOverlap[j].seqIdx == -1 )
					continue ;
				printf( "%d: %s %d %d; %d %d; %d %lf\n", j, refSet.GetSeqName( geneOverlap[j].seqIdx ), 
					geneOverlap[j].readStart, geneOverlap[j].readEnd,
					geneOverlap[j].seqStart, geneOverlap[j].seqEnd,
					geneOverlap[j].strand, geneOverlap[j].similarity ) ;
			}*/
			for ( j = 0 ; j < 4 ; ++j )
			{
				int l ;
				if ( geneOverlap[j].seqIdx == -1 )
					continue ;
				for ( l = j + 1 ; l < 4 ; ++l )
				{
					if ( geneOverlap[l].seqIdx == -1 )
						continue ;
					
					if ( geneOverlap[j].readEnd - 10 > geneOverlap[l].readStart ) 
					{
						filter = true ;	
						break ;
					}
				}
				if ( filter )
					break ;
			}

			if ( geneOverlap[3].seqIdx != -1 && geneOverlap[0].seqIdx == -1 && geneOverlap[2].seqIdx == -1 ) // From constant gene.
			{
				//if ( geneOverlap[3].readEnd - geneOverlap[3].readStart + 1 < sortedReads[i].len )
				//{
				if ( geneOverlap[3].seqStart >= constantGeneEnd )
					filter = true ;
				else if ( constantGeneEnd <= 200 && geneOverlap[3].seqStart >= 100 && ( geneOverlap[3].strand == 1 
					|| geneOverlap[3].readEnd - geneOverlap[3].readStart + 1 < sortedReads[i].len ) ) // This filter is for the case where short constant gene sequence is good enough. When constatnGeneEnd > 200, we assume user is not only interested in CDR3, so we filter fewer C gene reads.
					filter = true ;
				//}
				
			}
			/*for ( j = 0 ; j < 4 ; ++j )
				printf( "%d ", geneOverlap[j].seqIdx ) ;
			printf( "\n" ) ;*/
			if ( filter ) 
				addRet = -1 ;
			else
			{
				char name[5] ;
				name[0] = '\0' ;
				strand = 0 ;
				int ambiguousStrand = 0 ;
				for ( j = 0 ; j < 4 ; ++j )
					if ( geneOverlap[j].seqIdx != -1 )
					{
						char *s = refSet.GetSeqName( geneOverlap[j].seqIdx ) ;
						name[0] = s[0] ; name[1] = s[1] ; name[2] = s[2] ; name[3] = s[3] ;
						name[4] = '\0' ;
						if ( strand != 0 && strand != geneOverlap[j].strand )
							ambiguousStrand = 1 ;
						strand = geneOverlap[j].strand ;
					}
				if ( ambiguousStrand )
					strand = 0 ;

				double similarityThreshold = 0.9 ;
				if ( sortedReads[i].minCnt >= 20 )
					similarityThreshold = 0.97 ;
				else if ( sortedReads[i].minCnt >= 2
					|| ( sortedReads[i].minCnt >= 5 && firstReadLen > 200 ) )
				{
					//double tmp = 1.0 - 4.0 / sortedReads[i].len ;
					//similarityThreshold = tmp < 0.95 ? tmp : 0.95 ;
					//if ( similarityThreshold < 0.9 )
					//	similarityThreshold = 0.9 ;
					similarityThreshold = 0.95 ;
				}

				if ( name[0] == 'T' && similarityThreshold < 0.95 ) // TCR
					similarityThreshold = 0.95 ;
				
				// When using barcode, we could be more casual.
				if ( hasBarcode || trimLevel > 1 )
					similarityThreshold = 0.9 ;

				//if ( similarityThreshold > 0.9 && geneOverlap[0].seqIdx != -1 && geneOverlap[2].seqIdx != -1 
				//	&& geneOverlap[0].readEnd < geneOverlap[2].readStart )
				//	similarityThreshold = 1.0 ;

				addRet = seqSet.AddRead( sortedReads[i].read, name, strand, sortedReads[i].barcode, 
						hasBarcode ? (sortedReads[i].minCnt + sortedReads[i].barcodeMinCnt + 1) / 2 : sortedReads[i].minCnt, trimLevel > 1, similarityThreshold ) ;
				
				if ( addRet < 0 )
				{
					// Check the overall hit length.
					int matchCnt = 0 ;
					for ( j = 0 ; j < 4 ; ++j )
					{
						if ( geneOverlap[j].seqIdx != -1 )
							matchCnt += geneOverlap[j].matchCnt / 2 ;
					}
					filter = true ;
					if ( matchCnt >= 31 )
						filter = false ;
					else
					{
						// The anchor is very close to the end of V or J gene
						if ( geneOverlap[0].seqIdx != -1 && geneOverlap[2].seqIdx != -1 
							&& geneOverlap[0].readEnd < geneOverlap[2].readStart )
						{
							filter = false ;
						}
						else if ( geneOverlap[0].seqIdx != -1 )
						{
							if ( geneOverlap[0].seqEnd >= 
								refSet.GetSeqConsensusLen( geneOverlap[0].seqIdx ) - 17 )
								filter = false ;
						}
						else if ( geneOverlap[2].seqIdx != -1 )
						{
							if ( geneOverlap[2].seqStart <= 17 )
								filter = false ;
						}
					}
					for ( j = 0 ; j < 4 ; ++j )
						if ( geneOverlap[j].seqIdx != -1 )
						{
							break ;
						}
					
					if ( !filter )
					{
						addRet = seqSet.InputNovelRead( refSet.GetSeqName( geneOverlap[j].seqIdx ), 
							sortedReads[i].read, geneOverlap[j].strand, sortedReads[i].barcode ) ;

					}
					else if ( goodCandidate[i] )
					{
						// The mate is matched well, we use motif to check whether this could be
						//   from CDR3.
						// Here info is set during good candidate set.
						if ( seqSet.HasMotif( sortedReads[i].read, -sortedReads[ sortedReads[i].info ].strand ) )
							addRet = seqSet.InputNovelRead( "Novel", sortedReads[i].read, 
									-sortedReads[ sortedReads[i].info ].strand,
									sortedReads[i].barcode ) ;
					}
					//printf( "hello %d %d. %d %d %d %d\n", i, addRet, geneOverlap[0].seqIdx,
					//		geneOverlap[1].seqIdx, geneOverlap[2].seqIdx, geneOverlap[3].seqIdx ) ;
				}
				sortedReads[i].strand = strand ;
			}
		}
		else
		{
			//printf( "saved time\n" ) ;
			if ( prevAddRet != -1 && prevAddRet != -3 )
				addRet = seqSet.RepeatAddRead( sortedReads[i].read ) ;
			else if ( prevAddRet == -3 )
				addRet = -3 ;

			sortedReads[i].strand = sortedReads[i - 1].strand ;
		}
		
		if ( addRet == -2 )
			rescueReadIdx.push_back( i ) ;
		else if ( addRet >= 0 )
		{
			++assembledReadCnt ;
			assembledReadIdx.push_back( i ) ;

			if ( sortedReads[i].mateIdx > i ) // This handles the case mateIdx = -1
			{
				bool good = false ;
				bool maySpan = false ;
				if ( geneOverlap[0].seqIdx != -1 && geneOverlap[0].similarity >= 0.9 
						&& sortedReads[i].strand == 1 )
				{
					good = true ;
					if ( geneOverlap[2].seqIdx != -1 && geneOverlap[2].readStart > geneOverlap[0].readEnd )
						maySpan = true ;
					if ( geneOverlap[3].seqIdx != -1 && geneOverlap[3].readStart > geneOverlap[0].readEnd )
						maySpan = true ;
				}

				for ( j = 2 ; j <= 3 ; ++j )
				{
					if ( geneOverlap[j].seqIdx != -1 && geneOverlap[j].similarity >= 0.9 
							&& sortedReads[i].strand == -1 )
					{
						good = true ;
						if ( geneOverlap[0].seqIdx != -1 && geneOverlap[j].readStart > geneOverlap[0].readEnd )
							maySpan = true ;
					}

				}

				if ( maySpan ) // If the read can span V,J, then there is no need to force add its mate.
					good = false ;

				if (good && !goodCandidate[ sortedReads[i].mateIdx ]) 
				{
					// We allow it go beyond the restriction of barcode
					// so it can borrow information across cells.
					int tag = sortedReads[i].mateIdx ;
					for ( j = tag - 1 ; j > 0 ; --j )
					{
						if ( !strcmp(sortedReads[j].read, sortedReads[tag].read) )
						{
							goodCandidate[j] = good ;
							sortedReads[j].info = i;
						}
						else 
							break ;
					}

					for ( j = tag + 1 ; j < readCnt ; ++j )
					{
						if ( !strcmp(sortedReads[j].read, sortedReads[tag].read) )
						{
							goodCandidate[j] = good ;
							sortedReads[j].info = i;
						}
						else 
							break ;
					}
				}

				if (good)
				{
					goodCandidate[ sortedReads[i].mateIdx ] = good ;
					sortedReads[ sortedReads[i].mateIdx ].info = i;
				}
			} // end if for setting good candidate
			
			// Check whether the many barcode is already finished, so we can purge them from the index and posWeight array
			if (hasBarcode && !keepMissingBarcode && sortedReads[i].barcode != -1)
			{
				int barcode = sortedReads[i].barcode ;
				++barcodeReadCount[barcode] ;
				if (barcodeReadCount[barcode] >= barcodeTotalReadCount[barcode])
				{
					finishedBarcodes[barcode] = barcodeTotalReadCount[barcode] ;
					if (finishedBarcodes.size() >= 1) // clean up the memory once every barcode finishes. The sort makes sure the barocde is processed one by one, so it can be efficiently removed.
					{
						seqSet.ReleaseFinishedBarcodeSeq(finishedBarcodes, /*release_index*/true, contigMinCov, /*early_stop*/true) ;
						finishedBarcodes.clear() ;
					}
				}
			}
		}

		if ( assembledReadCnt > 0 && assembledReadCnt % 10000 == 0 
				&& !hasBarcode)
			seqSet.UpdateAllConsensus() ;

		if ( ( i + 1 ) % 100000 == 0 )
			PrintLog( "Processed %d reads (%d are used for assembly).", i + 1, assembledReadCnt ) ;
		
		prevAddRet = addRet ;
#ifdef DEBUG
		printf( "done\n" ) ;
#endif

		if ( seqSet.Size() > changeKmerLengthThreshold && indexKmerLength < 16 && !hasBarcode )
		{
			changeKmerLengthThreshold *= 4 ;
			indexKmerLength += 2 ;
			seqSet.ChangeKmerLength( indexKmerLength ) ;
		}
	}
	seqSet.UpdateAllConsensus() ;
	//printf("index=%lu reads=%lu\n", seqSet.GetSpace(), sortedReads.size() * sizeof(sortedReads[0])) ;
	PrintLog(  "Assembled %d reads.", assembledReadCnt ) ;

	// Go through the second round.
	// TODO: user-defined number of rounds.
	/*seqSet.ResetPosWeight() ;
	reads.Rewind() ;
	while ( reads.Next() )
	{
		//printf( "%s %s\n", reads.id, reads.seq ) ;
		//fflush( stdout ) ;
		seqSet.AddRead( reads.seq ) ;
		//printf( "done\n" ) ;
	}*/

	int rescueReadCnt = rescueReadIdx.size() ;
	if ( firstReadLen > 200 )
		rescueReadCnt = 0 ;
	PrintLog( "Try to rescue %d reads for assembly.", rescueReadCnt) ;
	assembledReadCnt = 0 ;
	for ( i = 0 ; i < rescueReadCnt ; ++i )
	{
#ifdef DEBUG
		printf( "%s %s %d %lf\n", sortedReads[ rescueReadIdx[i] ].id, 
			sortedReads[ rescueReadIdx[i] ].read, sortedReads[  rescueReadIdx[i]  ].medianCnt, 
			sortedReads[ rescueReadIdx[i] ].avgCnt ) ;
		fflush( stdout ) ;
#endif
		int addRet = -1 ;
		char name[2] = "" ;
		
		double similarityThreshold = 0.9 ;
		if ( sortedReads[ rescueReadIdx[i] ].minCnt >= 20 )
			similarityThreshold = 0.97 ;
		else if ( sortedReads[ rescueReadIdx[i] ].minCnt >= 2 )
		{
			//double tmp = 1.0 - 4.0 / sortedReads[i].len ;
			//similarityThreshold = tmp < 0.95 ? tmp : 0.95 ;
			//if ( similarityThreshold < 0.9 )
			//	similarityThreshold = 0.9 ;
			similarityThreshold = 0.95 ;
		}
		
		int strand = 0 ;
		addRet = seqSet.AddRead( sortedReads[ rescueReadIdx[i] ].read, name, strand, sortedReads[ rescueReadIdx[i] ].barcode, 
			1, trimLevel > 1 ? true : false, similarityThreshold ) ;
		sortedReads[ rescueReadIdx[i] ].strand = strand ;

		if ( addRet >= 0 )
		{
			++assembledReadCnt ;
			assembledReadIdx.push_back( rescueReadIdx[i] ) ;
		}
#ifdef DEBUG
		printf( "rescue done\n" ) ;
#endif
	}
	seqSet.UpdateAllConsensus() ;
	PrintLog( "Rescued %d reads.", assembledReadCnt ) ;
	//exit( 1 ) ;
	/*seqSet.Clean( true ) ;
	fprintf( stderr, "Found %d raw contigs.\nStart to merge raw contigs.\n", seqSet.GetSeqCnt() ) ;
	
	i = seqSet.Assemble() ;
	seqSet.UpdateAllConsensus() ;
	fprintf( stderr, "Reduce %d raw contigs.\n", i ) ;
	
	fprintf( stderr, "Annotate the contigs.\n" ) ;
	seqSet.Annotate( refSet ) ;*/

	if (contigMinCov > 0)
	{
		seqSet.ReleaseShallowContigs(contigMinCov) ;
	}

	// Output the preliminary assembly.
	//seqSet.Clean( true ) ;
	FILE *fp ;
	if ( outputPrefix[0] != '-' )
	{
		sprintf( buffer, "%s_raw.out", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
	}
	else
		fp = stdout ;
	
	if ( hasBarcode )
		seqSet.Output( fp, &barcodeIntToStr ) ;
	else
		seqSet.Output( fp, NULL ) ;
	fflush( fp ) ;
	
	if ( outputPrefix[0] != '-' )
		fclose( fp ) ;

	//sprintf( buffer, "%s_assembled_reads.fa", outputPrefix ) ;

	std::vector<struct _assignRead> assembledReads ;
	assembledReadCnt = assembledReadIdx.size() ;
	for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		struct _assignRead nr ;
		nr.id = sortedReads[ assembledReadIdx[i] ].id ;
		nr.read = sortedReads[ assembledReadIdx[i] ].read ;
		nr.barcode = sortedReads[ assembledReadIdx[i] ].barcode ;
		nr.umi = sortedReads[ assembledReadIdx[i] ].umi ;
		nr.info = assembledReadIdx[i] ;
		nr.overlap.seqIdx = -1 ;
		nr.overlap.strand = sortedReads[ assembledReadIdx[i] ].strand ;
		assembledReads.push_back( nr ) ;
	}
	//std::sort( assembledReads.begin(), assembledReads.end(), CompSortReadById ) ;	
	
	sprintf( buffer, "%s_assembled_reads.fa", outputPrefix ) ;
	fp = fopen( buffer, "w" ) ;
	for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		buffer[0] = '\0' ;
		if ( hasBarcode )
		{
			char *p = buffer + strlen( buffer ) ;
			sprintf( p, " barcode:%s", barcodeIntToStr[ assembledReads[i].barcode ].c_str() ) ;
		}
		if ( hasUmi )
		{
			char *p = buffer + strlen( buffer ) ;
			sprintf( p, " umi:%d", assembledReads[i].umi ) ;
		}

		fprintf( fp, ">%s %d %d %d%s\n%s\n", assembledReads[i].id, assembledReads[i].overlap.strand,
			sortedReads[ assembledReadIdx[i] ].minCnt, 
			sortedReads[ assembledReadIdx[i] ].medianCnt, buffer,
			assembledReads[i].read ) ;
	}
	fclose( fp ) ;

	if ( skipMateExtension || !hasMate || hasBarcode)
	{
		FILE *fp ;
		if ( outputPrefix[0] != '-' )
		{
			sprintf( buffer, "%s_final.out", outputPrefix ) ;
			fp = fopen( buffer, "w" ) ;
		}
		else
			fp = stdout ;

		if ( hasBarcode )
			seqSet.Output( fp, &barcodeIntToStr ) ;
		else
			seqSet.Output( fp, NULL ) ;
		fflush( fp ) ;

		if ( outputPrefix[0] != '-' )
			fclose( fp ) ;

		for ( i = 0 ; i < readCnt ; ++i )
		{
			free( sortedReads[i].id ) ;
			free( sortedReads[i].read ) ;
			if ( sortedReads[i].qual != NULL )
				free( sortedReads[i].qual ) ;
		}
		return 0 ;	
	}

	SeqSet extendedSeq( indexKmerLength > 17 ? indexKmerLength : 17 ) ;
	extendedSeq.InputSeqSet( seqSet, false ) ;
	struct _overlap assign ;
	//memset( &assign, -1, sizeof( assign ) ) ;

	/*fp = fopen( "assign_results.out", "r" ) ;
	for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		fread( &assembledReads[i].overlap, sizeof( assembledReads[i].overlap ), 1, fp ) ;
		//printf( "%d %d %d\n", assembledReads[i].overlap.seqIdx, assembledReads[i].overlap.readStart, 
		//	assembledReads[i].overlap.readEnd ) ;
	}*/
	//fp = fopen( "assign_results.out", "w" ) ;

	// Add the distant constant genes from the annotation
	/*int refSetSize = refSet.GetSeqCnt() ;
	int constantGeneAddCnt = 0 ;
	for ( i = 0 ; i < refSetSize ; ++i )
	{
		char *name = GetSeqName( i ) ;
		if ( refSet.GetGeneType( name ) != 3 || refSet.GetSeqConsensusLen( i ) <= constantGeneEnd )
			continue ;
		extendedSeq.InputRefSeq( name, refSet.GetSeqConsensus( i ) + constantGeneEnd ) ;
	}*/
	if ( firstReadLen > 200 )
		extendedSeq.SetIsLongSeqSet( true ) ;
	extendedSeq.SetNovelSeqSimilarity( 0.95 ) ;
	if ( threadCnt <= 1 )
	{
		for ( i = 0 ; i < assembledReadCnt ; ++i )
		{
			if ( i == 0 || strcmp( assembledReads[i].read, assembledReads[i - 1].read ) )
				extendedSeq.AssignRead( assembledReads[i].read, assembledReads[i].overlap.strand, 
						assembledReads[i].barcode, assign ) ;
			assembledReads[i].overlap = assign ;	
			if ( ( i + 1 ) % 100000 == 0 )
			{
				PrintLog( "Processed %d reads for extension.", i + 1 ) ;
			}
			//fprintf( fp, "%s %s %d\n", assembledReads[i].id, assembledReads[i].read, assign.seqIdx ) ;
			//fwrite( &assign, sizeof( assign ), 1, fp ) ;
		}
	}
	else
	{
		pthread_t *threads = new pthread_t[ threadCnt ] ;
		struct _assignReadsThreadArg *args = new struct _assignReadsThreadArg[threadCnt] ;
		pthread_attr_t attr ;
		
		pthread_attr_init( &attr ) ;
		pthread_attr_setdetachstate( &attr, PTHREAD_CREATE_JOINABLE ) ;
		
		for ( i = 0 ; i < threadCnt ; ++i )
		{
			args[i].tid = i ;
			args[i].threadCnt = threadCnt ;
			args[i].seqSet = &extendedSeq ;
			args[i].pAssembledReads = &assembledReads ;
			args[i].assembledReadCnt = assembledReadCnt ;
			pthread_create( &threads[i], &attr, AssignReads_Thread, (void *)( args + i ) ) ;
		}

		for ( i = 0 ; i < threadCnt ; ++i )
			pthread_join( threads[i], NULL ) ;
	

		delete[] threads ;
		delete[] args ;
	}
	extendedSeq.SetNovelSeqSimilarity( 0.9 ) ;
	extendedSeq.RecomputePosWeight( assembledReads ) ;
	//fclose( fp ) ;
	//exit( 1 ) ;

#ifdef DEBUG
	for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		printf( "%s %d %lf %d\n", assembledReads[i].id, assembledReads[i].overlap.seqIdx, assembledReads[i].overlap.similarity,
			assembledReads[i].overlap.strand ) ;
	}
#endif

	/*extendedSeq.BreakFalseAssembly( assembledReads ) ;
	
	if ( outputPrefix[0] != '-' )
	{
		sprintf( buffer, "%s_contig.out", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
	}
	else
		fp = stdout ;

	extendedSeq.Output( fp ) ;
	fflush( fp ) ;
	
	if ( outputPrefix[0] != '-' )
		fclose( fp ) ;*/

	int avgReadLen = 1 ;
	for ( i = 0 ; i < assembledReadCnt && i < 1000 ; ++i )
		avgReadLen += strlen( assembledReads[i].read ) ;
	if ( i > 0 )
		avgReadLen /= i ;
	
	PrintLog( "Extend assemblies by mate pair information." ) ;
	extendedSeq.ExtendSeqFromReads( assembledReads, 17, refSet ) ; //( avgReadLen / 30 < 31 ) ? 31 : ( avgReadLen / 3 )  ) ;
	extendedSeq.UpdateAllConsensus() ;
		
	/*if ( outputPrefix[0] != '-' )
	{
		sprintf( buffer, "%s_extended.out", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
	}
	else
		fp = stdout ;

	extendedSeq.Output( fp ) ;
	fflush( fp ) ;
	
	if ( outputPrefix[0] != '-' )
		fclose( fp ) ;*/
	
	// Now the assembled reads should be sorted by their read id.
	//PrintLog( "Rescue low frequency reads." ) ;
	SeqSet lowFreqSeqSet( indexKmerLength ) ;
	KmerCount lowFreqKmerCount( 31 ) ;
	std::vector<int> lowFreqReadsIdx ;
	for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		continue ;
		if ( assembledReads[i].overlap.seqIdx != -1 ) 
		//if ( assembledReads[i].overlap.seqIdx != -1 && assembledReads[i + 1].overlap.seqIdx != -1 
		//	&& assembledReads[i].overlap.seqIdx == assembledReads[i + 1].overlap.seqIdx )
		{
			continue ;
		}

		if ( sortedReads[ assembledReads[i].info ].minCnt < 4 )
			continue ;
		
		int nCnt = 0 ;
		for ( j = 0 ; assembledReads[i].read[j] ; ++j )
			if ( assembledReads[i].read[j] == 'N' )
				++nCnt ;
		if ( nCnt >= 1 )
		{
			continue ;
		}

		lowFreqKmerCount.AddCount( assembledReads[i].read ) ;
		lowFreqReadsIdx.push_back(i) ;
	}
	k = 0 ;
	int lowFreqCnt = lowFreqReadsIdx.size() ;
	std::vector<struct _sortRead> lowFreqReads ;
	for ( i = 0 ; i < lowFreqCnt ; ++i )
	{
		// At least one of the mate should have at least 2 kmer count.
		int a = lowFreqReadsIdx[i] ;
		int minCntA, tmp ;
		float tmpf ;
		lowFreqKmerCount.GetCountStatsAndTrim( assembledReads[a].read, NULL, 
			minCntA, tmp, tmpf ) ;

		// Then these two pairs must overlap with each other.
		struct _sortRead nr ;
		nr.id = NULL ;
		nr.read = assembledReads[a].read ;
		nr.qual = NULL ;
		nr.minCnt = minCntA ; 
		nr.medianCnt = tmp ;
		nr.avgCnt = tmpf ;
		nr.strand = assembledReads[a].overlap.strand ;
		lowFreqReads.push_back( nr ) ;
	}
	std::sort( lowFreqReads.begin(), lowFreqReads.end() ) ;
	/*for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		char *r = strdup( assembledReads[i].read ) ;
		lowFreqMergedReads.push_back( r ) ;
		++k ;
		
	}*/
	lowFreqCnt = lowFreqReads.size() ;
	//PrintLog( "Found %d low frequency reads.", lowFreqCnt ) ;
	for ( i = 0 ; i < lowFreqCnt ; ++i )
	{
		char name[10] = "" ;
		name[0] = '\0' ;
		//printf( "%s\n", lowFreqReads[i].read ) ;
		//fflush( stdout ) ;
		//printf( ">r%d\n%s\n", i, lowFreqReads[i].read ) ;
		extendedSeq.AssignRead( lowFreqReads[i].read, lowFreqReads[i].strand, lowFreqReads[i].barcode, assign ) ;
		if ( assign.seqIdx != -1 )
		{
			if ( assign.similarity < 1.0 )
				extendedSeq.AddAssignedRead( lowFreqReads[i].read, assign ) ;
		}
		//else if ( extendedSeq.AddRead( lowFreqReads[i].read, name, 1.0 ) == -1 )
		else
		{
			int strand ;
			if ( lowFreqReads[i].minCnt < 2 )
			{
				lowFreqSeqSet.AssignRead( lowFreqReads[i].read, lowFreqReads[i].strand, 
						lowFreqReads[i].barcode, assign ) ;
				if ( assign.seqIdx != -1 )
				{
					lowFreqSeqSet.AddAssignedRead( lowFreqReads[i].read, assign ) ;
				}
			}
			else if ( lowFreqSeqSet.AddRead( lowFreqReads[i].read, name, strand, lowFreqReads[i].barcode, 1, false, 0.97 ) < 0 )
			{
				struct _overlap geneOverlap[4] ;
				//buffer[0] = '\0' ;
				refSet.AnnotateRead( lowFreqReads[i].read, 0, geneOverlap, NULL, NULL ) ;

				int componentCnt = 0 ;
				int lastJ = -1 ;
				for ( j = 0 ; j < 4 ; ++j )
				{
					if ( geneOverlap[j].seqIdx != -1 )
					{
						++componentCnt ;
						lastJ = j ;
					}
				}
				if ( componentCnt >= 1 && lowFreqReads[i].minCnt >= 2 )
				{
					j = lastJ ;
					lowFreqSeqSet.InputNovelRead( refSet.GetSeqName( geneOverlap[j].seqIdx ), 
							lowFreqReads[i].read, geneOverlap[j].strand, lowFreqReads[i].barcode ) ;
				}

			}
		}
	}
	
	extendedSeq.UpdateAllConsensus() ;
	lowFreqSeqSet.UpdateAllConsensus() ;
	extendedSeq.InputSeqSet( lowFreqSeqSet, false ) ;	
	
	
	PrintLog( "Remove redundant assemblies." ) ;
	extendedSeq.ChangeKmerLength( 31 ) ;
	//i = extendedSeq.ExtendSeqFromSeqOverlap( ( avgReadLen / 2 < 31 ) ? 31 : ( avgReadLen / 2 ) ) ;
	i = extendedSeq.RemoveRedundantSeq() ;
	
	
	if ( outputPrefix[0] != '-' )
	{
		sprintf( buffer, "%s_final.out", outputPrefix ) ;
		fp = fopen( buffer, "w" ) ;
	}
	else
		fp = stdout ;
	
	if ( hasBarcode )
		extendedSeq.Output( fp, &barcodeIntToStr ) ;
	else
		extendedSeq.Output( fp, NULL ) ;
	fflush( fp ) ;
	
	if ( outputPrefix[0] != '-' )
		fclose( fp ) ;
	
	/*for ( i = 0 ; i < assembledReadCnt ; ++i )
	{
		fprintf( fp, ">%s\n%s\n", assembledReads[i].id, assembledReads[i].read ) ;
	}
	fclose( fp ) ;*/
	if ( readCnt == 0 )
		readCnt = assembledReadCnt ;	
	for ( i = 0 ; i < readCnt ; ++i )
	{
		free( sortedReads[i].id ) ;
		free( sortedReads[i].read ) ;
		if ( sortedReads[i].qual != NULL )
			free( sortedReads[i].qual ) ;
	}

	PrintLog( "Finish assembly." ) ;
	return 0 ;
}