mrfast: mrfast-2.1.0.4/MrFAST.c comparison

comparison mrfast-2.1.0.4/MrFAST.c @ 0:7b3dc85dc7fd

Uploaded mrfast source tarball

author	calkan
date	Tue, 21 Feb 2012 10:29:47 -0500
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:7b3dc85dc7fd
+/*
+* Copyright (c) <2008 - 2012>, University of Washington, Simon Fraser University
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without modification,
+* are permitted provided that the following conditions are met:
+*
+* Redistributions of source code must retain the above copyright notice, this list
+* of conditions and the following disclaimer.
+* - Redistributions in binary form must reproduce the above copyright notice, this
+*   list of conditions and the following disclaimer in the documentation and/or other
+*   materials provided with the distribution.
+* - Neither the names of the University of Washington, Simon Fraser University,
+*   nor the names of its contributors may be
+*   used to endorse or promote products derived from this software without specific
+*   prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+/*
+Authors:
+	Farhad Hormozdiari
+	Faraz Hach
+	Can Alkan
+Emails:
+	farhadh AT uw DOT edu
+	fhach AT cs DOT sfu DOT ca
+calkan AT uw DOT edu
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <dirent.h>
+#include <xmmintrin.h>
+#include <emmintrin.h>
+#include <mmintrin.h>
+#include "Common.h"
+#include "Reads.h"
+#include "HashTable.h"
+#include "Output.h"
+#include "MrFAST.h"
+#include "RefGenome.h"
+#define min(a,b) ((a)>(b)?(b):(a))
+#define min3(a,b,c) ((a)>(b)?(b>c?c:b):(a>c?c:a))
+#define CHARCODE(a) (a=='A' ? 0 : (a=='C' ? 1 : (a=='G' ? 2 : (a=='T' ? 3 : 4))))
+#define MAX_REF_SIZE	18
+float calculateScore(int index, char *seq, char *qual, char *md);
+unsigned char		mrFAST = 1;
+char				*versionNumberF="0.4";
+long long			verificationCnt = 0;
+long long 			mappingCnt = 0;
+long long			mappedSeqCnt = 0;
+long long			completedSeqCnt = 0;
+char				*mappingOutput;
+/**********************************************/
+char				*_msf_refGen = NULL;
+int				_msf_refGenLength = 0;
+int				_msf_refGenOffset = 0;
+char				*_msf_refGenName = NULL;
+int				_msf_refGenBeg;
+int				_msf_refGenEnd;
+IHashTable			*_msf_hashTable = NULL;
+int				*_msf_samplingLocs;
+int				*_msf_samplingLocsEnds;
+int				_msf_samplingLocsSize;
+Read				*_msf_seqList;
+int				_msf_seqListSize;
+Pair 				*_msf_sort_seqList = NULL;
+int 				*_msf_map_sort_seqList;
+ReadIndexTable			*_msf_rIndex = NULL;
+int				_msf_rIndexSize;
+int				_msf_rIndexMax;
+SAM				_msf_output;
+OPT_FIELDS			*_msf_optionalFields;
+char				*_msf_op;
+int				*_msf_verifiedLocs = NULL;
+char				_msf_numbers[200][3];
+char				_msf_cigar[5];
+MappingInfo			*_msf_mappingInfo;
+int				*_msf_seqHits;
+int				_msf_openFiles = 0;
+int				_msf_maxLSize=0;
+int				_msf_maxRSize=0;
+BestFullMappingInfo			*bestHitMappingInfo;
+/*************************/
+int				_msf_maxFile=0;
+char				_msf_fileName[4000][200][2][FILE_NAME_LENGTH];
+int				_msf_fileCount[4000];
+char				*_msf_readHasConcordantMapping; //boolean if a read has concordant mapping :D
+int 				*_msf_oeaMapping;
+int 				*_msf_discordantMapping;
+FILE *bestConcordantFILE;
+FILE *bestDiscordantFILE;
+int counter = 0;
+int scoreF[200][200];
+int scoreB[200][200];
+int score[200][200];
+int direction1[200][200];
+int direction2[200][200];
+__m128i MASK;
+int 			lookUpTable[15625][15625];
+/**************************************************Methods***************************************************/
+int smallEditDistanceF(char *a, int lena, char *b, int lenb)
+{
+int matrix[20][20];
+int i = 0;
+int j = 0;
+for(i = 0; i <= lena; i++)
+{
+matrix[0][i] = i;
+}
+for(i = 0; i <= lenb; i++)
+{
+matrix[i][0] = i;
+}
+for(i = 1; i <= lenb; i++)
+{
+for(j = 1; j <= lena; j++)
+{
+matrix[i][j] = min3(matrix[i-1][j-1]+ (a[j-1] != b[i-1]),matrix[i][j-1]+1 ,matrix[i-1][j]+1);
+}
+}
+return (matrix[lenb][lena]>errThreshold?-1:matrix[lenb][lena]);
+}
+int smallEditDistanceB(char *a, int lena, char *b, int lenb)
+{
+int matrix[20][20];
+int i = 0;
+int j = 0;
+for(i = 0; i <= lena; i++)
+{
+matrix[0][i] = i;
+}
+for(i = 0; i <= lenb; i++)
+{
+matrix[i][0] = i;
+}
+for(i = 1; i <= lenb; i++)
+{
+for(j = 1; j <= lena; j++)
+{
+matrix[i][j] = min3(matrix[i-1][j-1]+ (*(a-j+1) != *(b-i+1)),matrix[i][j-1]+1 ,matrix[i-1][j]+1);
+}
+}
+return (matrix[lenb][lena]>errThreshold?-1:matrix[lenb][lena]);
+}
+char fastEditDistance(int per1, int per2)
+{
+	int i = 0;
+	int j = 0;
+	char str1[7];
+	char str2[7];
+	int val1 = per1;
+	int val2 = per2;
+	int index = 0;
+	int mod = 0;
+	int matrix[7][7];
+	int min = 20;
+	while(index < 6)
+	{
+		mod = val1%5;
+		str1[5-index] = (mod==0 ? 'A':(mod==1?'C':mod==2?'G':(mod==3)?'T':'N'));
+		val1 = val1 /5;
+		index++;
+	}
+	str1[6] = '\0';
+	index = 0;
+	while(index < 6)
+	{
+		mod=val2%5;
+		str2[5-index] = (mod==0 ? 'A':(mod==1?'C':mod==2?'G':(mod==3)?'T':'N'));
+		val2 = val2 / 5;
+		index++;
+	}
+	str2[6] = '\0';
+	for(i = 0; i < 7; i++)
+	{
+		matrix[0][i] = i;
+		matrix[i][0] = i;
+	}
+	for(i = 1; i < 7; i++)
+	{
+		for(j = 1; j < 7; j++)
+		{
+			matrix[i][j] = min3(matrix[i-1][j-1]+ (str1[i-1] != str2[j-1]),matrix[i][j-1]+1 ,matrix[i-1][j]+1);
+		}
+	}
+	for(i = 0; i < 7; i++)
+	{
+		if(matrix[i][6] < min)
+			min = matrix[i][6];
+	}
+	for(i = 0; i < 7; i++)
+	{
+		if(matrix[6][i] < min)
+			min = matrix[6][i];
+	}
+	return min;
+}
+void initLookUpTable()
+{
+	int i = 0;
+	MASK = _mm_insert_epi16(MASK,1,0);
+	MASK = _mm_insert_epi16(MASK,1,1);
+	MASK = _mm_insert_epi16(MASK,1,2);
+	MASK = _mm_insert_epi16(MASK,1,3);
+	MASK = _mm_insert_epi16(MASK,1,4);
+	MASK = _mm_insert_epi16(MASK,0,5);
+	MASK = _mm_insert_epi16(MASK,0,6);
+	MASK = _mm_insert_epi16(MASK,0,7);
+	for(i = 0 ; i < errThreshold + 1; i++)
+	{
+		scoreF[0][i] = i;
+		scoreF[i][0] = i;
+	}
+	for(i = 0 ; i < errThreshold + 1; i++)
+	{
+		scoreB[0][i] = i;
+		scoreB[i][0] = i;
+	}
+}
+int backwardEditDistanceSSE2Odd(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int e = errThreshold;
+	char flag = 0;
+	int minError = 2*e;
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i Error;
+	__m128i tmp;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Error = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	if(lenb <= e)
+	{
+		return smallEditDistanceB(a,lena,b,lenb);
+	}
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,2*e,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,1,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,2);
+	R0 = _mm_insert_epi16(R0,0,0);
+	R1 = _mm_insert_epi16(R1,1,0);
+	R1 = _mm_insert_epi16(R1,1,1);
+	for(i=2; i <= e; i++)
+	{
+		//set side
+		Side1 = _mm_slli_si128(Side1,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_slli_si128(Down1,2);
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=1;j<=i-1;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(i/2-1+(i/2-j))) != *(a-(i/2-1-(i/2-j))),0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R0 = _mm_min_epi16(R1+Side1, _mm_slli_si128(R0,2)+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down1);
+		}
+		else
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=i/2-1;j>=-i/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-((i+1)/2+j-1)) != *(a-((i-1)/2-j-1)),0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down1);
+		}
+	}
+	Error = _mm_xor_si128(Error, Error);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Error = _mm_insert_epi16(Error,e,0);
+	Side1 = _mm_insert_epi16(Side2,2*e,0);
+	Side2 = _mm_insert_epi16(Side2,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Error = _mm_slli_si128(Error, 2);
+		Error = _mm_insert_epi16(Error, e, 0);
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,2*e,0);
+	for(; i <= 2*lenb-(e-1);i++)
+	{
+		flag = 0;
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			for(j=e/2;j>=-e/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(i/2-1+j)) != *(a-(i/2-1-j)),0);
+			}
+			R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+			R0 = _mm_min_epi16(R0, R1+Down1);
+			if(_mm_extract_epi16(R0,0) <= e)
+				flag = 1;
+			tmp = _mm_srli_si128(R0,2);
+			for(j=0; j <= e;j++)
+			{
+				if(_mm_extract_epi16(tmp,0) <= e)
+					flag = 1;
+				tmp = _mm_srli_si128(tmp,2);
+			}
+			if(flag == 0)
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			for(j=e/2;j>=-e/2-1;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-((i+1)/2+j-1)) != *(a-((i)/2-j-1)),0);
+			}
+			//	printf("@%d %d %d %d\n", _mm_extract_epi16(Diag,0), _mm_extract_epi16(Diag,1), _mm_extract_epi16(Diag,2),
+			//              _mm_extract_epi16(Diag,3));
+			R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+			//	printf("#~%d %d %d %d\n", _mm_extract_epi16(R1,0), _mm_extract_epi16(R1,1), _mm_extract_epi16(R1,2),
+			//              _mm_extract_epi16(R1,3));
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+			//	printf("$%d %d %d %d\n", _mm_extract_epi16(Side2,0), _mm_extract_epi16(Side2,1), _mm_extract_epi16(Side2,2),
+			//              _mm_extract_epi16(Side2,3));
+			//	printf("#%d %d %d %d\n", _mm_extract_epi16(R1,0), _mm_extract_epi16(R1,1), _mm_extract_epi16(R1,2),
+			//		_mm_extract_epi16(R1,3));
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	//first cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-3)) != *(a-lena), 0);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-2)) != *(a-(lena-1)), 1);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-1)) != *(a-(lena-2)), 2);
+	Diag = _mm_insert_epi16(Diag, 2*e, 3);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R1,2));
+	//second cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-2)) != *(a-(lena)), 0);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-1)) != *(a-(lena-1)), 1);
+	Diag = _mm_insert_epi16(Diag, 2*e, 2);
+	R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0, R1+Down1);
+	minError = min(minError, _mm_extract_epi16(R0,1));
+	//third cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-2)) != *(a-(lena+1)), 0);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-1)) != *(a-(lena)), 1);
+	Diag = _mm_insert_epi16(Diag, 2*e, 2);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	//forth
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-1)) != *(a-(lena+1)), 0);
+	Diag = _mm_insert_epi16(Diag, 2*e, 1);
+	R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0, R1+Down1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	//fifth
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, *(b-(lenb-1)) != *(a-(lena+2)), 0);
+	Diag = _mm_insert_epi16(Diag, 2*e, 1);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	if(minError > e)
+		return -1;
+	return minError;
+}
+int backwardEditDistanceSSE2G(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int e = errThreshold;
+	char flag = 0;
+	int minError = 2*e;
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i Error;
+	__m128i tmp;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Error = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	if(lenb <= e)
+	{
+		return smallEditDistanceB(a,lena,b,lenb);
+	}
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,2*e,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,1,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,2);
+	R0 = _mm_insert_epi16(R0,0,0);
+	R1 = _mm_insert_epi16(R1,1,0);
+	R1 = _mm_insert_epi16(R1,1,1);
+	for(i=2; i <= e; i++)
+	{
+		//set side
+		Side1 = _mm_slli_si128(Side1,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_slli_si128(Down1,2);
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=1;j<=i-1;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(i/2-1+(i/2-j))) != *(a-(i/2-1-(i/2-j))),0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R0 = _mm_min_epi16(R1+Side1, _mm_slli_si128(R0,2)+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down1);
+		}
+		else
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=i/2-1;j>=-i/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-((i+1)/2+j-1)) != *(a-((i-1)/2-j-1)),0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down1);
+		}
+	}
+	Error = _mm_xor_si128(Error, Error);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Error = _mm_insert_epi16(Error,e,0);
+	Side2 = _mm_insert_epi16(Side2,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Error = _mm_slli_si128(Error, 2);
+		Error = _mm_insert_epi16(Error, e, 0);
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,2*e,0);
+	for(; i <= 2*lenb-(e-1);i++)
+	{
+		flag = 0;
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			for(j=e/2;j>=-e/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(i/2-1+j)) != *(a-(i/2-1-j)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			if(_mm_extract_epi16(R0,0) <= e)
+				flag = 1;
+			tmp = _mm_srli_si128(R0,2);
+			for(j=0; j <= e;j++)
+			{
+				if(_mm_extract_epi16(tmp,0) <= e)
+					flag = 1;
+				tmp = _mm_srli_si128(tmp,2);
+			}
+			if(flag == 0)
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			for(j=-e/2+1;j<=e/2;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-((i+1)/2-j-1)) != *(a-((i-1)/2+j-1)),0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	j=0;
+	int tmpE = e;
+	for(;j<2*(e-2)+1;j++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		//set the first element
+		if(j==0)
+		{
+			for( k=0;k<=e-1;k++  )
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < e-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else if(j%2 == 0)
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < tmpE-1;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			tmp = _mm_srli_si128(R1,2);
+			for(k=0; k < tmpE-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		i++;
+	}
+	//Diag
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, 2*e, 0);
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-2)) != *(b-(lenb-1)), 1);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,1,1);
+	Down1 = _mm_insert_epi16(Down1, 2*e, 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 1);
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+	R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-1)) != *(b-(lenb-1)), 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 0);
+	R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	if(minError > e)
+		return -1;
+	return minError;
+}
+inline int backwardEditDistanceSSE2Extention(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int i0;
+	int i1;
+	int i2;
+	int i4;
+	int i5;
+	int e = 4;
+	int mismatch = errThreshold;
+	int minError = 2*errThreshold;
+	int index = 0;
+	int tmpValue = 0;
+	if(lenb <= e)
+	{
+		return smallEditDistanceB(a,lena,b,lenb);
+	}
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i tmp;
+	__m128i SeqA, SeqB;
+	__m128i Result;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	SeqA = _mm_setzero_si128 ();
+	SeqB = _mm_setzero_si128 ();
+	Result = _mm_setzero_si128 ();
+	/* end initialize */
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Diag = _mm_insert_epi16(Diag,minError,0);
+	i0 = (a[0] != b[0]);
+	i1 = min(i0, ( *(a-1)!=*b) )+1;
+	i2 = min(i0,( a[0] != *(b-1) ) )+1;
+	i0 = min3( i0+ ( *(a-1)!=*(b-1) ),i1+1,i2+1);
+	i4 = min(i1, ( *(a-2)!=b[0] )+1)+1;
+	i5 = min(i2, (a[0] != *(b-2))+1)+1;
+	R1 = _mm_insert_epi16(R1, 3, 0);
+	R1 = _mm_insert_epi16(R1, i1, 1);
+	R1 = _mm_insert_epi16(R1, i2, 2);
+	R1 = _mm_insert_epi16(R1, 3, 3);
+	R0 = _mm_insert_epi16(R0, 4, 0);
+	R0 = _mm_insert_epi16(R0, i4, 1);
+	R0 = _mm_insert_epi16(R0, i0, 2);
+	R0 = _mm_insert_epi16(R0, i5, 3);
+	R0 = _mm_insert_epi16(R0, 4, 4);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Side2 = _mm_insert_epi16(Side2,minError,0);
+	Down1 = _mm_insert_epi16(Down1,minError,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	index = 0;
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		SeqA = _mm_slli_si128(SeqA, 2);
+		SeqB = _mm_slli_si128(SeqB, 2);
+		SeqA = _mm_insert_epi16(SeqA,*(a-index),0);
+		SeqB = _mm_insert_epi16(SeqB,*(b-index),0);
+		index++;
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,minError,0);
+	index = 4;
+	i = 5;
+	int loopEnd = 2*lenb-(e-1);
+	for(; i <= loopEnd ;i++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			SeqA = _mm_slli_si128(SeqA, 2);
+			SeqB = _mm_slli_si128(SeqB, 2);
+			SeqA = _mm_insert_epi16(SeqA,*(a-(index)),0);
+			SeqB = _mm_insert_epi16(SeqB,*(b-(index)),0);
+			index++;
+			tmp = _mm_shufflelo_epi16(SeqB,27);
+			tmp = _mm_slli_si128(tmp, 2);
+			tmpValue = _mm_extract_epi16(tmp, 5);
+			tmp = _mm_insert_epi16(tmp, tmpValue, 0);
+			Result = _mm_cmpeq_epi16(SeqA, tmp);
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			if(_mm_extract_epi16(R0, 0) > errThreshold && _mm_extract_epi16(R0, 1) > errThreshold && _mm_extract_epi16(R0, 2) > errThreshold
+					&& _mm_extract_epi16(R0, 3) > errThreshold && _mm_extract_epi16(R0, 4) > errThreshold && _mm_extract_epi16(R1, 0) > errThreshold
+					&& _mm_extract_epi16(R1, 1) > errThreshold && _mm_extract_epi16(R1, 2) > errThreshold  && _mm_extract_epi16(R1, 3) > errThreshold)
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			Result = _mm_cmpeq_epi16(SeqA, _mm_shufflelo_epi16(SeqB,27));
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	j=0;
+	int tmpE = e;
+	for(;j<2*(e-2)+1;j++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		//set the first element
+		if(j==0)
+		{
+			for( k=0;k<=e-1;k++  )
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < e-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else if(j%2 == 0)
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < tmpE-1;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			tmp = _mm_srli_si128(R1,2);
+			for(k=0; k < tmpE-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		i++;
+	}
+	//Diag
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, 2*errThreshold, 0);
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-2)) != *(b-(lenb-1)), 1);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,1,1);
+	Down1 = _mm_insert_epi16(Down1, 2*errThreshold, 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 1);
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+	R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-1)) != *(b-(lenb-1)), 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 0);
+	R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	if(minError > mismatch)
+		return -1;
+	return minError;
+}
+int backwardEditDistance4SSE2(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int i0;
+	int i1;
+	int i2;
+	int i4;
+	int i5;
+	int e = errThreshold;
+	int minError = 2*e;
+	int index = 0;
+	int tmpValue = 0;
+	if(lenb <= e)
+	{
+		return smallEditDistanceB(a,lena,b,lenb);
+	}
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i tmp;
+	__m128i SeqA, SeqB;
+	__m128i Result;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	SeqA = _mm_setzero_si128 ();
+	SeqB = _mm_setzero_si128 ();
+	Result = _mm_setzero_si128 ();
+	/* end initialize */
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	i0 = (a[0] != b[0]);
+	i1 = min(i0, ( *(a-1)!=*b) )+1;
+	i2 = min(i0,( a[0] != *(b-1) ) )+1;
+	i0 = min3( i0+ ( *(a-1)!=*(b-1) ),i1+1,i2+1);
+	i4 = min(i1, ( *(a-2)!=b[0] )+1)+1;
+	i5 = min(i2, (a[0] != *(b-2))+1)+1;
+	R1 = _mm_insert_epi16(R1, 3, 0);
+	R1 = _mm_insert_epi16(R1, i1, 1);
+	R1 = _mm_insert_epi16(R1, i2, 2);
+	R1 = _mm_insert_epi16(R1, 3, 3);
+	R0 = _mm_insert_epi16(R0, 4, 0);
+	R0 = _mm_insert_epi16(R0, i4, 1);
+	R0 = _mm_insert_epi16(R0, i0, 2);
+	R0 = _mm_insert_epi16(R0, i5, 3);
+	R0 = _mm_insert_epi16(R0, 4, 4);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Side2 = _mm_insert_epi16(Side2,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	index = 0;
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		SeqA = _mm_slli_si128(SeqA, 2);
+		SeqB = _mm_slli_si128(SeqB, 2);
+		SeqA = _mm_insert_epi16(SeqA,*(a-index),0);
+		SeqB = _mm_insert_epi16(SeqB,*(b-index),0);
+		index++;
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,2*e,0);
+	index = 4;
+	i = 5;
+	int loopEnd = 2*lenb-(e-1);
+	for(; i <= loopEnd ;i++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			SeqA = _mm_slli_si128(SeqA, 2);
+			SeqB = _mm_slli_si128(SeqB, 2);
+			SeqA = _mm_insert_epi16(SeqA,*(a-(index)),0);
+			SeqB = _mm_insert_epi16(SeqB,*(b-(index)),0);
+			index++;
+			tmp = _mm_shufflelo_epi16(SeqB,27);
+			tmp = _mm_slli_si128(tmp, 2);
+			tmpValue = _mm_extract_epi16(tmp, 5);
+			tmp = _mm_insert_epi16(tmp, tmpValue, 0);
+			Result = _mm_cmpeq_epi16(SeqA, tmp);
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			//tmp = _mm_sub_epi16(Error, R0);
+			//i0 = _mm_movemask_epi8(tmp);
+			if( _mm_extract_epi16(R0, 0) > e && _mm_extract_epi16(R0, 1) > e && _mm_extract_epi16(R0, 2) > e
+					&& _mm_extract_epi16(R0, 3) > e && _mm_extract_epi16(R0, 4) > e && _mm_extract_epi16(R1, 0) > e &&
+					_mm_extract_epi16(R1, 1) > e && _mm_extract_epi16(R1, 2) > e && _mm_extract_epi16(R1, 3) > e  )
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			Result = _mm_cmpeq_epi16(SeqA, _mm_shufflelo_epi16(SeqB,27));
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	j=0;
+	int tmpE = e;
+	for(;j<2*(e-2)+1;j++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		//set the first element
+		if(j==0)
+		{
+			for( k=0;k<=e-1;k++  )
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < e-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else if(j%2 == 0)
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < tmpE-1;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, *(b-(lenb-1-k)) != *(a-((i-lenb)-1+k)),0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			tmp = _mm_srli_si128(R1,2);
+			for(k=0; k < tmpE-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		i++;
+	}
+	//Diag
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, 2*e, 0);
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-2)) != *(b-(lenb-1)), 1);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,1,1);
+	Down1 = _mm_insert_epi16(Down1, 2*e, 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 1);
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+	R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	Diag = _mm_insert_epi16(Diag, *(a-(lenb+e-1)) != *(b-(lenb-1)), 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 0);
+	R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	if(minError > e)
+		return -1;
+	return minError;
+}
+inline int forwardEditDistanceSSE2Extention(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int i0=0;
+	int i1=0;
+	int i2=0;
+	int i4=0;
+	int i5=0;
+	int mismatch = errThreshold;
+	int e = 4;
+	int minError = 4*mismatch+1;
+	int index = 0;
+	int tmpValue = 0;
+	if(lenb <= e)
+	{
+		return smallEditDistanceF(a,lena,b,lenb);
+	}
+	register __m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i tmp;
+	register __m128i SeqA, SeqB;
+	__m128i Result;
+	__m128i tmpSeqA;
+	__m128i tmpSeqB;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	SeqA = _mm_setzero_si128 ();
+	SeqB = _mm_setzero_si128 ();
+	Result = _mm_setzero_si128 ();
+	/* end initialize */
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Diag = _mm_insert_epi16(Diag,minError,0);
+	i0 = (a[0] != b[0]);
+	i1 = min(i0, (a[1]!=b[0]))+1;
+	i2 = min(i0,(a[0]!=b[1]))+1;
+	i0 = min3(i0+(a[1]!=b[1]),i1+1,i2+1);
+	i4 = min(i1, (a[2]!=b[0])+1)+1;
+	i5 = min(i2, (a[0]!=b[2])+1)+1;
+	R1 = _mm_insert_epi16(R1, 3, 0);
+	R1 = _mm_insert_epi16(R1, i1, 1);
+	R1 = _mm_insert_epi16(R1, i2, 2);
+	R1 = _mm_insert_epi16(R1, 3, 3);
+	R0 = _mm_insert_epi16(R0, 4, 0);
+	R0 = _mm_insert_epi16(R0, i4, 1);
+	R0 = _mm_insert_epi16(R0, i0, 2);
+	R0 = _mm_insert_epi16(R0, i5, 3);
+	R0 = _mm_insert_epi16(R0, 4, 4);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Side2 = _mm_insert_epi16(Side2,minError,0);
+	Down1 = _mm_insert_epi16(Down1,minError,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	index = 0;
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		SeqA = _mm_slli_si128(SeqA, 2);
+		SeqB = _mm_slli_si128(SeqB, 2);
+		SeqA = _mm_insert_epi16(SeqA,a[index],0);
+		SeqB = _mm_insert_epi16(SeqB,b[index],0);
+		index++;
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,minError,0);
+	index = 4;
+	i = 5;
+	int loopEnd = 2*lenb-(e-1);
+	for(; i <= loopEnd ;i++)
+	{
+		if( i%2 == 0)
+		{
+			tmpSeqA = _mm_slli_si128(SeqA, 2);
+			tmpSeqB = _mm_slli_si128(SeqB, 2);
+			SeqA = _mm_insert_epi16(tmpSeqA,a[index],0);
+			SeqB = _mm_insert_epi16(tmpSeqB,b[index],0);
+			index++;
+			tmp = _mm_shufflelo_epi16(SeqB,27);
+			tmp = _mm_slli_si128(tmp, 2);
+			tmpValue = _mm_extract_epi16(tmp, 5);
+			tmp = _mm_insert_epi16(tmp, tmpValue, 0);
+			Result = _mm_cmpeq_epi16(SeqA, tmp);
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			if(_mm_extract_epi16(R0, 0) > errThreshold && _mm_extract_epi16(R0, 1) > errThreshold && _mm_extract_epi16(R0, 2) > errThreshold
+					&& _mm_extract_epi16(R0, 3) > errThreshold && _mm_extract_epi16(R0, 4) > errThreshold &&
+					_mm_extract_epi16(R1, 0) > errThreshold && _mm_extract_epi16(R1, 1) > errThreshold &&
+					_mm_extract_epi16(R1, 2) > errThreshold && _mm_extract_epi16(R1, 3) > errThreshold)
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			Result = _mm_cmpeq_epi16(SeqA, _mm_shufflelo_epi16(SeqB,27));
+			Diag =  _mm_andnot_si128(Result, MASK);
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+j=0;
+int tmpE = e;
+for(;j<2*(e-2)+1;j++)
+{
+			Diag = _mm_xor_si128(Diag, Diag);
+			//set the first element
+			if(j==0)
+			{
+				for( k=0;k<=e-1;k++  )
+				{
+					Diag = _mm_slli_si128(Diag, 2);
+					Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+				}
+				R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+				R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+				tmpE--;
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+			else if(j%2 == 0)
+			{
+				for(k=0;k<tmpE;k++)
+				{
+					Diag = _mm_slli_si128(Diag, 2);
+					Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+				}
+				R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+				R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+				tmpE--;
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < tmpE-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+			else
+			{
+				for(k=0;k<tmpE;k++)
+				{
+					Diag = _mm_slli_si128(Diag, 2);
+					Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+				}
+				R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+				R1 = _mm_min_epi16(R1, R0+Down1);
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < tmpE-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+			i++;
+		}
+		//Diag
+Diag = _mm_xor_si128(Diag,Diag);
+Diag = _mm_insert_epi16(Diag, minError, 0);
+Diag = _mm_insert_epi16(Diag, a[lenb+e-2] != b[lenb-1], 1);
+Side1 = _mm_insert_epi16(Side1,1,0);
+Side1 = _mm_insert_epi16(Side1,1,1);
+Down1 = _mm_insert_epi16(Down1, minError, 0);
+Down1 = _mm_insert_epi16(Down1, 1, 1);
+R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+minError = min(minError, _mm_extract_epi16(R1,1));
+Diag = _mm_insert_epi16(Diag, a[lenb+e-1] != b[lenb-1], 0);
+Down1 = _mm_insert_epi16(Down1, 1, 0);
+R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+minError = min(minError, _mm_extract_epi16(R0,0));
+if(minError > mismatch)
+return -1;
+return minError;
+}
+int forwardEditDistance4SSE2(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+return 0;
+int i = 0;
+int j = 0;
+int k = 0;
+int i0=0;
+	int i1=0;
+	int i2=0;
+	int i4=0;
+	int i5=0;
+int e = errThreshold;
+int minError = 2*e;
+	int index = 0;
+	int tmpValue = 0;
+	if(lenb <= e)
+{
+return smallEditDistanceF(a,lena,b,lenb);
+}
+register __m128i R0, R1;
+__m128i Diag;
+__m128i Side1, Side2;
+__m128i Down1, Down2;
+__m128i tmp;
+	register __m128i SeqA, SeqB;
+__m128i Result;
+	__m128i tmpSeqA;
+	__m128i tmpSeqB;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	SeqA = _mm_setzero_si128 ();
+	SeqB = _mm_setzero_si128 ();
+	Result = _mm_setzero_si128 ();
+	/* end initialize */
+R1 = _mm_xor_si128(R1, R1);
+R0 = _mm_xor_si128(R0, R0);
+Diag = _mm_xor_si128(Diag, Diag);
+Diag = _mm_insert_epi16(Diag,2*e,0);
+	i0 = (a[0] != b[0]);
+	i1 = min(i0, (a[1]!=b[0]))+1;
+	i2 = min(i0,(a[0]!=b[1]))+1;
+	i0 = min3(i0+(a[1]!=b[1]),i1+1,i2+1);
+	i4 = min(i1, (a[2]!=b[0])+1)+1;
+	i5 = min(i2, (a[0]!=b[2])+1)+1;
+	R1 = _mm_insert_epi16(R1, 3, 0);
+R1 = _mm_insert_epi16(R1, i1, 1);
+R1 = _mm_insert_epi16(R1, i2, 2);
+	R1 = _mm_insert_epi16(R1, 3, 3);
+	R0 = _mm_insert_epi16(R0, 4, 0);
+	R0 = _mm_insert_epi16(R0, i4, 1);
+	R0 = _mm_insert_epi16(R0, i0, 2);
+	R0 = _mm_insert_epi16(R0, i5, 3);
+	R0 = _mm_insert_epi16(R0, 4, 4);
+Side2 = _mm_xor_si128(Side2, Side2);
+Down2 = _mm_xor_si128(Down2, Down2);
+Down1 = _mm_xor_si128(Down1, Down1);
+	Side1 = _mm_xor_si128(Side1, Side1);
+Side2 = _mm_insert_epi16(Side2,2*e,0);
+Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	index = 0;
+	for(j=0; j < e; j++)
+{
+Side2 = _mm_slli_si128(Side2, 2);
+Side2 = _mm_insert_epi16(Side2,1,0);
+Down1 = _mm_slli_si128(Down1, 2);
+Down1 = _mm_insert_epi16(Down1,1,0);
+Down2 = _mm_slli_si128(Down2, 2);
+Down2 = _mm_insert_epi16(Down2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+	        Side1 = _mm_insert_epi16(Side1,1,0);
+		SeqA = _mm_slli_si128(SeqA, 2);
+SeqB = _mm_slli_si128(SeqB, 2);
+SeqA = _mm_insert_epi16(SeqA,a[index],0);
+SeqB = _mm_insert_epi16(SeqB,b[index],0);
+index++;
+}
+Down2= _mm_slli_si128(Down2, 2);
+Down2 = _mm_insert_epi16(Down2,2*e,0);
+	index = 4;
+	i = 5;
+	int loopEnd = 2*lenb-(e-1);
+for(; i <= loopEnd ;i++)
+{
+//Diag = _mm_xor_si128(Diag, Diag);
+if( i%2 == 0)
+{
+			tmpSeqA = _mm_slli_si128(SeqA, 2);
+tmpSeqB = _mm_slli_si128(SeqB, 2);
+SeqA = _mm_insert_epi16(tmpSeqA,a[index],0);
+SeqB = _mm_insert_epi16(tmpSeqB,b[index],0);
+index++;
+			tmp = _mm_shufflelo_epi16(SeqB,27);
+			tmp = _mm_slli_si128(tmp, 2);
+			tmpValue = _mm_extract_epi16(tmp, 5);
+			tmp = _mm_insert_epi16(tmp, tmpValue, 0);
+Result = _mm_cmpeq_epi16(SeqA, tmp);
+Diag =  _mm_andnot_si128(Result, MASK);
+R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+if(_mm_extract_epi16(R0, 0) > e && _mm_extract_epi16(R0, 1) > e && _mm_extract_epi16(R0, 2) > e
+				&& _mm_extract_epi16(R0, 3) > e && _mm_extract_epi16(R0, 4) > e && _mm_extract_epi16(R1, 0) > e &&
+				  _mm_extract_epi16(R1, 1) > e && _mm_extract_epi16(R1, 2) > e && _mm_extract_epi16(R1, 3) > e)
+return -1;
+	               if(i == 2*lenb-e)
+{
+tmp = _mm_srli_si128(R0,2);
+for(k=0; k < e-1;k++)
+tmp = _mm_srli_si128(tmp,2);
+minError = _mm_extract_epi16(tmp,0);
+}
+}
+else
+{
+Result = _mm_cmpeq_epi16(SeqA, _mm_shufflelo_epi16(SeqB,27));
+Diag =  _mm_andnot_si128(Result, MASK);
+R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+R1 = _mm_min_epi16(R1, R0+Down1);
+if(i >= 2*lenb-e)
+{
+tmp = _mm_srli_si128(R1,2);
+for(k=0; k < e-2;k++)
+tmp = _mm_srli_si128(tmp,2);
+minError = min(minError, _mm_extract_epi16(tmp,0));
+}
+}
+	}
+j=0;
+int tmpE = e;
+for(;j<2*(e-2)+1;j++)
+{
+Diag = _mm_xor_si128(Diag, Diag);
+//set the first element
+if(j==0)
+{
+for( k=0;k<=e-1;k++  )
+{
+Diag = _mm_slli_si128(Diag, 2);
+Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+}
+R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+tmpE--;
+tmp = _mm_srli_si128(R0,2);
+for(k=0; k < e-2;k++)
+tmp = _mm_srli_si128(tmp,2);
+minError = min(minError, _mm_extract_epi16(tmp,0));
+}
+		else if(j%2 == 0)
+{
+for(k=0;k<tmpE;k++)
+{
+Diag = _mm_slli_si128(Diag, 2);
+Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+}
+R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+tmpE--;
+tmp = _mm_srli_si128(R0,2);
+for(k=0; k < tmpE-1;k++)
+tmp = _mm_srli_si128(tmp,2);
+minError = min(minError, _mm_extract_epi16(tmp,0));
+}
+		else
+{
+for(k=0;k<tmpE;k++)
+{
+Diag = _mm_slli_si128(Diag, 2);
+Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+}
+R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+R1 = _mm_min_epi16(R1, R0+Down1);
+tmp = _mm_srli_si128(R1,2);
+for(k=0; k < tmpE-2;k++)
+tmp = _mm_srli_si128(tmp,2);
+minError = min(minError, _mm_extract_epi16(tmp,0));
+}
+i++;
+}
+	//Diag
+Diag = _mm_xor_si128(Diag,Diag);
+Diag = _mm_insert_epi16(Diag, 2*e, 0);
+Diag = _mm_insert_epi16(Diag, a[lenb+e-2] != b[lenb-1], 1);
+Side1 = _mm_insert_epi16(Side1,1,0);
+Side1 = _mm_insert_epi16(Side1,1,1);
+Down1 = _mm_insert_epi16(Down1, 2*e, 0);
+Down1 = _mm_insert_epi16(Down1, 1, 1);
+R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+minError = min(minError, _mm_extract_epi16(R1,1));
+Diag = _mm_insert_epi16(Diag, a[lenb+e-1] != b[lenb-1], 0);
+Down1 = _mm_insert_epi16(Down1, 1, 0);
+R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+minError = min(minError, _mm_extract_epi16(R0,0));
+if(minError > e)
+return -1;
+return minError;
+}
+int forwardEditDistanceSSE2Odd(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int e = errThreshold;
+	int minError = 2*e;
+	char flag = 0;
+	if(lenb <= e)
+	{
+		return smallEditDistanceF(a,lena,b,lenb);
+	}
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i Error;
+	__m128i tmp;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Error = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,2*e,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,1,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,2);
+	R0 = _mm_insert_epi16(R0,0,0);
+	R1 = _mm_insert_epi16(R1,1,0);
+	R1 = _mm_insert_epi16(R1,1,1);
+	for(i=2; i <= e; i++)
+	{
+		//set side
+		Side1 = _mm_slli_si128(Side1,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_slli_si128(Down1,2);
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=1;j<=i-1;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[i/2-1+(i/2-j)] != a[i/2-1-(i/2-j)],0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R0 = _mm_min_epi16(R1+Side1, _mm_slli_si128(R0,2)+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down1);
+		}
+		else
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=i/2-1;j>=-i/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[(i+1)/2+j-1] != a[(i-1)/2-j-1],0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down1);
+		}
+	}
+	Error = _mm_xor_si128(Error, Error);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Error = _mm_insert_epi16(Error,e,0);
+	Side2 = _mm_insert_epi16(Side2,2*e,0);
+	Side1 = _mm_insert_epi16(Side2,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Side1 = _mm_slli_si128(Side1, 2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Error = _mm_slli_si128(Error, 2);
+		Error = _mm_insert_epi16(Error, e, 0);
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,2*e,0);
+	for(; i <= 2*lenb-(e-1);i++)
+	{
+		flag = 0;
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			for(j=e/2;j>=-e/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[i/2-1+j] != a[i/2-1-j],0);
+			}
+			R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+			R0 = _mm_min_epi16(R0, R1+Down1);
+			if(_mm_extract_epi16(R0,0) <= e)
+				flag = 1;
+			tmp = _mm_srli_si128(R0,2);
+			for(j=0; j < e-1;j++)
+			{
+				if(_mm_extract_epi16(tmp,0) <= e)
+					flag = 1;
+				tmp = _mm_srli_si128(tmp,2);
+			}
+			//			printf("#%d %d %d\n", _mm_extract_epi16(R0,0), _mm_extract_epi16(R0,1), _mm_extract_epi16(R0,2));
+			if(flag == 0)
+				return -1;
+			if(i == 2*lenb-(e-1))
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			for(j=e/2;j>=-e/2-1;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[(i+1)/2+j-1] != a[(i)/2-j-1],0);
+			}
+			R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+			//printf("#%d %d %d %d\n", _mm_extract_epi16(R1,0), _mm_extract_epi16(R1,1), _mm_extract_epi16(R1,2),
+			//	_mm_extract_epi16(R1,3));
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	//first cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, b[lenb-3] != a[lena], 0);
+	Diag = _mm_insert_epi16(Diag, b[lenb-2] != a[lena-1], 1);
+	Diag = _mm_insert_epi16(Diag, b[lenb-1] != a[lena-2], 2);
+	Diag = _mm_insert_epi16(Diag, 2*e, 3);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R1,2));
+	//second cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, b[lenb-2] != a[lena], 0);
+	Diag = _mm_insert_epi16(Diag, b[lenb-1] != a[lena-1], 1);
+	Diag = _mm_insert_epi16(Diag, 2*e, 2);
+	R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0, R1+Down1);
+	minError = min(minError, _mm_extract_epi16(R0,1));
+	//third cell
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, b[lenb-2] != a[lena+1], 0);
+	Diag = _mm_insert_epi16(Diag, b[lenb-1] != a[lena], 1);
+	Diag = _mm_insert_epi16(Diag, 2*e, 2);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	//forth
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, b[lenb-1] != a[lena+1], 0);
+	Diag = _mm_insert_epi16(Diag, 2*e, 1);
+	R0 = _mm_min_epi16(_mm_srli_si128(R1,2)+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0, R1+Down1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	//fifth
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, b[lenb-1] != a[lena+2], 0);
+	Diag = _mm_insert_epi16(Diag, 2*e, 1);
+	R1 = _mm_min_epi16(R0+Side2, R1+Diag);
+	R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down2);
+	minError = min(minError, _mm_extract_epi16(R1,0));
+	if(minError > e)
+		return -1;
+	return minError;
+}
+int forwardEditDistanceSSE2G(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+		return 0;
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int e = errThreshold;
+	int minError = 2*e;
+	char flag = 0;
+	if(lenb <= e)
+	{
+		return smallEditDistanceF(a,lena,b,lenb);
+	}
+	__m128i R0, R1;
+	__m128i Diag;
+	__m128i Side1, Side2;
+	__m128i Down1, Down2;
+	__m128i Error;
+	__m128i tmp;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Error = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	R1 = _mm_xor_si128(R1, R1);
+	R0 = _mm_xor_si128(R0, R0);
+	Diag = _mm_xor_si128(Diag, Diag);
+	Side1 = _mm_xor_si128(Side1, Side1);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,2*e,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,1,1);
+	Down1 = _mm_insert_epi16(Down1,2*e,2);
+	R0 = _mm_insert_epi16(R0,0,0);
+	R1 = _mm_insert_epi16(R1,1,0);
+	R1 = _mm_insert_epi16(R1,1,1);
+	for(i=2; i <= e; i++)
+	{
+		//set side
+		Side1 = _mm_slli_si128(Side1,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_slli_si128(Down1,2);
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=1;j<=i-1;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[i/2-1+(i/2-j)] != a[i/2-1-(i/2-j)],0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R0 = _mm_min_epi16(R1+Side1, _mm_slli_si128(R0,2)+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down1);
+		}
+		else
+		{
+			Diag = _mm_insert_epi16(Diag,2*e,0);
+			for(j=i/2-1;j>=-i/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[(i+1)/2+j-1] != a[(i-1)/2-j-1],0);
+			}
+			Diag = _mm_slli_si128(Diag, 2);
+			Diag = _mm_insert_epi16(Diag, 2*e,0);
+			R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+			R1 = _mm_min_epi16(R1,  _mm_slli_si128(R0,2)+Down1);
+		}
+	}
+	Error = _mm_xor_si128(Error, Error);
+	Side2 = _mm_xor_si128(Side2, Side2);
+	Down2 = _mm_xor_si128(Down2, Down2);
+	Down1 = _mm_xor_si128(Down1, Down1);
+	Error = _mm_insert_epi16(Error,e,0);
+	Side2 = _mm_insert_epi16(Side2,2*e,0);
+	Down1 = _mm_insert_epi16(Down1,2*e,0);
+	for(j=0; j < e; j++)
+	{
+		Side2 = _mm_slli_si128(Side2, 2);
+		Side2 = _mm_insert_epi16(Side2,1,0);
+		Down1 = _mm_slli_si128(Down1, 2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down2 = _mm_slli_si128(Down2, 2);
+		Down2 = _mm_insert_epi16(Down2,1,0);
+		Error = _mm_slli_si128(Error, 2);
+		Error = _mm_insert_epi16(Error, e, 0);
+	}
+	Down2= _mm_slli_si128(Down2, 2);
+	Down2 = _mm_insert_epi16(Down2,2*e,0);
+	for(; i <= 2*lenb-(e-1);i++)
+	{
+		flag = 0;
+		Diag = _mm_xor_si128(Diag, Diag);
+		if( i%2 == 0)
+		{
+			for(j=e/2;j>=-e/2;j--)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[i/2-1+j] != a[i/2-1-j],0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			if(_mm_extract_epi16(R0,0) <= e)
+				flag = 1;
+			tmp = _mm_srli_si128(R0,2);
+			for(j=0; j < e-1;j++)
+			{
+				if(_mm_extract_epi16(tmp,0) <= e)
+					flag = 1;
+				tmp = _mm_srli_si128(tmp,2);
+			}
+			if(flag == 0)
+				return -1;
+			if(i == 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R0,2);
+				for(k=0; k < e-1;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = _mm_extract_epi16(tmp,0);
+			}
+		}
+		else
+		{
+			for(j=-e/2+1;j<=e/2;j++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[(i+1)/2-j-1] != a[(i-1)/2+j-1],0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			if(i >= 2*lenb-e)
+			{
+				tmp = _mm_srli_si128(R1,2);
+				for(k=0; k < e-2;k++)
+					tmp = _mm_srli_si128(tmp,2);
+				minError = min(minError, _mm_extract_epi16(tmp,0));
+			}
+		}
+	}
+	j=0;
+	int tmpE = e;
+	for(;j<2*(e-2)+1;j++)
+	{
+		Diag = _mm_xor_si128(Diag, Diag);
+		//set the first element
+		if(j==0)
+		{
+			for( k=0;k<=e-1;k++  )
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < e-2;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else if(j%2 == 0)
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+			}
+			R0 = _mm_min_epi16(R1+Side2, R0+Diag);
+			R0 = _mm_min_epi16(R0, _mm_slli_si128(R1,2)+Down2);
+			tmpE--;
+			tmp = _mm_srli_si128(R0,2);
+			for(k=0; k < tmpE-1;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		else
+		{
+			for(k=0;k<tmpE;k++)
+			{
+				Diag = _mm_slli_si128(Diag, 2);
+				Diag = _mm_insert_epi16(Diag, b[lenb-1-k] != a[(i-lenb)-1+k],0);
+			}
+			R1 = _mm_min_epi16(_mm_srli_si128(R0,2)+Side1, R1+Diag);
+			R1 = _mm_min_epi16(R1, R0+Down1);
+			tmp = _mm_srli_si128(R1,2);
+			for(k=0; k < tmpE-1;k++)
+				tmp = _mm_srli_si128(tmp,2);
+			minError = min(minError, _mm_extract_epi16(tmp,0));
+		}
+		i++;
+	}
+	//Diag
+	Diag = _mm_xor_si128(Diag,Diag);
+	Diag = _mm_insert_epi16(Diag, 2*e, 0);
+	Diag = _mm_insert_epi16(Diag, a[lenb+e-2] != b[lenb-1], 1);
+	Side1 = _mm_insert_epi16(Side1,1,0);
+	Side1 = _mm_insert_epi16(Side1,1,1);
+	Down1 = _mm_insert_epi16(Down1, 2*e, 0);
+	Down1 = _mm_insert_epi16(Down1, 1, 1);
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+	R1 = _mm_min_epi16(R1, _mm_slli_si128(R0,2)+Down1);
+	minError = min(minError, _mm_extract_epi16(R1,1));
+	Diag = _mm_insert_epi16(Diag, a[lenb+e-1] != b[lenb-1], 1);
+	Down1 = _mm_insert_epi16(Down1, 1, 0);
+	R0 = _mm_min_epi16(R1+Down1,R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_srli_si128(R1,2)+Side1);
+	minError = min(minError, _mm_extract_epi16(R0,0));
+	if(minError > e)
+		return -1;
+	return minError;
+}
+int forwardEditDistance2SSE2(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+return 0;
+int i0 = 0;
+int i1 = 0;
+int error;                      //0: if the two character are equal 1: if not
+int i = 0;                      //loop index
+int e = 2;                      //error bound
+	int totalError = 0;
+__m128i R0;
+__m128i R1;
+__m128i Side1, Side2,Side;                    //side matrix
+__m128i Down1, Down2,Down;                    //down matrix
+__m128i Diag;
+__m128i tmp;
+	__m128i ERROR_REACH;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Side = _mm_setzero_si128 ();
+Down = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+ERROR_REACH = _mm_setzero_si128 ();
+	/* end initialize */
+	if(lenb <= e)
+{
+return smallEditDistanceF(a,lena,b,lenb);
+}
+	ERROR_REACH = _mm_set_epi16(0,0,0,0,0,e,e,e);
+	R0 = _mm_insert_epi16(R0,0,0);
+R1 = _mm_insert_epi16(R1,1,0);
+R1 = _mm_insert_epi16(R1,1,1);
+//        error = ((a[0]) != (b[0]));
+	Diag  = _mm_set_epi16(0,0,0,0,0,2*e,((a[0]) != (b[0])),2*e);
+	Side1 = _mm_set_epi16(0,0,0,0,0,2*e,1,1);
+	Side2 = _mm_set_epi16(0,0,0,0,0,1,1,2*e);
+	Down1 = _mm_set_epi16(0,0,0,0,0,2*e,1,1);
+	Down2 = _mm_set_epi16(0,0,0,0,0,1,1,2*e);
+tmp = _mm_slli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+R0 = _mm_min_epi16(R0,tmp+Down2);
+	for (i = 3; i  < 2*lena; i++)
+{
+if(i % 2 ==1)
+{
+			Diag = _mm_xor_si128(Diag, Diag);
+error = ((a[(i+1)/2-1]) != (b[(i-1)/2-1]));
+Diag = _mm_insert_epi16(Diag,error,0);
+error = ((a[(i-1)/2-1]) != (b[(i+1)/2-1]));
+Diag = _mm_insert_epi16(Diag,error,1);
+//		        Diag  = _mm_set_epi16(0, 0, 0, 0, 0, 0, ((a[(i-1)/2-1]) != (b[(i+1)/2-1])) ,((a[(i+1)/2-1]) != (b[(i-1)/2-1])));
+tmp = _mm_srli_si128(R0,2);
+R1 = _mm_min_epi16(tmp+Side1, R1+Diag);
+R1 = _mm_min_epi16(R1,R0+Down1);
+if(i > 2 * lenb - 2)
+			{
+	                        i1 = _mm_extract_epi16(R1, 1);
+totalError = min(totalError, i1);
+			}
+}
+		else if(i % 2 == 0)
+{
+error = ((a[i/2]) != (b[i/2-2]));
+Diag = _mm_insert_epi16(Diag,error,0);
+error = ((a[i/2-1]) != (b[i/2-1]));
+Diag = _mm_insert_epi16(Diag,error,1);
+error = ((a[i/2-2]) != (b[i/2]));
+Diag = _mm_insert_epi16(Diag,error,2);
+//                Diag  = _mm_set_epi16(0, 0, 0, 0, 0, ((a[i/2-2]) != (b[i/2])) , ((a[i/2-1]) != (b[i/2-1]))  , ((a[i/2]) != (b[i/2-2])) );
+tmp = _mm_slli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+R0 = _mm_min_epi16(R0,tmp+Down2);
+			tmp = _mm_sub_epi16(ERROR_REACH, R0);
+i0 = _mm_movemask_epi8(tmp);
+if(i0 == 63 &&  _mm_extract_epi16(R1,0) > errThreshold  && _mm_extract_epi16(R1,1) > errThreshold && i < 2 * lenb - 2)
+return -1;
+if(i == 2 * lenb - 2) {
+				totalError = _mm_extract_epi16(R0, 2);
+			}
+}
+}
+Down1 = _mm_insert_epi16(Down1,2*e,0);
+	//fill the first part of the error
+	error = ((a[i/2]) != (b[i/2-2]));
+Diag = _mm_insert_epi16(Diag,error,0);
+	error = ((a[i/2-1]) != (b[i/2-1]));
+	Diag = _mm_insert_epi16(Diag,error,1);
+	Diag = _mm_insert_epi16(Diag,2*e,2);
+//        Diag  = _mm_set_epi16(0, 0, 0, 0, 0, 2*e , ((a[i/2-1]) != (b[i/2-1]))  , ((a[i/2]) != (b[i/2-2])) );
+	R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_slli_si128(R1,2)+Down1);
+//	i0 = _mm_extract_epi16(R0, 0);
+	i1 = _mm_extract_epi16(R0, 1);
+	totalError = min(totalError, i1);
+	//fill the second part of the error
+	i++;
+	Diag = _mm_xor_si128(Diag, Diag);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	error = ((a[i/2]) != (b[lenb-1]));
+Diag = _mm_insert_epi16(Diag,error,1);
+	Diag = _mm_insert_epi16(Diag,2*e,2);
+//        Diag  = _mm_set_epi16(0, 0, 0, 0, 0, 2*e , ((a[i/2]) != (b[lenb-1]))  , 2*e );
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+R1 = _mm_min_epi16(R1,_mm_slli_si128(R0,2)+Down1);
+//        i0 = _mm_extract_epi16(R1, 0);
+i1 = _mm_extract_epi16(R1, 1);
+totalError = min(totalError, i1);
+	//fill the last the last element of the matrix
+	i++;
+	Diag = _mm_xor_si128(Diag, Diag);
+	error = ((a[i/2]) != (b[lenb-1]));
+Diag = _mm_insert_epi16(Diag,error,0);
+//        Diag  = _mm_set_epi16(0, 0, 0, 0, 0, 0 , 0  , ((a[i/2]) != (b[lenb-1])) );
+Down = _mm_insert_epi16(Down,1,0);
+Side = _mm_insert_epi16(Side,1,0);
+tmp = _mm_srli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Down, _mm_srli_si128(R0,2)+Diag);
+R0 = _mm_min_epi16(R0,tmp+Side);
+i0 = _mm_extract_epi16(R0, 0);
+totalError = min(totalError, i0);
+	if(totalError > e)
+		return -1;
+	return totalError;
+}
+int backwardEditDistance2SSE2(char *a, int lena, char *b,int lenb)
+{
+	if(lenb == 0 || lena == 0)
+return 0;
+int i0 = 0;
+int i1 = 0;
+int error;                      //0: if the two character are equal 1: if not
+int i = 0;                      //loop index
+int e = 2;                      //error bound
+	int totalError = 0;
+__m128i R0;
+__m128i R1;
+__m128i Side1, Side2,Side;                    //side matrix
+__m128i Down1, Down2,Down;                    //down matrix
+__m128i Diag;                   	      //diag matrix
+__m128i tmp;
+	__m128i ERROR_REACH;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Side = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+	Down = _mm_setzero_si128 ();
+	ERROR_REACH = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	if(lenb <= e)
+{
+return smallEditDistanceB(a,lena,b,lenb);
+}
+	ERROR_REACH = _mm_set_epi16(0,0,0,0,0,e,e,e);
+	R0 = _mm_insert_epi16(R0,0,0);
+R1 = _mm_insert_epi16(R1,1,0);
+R1 = _mm_insert_epi16(R1,1,1);
+error = ((a[0]) != (b[0]));
+Diag = _mm_insert_epi16(Diag,2*e,0);
+Diag = _mm_insert_epi16(Diag,error,1);
+Diag = _mm_insert_epi16(Diag,2*e,2);
+Side1 = _mm_insert_epi16(Side1,1,0);
+Side1 = _mm_insert_epi16(Side1,1,1);
+Side1 = _mm_insert_epi16(Side1,2*e,2);
+Side2 = _mm_insert_epi16(Side2,2*e,0);
+Side2 = _mm_insert_epi16(Side2,1,1);
+Side2 = _mm_insert_epi16(Side2,1,2);
+	Down1 = _mm_insert_epi16(Down1,1,0);
+Down1 = _mm_insert_epi16(Down1,1,1);
+Down1 = _mm_insert_epi16(Down1,2*e,2);
+Down2 = _mm_insert_epi16(Down2,2*e,0);
+Down2 = _mm_insert_epi16(Down2,1,1);
+Down2 = _mm_insert_epi16(Down2,1,2);
+tmp = _mm_slli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+R0 = _mm_min_epi16(R0,tmp+Down2);
+//        printf("%d %d %d\n", _mm_extract_epi16(R0,0), _mm_extract_epi16(R0,1), _mm_extract_epi16(R0,2));
+	for (i = 3; i  < 2*lena; i++)
+{
+if(i % 2 ==1)
+{
+Diag = _mm_sub_epi8(Diag, Diag);
+error = ( *(a-((i+1)/2-1)) != *(b-((i-1)/2-1)) );
+Diag = _mm_insert_epi16(Diag,error,0);
+error = ( *(a-((i-1)/2-1)) != *(b-((i+1)/2-1)) );
+Diag = _mm_insert_epi16(Diag,error,1);
+			//printf("#%d #%d\n", _mm_extract_epi16(Diag,0), _mm_extract_epi16(Diag,1));
+tmp = _mm_srli_si128(R0,2);
+R1 = _mm_min_epi16(tmp+Side1, R1+Diag);
+R1 = _mm_min_epi16(R1,R0+Down1);
+if(i > 2 * lenb - 2) {
+	                        i1 = _mm_extract_epi16(R1, 1);
+totalError = min(totalError, i1);
+			}
+//			printf("%d %d\n", _mm_extract_epi16(R1,0), _mm_extract_epi16(R1,1));
+}
+		else if(i % 2 == 0)
+{
+error = ( *(a-(i/2))   != *(b-(i/2-2)) );
+Diag = _mm_insert_epi16(Diag,error,0);
+error = ( *(a-(i/2-1)) != *(b-(i/2-1)) );
+Diag = _mm_insert_epi16(Diag,error,1);
+error = ( *(a-(i/2-2))    != *(b-(i/2)));
+Diag = _mm_insert_epi16(Diag,error,2);
+tmp = _mm_slli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+R0 = _mm_min_epi16(R0,tmp+Down2);
+			tmp = _mm_sub_epi16(ERROR_REACH, R0);
+			i0 = _mm_movemask_epi8(tmp);
+if(i0 == 63 && _mm_extract_epi16(R1,0) > errThreshold && _mm_extract_epi16(R1,1) > errThreshold && i < 2 * lenb - 2)
+return -1;
+if(i == 2 * lenb - 2) {
+				totalError = _mm_extract_epi16(R0, 2);
+			}
+}
+}
+Down1 = _mm_insert_epi16(Down1,2*e,0);
+	//fill the first part of the error
+	error = ( *(a-(i/2))  != *(b-(i/2-2)) );
+Diag = _mm_insert_epi16(Diag,error,0);
+	error = ( *(a-(i/2-1)) != *(b-(i/2-1)) );
+	Diag = _mm_insert_epi16(Diag,error,1);
+	Diag = _mm_insert_epi16(Diag,2*e,2);
+	R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+	R0 = _mm_min_epi16(R0,_mm_slli_si128(R1,2)+Down1);
+	i0 = _mm_extract_epi16(R0, 0);
+	i1 = _mm_extract_epi16(R0, 1);
+	totalError = min(totalError, i1);
+	//fill the second part of the error
+	i++;
+	Diag = _mm_sub_epi8(Diag, Diag);
+	Diag = _mm_insert_epi16(Diag,2*e,0);
+	error = ( *(a-(i/2)) != *(b-(lenb-1)) );
+Diag = _mm_insert_epi16(Diag,error,1);
+	Diag = _mm_insert_epi16(Diag,2*e,2);
+	R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+R1 = _mm_min_epi16(R1,_mm_slli_si128(R0,2)+Down1);
+i0 = _mm_extract_epi16(R1, 0);
+i1 = _mm_extract_epi16(R1, 1);
+totalError = min(totalError, i1);
+	//fill the last the last element of the matrix
+	i++;
+	Diag = _mm_sub_epi8(Diag, Diag);
+	error = ( *(a-(i/2)) != *(b-(lenb-1)) );
+Diag = _mm_insert_epi16(Diag,error,0);
+Down = _mm_insert_epi16(Down,1,0);
+Side = _mm_insert_epi16(Side,1,0);
+tmp = _mm_srli_si128(R1,2);
+R0 = _mm_min_epi16(R1+Down, _mm_srli_si128(R0,2)+Diag);
+R0 = _mm_min_epi16(R0,tmp+Side);
+i0 = _mm_extract_epi16(R0, 0);
+totalError = min(totalError, i0);
+	if(totalError > e || totalError == 0)
+		return -1;
+return totalError;
+}
+void initBestMapping(int totalReadNumber)
+{
+	int i = 0;
+	bestHitMappingInfo = getMem(totalReadNumber * sizeof(BestFullMappingInfo));
+	for(i = 0; i < totalReadNumber; i++) {
+		bestHitMappingInfo[i].loc = -1;
+	}
+}
+void finalizeBestSingleMapping()
+{
+	int i = 0;
+	char *_tmpQual, *_tmpSeq;
+	char rqual[SEQ_LENGTH + 1];
+	rqual[SEQ_LENGTH]='\0';
+	for(i = 0; i < _msf_seqListSize; i++)
+	{
+		if(_msf_seqList[i].hits[0] != 0)
+		{
+			if (bestHitMappingInfo[i].dir)
+			{
+				reverse(_msf_seqList[i].qual, rqual, SEQ_LENGTH);
+				_tmpQual = rqual;
+				_tmpSeq = _msf_seqList[i].rseq;
+			}
+			else
+			{
+				_tmpQual = _msf_seqList[i].qual;
+				_tmpSeq = _msf_seqList[i].seq;
+			}
+			_msf_output.QNAME               = _msf_seqList[i].name;
+			_msf_output.FLAG                = 16 * bestHitMappingInfo[i].dir;
+			_msf_output.RNAME               = bestHitMappingInfo[i].chr;
+			_msf_output.POS                 = bestHitMappingInfo[i].loc;
+			_msf_output.MAPQ                = 255;
+			_msf_output.CIGAR               = bestHitMappingInfo[i].cigar ;
+			_msf_output.MRNAME              = "*";
+			_msf_output.MPOS                = 0;
+			_msf_output.ISIZE               = 0;
+			_msf_output.SEQ                 = _tmpSeq;
+			_msf_output.QUAL                = _tmpQual;
+			_msf_output.optSize             = 2;
+			_msf_output.optFields   = _msf_optionalFields;
+			_msf_optionalFields[0].tag = "NM";
+			_msf_optionalFields[0].type = 'i';
+			_msf_optionalFields[0].iVal = bestHitMappingInfo[i].err;
+			_msf_optionalFields[1].tag = "MD";
+			_msf_optionalFields[1].type = 'Z';
+			_msf_optionalFields[1].sVal = bestHitMappingInfo[i].md;
+			output(_msf_output);
+		}
+	}
+	freeMem(bestHitMappingInfo, _msf_seqListSize * sizeof(FullMappingInfo));
+}
+/**********************************************/
+int compare (const void *a, const void *b)
+{
+	return ((Pair *)a)->hv - ((Pair *)b)->hv;
+	/*char *s1 = ((Pair *)a)->hv;
+	char *s2 = ((Pair *)b)->hv;
+	int i = 0;
+	int diff = 0;
+	int sign = 0;
+	for(i = 0; i < SEQ_LENGTH; i++)
+	{
+		diff += (s1[i] != s2[i]);
+		if(s1[i] > s2[i])
+			sign++;
+		else if(s1[i] < s2[i])
+			sign--;
+	}
+	return diff*sign;*/
+//	return strncmp(s1, s2,SEQ_LENGTH);
+}
+/**********************************************/
+void preProcessReads()
+{
+	int i = 0;
+	_msf_sort_seqList = getMem(_msf_seqListSize * sizeof(Pair));
+	for(i = 0; i < _msf_seqListSize; i++)
+	{
+		_msf_sort_seqList[i].hv = hashVal(_msf_seqList[i].seq);
+		_msf_sort_seqList[i].readNumber = i;
+	}
+	qsort(_msf_sort_seqList, _msf_seqListSize, sizeof(Pair), compare);
+	/*
+	for(i = 0; i < _msf_seqListSize; i++)
+{
+		//printf("%s\n", _msf_sort_seqList[i].hv);
+	}
+	*/
+	_msf_map_sort_seqList = getMem(_msf_seqListSize * sizeof(int));
+	for(i = 0; i < _msf_seqListSize; i++)
+		_msf_map_sort_seqList[_msf_seqList[i].readNumber] = i;
+}
+/**********************************************/
+int verifySingleEnd(int index, char* seq, int offset)
+{
+	int curOff = 0;
+	int i;
+	char *ref;
+	int err;
+	int errCnt =0;
+	int errCntOff = 0;
+	int NCntOff = 0;
+	ref = _msf_refGen + index - 1;
+	verificationCnt++;
+	for (i = 0; i < SEQ_LENGTH; i++)
+	{
+		err	= *ref != *seq;
+		errCnt += err;
+		if (errCnt > errThreshold)
+		{
+			return -1;
+		}
+		if (i >= _msf_samplingLocs[curOff] && i <= _msf_samplingLocsEnds[curOff])
+		{
+			errCntOff +=  err;
+			NCntOff += (*seq == 'N');
+		}
+		else if (curOff < _msf_samplingLocsSize && i>=_msf_samplingLocs[curOff+1])
+		{
+			if (errCntOff == 0 && NCntOff == 0 && offset > curOff)
+			{
+				return -1;
+			}
+			errCntOff = 0;
+			NCntOff = 0;
+			curOff++;
+			if ( i >= _msf_samplingLocs[curOff])
+			{
+				errCntOff += err;
+				NCntOff += (*seq == 'N');
+			}
+		}
+		ref++;
+		seq++;
+	}
+	return errCnt;
+}
+/*********************************************/
+void initFAST(Read *seqList, int seqListSize, int *samplingLocs, int samplingLocsSize, char *genFileName)
+{
+	int i;
+	if (_msf_optionalFields == NULL)
+	{
+		_msf_op = getMem(SEQ_LENGTH);
+		if (pairedEndMode)
+		{
+			_msf_optionalFields = getMem(8*sizeof(OPT_FIELDS));
+		}
+		else
+		{
+			_msf_optionalFields = getMem(2*sizeof(OPT_FIELDS));
+		}
+		for (i=0; i<200;i++)
+		{
+			sprintf(_msf_numbers[i],"%d%c",i, '\0');
+		}
+		sprintf(_msf_cigar, "%dM", SEQ_LENGTH);
+	}
+	if (_msf_samplingLocsEnds == NULL)
+	{
+		_msf_samplingLocs = samplingLocs;
+		_msf_samplingLocsSize = samplingLocsSize;
+		_msf_samplingLocsEnds = getMem(sizeof(int)*_msf_samplingLocsSize);
+		for (i=0; i<_msf_samplingLocsSize; i++)
+		{
+			_msf_samplingLocsEnds[i]=_msf_samplingLocs[i]+WINDOW_SIZE-1;
+		}
+		_msf_seqList = seqList;
+		_msf_seqListSize = seqListSize;
+		preProcessReads();
+		_msf_oeaMapping = getMem(_msf_seqListSize * sizeof(int));
+		for(i = 0; i < _msf_seqListSize; i++)
+		{
+			_msf_oeaMapping[i] = 0;
+		}
+		_msf_discordantMapping = getMem(_msf_seqListSize * sizeof(int));
+		for(i = 0; i < _msf_seqListSize; i++)
+		{
+			_msf_discordantMapping[i] = 0;
+		}
+	}
+	if (_msf_refGenName == NULL)
+	{
+		_msf_refGenName = getMem(4*SEQ_LENGTH);
+	}
+	_msf_refGen =  getRefGenome();
+	_msf_refGenLength = strlen(_msf_refGen);
+	_msf_refGenOffset = getRefGenomeOffset();
+	snprintf(_msf_refGenName, 4*SEQ_LENGTH,"%s%c", getRefGenomeName(), '\0');
+	_msf_refGenName[strlen(getRefGenomeName())] = '\0';
+	if (_msf_verifiedLocs != NULL){
+		freeMem(_msf_verifiedLocs, sizeof(int) * (_msf_refGenLength+1));
+	}
+	_msf_verifiedLocs = (int *) getMem(sizeof(int)*(_msf_refGenLength+1));
+	for (i=0; i<=_msf_refGenLength; i++)
+		_msf_verifiedLocs[i] = _msf_seqListSize*10+1;
+	if (pairedEndMode && _msf_seqHits == NULL)
+	{
+		_msf_mappingInfo  = getMem(seqListSize * sizeof (MappingInfo));
+		for (i=0; i<seqListSize; i++)
+		{
+			//_msf_mappingInfo[i].next = getMem(sizeof(MappingLocations));
+			_msf_mappingInfo[i].next = NULL;
+			_msf_mappingInfo[i].size = 0;
+		}
+		_msf_seqHits = getMem((_msf_seqListSize) * sizeof(int));
+		for (i=0; i<_msf_seqListSize; i++)
+		{
+			_msf_seqHits[i] = 0;
+		}
+		_msf_readHasConcordantMapping = getMem(_msf_seqListSize / 2 * sizeof(char));
+		for(i = 0; i < _msf_seqListSize/2; i++)
+		{
+			_msf_readHasConcordantMapping[i] = 0;
+		}
+		initLoadingRefGenome(genFileName);
+	}
+	if (_msf_refGenOffset == 0)
+	{
+		_msf_refGenBeg = 1;
+	}
+	else
+	{
+		_msf_refGenBeg = CONTIG_OVERLAP - SEQ_LENGTH + 2;
+	}
+	_msf_refGenEnd = _msf_refGenLength - SEQ_LENGTH + 1;
+}
+/**********************************************/
+void finalizeFAST()
+{
+	freeMem(_msf_seqHits, (_msf_seqListSize) * sizeof(int));
+	freeMem(_msf_refGenName, 4*SEQ_LENGTH);
+	/*
+	  int i;
+	  for (i=0; i<_msf_rIndexSize; i++)
+	{
+		freeMem(_msf_rIndex[i].seqInfo, _msf_rIndex[i].seqInfo[0]+1);
+	}
+	freeMem(_msf_rIndex, _msf_rIndexSize);*/
+	freeMem(_msf_map_sort_seqList, sizeof(Pair) * _msf_seqListSize);
+	freeMem(_msf_sort_seqList, sizeof(int) * _msf_seqListSize);
+}
+/*
+Will apply the Levenshtein Dynamic programming.
+Different from verifySingleEndEditDistance fucntion
+	as in this fucntion only one dynamic table is made while
+	in verifySingleEndEditDistance two dynamic table is made
+	for each right and left string
+*/
+int editDistance(int refIndex, char *seq, int seqLength, char *matrix)
+{
+	int i = 0;
+	int size = 0;
+	int error = 0;
+	int rIndex = 0;
+int directionIndex = 0;
+	int min = 0;
+	int minIndex =0;
+	int tempUp = 0;
+	int tempDown = 0;
+	char *ref;
+	int errorString = 0;
+/*
+1: Up
+2: Side
+3: Diagnoal Match
+4: Diagnoal Mismatch
+*/
+	int upValue;
+	int diagValue;
+	int sideValue;
+	ref = _msf_refGen + refIndex - 1;
+	rIndex = 1;
+	for(i=0; i <= errThreshold; i++)
+	{
+		score[0][i] = i;
+		score[i][0] = i;
+	}
+	while(rIndex <= seqLength +errThreshold)
+	{
+		tempUp   = ((rIndex - errThreshold) > 0 ? ((rIndex > seqLength) ? seqLength - errThreshold :rIndex - errThreshold) : 1 );
+		tempDown =  ((rIndex >= seqLength-errThreshold ) ? seqLength+1  :rIndex + errThreshold + 1);
+		for(i = tempUp ; i < tempDown ; i++)
+		{
+			errorString = (*(ref+rIndex-1) == *(seq+i-1));
+			upValue   = score[i-1][rIndex]+1;
+			diagValue = score[i-1][rIndex-1]+ !errorString;
+			sideValue = score[i][rIndex-1]+1;
+			if(i != tempUp &&  i != tempDown-1)
+				score[i][rIndex] = min3(sideValue, diagValue , upValue);
+			else if(  (i == ((rIndex - errThreshold) > 0 ? rIndex - errThreshold : 1)) && rIndex <= seqLength  )
+				score[i][rIndex] = min(sideValue, diagValue);
+			else if(rIndex > seqLength &&  (i == seqLength - errThreshold) )
+				 score[i][rIndex] = sideValue;
+			else
+				score[i][rIndex] = min(diagValue , upValue);
+			if(i == tempUp)
+				error = score[i][rIndex];
+			else if(error > score[i][rIndex])
+error = score[i][rIndex];
+		}
+		rIndex++;
+	}
+	min  = score[seqLength][seqLength+errThreshold];
+	minIndex = seqLength + errThreshold;
+	// Find the Best error for all the possible ways.
+	for(i = 1; i <= 2*errThreshold; i++)
+	{
+		if(min >= score[seqLength][seqLength+errThreshold-i] && seqLength+errThreshold-i > 0)
+		{
+			min = score[seqLength][seqLength+errThreshold-i];
+			minIndex = seqLength+errThreshold-i;
+		}
+	}
+	error = score[seqLength][minIndex];
+	directionIndex = seqLength;
+	rIndex = minIndex;
+	while(directionIndex != 0 || rIndex != 0)
+{
+		if(rIndex == 0)
+		{
+			if(score[directionIndex][rIndex] - score[directionIndex-1][rIndex] == 1)
+{
+matrix[size] = *(seq+directionIndex-1);
+size++;
+matrix[size] = 'I';
+directionIndex--;
+}
+		}
+		else if(directionIndex == 0)
+		{
+			if(score[directionIndex][rIndex] - score[directionIndex][rIndex-1] == 1)
+{
+matrix[size] = *(ref+rIndex-1);
+size++;
+matrix[size] = 'D';
+rIndex--;
+}
+		}
+		else if(directionIndex-rIndex == errThreshold)
+		{
+			if(score[directionIndex][rIndex] - score[directionIndex-1][rIndex] == 1)
+			{
+				matrix[size] = *(seq+directionIndex-1);
+size++;
+matrix[size] = 'I';
+directionIndex--;
+			}
+			else if( score[directionIndex][rIndex] - score[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrix[size] = *(ref+rIndex-1);
+	                        rIndex--;
+	                        directionIndex--;
+			}
+			else
+			{
+				matrix[size] = 'M';
+	                        rIndex--;
+	                directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+		        if(score[directionIndex][rIndex] - score[directionIndex][rIndex-1] == 1)
+{
+				matrix[size] = *(ref+rIndex-1);
+size++;
+matrix[size] = 'D';
+rIndex--;
+}
+else if( score[directionIndex][rIndex] - score[directionIndex-1][rIndex-1] == 1 )
+{
+matrix[size] = *(ref+rIndex-1);
+rIndex--;
+directionIndex--;
+}
+else
+{
+matrix[size] = 'M';
+rIndex--;
+directionIndex--;
+}
+		}
+		else
+		{
+			if(score[directionIndex][rIndex] - score[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+{
+				matrix[size] = *(seq+directionIndex-1);
+size++;
+matrix[size] = 'I';
+directionIndex--;
+}
+			else if(score[directionIndex][rIndex] - score[directionIndex][rIndex-1] == 1 && rIndex != 0)
+{
+matrix[size] = *(ref+rIndex-1);
+size++;
+matrix[size] = 'D';
+rIndex--;
+}
+else if( score[directionIndex][rIndex] - score[directionIndex-1][rIndex-1] == 1 )
+{
+matrix[size] = *(ref+rIndex-1);
+rIndex--;
+directionIndex--;
+}
+else
+{
+matrix[size] = 'M';
+rIndex--;
+directionIndex--;
+}
+		}
+size++;
+}
+matrix[size] = '\0';
+	char returnString[200];
+	returnString[0] = '\0';
+	reverse(matrix, returnString, size);
+	sprintf(matrix, "%s", returnString);
+	return error;
+}
+/*
+	Will apply the Levenshtein Dynamic programming.
+	in both right and left direction as long as the
+	threshould error is reached or end of string length
+*/
+int msfHashVal(char *seq)
+{
+int i=0;
+int val=0, numericVal=0;
+while(i<6)
+{
+switch (seq[i])
+{
+case 'A':
+numericVal = 0;
+break;
+case 'C':
+numericVal = 1;
+break;
+case 'G' :
+numericVal = 2;
+break;
+case 'T':
+numericVal = 3;
+break;
+default:
+return -1;
+break;
+}
+val = (val << 2)|numericVal;
+i++;
+}
+return val;
+}
+int verifySingleEndEditDistance2(int refIndex, char *lSeq, int lSeqLength, char *rSeq, int rSeqLength, int segLength, char *matrix, int *map_location,  short *seqHashValue)
+{
+	int i = 0;
+	char *  ref;
+	char *  tempref;
+	int rIndex = 0;				//reference Index
+	int e = errThreshold;
+	int error        = 0;
+	int error1       = 0;
+	int totalError   = 0;
+	/*
+		1: Up
+		2: Side
+		3: Diagnoal Match
+		4: Diagnoal Mismatch
+	*/
+	int minIndex1 = 0;
+	int minIndex2 = 0;
+	int directionIndex = 0;
+	int size = 0;
+	int startIndex1 = 0;
+	rIndex = 1;
+	char matrixR[200];
+	char matrixL[200];
+	ref = _msf_refGen + refIndex - 1;
+	tempref = _msf_refGen + refIndex - 1;
+	int jumpIndex = 0;
+	if(rSeqLength != 0)
+	{
+		error1 = forwardEditDistance2SSE2(ref+segLength+jumpIndex, rSeqLength-jumpIndex, rSeq+jumpIndex, rSeqLength-jumpIndex);
+		if(error1 == -1)
+			return -1;
+	}
+	if(lSeqLength != 0)
+	{
+		error = backwardEditDistance2SSE2(ref-1, lSeqLength, lSeq+lSeqLength-1, lSeqLength);
+		if(error == -1)
+		{
+			return -1;
+		}
+	}
+	matrixL[0] = '\0';
+	matrixR[0] = '\0';
+	ref = _msf_refGen + refIndex - 1;
+	if(error1+error > errThreshold)
+		return -1;
+	ref = _msf_refGen + refIndex - 1;
+	rIndex = startIndex1+1;
+	int i0 = 0;
+	int i1 = 0;
+	int i2 = 0;
+	__m128i R0;
+	__m128i R1;
+	__m128i Side1, Side2,Side;                    //side matrix
+	__m128i Down1, Down2,Down;                    //down matrix
+	__m128i Diag;                   //
+	__m128i tmp;
+	/* initialize */
+	R0 = _mm_setzero_si128 ();
+	R1 = _mm_setzero_si128 ();
+	Diag = _mm_setzero_si128 ();
+	Side1 = _mm_setzero_si128 ();
+	Side2 = _mm_setzero_si128 ();
+	Down1 = _mm_setzero_si128 ();
+	Down2 = _mm_setzero_si128 ();
+Down = _mm_setzero_si128 ();
+Side = _mm_setzero_si128 ();
+	tmp = _mm_setzero_si128 ();
+	/* end initialize */
+	int mismatch[3] = {0,0,0};
+	if(lSeqLength != 0)
+	{
+		char *a;
+		char *b;
+		a = ref-1;
+		b = lSeq+lSeqLength-1;
+		R0 = _mm_insert_epi16(R0,0,0);
+		score[0][0] = 0;
+		R1 = _mm_insert_epi16(R1,1,0);
+		R1 = _mm_insert_epi16(R1,1,1);
+		score[1][0] = 1;
+		direction1[1][0] = 1;
+		score[0][1] = 1;
+		direction1[0][1] = 2;
+		mismatch[0] = ((a[0]) != (b[0]));
+		Diag = _mm_insert_epi16(Diag,2*e,0);
+		Diag = _mm_insert_epi16(Diag,mismatch[0],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Side1 = _mm_insert_epi16(Side1,1,1);
+		Side1 = _mm_insert_epi16(Side1,2*e,2);
+		Side2 = _mm_insert_epi16(Side2,2*e,0);
+		Side2 = _mm_insert_epi16(Side2,1,1);
+		Side2 = _mm_insert_epi16(Side2,1,2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,1);
+		Down1 = _mm_insert_epi16(Down1,2*e,2);
+		Down2 = _mm_insert_epi16(Down2,2*e,0);
+		Down2 = _mm_insert_epi16(Down2,1,1);
+		Down2 = _mm_insert_epi16(Down2,1,2);
+		tmp = _mm_slli_si128(R1,2);
+		R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+		R0 = _mm_min_epi16(R0,tmp+Down2);
+		i0 = _mm_extract_epi16(R0, 0);
+		i1 = _mm_extract_epi16(R0, 1);
+		i2 = _mm_extract_epi16(R0, 2);
+		score[0][2] = i0;
+		score[1][1] = i1;
+		score[2][0] = i2;
+		direction1[0][2] = 2;
+		direction1[1][1] = ((mismatch[0] == 0)? 3 : 4);
+		direction1[2][0] = 1;
+		for (i = 3; i  < 2*lSeqLength; i++)
+		{
+			if(i % 2 ==1)
+			{
+				Diag = _mm_sub_epi8(Diag, Diag);
+				mismatch[0] = ( *(a-((i+1)/2-1)) != *(b-((i-1)/2-1)) );
+				Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+				mismatch[1] = ( *(a-((i-1)/2-1)) != *(b-((i+1)/2-1)) );
+				Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+				tmp = _mm_srli_si128(R0,2);
+				R1 = _mm_min_epi16(tmp+Side1, R1+Diag);
+				R1 = _mm_min_epi16(R1,R0+Down1);
+				i0 = _mm_extract_epi16(R1, 0);
+				i1 = _mm_extract_epi16(R1, 1);
+				score[i/2][i/2+1] = i0;
+				score[i/2+1][i/2] = i1;
+				direction1[i/2][i/2+1] = (score[i/2][i/2+1]==score[i/2-1][i/2] && mismatch[0] == 0)  ? 3 :
+							 (score[i/2][i/2+1]-score[i/2-1][i/2+1]==1)? 1 :
+							 (score[i/2][i/2+1]-score[i/2][i/2]==1)    ? 2 : 4;
+				direction1[i/2+1][i/2] = (score[i/2+1][i/2]==score[i/2][i/2-1] && mismatch[1] == 0)  ? 3 :
+(score[i/2+1][i/2]-score[i/2][i/2]==1)    ? 1 :
+(score[i/2+1][i/2]-score[i/2+1][i/2-1]==1)? 2 : 4;
+				if(i > 2 * lSeqLength - 2)
+				{
+					error = min(error, i1);
+					if(error == i1)
+						minIndex1 = i-lSeqLength;
+				}
+			}
+			else if(i % 2 == 0)
+			{
+				mismatch[0] = ( *(a-(i/2))   != *(b-(i/2-2)) );
+				Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+				mismatch[1] = ( *(a-(i/2-1)) != *(b-(i/2-1)) );
+				Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+				mismatch[2] = ( *(a-(i/2-2)) != *(b-(i/2))   );
+				Diag = _mm_insert_epi16(Diag,mismatch[2],2);
+				tmp = _mm_slli_si128(R1,2);
+				R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+				R0 = _mm_min_epi16(R0,tmp+Down2);
+				i0 = _mm_extract_epi16(R0, 0);
+				i1 = _mm_extract_epi16(R0, 1);
+				i2 = _mm_extract_epi16(R0, 2);
+				score[i/2-1][i/2+1] = i0;
+				score[i/2][i/2] = i1;
+				score[i/2+1][i/2-1] = i2;
+				direction1[i/2-1][i/2+1] = (score[i/2-1][i/2+1]==score[i/2-2][i/2] && mismatch[0] == 0)  ? 3 : (score[i/2-1][i/2+1]-score[i/2-1][i/2]==1)  ? 2 : 4;
+direction1[i/2][i/2] =   (score[i/2][i/2]==score[i/2-1][i/2-1] && mismatch[1] == 0)  ? 3 :
+(score[i/2][i/2]-score[i/2-1][i/2]==1)    ? 1 :
+(score[i/2][i/2]-score[i/2][i/2-1]==1)    ? 2 : 4;
+				direction1[i/2+1][i/2-1] =  (score[i/2+1][i/2-1]==score[i/2][i/2-2] && mismatch[2] == 0)  ? 3 :
+(score[i/2+1][i/2-1]-score[i/2][i/2-1]==1) ? 1 : 4;
+				if( (i/2) % segLength == 0 && i1 == 0)	// the segment has been processed no need to process it again
+				{
+					return -1;
+				}
+				if(i == 2 * lSeqLength - 2)
+				{
+					error = i2;
+					minIndex1 = i-lSeqLength;
+				}
+			}
+		}
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		//fill the first part of the error
+		mismatch[0] = ( *(a-(i/2))   != *(b-(i/2-2)) );
+		Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+		mismatch[1] = ( *(a-(i/2-1)) !=*(b-(i/2-1)) );
+		Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+		R0 = _mm_min_epi16(R0,_mm_slli_si128(R1,2)+Down1);
+		i0 = _mm_extract_epi16(R0, 0);
+		i1 = _mm_extract_epi16(R0, 1);
+		error = min(error, i1);
+		if(error == i1)
+			minIndex1 = i-lSeqLength;
+		score[i/2-1][i/2+1] = i0;
+score[i/2][i/2] = i1;
+direction1[i/2-1][i/2+1] = (score[i/2-1][i/2+1]==score[i/2-2][i/2] && mismatch[0] == 0)  ? 3 :
+(score[i/2-1][i/2+1]-score[i/2-1][i/2])  ? 2 : 4;
+direction1[i/2][i/2] = (score[i/2][i/2]==score[i/2-1][i/2-1] && mismatch[1] == 0)  ? 3 :
+(score[i/2][i/2]-score[i/2-1][i/2]==1)    ? 1 :
+(score[i/2][i/2]-score[i/2][i/2-1]==1)? 2 : 4;
+		//fill the second part of the error
+		i++;
+		Diag = _mm_sub_epi8(Diag, Diag);
+		Diag = _mm_insert_epi16(Diag,2*e,0);
+		mismatch[0] = ( *(a-(i/2)) != *(b-(lSeqLength-1)) );
+		Diag = _mm_insert_epi16(Diag,mismatch[0],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+		R1 = _mm_min_epi16(R1,_mm_slli_si128(R0,2)+Down1);
+		i0 = _mm_extract_epi16(R1, 0);
+		i1 = _mm_extract_epi16(R1, 1);
+		error = min(error, i1);
+		if(error == i1)
+			minIndex1 = i-lSeqLength;
+		score[i/2-1][i/2+2] = i0;
+		score[i/2][i/2+1] = i1;
+direction1[i/2-1][i/2+2] = (score[i/2-1][i/2+2]==score[i/2-2][i/2+1] && mismatch[0] == 0)  ? 3 :
+(score[i/2-1][i/2+2]-score[i/2-1][i/2+1]==1)    ? 2 : 4;
+direction1[i/2][i/2+1] = (score[i/2][i/2+1]==score[i/2-1][i/2])  ? 3 :
+(score[i/2][i/2+1]-score[i/2-1][i/2+1]==1)? 1 :
+(score[i/2][i/2+1]-score[i/2][i/2]==1)? 2 : 4;
+		//fill the last the last element of the matrix
+		i++;
+		Diag = _mm_sub_epi8(Diag, Diag);
+		mismatch[0] = ( *(a-(i/2)) != *(b-(lSeqLength-1)) );
+		Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+		Down = _mm_insert_epi16(Down,1,0);
+		Side = _mm_insert_epi16(Side,1,0);
+		tmp = _mm_srli_si128(R1,2);
+		R0 = _mm_min_epi16(R1+Down, R0+Diag);
+		R0 = _mm_min_epi16(R0,tmp+Side);
+		i0 = _mm_extract_epi16(R0, 0);
+		error = min(error, i0);
+		if(error == 0)
+			return -1;
+		if(error == i0)
+			minIndex1 = i-lSeqLength;
+		if(mismatch[0] == 0)
+			direction1[lSeqLength][lSeqLength+errThreshold] = 3;
+		else
+		{
+			if(score[lSeqLength][lSeqLength+errThreshold] - score[lSeqLength][lSeqLength+errThreshold-1] == 1)
+				direction1[lSeqLength][lSeqLength+errThreshold] = 2;
+			else if(score[lSeqLength][lSeqLength+errThreshold] - score[lSeqLength-1][lSeqLength+errThreshold] == 1)
+				direction1[lSeqLength][lSeqLength+errThreshold] = 1;
+			else
+				direction1[lSeqLength][lSeqLength+errThreshold] = 4;
+		}
+	}
+	error1   = error;
+	error    = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	*map_location = ((lSeqLength == 0) ? refIndex : refIndex - rIndex) ;
+	ref = ref + segLength;
+	if(rSeqLength <= e)
+	{
+		char *a;
+		char *b;
+		int tmp_index = 0;
+		a = ref;
+		b = rSeq;
+		for(tmp_index = 0; tmp_index < rSeqLength; tmp_index++)
+		{
+			matrixR[tmp_index] = (a[tmp_index]==b[tmp_index]) ? 'M' : a[tmp_index]  ;
+		}
+		matrixR[tmp_index] = '\0';
+	}
+	else if(rSeqLength != 0 && rSeqLength >= e)
+	{
+		char *a;
+		char *b;
+		a = ref;
+		b = rSeq;
+		R0 = _mm_sub_epi8(R0, R0);
+		R1 = _mm_sub_epi8(R1, R1);
+		R0 = _mm_insert_epi16(R0,0,0);
+		score[0][0] = 0;
+		R1 = _mm_insert_epi16(R1,1,0);
+		R1 = _mm_insert_epi16(R1,1,1);
+		score[1][0] = 1;
+		direction2[1][0] = 1;
+		score[0][1] = 1;
+		direction2[0][1] = 2;
+		mismatch[0] = ((a[0]) != (b[0]));
+		Diag = _mm_insert_epi16(Diag,2*e,0);
+		Diag = _mm_insert_epi16(Diag,mismatch[0],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		Side1 = _mm_insert_epi16(Side1,1,0);
+		Side1 = _mm_insert_epi16(Side1,1,1);
+		Side1 = _mm_insert_epi16(Side1,2*e,2);
+		Side2 = _mm_insert_epi16(Side2,2*e,0);
+		Side2 = _mm_insert_epi16(Side2,1,1);
+		Side2 = _mm_insert_epi16(Side2,1,2);
+		Down1 = _mm_insert_epi16(Down1,1,0);
+		Down1 = _mm_insert_epi16(Down1,1,1);
+		Down1 = _mm_insert_epi16(Down1,2*e,2);
+		Down2 = _mm_insert_epi16(Down2,2*e,0);
+		Down2 = _mm_insert_epi16(Down2,1,1);
+		Down2 = _mm_insert_epi16(Down2,1,2);
+		tmp = _mm_slli_si128(R1,2);
+		R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+		R0 = _mm_min_epi16(R0,tmp+Down2);
+		i0 = _mm_extract_epi16(R0, 0);
+		i1 = _mm_extract_epi16(R0, 1);
+		i2 = _mm_extract_epi16(R0, 2);
+		score[0][2] = i0;
+		score[1][1] = i1;
+		score[2][0] = i2;
+		direction2[0][2] = 2;
+		direction2[1][1] = ((mismatch[0] == 0)? 3 : 4);
+		direction2[2][0] = 1;
+		for (i = 3; i  < 2*rSeqLength; i++)
+		{
+			if(i % 2 ==1)
+			{
+				Diag = _mm_sub_epi8(Diag, Diag);
+				mismatch[0] = ((a[(i+1)/2-1]) != (b[(i-1)/2-1]));
+				Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+				mismatch[1] = ((a[(i-1)/2-1]) != (b[(i+1)/2-1]));
+				Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+				tmp = _mm_srli_si128(R0,2);
+				R1 = _mm_min_epi16(tmp+Side1, R1+Diag);
+				R1 = _mm_min_epi16(R1,R0+Down1);
+				i0 = _mm_extract_epi16(R1, 0);
+				i1 = _mm_extract_epi16(R1, 1);
+				score[i/2][i/2+1] = i0;
+				score[i/2+1][i/2] = i1;
+				direction2[i/2][i/2+1] = (score[i/2][i/2+1]==score[i/2-1][i/2] && mismatch[0] == 0)  ? 3 :
+										 (score[i/2][i/2+1]-score[i/2-1][i/2+1]==1)? 1 :
+										 (score[i/2][i/2+1]-score[i/2][i/2]==1)    ? 2 : 4;
+				direction2[i/2+1][i/2] = (score[i/2+1][i/2]==score[i/2][i/2-1] && mismatch[1] == 0)  ? 3 :
+										 (score[i/2+1][i/2]-score[i/2][i/2]==1)    ? 1 :
+										 (score[i/2+1][i/2]-score[i/2+1][i/2-1]==1)? 2 : 4;
+				if(i > 2 * rSeqLength - 2)
+				{
+					error = min(error, i1);
+					if(error == i1)
+						minIndex2 = i-rSeqLength;
+				}
+			}
+			else if(i % 2 == 0)
+			{
+				mismatch[0] = ((a[i/2]) != (b[i/2-2]));
+				Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+				mismatch[1] = ((a[i/2-1]) != (b[i/2-1]));
+				Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+				mismatch[2] = ((a[i/2-2]) != (b[i/2]));
+				Diag = _mm_insert_epi16(Diag,mismatch[2],2);
+				tmp = _mm_slli_si128(R1,2);
+				R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+				R0 = _mm_min_epi16(R0,tmp+Down2);
+				i0 = _mm_extract_epi16(R0, 0);
+				i1 = _mm_extract_epi16(R0, 1);
+				i2 = _mm_extract_epi16(R0, 2);
+				score[i/2-1][i/2+1] = i0;
+				score[i/2][i/2] = i1;
+				score[i/2+1][i/2-1] = i2;
+				direction2[i/2-1][i/2+1] = (score[i/2-1][i/2+1]==score[i/2-2][i/2] && mismatch[0] == 0)  ? 3 :
+										   (score[i/2-1][i/2+1]-score[i/2-1][i/2]==1)  ? 2 : 4;
+				direction2[i/2][i/2] =   (score[i/2][i/2]==score[i/2-1][i/2-1] && mismatch[1] == 0)  ? 3 :
+										 (score[i/2][i/2]-score[i/2-1][i/2]==1)    ? 1 :
+										 (score[i/2][i/2]-score[i/2][i/2-1]==1)    ? 2 : 4;
+				direction2[i/2+1][i/2-1] =  (score[i/2+1][i/2-1]==score[i/2][i/2-2] && mismatch[2]==0)  ? 3 :
+											(score[i/2+1][i/2-1]-score[i/2][i/2-1]==1) ? 1 : 4;
+				if(i == 2 * rSeqLength - 2)
+				{
+					error = i2;
+					minIndex2 = i-rSeqLength;
+				}
+			}
+		}
+		Down1 = _mm_insert_epi16(Down1,2*e,0);
+		//fill the first part of the error
+		mismatch[0] = ((a[i/2]) != (b[i/2-2]));
+		Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+		mismatch[1] = ((a[i/2-1]) != (b[i/2-1]));
+		Diag = _mm_insert_epi16(Diag,mismatch[1],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		R0 = _mm_min_epi16(R1+Side1, R0+Diag);
+		R0 = _mm_min_epi16(R0,_mm_slli_si128(R1,2)+Down1);
+		i0 = _mm_extract_epi16(R0, 0);
+		i1 = _mm_extract_epi16(R0, 1);
+		error = min(error, i1);
+		if(error == i1)
+			minIndex2 = i-rSeqLength;
+		score[i/2-1][i/2+1] = i0;
+		score[i/2][i/2] = i1;
+direction2[i/2-1][i/2+1] = (score[i/2-1][i/2+1]==score[i/2-2][i/2] && mismatch[0] == 0)  ? 3 :
+										   (score[i/2-1][i/2+1]-score[i/2-1][i/2]==1)    ? 2 : 4;
+direction2[i/2][i/2] = (score[i/2][i/2]==score[i/2-1][i/2-1] && mismatch[1] == 0)  ? 3 :
+(score[i/2][i/2]-score[i/2-1][i/2]==1)    ? 1 :
+(score[i/2][i/2]-score[i/2][i/2-1]==1)? 2 : 4;
+		//fill the second part of the error
+		i++;
+		Diag = _mm_sub_epi8(Diag, Diag);
+		Diag = _mm_insert_epi16(Diag,2*e,0);
+		mismatch[0] = ((a[i/2]) != (b[rSeqLength-1]));
+		Diag = _mm_insert_epi16(Diag,mismatch[0],1);
+		Diag = _mm_insert_epi16(Diag,2*e,2);
+		R1 = _mm_min_epi16(R0+Side1, _mm_slli_si128(R1,2)+Diag);
+		R1 = _mm_min_epi16(R1,_mm_slli_si128(R0,2)+Down1);
+		i0 = _mm_extract_epi16(R1, 0);
+		i1 = _mm_extract_epi16(R1, 1);
+		error = min(error, i1);
+		if(error == i1)
+			minIndex2 = i-rSeqLength;
+		score[i/2-1][i/2+2] = i0;
+		score[i/2][i/2+1] = i1;
+direction2[i/2-1][i/2+2] = (score[i/2-1][i/2+2]==score[i/2-2][i/2+1] && mismatch[0] == 0)  ? 3 :
+(score[i/2-1][i/2+2]-score[i/2-1][i/2+1]==1)  ? 2 : 3;
+direction2[i/2][i/2+1] = (score[i/2][i/2+1]==score[i/2-1][i/2] && mismatch[0] == 0)  ? 3 :
+(score[i/2][i/2+1]-score[i/2-1][i/2+1]==1)? 1 :
+(score[i/2][i/2+1]-score[i/2][i/2]==1)? 2 : 4;
+		//fill the last the last element of the matrix
+		i++;
+		Diag = _mm_sub_epi8(Diag, Diag);
+		mismatch[0] = ((a[i/2]) != (b[rSeqLength-1]));
+		Diag = _mm_insert_epi16(Diag,mismatch[0],0);
+		Down = _mm_sub_epi8(Down, Down);
+		Down = _mm_insert_epi16(Down,1,0);
+		Side = _mm_sub_epi8(Side, Side);
+		Side = _mm_insert_epi16(Side,1,0);
+		tmp = _mm_srli_si128(R1,2);
+		R0 = _mm_min_epi16(R1+Down, R0+Diag);
+		R0 = _mm_min_epi16(R0,tmp+Side);
+		i0 = _mm_extract_epi16(R0, 0);
+		error = min(error, i0);
+		if(error == i0)
+			minIndex2 = i-rSeqLength;
+		if(mismatch[0] == 0)
+			direction2[rSeqLength][rSeqLength+errThreshold] = 3;
+		else
+		{
+			if(score[rSeqLength][rSeqLength+errThreshold] - score[rSeqLength][rSeqLength+errThreshold-1] == 1)
+				direction2[lSeqLength][lSeqLength+errThreshold] = 2;
+			else if(score[rSeqLength][rSeqLength+errThreshold] - score[rSeqLength-1][rSeqLength+errThreshold] == 1)
+				direction2[rSeqLength][rSeqLength+errThreshold] = 1;
+			else
+				direction2[rSeqLength][rSeqLength+errThreshold] = 4;
+		}
+	}
+	totalError = error1 + error;
+	size = 0;
+	directionIndex = rSeqLength;
+	rIndex = minIndex2;
+	if(rSeqLength > e)
+	{
+		while(directionIndex != 0 || rIndex != 0)
+		{
+			if(direction2[directionIndex][rIndex] == 3)
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+			else if(direction2[directionIndex][rIndex] == 4)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else if(direction2[directionIndex][rIndex] == 2)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			size++;
+		}
+		matrixR[size] = '\0';
+	}
+	size = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(direction1[directionIndex][rIndex] == 3)
+		{
+			matrixL[size] = 'M';
+			rIndex--;
+			directionIndex--;
+		}
+		else if(direction1[directionIndex][rIndex] == 4)
+		{
+			matrixL[size] = *(tempref-rIndex);
+			rIndex--;
+			directionIndex--;
+		}
+		else if(direction1[directionIndex][rIndex] == 2)
+		{
+			matrixL[size] = 'D';
+			size++;
+			matrixL[size] = *(tempref-rIndex);
+			rIndex--;
+		}
+		else
+		{
+			matrixL[size] = 'I';
+			size++;
+			matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+			directionIndex--;
+		}
+		size++;
+	}
+	matrixL[size] = '\0';
+	char middle[200];
+	middle[0] = '\0';
+	for(i = 0; i < segLength; i++)
+		middle[i] = 'M';
+	middle[segLength] = '\0';
+	char rmatrixR[200];
+	reverse(matrixR, rmatrixR, strlen(matrixR));
+	sprintf(matrix, "%s%s%s", matrixL, middle, rmatrixR);
+	return totalError;
+}
+int verifySingleEndEditDistance4(int refIndex, char *lSeq, int lSeqLength, char *rSeq, int rSeqLength, int segLength, char *matrix, int *map_location,  short *seqHashValue)
+{
+	int i = 0;
+	char *  ref;
+	char *  tempref;
+	int rIndex = 0;				//reference Index
+	int error        = 0;
+	int error1       = 0;
+	int error2   	 = 0;
+	int error3 	 = 0;
+	int totalError   = 0;
+	int errorSegment = 0;
+	int ERROR_BOUND = errThreshold;
+	/*
+		1: Up
+		2: Side
+		3: Diagnoal Match
+		4: Diagnoal Mismatch
+	*/
+	int min = 0;
+	int minIndex1 = 0;
+	int minIndex2 = 0;
+	int directionIndex = 0;
+	int size = 0;
+	ref = _msf_refGen + refIndex - 1;
+	tempref = _msf_refGen + refIndex - 1;
+	if(lSeqLength != 0)
+	{
+		error3 = backwardEditDistance4SSE2(ref-1, lSeqLength, lSeq+lSeqLength-1, lSeqLength);
+		if(error3 == -1 || error3 == 0){
+			return -1;
+		}
+	}
+	if(rSeqLength != 0)
+	{
+		error2 = forwardEditDistance4SSE2(ref+segLength, rSeqLength, rSeq, rSeqLength);
+		if(error2 == -1)
+			return -1;
+	}
+	if(error2 + error3 > errThreshold)
+		return -1;
+	rIndex = 1;
+	int prevError = 0;
+	int tempUp = 0;
+	int tempDown = 0;
+	int errorString = 0;
+	int upValue;
+	int diagValue;
+	int sideValue;
+	while(rIndex <= lSeqLength+errThreshold  && lSeqLength != 0)
+	{
+		tempUp   = ((rIndex - ERROR_BOUND) > 0 ? ((rIndex > lSeqLength) ? lSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1 );
+		tempDown =  ((rIndex >= lSeqLength-ERROR_BOUND ) ? lSeqLength+1  :rIndex + ERROR_BOUND + 1);
+		for(i = tempUp ; i < tempDown ; i++)
+		{
+			errorString = (*(ref-rIndex) == *(lSeq+lSeqLength-i));
+			upValue   = scoreB[i-1][rIndex]+1;
+			diagValue = scoreB[i-1][rIndex-1]+ !errorString;
+			sideValue = scoreB[i][rIndex-1]+1;
+			if(i != tempUp &&  i != tempDown-1)
+				scoreB[i][rIndex] = min3(sideValue, diagValue , upValue);
+			else if(  (i == ((rIndex - ERROR_BOUND) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= lSeqLength  )
+				scoreB[i][rIndex] = min(sideValue, diagValue);
+			else if(rIndex > lSeqLength &&  (i == lSeqLength - ERROR_BOUND) )
+				 scoreB[i][rIndex] = sideValue;
+			else
+				scoreB[i][rIndex] = min(diagValue , upValue);
+			if(i == tempUp)
+				error = scoreB[i][rIndex];
+			else if(error > scoreB[i][rIndex])
+error = scoreB[i][rIndex];
+		}
+		if(rIndex <= lSeqLength)
+		{
+			errorSegment = error-prevError;
+		}
+		rIndex++;
+	}
+	if(lSeqLength != 0)
+	{
+		min  = scoreB[lSeqLength][lSeqLength+errThreshold];
+		minIndex1 = lSeqLength + errThreshold;
+		// Find the Best error for all the possible ways.
+		for(i = 1; i <= 2*errThreshold; i++)
+		{
+			if(min >= scoreB[lSeqLength][lSeqLength+errThreshold-i] && lSeqLength+errThreshold-i > 0)
+			{
+				min = scoreB[lSeqLength][lSeqLength+errThreshold-i];
+				minIndex1 = lSeqLength+errThreshold-i;
+			}
+		}
+		error = scoreB[lSeqLength][minIndex1];
+	}
+	error1   = error;
+	error    = 0;
+	errorSegment = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	*map_location = ((lSeqLength == 0) ? refIndex : refIndex - rIndex) ;
+	ref = ref + segLength;
+	if(rSeqLength != 0)
+	{
+		rIndex = 1;
+		while(rIndex <= rSeqLength+errThreshold-error1)
+		{
+			tempUp   = (rIndex - ERROR_BOUND) > 0 ? ((rIndex > rSeqLength) ? rSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1;
+			tempDown = ((rIndex >= rSeqLength- ERROR_BOUND ) ? rSeqLength+1  :rIndex + ERROR_BOUND + 1);
+			for(i = tempUp; i < tempDown ; i++)
+			{
+				errorString = (*(ref+rIndex-1) == *(rSeq+i-1));
+				upValue   = scoreF[i-1][rIndex]+1;
+				diagValue = scoreF[i-1][rIndex-1]+ !errorString;
+				sideValue = scoreF[i][rIndex-1]+1;
+				if(i != tempUp &&  i != tempDown-1)
+					scoreF[i][rIndex] = min3(sideValue, diagValue , upValue);
+				else if(  (i == ((rIndex - ERROR_BOUND ) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= rSeqLength  )
+					scoreF[i][rIndex] = min(sideValue, diagValue);
+				else if(rIndex > rSeqLength &&  (i == rSeqLength - ERROR_BOUND) )
+					scoreF[i][rIndex] = sideValue;
+				else
+					scoreF[i][rIndex] = min(diagValue , upValue);
+				if(i == tempUp)
+					error = scoreF[i][rIndex];
+				if(error > scoreF[i][rIndex])
+					error = scoreF[i][rIndex];
+			}
+			if(rIndex <= rSeqLength)
+			{
+				errorSegment = error;
+			}
+			rIndex++;
+		}
+		min  = scoreF[rSeqLength][rSeqLength+errThreshold-error1];
+		minIndex2 = rSeqLength + errThreshold-error1;
+		// Find the Best error for all the possible ways.
+		for(i = 1; i <= 2*(errThreshold-error1); i++)
+		{
+			if(min > scoreF[rSeqLength][rSeqLength+errThreshold-error1-i] && rSeqLength+errThreshold-error1-i > 0)
+			{
+				min = scoreF[rSeqLength][rSeqLength+errThreshold-error1-i];
+				minIndex2 = rSeqLength+errThreshold-error1-i;
+			}
+		}
+		error = scoreF[rSeqLength][minIndex2];
+	}
+	totalError = error + error1;
+	if(errThreshold > 4)
+		printf("ERROR in errorThreshold.\n");
+	if(totalError != error2 + error3 && totalError > errThreshold)
+	{
+		printf("ErrorF=%d, ErrorB=%d Error=%d Error=%d\n", error2,error3,error1,error);
+		scanf("%d", &i);
+	}
+	char matrixR[200];
+	char matrixL[200];
+	matrixR[0] = '\0';
+	matrixL[0] = '\0';
+	size = 0;
+	directionIndex = rSeqLength;
+	rIndex = minIndex2;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixR[size] = '\0';
+	size = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixL[size] = '\0';
+	char middle[200];
+	middle[0] = '\0';
+	for(i = 0; i < segLength; i++)
+		middle[i] = 'M';
+	middle[segLength] = '\0';
+	char rmatrixR[200];
+	reverse(matrixR, rmatrixR, strlen(matrixR));
+	sprintf(matrix, "%s%s%s", matrixL, middle, rmatrixR);
+	return totalError;
+}
+int verifySingleEndEditDistanceExtention(int refIndex, char *lSeq, int lSeqLength, char *rSeq, int rSeqLength, int segLength,
+					char *matrix, int *map_location, short *seqHashValue)
+{
+	int i = 0;
+	char *  ref;
+	char *  tempref;
+	int rIndex = 0;				//reference Index
+	int error        = 0;
+	int error1       = 0;
+	int error2   	 = 0;
+	int error3 	 = 0;
+	int totalError   = 0;
+	int errorSegment = 0;
+	int ERROR_BOUND = min(4, errThreshold);
+	/*
+		1: Up
+		2: Side
+		3: Diagnoal Match
+		4: Diagnoal Mismatch
+	*/
+	int min = 0;
+	int minIndex1 = 0;
+	int minIndex2 = 0;
+	int directionIndex = 0;
+	int size = 0;
+	ref = _msf_refGen + refIndex - 1;
+	tempref = _msf_refGen + refIndex - 1;
+	if(lSeqLength != 0)
+	{
+		error3 = backwardEditDistanceSSE2Extention(ref-1, lSeqLength, lSeq+lSeqLength-1, lSeqLength);
+		if(error3 == -1){
+			return -1;
+		}
+	}
+	if(rSeqLength != 0)
+	{
+		error2 = forwardEditDistanceSSE2Extention(ref+segLength, rSeqLength, rSeq, rSeqLength);
+		if(error2 == -1)
+			return -1;
+	}
+	if(error2 + error3 > errThreshold)
+		return -1;
+	rIndex = 1;
+	int prevError = 0;
+	int tempUp = 0;
+	int tempDown = 0;
+	int errorString = 0;
+	int upValue;
+	int diagValue;
+	int sideValue;
+	if(lSeqLength > ERROR_BOUND)
+	{
+		while(rIndex <= lSeqLength+ERROR_BOUND  && lSeqLength != 0)
+		{
+			tempUp   = ((rIndex - ERROR_BOUND) > 0 ? ((rIndex > lSeqLength) ? lSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1 );
+			tempDown =  ((rIndex >= lSeqLength-ERROR_BOUND ) ? lSeqLength+1  :rIndex + ERROR_BOUND + 1);
+			for(i = tempUp ; i < tempDown ; i++)
+			{
+				errorString = (*(ref-rIndex) == *(lSeq+lSeqLength-i));
+				upValue   = scoreB[i-1][rIndex]+1;
+				diagValue = scoreB[i-1][rIndex-1]+ !errorString;
+				sideValue = scoreB[i][rIndex-1]+1;
+				if(i != tempUp &&  i != tempDown-1)
+					scoreB[i][rIndex] = min3(sideValue, diagValue , upValue);
+				else if(  (i == ((rIndex - ERROR_BOUND) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= lSeqLength  )
+					scoreB[i][rIndex] = min(sideValue, diagValue);
+				else if(rIndex > lSeqLength &&  (i == lSeqLength - ERROR_BOUND) )
+					scoreB[i][rIndex] = sideValue;
+				else
+					scoreB[i][rIndex] = min(diagValue , upValue);
+				if(i == tempUp)
+					error = scoreB[i][rIndex];
+				else if(error > scoreB[i][rIndex])
+					error = scoreB[i][rIndex];
+			}
+			if(rIndex <= lSeqLength)
+			{
+				errorSegment = error-prevError;
+			}
+			rIndex++;
+		}
+		if(lSeqLength != 0)
+		{
+			min  = scoreB[lSeqLength][lSeqLength+ERROR_BOUND];
+			minIndex1 = lSeqLength + ERROR_BOUND;
+			// Find the Best error for all the possible ways.
+			for(i = 1; i <= 2*ERROR_BOUND; i++)
+			{
+				if(min >= scoreB[lSeqLength][lSeqLength+ERROR_BOUND-i] && lSeqLength+ERROR_BOUND-i > 0)
+				{
+					min = scoreB[lSeqLength][lSeqLength+ERROR_BOUND-i];
+					minIndex1 = lSeqLength+ERROR_BOUND-i;
+				}
+			}
+			error = scoreB[lSeqLength][minIndex1];
+		}
+	}
+	else
+	{
+		int j = 0;
+		for(i = 1; i <= lSeqLength; i++)
+		{
+			for(j = 1; j <= lSeqLength; j++)
+			{
+				scoreB[i][j] = min3(scoreB[i-1][j-1]+ (*(ref-j) != *(lSeq+lSeqLength-i) ),scoreB[i][j-1]+1 ,scoreB[i-1][j]+1);
+			}
+		}
+		error = scoreB[lSeqLength][lSeqLength];
+		minIndex1 = lSeqLength;
+	}
+	error1   = error;
+	error    = 0;
+	errorSegment = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	*map_location = ((lSeqLength == 0) ? refIndex : refIndex - rIndex) ;
+	ref = ref + segLength;
+	if(rSeqLength != 0 && rSeqLength > ERROR_BOUND)
+	{
+		ERROR_BOUND = min(ERROR_BOUND, rSeqLength);
+		if(rSeqLength == ERROR_BOUND)
+		{
+			for(i=0; i < 2*ERROR_BOUND; i++)
+				scoreF[0][i] = i;
+		}
+		rIndex = 1;
+		while(rIndex <= rSeqLength+ERROR_BOUND)
+		{
+			tempUp   = (rIndex - ERROR_BOUND) > 0 ? ((rIndex > rSeqLength) ? rSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1;
+			tempDown = ((rIndex >= rSeqLength- ERROR_BOUND ) ? rSeqLength+1  :rIndex + ERROR_BOUND + 1);
+			for(i = tempUp; i < tempDown ; i++)
+			{
+				errorString = (*(ref+rIndex-1) == *(rSeq+i-1));
+				upValue   = scoreF[i-1][rIndex]+1;
+				diagValue = scoreF[i-1][rIndex-1]+ !errorString;
+				sideValue = scoreF[i][rIndex-1]+1;
+				if(i != tempUp &&  i != tempDown-1)
+					scoreF[i][rIndex] = min3(sideValue, diagValue , upValue);
+				else if(  (i == ((rIndex - ERROR_BOUND ) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= rSeqLength  )
+					scoreF[i][rIndex] = min(sideValue, diagValue);
+				else if(rIndex > rSeqLength &&  (i == rSeqLength - ERROR_BOUND) )
+					scoreF[i][rIndex] = sideValue;
+				else
+					scoreF[i][rIndex] = min(diagValue , upValue);
+				if(i == tempUp)
+					error = scoreF[i][rIndex];
+				if(error > scoreF[i][rIndex])
+					error = scoreF[i][rIndex];
+			}
+			if(rIndex <= rSeqLength)
+			{
+				errorSegment = error;
+			}
+			rIndex++;
+		}
+		min  = scoreF[rSeqLength][rSeqLength+ERROR_BOUND];
+		minIndex2 = rSeqLength + ERROR_BOUND;
+		// Find the Best error for all the possible ways.
+		for(i = 1; i <= 2*ERROR_BOUND; i++)
+		{
+			if(min > scoreF[rSeqLength][rSeqLength+ERROR_BOUND-i] && rSeqLength+ERROR_BOUND-i > 0)
+			{
+				min = scoreF[rSeqLength][rSeqLength+ERROR_BOUND-i];
+				minIndex2 = rSeqLength+ERROR_BOUND-i;
+			}
+		}
+		error = scoreF[rSeqLength][minIndex2];
+	}
+	else
+	{
+		int j = 0;
+		for(i = 1; i <= rSeqLength; i++)
+		{
+			for(j = 1; j <= rSeqLength; j++)
+			{
+				scoreF[i][j] = min3(scoreF[i-1][j-1]+ (*(ref+j-1) != *(rSeq+i-1) ),scoreF[i][j-1]+1 ,scoreF[i-1][j]+1);
+			}
+		}
+		error = scoreF[rSeqLength][rSeqLength];
+		minIndex2 = rSeqLength;
+	}
+totalError = error + error1;
+	if(totalError  != error2+error3)
+	{
+		for(i = 0; i < lSeqLength; i++)
+			printf("%c", *(tempref-1-i));
+		printf("\n");
+		for(i = 0; i < lSeqLength; i++)
+			printf("%c", *(lSeq+i));
+		printf("\n");
+		for(i = 0; i < rSeqLength; i++)
+			printf("%c", *(tempref+segLength+i));
+		printf("\n");
+		for(i = 0; i < rSeqLength; i++)
+			printf("%c", *(rSeq+i));
+		printf("\n");
+		printf("ERROR=%d\n", totalError);
+		printf("ERROR_SSE=%d\n", error3+error2);
+		printf("ERROR_SSE_back=%d E_SSE_forw=%d\n", error3, error2);
+		printf("ERROR_back=%d E_forw=%d\n", error1, error);
+	}
+	char matrixR[200];
+	char matrixL[200];
+	matrixR[0] = '\0';
+	matrixL[0] = '\0';
+	size = 0;
+	directionIndex = rSeqLength;
+	rIndex = minIndex2;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixR[size] = '\0';
+	size = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixL[size] = '\0';
+	char middle[200];
+	middle[0] = '\0';
+	for(i = 0; i < segLength; i++)
+		middle[i] = 'M';
+	middle[segLength] = '\0';
+	char rmatrixR[200];
+	reverse(matrixR, rmatrixR, strlen(matrixR));
+	sprintf(matrix, "%s%s%s", matrixL, middle, rmatrixR);
+	return totalError;
+}
+int verifySingleEndEditDistance(int refIndex, char *lSeq, int lSeqLength, char *rSeq, int rSeqLength, int segLength, char *matrix, int *map_location, short *seqHashValue)
+{
+	int i = 0;
+	char *  ref;
+	char *  tempref;
+	int rIndex = 0;				//reference Index
+	int error        = 0;
+	int error1       = 0;
+	int error2 	 = 0;
+	int error3	 = 0;
+	int totalError   = 0;
+	int errorSegment = 0;
+	int ERROR_BOUND = errThreshold;
+	/*
+		1: Up
+		2: Side
+		3: Diagnoal Match
+		4: Diagnoal Mismatch
+	*/
+	int min = 0;
+	int minIndex1 = 0;
+	int minIndex2 = 0;
+	int directionIndex = 0;
+	int size = 0;
+	ref = _msf_refGen + refIndex - 1;
+	tempref = _msf_refGen + refIndex - 1;
+	if(rSeqLength != 0)
+	{
+		if(errThreshold %2 == 1)
+			error2 = forwardEditDistanceSSE2Odd(ref+segLength, rSeqLength, rSeq, rSeqLength);
+		else
+			error2 = forwardEditDistanceSSE2G(ref+segLength, rSeqLength, rSeq, rSeqLength);
+		if(error2 == -1)
+			return -1;
+	}
+	if(lSeqLength != 0)
+	{
+		if(errThreshold % 2 == 1)
+			error3 = backwardEditDistanceSSE2Odd(ref-1, lSeqLength, lSeq+lSeqLength-1, lSeqLength);
+		else
+			error3 = backwardEditDistanceSSE2G(ref-1, lSeqLength, lSeq+lSeqLength-1, lSeqLength);
+		if(error3 == -1 || error3 == 0){
+			return -1;
+		}
+	}
+	if(error3 + error2 > errThreshold)
+		return -1;
+	for(i = 0 ; i < errThreshold + 1; i++)
+	{
+		scoreB[0][i] = i;
+		scoreB[i][0] = i;
+	}
+	rIndex = 1;
+	int prevError = 0;
+	int tempUp = 0;
+	int tempDown = 0;
+	int errorString = 0;
+	int upValue;
+	int diagValue;
+	int sideValue;
+	while(rIndex <= lSeqLength+errThreshold  && lSeqLength != 0)
+	{
+		tempUp   = ((rIndex - ERROR_BOUND) > 0 ? ((rIndex > lSeqLength) ? lSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1 );
+		tempDown =  ((rIndex >= lSeqLength-ERROR_BOUND ) ? lSeqLength+1  :rIndex + ERROR_BOUND + 1);
+		for(i = tempUp ; i < tempDown ; i++)
+		{
+			errorString = (*(ref-rIndex) == *(lSeq+lSeqLength-i));
+			upValue   = scoreB[i-1][rIndex]+1;
+			diagValue = scoreB[i-1][rIndex-1]+ !errorString;
+			sideValue = scoreB[i][rIndex-1]+1;
+			if(i != tempUp &&  i != tempDown-1)
+				scoreB[i][rIndex] = min3(sideValue, diagValue , upValue);
+			else if(  (i == ((rIndex - ERROR_BOUND) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= lSeqLength  )
+				scoreB[i][rIndex] = min(sideValue, diagValue);
+			else if(rIndex > lSeqLength &&  (i == lSeqLength - ERROR_BOUND) )
+				 scoreB[i][rIndex] = sideValue;
+			else
+				scoreB[i][rIndex] = min(diagValue , upValue);
+			if(i == tempUp)
+				error = scoreB[i][rIndex];
+			else if(error > scoreB[i][rIndex])
+error = scoreB[i][rIndex];
+		}
+		if(rIndex <= lSeqLength)
+		{
+			errorSegment = error-prevError;
+		}
+		rIndex++;
+	}
+	if(lSeqLength != 0)
+	{
+		min  = scoreB[lSeqLength][lSeqLength+errThreshold];
+		minIndex1 = lSeqLength + errThreshold;
+		// Find the Best error for all the possible ways.
+		for(i = 1; i <= 2*errThreshold; i++)
+		{
+			if(min >= scoreB[lSeqLength][lSeqLength+errThreshold-i] && lSeqLength+errThreshold-i > 0)
+			{
+				min = scoreB[lSeqLength][lSeqLength+errThreshold-i];
+				minIndex1 = lSeqLength+errThreshold-i;
+			}
+		}
+		error = scoreB[lSeqLength][minIndex1];
+	}
+	error1   = error;
+	error    = 0;
+	errorSegment = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	*map_location = ((lSeqLength == 0) ? refIndex : refIndex - rIndex) ;
+	ref = ref + segLength;
+	if(rSeqLength != 0)
+	{
+		for(i = 0 ; i < errThreshold + 1; i++)
+		{
+			scoreF[0][i] = i;
+			scoreF[i][0] = i;
+		}
+		rIndex = 1;
+		while(rIndex <= rSeqLength+errThreshold-error1)
+		{
+			tempUp   = (rIndex - ERROR_BOUND) > 0 ? ((rIndex > rSeqLength) ? rSeqLength - ERROR_BOUND :rIndex - ERROR_BOUND) : 1;
+			tempDown = ((rIndex >= rSeqLength- ERROR_BOUND ) ? rSeqLength+1  :rIndex + ERROR_BOUND + 1);
+			for(i = tempUp; i < tempDown ; i++)
+			{
+				errorString = (*(ref+rIndex-1) == *(rSeq+i-1));
+				upValue   = scoreF[i-1][rIndex]+1;
+				diagValue = scoreF[i-1][rIndex-1]+ !errorString;
+				sideValue = scoreF[i][rIndex-1]+1;
+				if(i != tempUp &&  i != tempDown-1)
+					scoreF[i][rIndex] = min3(sideValue, diagValue , upValue);
+				else if(  (i == ((rIndex - ERROR_BOUND ) > 0 ? rIndex - ERROR_BOUND : 1)) && rIndex <= rSeqLength  )
+					scoreF[i][rIndex] = min(sideValue, diagValue);
+				else if(rIndex > rSeqLength &&  (i == rSeqLength - ERROR_BOUND) )
+					scoreF[i][rIndex] = sideValue;
+				else
+					scoreF[i][rIndex] = min(diagValue , upValue);
+				if(i == tempUp)
+					error = scoreF[i][rIndex];
+				if(error > scoreF[i][rIndex])
+					error = scoreF[i][rIndex];
+			}
+			if(rIndex <= rSeqLength)
+			{
+				errorSegment = error;
+			}
+			rIndex++;
+		}
+		min  = scoreF[rSeqLength][rSeqLength+errThreshold-error1];
+		minIndex2 = rSeqLength + errThreshold-error1;
+		// Find the Best error for all the possible ways.
+		for(i = 1; i <= 2*(errThreshold-error1); i++)
+		{
+			if(min > scoreF[rSeqLength][rSeqLength+errThreshold-error1-i] && rSeqLength+errThreshold-error1-i > 0)
+			{
+				min = scoreF[rSeqLength][rSeqLength+errThreshold-error1-i];
+				minIndex2 = rSeqLength+errThreshold-error1-i;
+			}
+		}
+		error = scoreF[rSeqLength][minIndex2];
+	}
+	totalError = error + error1;
+	if(totalError != error2 + error3 && totalError > errThreshold)
+	{
+		for(i = 0; i < lSeqLength; i++)
+printf("%c", *(tempref-1-i));
+printf("\n");
+for(i = 0; i < lSeqLength; i++)
+printf("%c", *(lSeq+i));
+printf("\n");
+for(i = 0; i < rSeqLength; i++)
+printf("%c", *(tempref+segLength+i));
+printf("\n");
+for(i = 0; i < rSeqLength; i++)
+printf("%c", *(rSeq+i));
+printf("\n");
+		printf("SSEF=%d SSEB%d\n", error2, error3);
+		printf("F=%d B=%d\n", error, error1);
+		scanf("%d", &i);
+	}
+	char matrixR[200];
+	char matrixL[200];
+	matrixR[0] = '\0';
+	matrixL[0] = '\0';
+	size = 0;
+	directionIndex = rSeqLength;
+	rIndex = minIndex2;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixR[size] = *(rSeq+directionIndex-1);
+				size++;
+				matrixR[size] = 'I';
+				directionIndex--;
+			}
+			else if(scoreF[directionIndex][rIndex] - scoreF[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				size++;
+				matrixR[size] = 'D';
+				rIndex--;
+			}
+			else if( scoreF[directionIndex][rIndex] - scoreF[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixR[size] = *(ref+rIndex-1);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixR[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixR[size] = '\0';
+	size = 0;
+	directionIndex = lSeqLength;
+	rIndex = minIndex1;
+	while(directionIndex != 0 || rIndex != 0)
+	{
+		if(directionIndex-rIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else if(rIndex - directionIndex == errThreshold)
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		else
+		{
+			if(scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex] == 1 && directionIndex != 0)
+			{
+				matrixL[size] = 'I';
+				size++;
+				matrixL[size] = *(lSeq+lSeqLength-directionIndex);
+				directionIndex--;
+			}
+			else if(scoreB[directionIndex][rIndex] - scoreB[directionIndex][rIndex-1] == 1 && rIndex != 0)
+			{
+				matrixL[size] = 'D';
+				size++;
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+			}
+			else if( scoreB[directionIndex][rIndex] - scoreB[directionIndex-1][rIndex-1] == 1 )
+			{
+				matrixL[size] = *(tempref-rIndex);
+				rIndex--;
+				directionIndex--;
+			}
+			else
+			{
+				matrixL[size] = 'M';
+				rIndex--;
+				directionIndex--;
+			}
+		}
+		size++;
+	}
+	matrixL[size] = '\0';
+	char middle[200];
+	middle[0] = '\0';
+	for(i = 0; i < segLength; i++)
+		middle[i] = 'M';
+	middle[segLength] = '\0';
+	char rmatrixR[200];
+	reverse(matrixR, rmatrixR, strlen(matrixR));
+	sprintf(matrix, "%s%s%s", matrixL, middle, rmatrixR);
+	return totalError;
+}
+int addCigarSize(int cnt){
+if (cnt<10) return 1;
+else if (cnt < 100) return 2;
+return 3;
+}
+/*
+	Generate Cigar from the back tracking matrix
+*/
+void generateCigar(char *matrix, int matrixLength, char *cigar)
+{
+	int i = 0;
+	int counterM=0;
+	int counterI=0;
+	int counterD=0;
+	int cigarSize = 0;
+	cigar[0] = '\0';
+	while(i < matrixLength)
+{
+if(matrix[i]=='M')
+{
+counterM++;
+if(counterI != 0)
+{
+			        sprintf(cigar, "%s%dI", cigar, counterI);
+				cigarSize += addCigarSize(counterI) + 1;
+				cigar[cigarSize] = '\0';
+				counterI=0;
+}
+			else if(counterD != 0)
+			{
+				sprintf(cigar, "%s%dD", cigar, counterD);
+				cigarSize += addCigarSize(counterD) + 1;
+				cigar[cigarSize] = '\0';
+counterD=0;
+			}
+}
+else if(matrix[i] == 'I')
+{
+if(counterM != 0)
+{
+sprintf(cigar, "%s%dM", cigar, counterM);
+				cigarSize += addCigarSize(counterM) + 1;
+				cigar[cigarSize] = '\0';
+counterM = 0;
+}
+else if(counterD != 0)
+{
+	sprintf(cigar, "%s%dD", cigar, counterD);
+				cigarSize += addCigarSize(counterD) + 1;
+				cigar[cigarSize] = '\0';
+counterD=0;
+}
+			counterI++;
+			i++;
+}
+		else if (matrix[i] == 'D')
+		{
+			if(counterM != 0)
+{
+sprintf(cigar, "%s%dM", cigar, counterM);
+				cigarSize += addCigarSize(counterM) + 1;
+				cigar[cigarSize] = '\0';
+counterM = 0;
+}
+else if(counterI != 0)
+{
+sprintf(cigar, "%s%dI", cigar, counterI);
+				cigarSize += addCigarSize(counterI) + 1;
+				cigar[cigarSize] = '\0';
+counterI=0;
+}
+			counterD++;
+			i++;
+		}
+		else
+		{
+			counterM++;
+if(counterI != 0)
+{
+sprintf(cigar, "%s%dI", cigar, counterI);
+				cigarSize += addCigarSize(counterI) + 1;
+				cigar[cigarSize] = '\0';
+counterI=0;
+}
+else if(counterD != 0)
+{
+sprintf(cigar, "%s%dD", cigar, counterD);
+				cigarSize += addCigarSize(counterD) + 1;
+				cigar[cigarSize] = '\0';
+counterD=0;
+}
+		}
+		i++;
+	}
+	if(counterM != 0)
+{
+sprintf(cigar, "%s%dM", cigar, counterM);
+		cigarSize += addCigarSize(counterM) + 1;
+		cigar[cigarSize] = '\0';
+counterM = 0;
+}
+else if(counterI != 0)
+{
+sprintf(cigar, "%s%dI", cigar, counterI);
+		cigarSize += addCigarSize(counterI) + 1;
+		cigar[cigarSize] = '\0';
+counterI = 0;
+}
+	else if(counterD != 0)
+	{
+		sprintf(cigar, "%s%dD", cigar, counterD);
+		cigarSize += addCigarSize(counterD) + 1;
+		cigar[cigarSize] = '\0';
+counterD = 0;
+	}
+	cigar[cigarSize] = '\0';
+}
+/*
+	Creates the Cigar output from the mismatching positions format  [0-9]+(([ACTGN]|\^[ACTGN]+)[0-9]+)*
+*/
+void generateCigarFromMD(char *mismatch, int mismatchLength, char *cigar)
+{
+	int i = 0;
+	int j = 0;
+	int start = 0;
+	int cigarSize = 0;
+cigar[0] = '\0';
+	while(i < mismatchLength)
+	{
+		if(mismatch[i] >= '0' && mismatch[i] <= '9')
+		{
+			start = i;
+			while(mismatch[i] >= '0' && mismatch[i] <= '9' && i < mismatchLength)
+				i++;
+			int value = atoi(mismatch+start);
+			for(j = 0;  j < value-1; j++)
+			{
+				cigar[cigarSize] = 'M';
+				cigarSize++;
+			}
+			cigar[cigarSize] = 'M';
+		}
+		else if(mismatch[i] == '^')
+		{
+			cigar[cigarSize] = 'I';
+			i++;
+		}
+		else if(mismatch[i] == '\'')
+		{
+			cigar[cigarSize] = 'D';
+			i++;
+		}
+		else
+		{
+			 cigar[cigarSize] = 'M';
+cigarSize++;
+		}
+		cigarSize++;
+		i++;
+	}
+	cigar[cigarSize] = '\0';
+}
+void generateSNPSAM(char *matrix, int matrixLength, char *outputSNP)
+{
+	int i = 0;
+	int counterM = 0;
+	int counterD = 0;
+	char delete[100];
+	int snpSize = 0;
+	outputSNP[0] = '\0';
+	delete[0] = '\0';
+	while(i < matrixLength)
+	{
+		if(matrix[i]=='M')
+		{
+			counterM++;
+			if(counterD != 0)
+			{
+				delete[counterD] = '\0';
+				counterD=0;
+sprintf(outputSNP, "%s^%s", outputSNP,delete);
+				snpSize += strlen(delete) + 1;
+				outputSNP[snpSize] = '\0';
+				delete[0] = '\0';
+			}
+		}
+		else if(matrix[i] == 'D')
+		{
+			if(counterM != 0)
+			{
+				sprintf(outputSNP, "%s%d", outputSNP, counterM);
+				snpSize += addCigarSize(counterM);
+				outputSNP[snpSize] = '\0';
+				counterM=0;
+				delete[counterD] = matrix[i+1];
+				i++;
+				counterD++;
+			}
+			else if(counterD != 0)
+			{
+				delete[counterD] = matrix[i+1];
+				counterD++;
+				i++;
+			}
+			else
+			{
+				delete[counterD] = matrix[i+1];
+				counterD++;
+i++;
+			}
+		}
+		else if(matrix[i] == 'I')
+		{
+			if(counterM != 0)
+{
+			  // sprintf(outputSNP, "%s%d\0", outputSNP, counterM);
+			  //counterM++;
+}
+else if(counterD != 0)
+{
+				delete[counterD] = '\0';
+				sprintf(outputSNP, "%s^%s", outputSNP, delete);
+				snpSize += strlen(delete) + 1;
+				outputSNP[snpSize] = '\0';
+				counterD = 0;
+				delete[0] = '\0';
+}
+			i++;
+		}
+		else
+		{
+			if(counterM != 0)
+			{
+				sprintf(outputSNP, "%s%d", outputSNP, counterM);
+				snpSize += addCigarSize(counterM);
+				outputSNP[snpSize] = '\0';
+counterM = 0;
+			}
+if(counterD != 0)
+{
+delete[counterD] = '\0';
+counterD=0;
+sprintf(outputSNP, "%s^%s", outputSNP, delete);
+				snpSize += strlen(delete) + 1;
+				outputSNP[snpSize] = '\0';
+				delete[0] = '\0';
+}
+			sprintf(outputSNP,"%s%c",outputSNP,matrix[i]);
+			snpSize += 1;
+			outputSNP[snpSize] = '\0';
+		}
+		i++;
+	}
+	if(counterM != 0)
+{
+	        sprintf(outputSNP, "%s%d", outputSNP, counterM);
+		snpSize += addCigarSize(counterM);
+		outputSNP[snpSize] = '\0';
+counterM = 0;
+}
+else if(counterD != 0)
+{
+	delete[counterD] = '\0';
+sprintf(outputSNP, "%s^%s", outputSNP, delete);
+		snpSize += strlen(delete) + 1;
+		outputSNP[snpSize] = '\0';
+counterD = 0;
+}
+	outputSNP[snpSize] = '\0';
+}
+/**********************************************/
+/*
+	direction = 0 forward
+		    1 backward
+*/
+void mapSingleEndSeq(unsigned int *l1, int s1, int readNumber, int readSegment, int direction)
+{
+		int j = 0;
+		int z = 0;
+		int *locs = (int *) l1;
+		char *_tmpSeq, *_tmpQual;
+		char rqual[SEQ_LENGTH+1];
+		rqual[SEQ_LENGTH]='\0';
+		int genLoc = 0;
+		int leftSeqLength = 0;
+		int rightSeqLength = 0;
+		int middleSeqLength = 0;
+		char matrix[200];
+		char editString[200];
+		char cigar[MAX_CIGAR_SIZE];
+		short *_tmpHashValue;
+		if (direction)
+		{
+			reverse(_msf_seqList[readNumber].qual, rqual, SEQ_LENGTH);
+			_tmpQual = rqual;
+			_tmpSeq = _msf_seqList[readNumber].rseq;
+			_tmpHashValue = _msf_seqList[readNumber].rhashValue;
+		}
+		else
+		{
+			_tmpQual = _msf_seqList[readNumber].qual;
+			_tmpSeq = _msf_seqList[readNumber].seq;
+_tmpHashValue = _msf_seqList[readNumber].hashValue;
+		}
+		int readId = 2*readNumber+direction;
+		for (z=0; z<s1; z++)
+		{
+			int map_location = 0;
+			int a = 0;
+			int o = readSegment;
+			genLoc = locs[z];//-_msf_samplingLocs[o];
+			if ( genLoc-_msf_samplingLocs[o] < _msf_refGenBeg ||
+			     genLoc-_msf_samplingLocs[o] > _msf_refGenEnd ||
+			     _msf_verifiedLocs[genLoc-_msf_samplingLocs[o]] == readId ||
+			     _msf_verifiedLocs[genLoc-_msf_samplingLocs[o]] == -readId
+			   )
+				continue;
+			int err = -1;
+			map_location = 0;
+			leftSeqLength = _msf_samplingLocs[o];
+			middleSeqLength = WINDOW_SIZE;
+			a = leftSeqLength + middleSeqLength;
+			rightSeqLength = SEQ_LENGTH - a;
+			if(errThreshold == 2)
+				 err = verifySingleEndEditDistance2(genLoc, _tmpSeq, leftSeqLength,
+								   _tmpSeq + a, rightSeqLength,
+								  middleSeqLength, matrix, &map_location, _tmpHashValue);
+			else if(errThreshold == 4)
+				err = verifySingleEndEditDistance4(genLoc, _tmpSeq, leftSeqLength,
+								   _tmpSeq + a, rightSeqLength,
+								  middleSeqLength, matrix, &map_location, _tmpHashValue);
+			else if(errThreshold ==3)
+				err = verifySingleEndEditDistance(genLoc, _tmpSeq, leftSeqLength,
+_tmpSeq + a, rightSeqLength,
+middleSeqLength, matrix, &map_location, _tmpHashValue);
+			/*else if(errThreshold == 6)
+				err = verifySingleEndEditDistance(genLoc, _tmpSeq, leftSeqLength,
+								   _tmpSeq + a, rightSeqLength,
+								  middleSeqLength, matrix, &map_location, _tmpHashValue);
+			*/
+			else
+				err = verifySingleEndEditDistanceExtention(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);
+			if(err != -1)
+			{
+				generateSNPSAM(matrix, strlen(matrix), editString);
+				generateCigar(matrix, strlen(matrix), cigar);
+			}
+			if(err != -1 && !bestMode)
+			{
+				mappingCnt++;
+				int j = 0;
+				int k = 0;
+				for(k = 0; k < readSegment+1; k++)
+				{
+					for(j = -errThreshold ; j <= errThreshold; j++)
+					{
+						if(genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j >= _msf_refGenBeg &&
+						   genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j <= _msf_refGenEnd)
+							_msf_verifiedLocs[genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j] = readId;
+					}
+				}
+				_msf_seqList[readNumber].hits[0]++;
+				_msf_output.QNAME		= _msf_seqList[readNumber].name;
+				_msf_output.FLAG		= 16 * direction;
+				_msf_output.RNAME		= _msf_refGenName;
+				_msf_output.POS			= map_location + _msf_refGenOffset;
+				_msf_output.MAPQ		= 255;
+				_msf_output.CIGAR		= cigar;
+				_msf_output.MRNAME		= "*";
+				_msf_output.MPOS		= 0;
+				_msf_output.ISIZE		= 0;
+				_msf_output.SEQ			= _tmpSeq;
+				_msf_output.QUAL		= _tmpQual;
+				_msf_output.optSize		= 2;
+				_msf_output.optFields	= _msf_optionalFields;
+				_msf_optionalFields[0].tag = "NM";
+				_msf_optionalFields[0].type = 'i';
+				_msf_optionalFields[0].iVal = err;
+				_msf_optionalFields[1].tag = "MD";
+				_msf_optionalFields[1].type = 'Z';
+				_msf_optionalFields[1].sVal = editString;
+				output(_msf_output);
+				if (_msf_seqList[readNumber].hits[0] == 1)
+				{
+					mappedSeqCnt++;
+				}
+				if ( maxHits == 0 )
+				{
+					_msf_seqList[readNumber].hits[0] = 2;
+				}
+				if ( maxHits!=0 && _msf_seqList[readNumber].hits[0] == maxHits)
+				{
+					completedSeqCnt++;
+					break;
+				}
+			}
+			else if(err != -1 && bestMode)
+			{
+				mappingCnt++;
+				_msf_seqList[readNumber].hits[0]++;
+				if (_msf_seqList[readNumber].hits[0] == 1)
+				{
+					mappedSeqCnt++;
+				}
+				if ( maxHits == 0 )
+				{
+					_msf_seqList[readNumber].hits[0] = 2;
+				}
+				if(err  < bestHitMappingInfo[readNumber].err || bestHitMappingInfo[readNumber].loc == -1)
+				{
+					setFullMappingInfo(readNumber, map_location + _msf_refGenOffset, direction, err, 0, editString, _msf_refGenName, cigar );
+				}
+			}
+			else
+			{
+				for(j = -errThreshold ; j <= errThreshold; j++)
+				{
+					if(genLoc+j > _msf_refGenBeg &&
+					   genLoc+j < _msf_refGenEnd)
+						_msf_verifiedLocs[genLoc+j] = -readId;
+				}
+			}
+		}
+}
+int mapAllSingleEndSeq()
+{
+	int i = 0;
+	int j = 0;
+	int k  = 0;
+	unsigned int *locs = NULL;
+	int prev_hash = 0;
+	for(i = 0; i < _msf_seqListSize; i++)
+	{
+		for(j = 0; j < _msf_samplingLocsSize; j++)
+		{
+			k = _msf_sort_seqList[i].readNumber;
+//			if(j != 0)
+//				if(strncmp(_msf_seqList[k].seq+_msf_samplingLocs[j], _msf_seqList[k].seq+_msf_samplingLocs[j-1], segSize) == 0)
+//					continue;
+//			if(prev_hash ==  hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]))
+//				continue;
+			locs = getCandidates ( hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]));
+			if ( locs != NULL)
+			{
+				mapSingleEndSeq(locs+1, locs[0],k ,j, 0);
+			}
+			prev_hash = hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]);
+		}
+	}
+	i = 0;
+	for(i = 0; i < _msf_seqListSize; i++)
+	{
+		for(j = 0; j < _msf_samplingLocsSize; j++)
+		{
+			k = _msf_sort_seqList[i].readNumber;
+//			if(j != 0)
+//				if(strncmp(_msf_seqList[k].rseq+_msf_samplingLocs[j], _msf_seqList[k].rseq+_msf_samplingLocs[j-1], segSize) == 0)
+//					continue;
+//			if(prev_hash ==  hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]))
+//				continue;
+			locs = getCandidates ( hashVal(_msf_seqList[k].rseq+_msf_samplingLocs[j]));
+			if ( locs != NULL)
+			{
+				mapSingleEndSeq(locs+1, locs[0],k ,j, 1);
+			}
+			prev_hash = hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]);
+		}
+	}
+	return 1;
+}
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+int compareOut (const void *a, const void *b)
+{
+	FullMappingInfo *aInfo = (FullMappingInfo *)a;
+	FullMappingInfo *bInfo = (FullMappingInfo *)b;
+	return aInfo->loc - bInfo->loc;
+}
+/**********************************************/
+/*
+	direction 0: Forward
+		  1: Reverse
+*/
+void mapPairEndSeqList(unsigned int *l1, int s1, int readNumber, int readSegment, int direction)
+{
+		int z = 0;
+		int *locs = (int *) l1;
+		char *_tmpSeq;
+		char rqual[SEQ_LENGTH+1];
+		char matrix[200];
+		char editString[200];
+		char cigar[MAX_CIGAR_SIZE];
+		short *_tmpHashValue;
+		int leftSeqLength = 0;
+		int middleSeqLength = 0;
+		int rightSeqLength =0;
+		int a = 0;
+		rqual[SEQ_LENGTH]='\0';
+		int r = readNumber;
+		char d = (direction==1)?-1:1;
+		if (d==-1)
+		{
+			_tmpSeq = _msf_seqList[readNumber].rseq;
+			_tmpHashValue = _msf_seqList[r].rhashValue;
+		}
+		else
+		{
+			_tmpSeq = _msf_seqList[readNumber].seq;
+			_tmpHashValue = _msf_seqList[r].hashValue;
+		}
+		int readId = 2*readNumber+direction;
+		for (z=0; z<s1; z++)
+		{
+			int genLoc = locs[z];//-_msf_samplingLocs[o];
+			int err = -1;
+			int map_location = 0;
+			int o = readSegment;
+			leftSeqLength = _msf_samplingLocs[o];
+			middleSeqLength = WINDOW_SIZE;
+			a = leftSeqLength + middleSeqLength;
+			rightSeqLength = SEQ_LENGTH - a;
+			if(genLoc - leftSeqLength < _msf_refGenBeg || genLoc + rightSeqLength + middleSeqLength > _msf_refGenEnd ||
+			   _msf_verifiedLocs[genLoc-_msf_samplingLocs[o]] == readId || _msf_verifiedLocs[genLoc-_msf_samplingLocs[o]] == -readId)
+				continue;
+			if(errThreshold == 2)
+				err = verifySingleEndEditDistance2(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);
+			else if(errThreshold == 4)
+				err = verifySingleEndEditDistance4(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);
+			else if(errThreshold ==3)
+				err = verifySingleEndEditDistance(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);
+			/*else if(errThreshold == 6)
+				err = verifySingleEndEditDistance(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);*/
+			else
+				err = verifySingleEndEditDistanceExtention(genLoc, _tmpSeq, leftSeqLength,
+						_tmpSeq + a, rightSeqLength,
+						middleSeqLength, matrix, &map_location, _tmpHashValue);
+			if (err != -1)
+			{
+				int j = 0;
+				int k = 0;
+				for(k = 0; k < readSegment+1; k++)
+				{
+					for(j = -errThreshold ; j <= errThreshold; j++)
+					{
+						if(genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j >= _msf_refGenBeg &&
+								genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j <= _msf_refGenEnd)
+							_msf_verifiedLocs[genLoc-(k*(_msf_samplingLocs[1]-_msf_samplingLocs[0]))+j] = readId;
+					}
+				}
+				generateSNPSAM(matrix, strlen(matrix), editString);
+				generateCigar(matrix, strlen(matrix), cigar);
+				MappingLocations *parent = NULL;
+				MappingLocations *child = _msf_mappingInfo[r].next;
+				genLoc = map_location + _msf_refGenOffset;
+				int i = 0;
+				for (i=0; i<(_msf_mappingInfo[r].size/MAP_CHUNKS); i++)
+				{
+					parent = child;
+					child = child->next;
+				}
+				if (child==NULL)
+				{
+					MappingLocations *tmp = getMem(sizeof(MappingLocations));
+					tmp->next = NULL;
+					tmp->loc[0]=genLoc * d;
+					tmp->err[0]=err;
+					tmp->cigarSize[0] = strlen(cigar);
+					sprintf(tmp->cigar[0],"%s", cigar);
+					tmp->mdSize[0] = strlen(editString);
+					sprintf(tmp->md[0],"%s", editString);
+					if (parent == NULL)
+						_msf_mappingInfo[r].next = tmp;
+					else
+						parent->next = tmp;
+				}
+				else
+				{
+					if(strlen(cigar) > SEQ_LENGTH || strlen(editString) > SEQ_LENGTH)
+					{
+					  printf("ERROR in %d read size(After mapping) exceedes cigar=%d md =%d cigar=%s md =%s\n", r,  (int)strlen(cigar), (int)strlen(editString), cigar, editString);
+					}
+					child->loc[_msf_mappingInfo[r].size % MAP_CHUNKS] = genLoc * d;
+					child->err[_msf_mappingInfo[r].size % MAP_CHUNKS] = err;
+					child->cigarSize[_msf_mappingInfo[r].size % MAP_CHUNKS] = strlen(cigar);
+					sprintf(child->cigar[_msf_mappingInfo[r].size % MAP_CHUNKS],"%s",cigar);
+					child->mdSize[_msf_mappingInfo[r].size % MAP_CHUNKS] = strlen(editString);
+					sprintf(child->md[_msf_mappingInfo[r].size % MAP_CHUNKS],"%s",editString);
+				}
+				_msf_mappingInfo[r].size++;
+			}
+			else
+			{
+				_msf_verifiedLocs[genLoc] = -readId;
+			}
+		}
+}
+/**********************************************/
+void mapPairedEndSeq()
+{
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	unsigned int *locs = NULL;
+	while ( i < _msf_seqListSize )
+	{
+		for(j = 0; j < _msf_samplingLocsSize; j++)
+		{
+			k = _msf_sort_seqList[i].readNumber;
+			locs = getCandidates ( hashVal(_msf_seqList[k].seq+_msf_samplingLocs[j]));
+			if ( locs != NULL)
+			{
+				mapPairEndSeqList(locs+1, locs[0],k ,j, 0);
+			}
+		}
+		i++;
+	}
+	i = 0;
+	while ( i < _msf_seqListSize )
+	{
+		for(j = 0; j < _msf_samplingLocsSize; j++)
+		{
+			k = _msf_sort_seqList[i].readNumber;
+			locs = getCandidates ( hashVal(_msf_seqList[k].rseq+_msf_samplingLocs[j]));
+			if ( locs != NULL)
+			{
+				mapPairEndSeqList(locs+1, locs[0],k ,j, 1);
+			}
+		}
+		i++;
+	}
+	char fname1[FILE_NAME_LENGTH];
+	char fname2[FILE_NAME_LENGTH];
+	MappingLocations *cur;
+	int tmpOut;
+	int lmax=0, rmax=0;
+	sprintf(fname1, "%s__%s__%s__%d__1.tmp",mappingOutputPath, _msf_refGenName, mappingOutput, _msf_openFiles);
+	sprintf(fname2, "%s__%s__%s__%d__2.tmp",mappingOutputPath, _msf_refGenName, mappingOutput, _msf_openFiles);
+	FILE* out;
+	FILE* out1 = fileOpen(fname1, "w");
+	FILE* out2 = fileOpen(fname2, "w");
+	_msf_openFiles++;
+	for (i=0; i<_msf_seqListSize; i++)
+	{
+		if (i%2==0)
+		{
+			out = out1;
+			if (lmax <  _msf_mappingInfo[i].size)
+			{
+				lmax = _msf_mappingInfo[i].size;
+			}
+		}
+		else
+		{
+			out = out2;
+			if (rmax < _msf_mappingInfo[i].size)
+			{
+				rmax = _msf_mappingInfo[i].size;
+			}
+		}
+		tmpOut = fwrite(&(_msf_mappingInfo[i].size), sizeof(int), 1, out);
+		if (_msf_mappingInfo[i].size > 0)
+		{
+			cur = _msf_mappingInfo[i].next;
+			for (j=0; j < _msf_mappingInfo[i].size; j++)
+			{
+				if ( j>0  && j%MAP_CHUNKS==0)
+				{
+					cur = cur->next;
+				}
+				if(cur->cigarSize[j % MAP_CHUNKS] > SEQ_LENGTH || cur->mdSize[j % MAP_CHUNKS] > SEQ_LENGTH)
+				{
+					printf("ERROR in %d read size exceeds cigar=%d md =%d cigar=%s md =%s\n", i,  cur->cigarSize[j % MAP_CHUNKS], cur->mdSize[j % MAP_CHUNKS], cur->cigar[j % MAP_CHUNKS], cur->md[j % MAP_CHUNKS]);
+				}
+				tmpOut = fwrite(&(cur->loc[j % MAP_CHUNKS]), sizeof(int), 1, out);
+				tmpOut = fwrite(&(cur->err[j % MAP_CHUNKS]), sizeof(int), 1, out);
+				tmpOut = fwrite(&(cur->cigarSize[j % MAP_CHUNKS]), sizeof(int), 1, out);
+				tmpOut = fwrite((cur->cigar[j % MAP_CHUNKS]), sizeof(char), (cur->cigarSize[j % MAP_CHUNKS]), out);
+				tmpOut = fwrite(&(cur->mdSize[j % MAP_CHUNKS]), sizeof(int), 1, out);
+				tmpOut = fwrite((cur->md[j % MAP_CHUNKS]), sizeof(char), (cur->mdSize[j % MAP_CHUNKS]), out);
+			}
+			_msf_mappingInfo[i].size = 0;
+			//_msf_mappingInfo[i].next = NULL;
+		}
+	}
+	_msf_maxLSize += lmax;
+	_msf_maxRSize += rmax;
+	fclose(out1);
+	fclose(out2);
+}
+void outputPairFullMappingInfo(FILE *fp, int readNumber)
+{
+	char *seq1, *seq2, *rseq1, *rseq2, *qual1, *qual2;
+	char rqual1[SEQ_LENGTH+1], rqual2[SEQ_LENGTH+1];
+	rqual1[SEQ_LENGTH] = rqual2[SEQ_LENGTH] = '\0';
+	seq1 = _msf_seqList[readNumber*2].seq;
+	rseq1 = _msf_seqList[readNumber*2].rseq;
+	qual1 = _msf_seqList[readNumber*2].qual;
+	reverse(_msf_seqList[readNumber*2].qual, rqual1, SEQ_LENGTH);
+	seq2 = _msf_seqList[readNumber*2+1].seq;
+	rseq2 = _msf_seqList[readNumber*2+1].rseq;
+	qual2 = _msf_seqList[readNumber*2+1].qual;
+	reverse(_msf_seqList[readNumber*2+1].qual, rqual2, SEQ_LENGTH);
+	if(bestHitMappingInfo[readNumber*2].loc == -1 && bestHitMappingInfo[readNumber*2+1].loc == -1)
+		return;
+	else
+	{
+		char *seq;
+		char *qual;
+		char d1;
+		char d2;
+		int isize;
+		int proper=0;
+		// ISIZE CALCULATION
+		// The distance between outer edges
+		isize = abs(bestHitMappingInfo[readNumber*2].loc - bestHitMappingInfo[readNumber*2+1].loc)+SEQ_LENGTH - 2;
+		if (bestHitMappingInfo[readNumber*2].loc - bestHitMappingInfo[readNumber*2+1].loc > 0)
+		{
+			isize *= -1;
+		}
+		d1 = (bestHitMappingInfo[readNumber*2].dir == -1)?1:0;
+		d2 = (bestHitMappingInfo[readNumber*2+1].dir == -1)?1:0;
+		if ( d1 )
+		{
+			seq = rseq1;
+			qual = rqual1;
+		}
+		else
+		{
+			seq = seq1;
+			qual = qual1;
+		}
+		if ( (bestHitMappingInfo[readNumber*2].loc < bestHitMappingInfo[readNumber*2+1].loc && !d1 && d2) ||
+			 (bestHitMappingInfo[readNumber*2].loc > bestHitMappingInfo[readNumber*2+1].loc && d1 && !d2) )
+		{
+			proper = 2;
+		}
+		else
+		{
+			proper = 0;
+		}
+		_msf_output.POS                 = bestHitMappingInfo[readNumber*2].loc;
+		_msf_output.MPOS                = bestHitMappingInfo[readNumber*2+1].loc;
+		_msf_output.FLAG                = 1+proper+16*d1+32*d2+64;
+		_msf_output.ISIZE               = isize;
+		_msf_output.SEQ                 = seq,
+		_msf_output.QUAL                = qual;
+		_msf_output.QNAME               = _msf_seqList[readNumber*2].name;
+		_msf_output.RNAME               = bestHitMappingInfo[readNumber*2].chr;
+		_msf_output.MAPQ                = 255;
+		_msf_output.CIGAR               = bestHitMappingInfo[readNumber*2].cigar;
+		_msf_output.MRNAME              = "=";
+		_msf_output.optSize     = 2;
+		_msf_output.optFields   = _msf_optionalFields;
+		_msf_optionalFields[0].tag = "NM";
+		_msf_optionalFields[0].type = 'i';
+		_msf_optionalFields[0].iVal = bestHitMappingInfo[readNumber*2].err;
+		_msf_optionalFields[1].tag = "MD";
+		_msf_optionalFields[1].type = 'Z';
+		_msf_optionalFields[1].sVal = bestHitMappingInfo[readNumber*2].md;
+		outputSAM(fp, _msf_output);
+		output(_msf_output);
+		if ( d2 )
+		{
+			seq = rseq2;
+			qual = rqual2;
+		}
+		else
+		{
+			seq = seq2;
+			qual = qual2;
+		}
+		_msf_output.POS                 = bestHitMappingInfo[readNumber*2+1].loc;
+		_msf_output.MPOS                = bestHitMappingInfo[readNumber*2].loc;
+		_msf_output.FLAG                = 1+proper+16*d2+32*d1+128;
+		_msf_output.ISIZE               = -isize;
+		_msf_output.SEQ                 = seq,
+		_msf_output.QUAL                = qual;
+		_msf_output.QNAME               = _msf_seqList[readNumber*2].name;
+		_msf_output.RNAME               = bestHitMappingInfo[readNumber*2].chr;
+		_msf_output.MAPQ                = 255;
+		_msf_output.CIGAR               = bestHitMappingInfo[readNumber*2+1].cigar;
+		_msf_output.MRNAME              = "=";
+		_msf_output.optSize     = 2;
+		_msf_output.optFields   = _msf_optionalFields;
+		_msf_optionalFields[0].tag = "NM";
+		_msf_optionalFields[0].type = 'i';
+		_msf_optionalFields[0].iVal = bestHitMappingInfo[readNumber*2+1].err;
+		_msf_optionalFields[1].tag = "MD";
+		_msf_optionalFields[1].type = 'Z';
+		_msf_optionalFields[1].sVal = bestHitMappingInfo[readNumber*2+1].md;
+		outputSAM(fp, _msf_output);
+		output(_msf_output);
+	}
+}
+/*
+Find the closet one to the c
+@return 0: if the x1 is closer to c
+	 1: if the x2 is closer to c
+	 2: if both distance are equal
+	 -1: if error
+*/
+int findNearest(int x1, int x2, int c)
+{
+	if (abs(x1 - c) > abs(x2 - c) )
+		return 0;
+	else if ( abs(x1 - c) <  abs(x2 - c) )
+		return 1;
+	else if ( abs(x1 - c) == abs(x2 - c) )
+		return 2;
+	else
+		return -1;
+}
+void initBestConcordantDiscordant(int readNumber)
+{
+char bestConcordantFileName[FILE_NAME_LENGTH];
+char bestDiscordantFileName[FILE_NAME_LENGTH];
+	//OPEN THE BEST CONCORDANT FILE
+//BEGIN{Farhad Hormozdiari}
+sprintf(bestConcordantFileName, "%s%s__BEST.CONCORDANT", mappingOutputPath, mappingOutput);
+bestConcordantFILE = fileOpen(bestConcordantFileName, "w");
+//END{Farhad Hormozdiari}
+//OPEN THE BEST DISCORDANT FILE
+//BEGIN{Farhad Hormozdiari}
+sprintf(bestDiscordantFileName, "%s%s__BEST.DISCORDANT", mappingOutputPath, mappingOutput);
+bestDiscordantFILE = fileOpen(bestDiscordantFileName, "w");
+//END{Farhad Hormozdiari}
+	initBestMapping(readNumber);
+}
+void finalizeBestConcordantDiscordant()
+{
+	int i = 0;
+	for(i = 0;  i<_msf_seqListSize/2; i++)
+	{
+		if(_msf_readHasConcordantMapping[i]==1)
+			outputPairFullMappingInfo(bestConcordantFILE, i);
+		else
+			outputPairFullMappingInfo(bestDiscordantFILE, i);
+	}
+	fclose(bestConcordantFILE);
+	fclose(bestDiscordantFILE);
+	freeMem(bestHitMappingInfo, _msf_seqListSize * sizeof(FullMappingInfo));
+}
+void setFullMappingInfo(int readNumber, int loc, int dir, int err, int score, char *md, char * refName, char *cigar)
+{
+	bestHitMappingInfo[readNumber].loc   = loc;
+	bestHitMappingInfo[readNumber].dir   = dir;
+	bestHitMappingInfo[readNumber].err   = err;
+	bestHitMappingInfo[readNumber].score = score;
+	strncpy(bestHitMappingInfo[readNumber].md,  md, strlen(md));
+	strncpy(bestHitMappingInfo[readNumber].chr, refName, strlen(refName));
+	strncpy(bestHitMappingInfo[readNumber].cigar, cigar, strlen(cigar));
+}
+void setPairFullMappingInfo(int readNumber, FullMappingInfo mi1, FullMappingInfo mi2)
+{
+	bestHitMappingInfo[readNumber*2].loc   = mi1.loc;
+	bestHitMappingInfo[readNumber*2].dir   = mi1.dir;
+	bestHitMappingInfo[readNumber*2].err   = mi1.err;
+	bestHitMappingInfo[readNumber*2].score = mi1.score;
+	snprintf(bestHitMappingInfo[readNumber*2].chr, MAX_REF_SIZE, "%s", _msf_refGenName);
+	strncpy(bestHitMappingInfo[readNumber*2].md, mi1.md, strlen(mi1.md));
+	strncpy(bestHitMappingInfo[readNumber*2].cigar, mi1.cigar, strlen(mi1.cigar));
+	/*
+	sprintf(bestHitMappingInfo[readNumber*2].md, "%s\0",  mi1.md);
+	sprintf(bestHitMappingInfo[readNumber*2].cigar, "%s\0", mi1.cigar);
+	*/
+	bestHitMappingInfo[readNumber*2+1].loc   = mi2.loc;
+	bestHitMappingInfo[readNumber*2+1].dir   = mi2.dir;
+	bestHitMappingInfo[readNumber*2+1].err   = mi2.err;
+	bestHitMappingInfo[readNumber*2+1].score = mi2.score;
+	snprintf(bestHitMappingInfo[readNumber*2+1].chr, MAX_REF_SIZE, "%s", _msf_refGenName);
+	/*
+	sprintf(bestHitMappingInfo[readNumber*2+1].md, "%s\0", mi2.md);
+	sprintf(bestHitMappingInfo[readNumber*2+1].cigar, "%s\0", mi2.cigar);
+	*/
+	strncpy(bestHitMappingInfo[readNumber*2+1].md, mi2.md, strlen(mi2.md));
+	strncpy(bestHitMappingInfo[readNumber*2+1].cigar, mi2.cigar, strlen(mi2.cigar));
+}
+/**********************************************/
+void outputPairedEnd()
+{
+	int i = 0;
+	char cigar[MAX_CIGAR_SIZE];
+	int tmpOut;
+	loadRefGenome(&_msf_refGen, &_msf_refGenName, &tmpOut);
+	FILE* in1[_msf_openFiles];
+	FILE* in2[_msf_openFiles];
+	char fname1[_msf_openFiles][FILE_NAME_LENGTH];
+	char fname2[_msf_openFiles][FILE_NAME_LENGTH];
+	// discordant
+	FILE *out=NULL, *out1=NULL;
+	char fname3[FILE_NAME_LENGTH];
+	char fname4[FILE_NAME_LENGTH];
+	int meanDistanceMapping = 0;
+	char *rqual1;
+	char *rqual2;
+	rqual1 = getMem((SEQ_LENGTH+1)*sizeof(char));
+	rqual2 = getMem((SEQ_LENGTH+1)*sizeof(char));
+	if (pairedEndDiscordantMode)
+	{
+		sprintf(fname3, "%s__%s__disc", mappingOutputPath, mappingOutput);
+		sprintf(fname4, "%s__%s__oea", mappingOutputPath, mappingOutput);
+		out = fileOpen(fname3, "a");
+		out1 = fileOpen(fname4, "a");
+	}
+	FullMappingInfo *mi1 = getMem(sizeof(FullMappingInfo) * _msf_maxLSize);
+	FullMappingInfo *mi2 = getMem(sizeof(FullMappingInfo) * _msf_maxRSize);
+	_msf_fileCount[_msf_maxFile] = 0;
+	for (i=0; i<_msf_openFiles; i++)
+	{
+		sprintf(fname1[i], "%s__%s__%s__%d__1.tmp", mappingOutputPath, _msf_refGenName, mappingOutput, i);
+		sprintf(_msf_fileName[_msf_maxFile][_msf_fileCount[_msf_maxFile]][0], "%s", fname1[i]);
+		sprintf(fname2[i], "%s__%s__%s__%d__2.tmp", mappingOutputPath, _msf_refGenName, mappingOutput, i);
+		sprintf(_msf_fileName[_msf_maxFile][_msf_fileCount[_msf_maxFile]][1], "%s", fname2[i]);
+		in1[i] = fileOpen(fname1[i], "r");
+		in2[i] = fileOpen(fname2[i], "r");
+		_msf_fileCount[_msf_maxFile]++;
+	}
+	_msf_maxFile++;
+	int size;
+	int j, k;
+	int size1, size2;
+	meanDistanceMapping = (pairedEndDiscordantMode==1)? (minPairEndedDiscordantDistance+maxPairEndedDiscordantDistance)/2 + SEQ_LENGTH
+																	  : (minPairEndedDistance + maxPairEndedDistance) / 2 + SEQ_LENGTH;
+	for (i=0; i<_msf_seqListSize/2; i++)
+	{
+		size1 = size2 = 0;
+		for (j=0; j<_msf_openFiles; j++)
+		{
+			tmpOut = fread(&size, sizeof(int), 1, in1[j]);
+			if ( size > 0 )
+			{
+				for (k=0; k<size; k++)
+				{
+					mi1[size1+k].dir = 1;
+					tmpOut = fread (&(mi1[size1+k].loc), sizeof(int), 1, in1[j]);
+					tmpOut = fread (&(mi1[size1+k].err), sizeof(int), 1, in1[j]);
+					tmpOut = fread (&(mi1[size1+k].cigarSize), sizeof(int), 1, in1[j]);
+					tmpOut = fread ((mi1[size1+k].cigar), sizeof(char), mi1[size1+k].cigarSize, in1[j]);
+					mi1[size1+k].cigar[mi1[size1+k].cigarSize] = '\0';
+					tmpOut = fread (&(mi1[size1+k].mdSize), sizeof(int), 1, in1[j]);
+					tmpOut = fread ((mi1[size1+k].md), sizeof(char), (mi1[size1+k].mdSize), in1[j]);
+					mi1[size1+k].md[mi1[size1+k].mdSize] = '\0';
+					if (mi1[size1+k].loc<1)
+					{
+						mi1[size1+k].loc *= -1;
+						mi1[size1+k].dir = -1;
+					}
+				}
+				qsort(mi1+size1, size, sizeof(FullMappingInfo), compareOut);
+				size1+=size;
+			}
+		}
+		for (j=0; j<_msf_openFiles; j++)
+		{
+			tmpOut = fread(&size, sizeof(int), 1, in2[j]);
+			if ( size > 0 )
+			{
+				for (k=0; k<size; k++)
+				{
+					mi2[size2+k].dir = 1;
+					tmpOut = fread (&(mi2[size2+k].loc), sizeof(int), 1, in2[j]);
+					tmpOut = fread (&(mi2[size2+k].err), sizeof(int), 1, in2[j]);
+					tmpOut = fread (&(mi2[size2+k].cigarSize), sizeof(int), 1, in2[j]);
+					tmpOut = fread ((mi2[size2+k].cigar), sizeof(char), mi2[size2+k].cigarSize, in2[j]);
+					mi2[size2+k].cigar[mi2[size2+k].cigarSize] = '\0';
+					tmpOut = fread (&(mi2[size2+k].mdSize), sizeof(int), 1, in2[j]);
+					tmpOut = fread ((mi2[size2+k].md), sizeof(char), mi2[size2+k].mdSize, in2[j]);
+					mi2[size2+k].md[mi2[size2+k].mdSize] = '\0';
+					if (mi2[size2+k].loc<1)
+					{
+						mi2[size2+k].loc *= -1;
+						mi2[size2+k].dir = -1;
+					}
+				}
+				qsort(mi2+size2, size, sizeof(FullMappingInfo), compareOut);
+				size2+=size;
+			}
+		}
+		int lm, ll, rl, rm;
+		int pos = 0;
+		if (pairedEndDiscordantMode)
+		{
+			for (j=0; j<size1; j++)
+			{
+				lm = mi1[j].loc - maxPairEndedDiscordantDistance + 1;
+				ll = mi1[j].loc - minPairEndedDiscordantDistance + 1;
+				rl = mi1[j].loc + minPairEndedDiscordantDistance - 1;
+				rm = mi1[j].loc + maxPairEndedDiscordantDistance - 1;
+				while (pos<size2 && mi2[pos].loc < lm)
+				{
+					pos++;
+				}
+				k = pos;
+				while (k<size2 && mi2[k].loc<=rm)
+				{
+					if ( mi2[k].loc <= ll || mi2[k].loc >= rl)
+					{
+						if ( (mi1[j].loc < mi2[k].loc && mi1[j].dir==1 && mi2[k].dir == -1)  ||
+							 (mi1[j].loc > mi2[k].loc && mi1[j].dir==-1 && mi2[k].dir == 1) )
+						{
+							_msf_seqList[i*2].hits[0]=1;
+							_msf_seqList[i*2+1].hits[0]=1;
+							if(nosamMode != 0)
+							{
+								size1=0;
+								size2=0;
+							}
+							break;
+						}
+					}
+					k++;
+				}
+			}
+			_msf_seqHits[i*2] += size1;
+			_msf_seqHits[i*2+1] += size2;
+			if (_msf_seqHits[i*2+1] * _msf_seqHits[i*2] > DISCORDANT_CUT_OFF && nosamMode != 0)
+			{
+				_msf_seqList[i*2].hits[0]=1;
+				_msf_seqList[i*2+1].hits[0]=1;
+				size1=0;
+				size2=0;
+			}
+			int tmp = 0;
+			int rNo = 0;
+			int loc = 0;
+			int err = 0;
+			float sc = 0;
+			char l = 0;
+			//write the OEA data
+			if(_msf_seqHits[i*2] == 0 )
+			{
+				for(k = 0;k < size2 && _msf_oeaMapping[i*2+1] < maxOEAOutput ;k++)
+				{
+					rNo = i*2+1;
+					loc = mi2[k].loc*mi2[k].dir;
+					err = mi2[k].err;
+					sc = mi2[k].score;
+					l = strlen(_msf_refGenName);
+					tmp = fwrite(&rNo, sizeof(int), 1, out1);
+					tmp = fwrite(&l, sizeof(char), 1, out1);
+					tmp = fwrite(_msf_refGenName, sizeof(char), l, out1);
+					tmp = fwrite(&loc, sizeof(int), 1, out1);
+					tmp = fwrite(&err, sizeof(int), 1, out1);
+					tmp = fwrite(&sc, sizeof(float), 1, out1);
+					if(mi2[k].cigarSize > SEQ_LENGTH || mi2[k].cigarSize <= 0)
+						printf("ERROR  CIGAR size=%d %s\n", mi2[k].cigarSize, _msf_seqList[i*2+1].seq);
+					tmp = fwrite (&(mi2[k].cigarSize), sizeof(int), 1, out1);
+					tmp = fwrite ((mi2[k].cigar), sizeof(char), mi2[k].cigarSize, out1);
+					tmp = fwrite (&(mi2[k].mdSize), sizeof(int), 1, out1);
+					tmp = fwrite ((mi2[k].md), sizeof(char), mi2[k].mdSize, out1);
+					_msf_oeaMapping[i*2+1]++;
+				}
+			}
+			if(_msf_seqHits[i*2+1] == 0)
+			{
+				for(j = 0;j < size1 && _msf_oeaMapping[i*2] < maxOEAOutput;j++)
+				{
+					rNo = i*2;
+					loc = mi1[j].loc*mi1[j].dir;
+					err = mi1[j].err;
+					sc = mi1[j].score;
+					l = strlen(_msf_refGenName);
+					tmp = fwrite(&rNo, sizeof(int), 1, out1);
+					tmp = fwrite(&l, sizeof(char), 1, out1);
+					tmp = fwrite(_msf_refGenName, sizeof(char), l, out1);
+					tmp = fwrite(&loc, sizeof(int), 1, out1);
+					tmp = fwrite(&err, sizeof(int), 1, out1);
+					tmp = fwrite(&sc, sizeof(float), 1, out1);
+					if(mi1[j].cigarSize > SEQ_LENGTH || mi1[j].cigarSize <= 0 )
+						printf("ERROR %d %s\n", mi1[j].cigarSize, _msf_seqList[i*2+1].seq);
+					tmp = fwrite (&(mi1[j].cigarSize), sizeof(int), 1, out1);
+					tmp = fwrite ((mi1[j].cigar), sizeof(char), mi1[j].cigarSize, out1);
+					tmp = fwrite (&(mi1[j].mdSize), sizeof(int), 1, out1);
+					tmp = fwrite ((mi1[j].md), sizeof(char), mi1[j].mdSize, out1);
+					_msf_oeaMapping[i*2]++;
+				}
+			}
+		}
+		char *seq1, *seq2, *rseq1, *rseq2, *qual1, *qual2;
+		rqual1[SEQ_LENGTH] = '\0';
+		rqual2[SEQ_LENGTH] = '\0';
+		rqual1[0] = '\0';
+		rqual2[0] = '\0';
+		seq1 = _msf_seqList[i*2].seq;
+		rseq1 = _msf_seqList[i*2].rseq;
+		qual1 = _msf_seqList[i*2].qual;
+		strncpy(rqual1, _msf_seqList[i*2].qual, SEQ_LENGTH);
+		seq2 = _msf_seqList[i*2+1].seq;
+		rseq2 = _msf_seqList[i*2+1].rseq;
+		qual2 = _msf_seqList[i*2+1].qual;
+		strncpy(rqual2, _msf_seqList[i*2+1].qual, SEQ_LENGTH);
+		if (pairedEndDiscordantMode)
+		{
+			for (k=0; k<size1; k++)
+			{
+				mi1[k].score = calculateScore(mi1[k].loc, (mi1[k].dir==-1)?rseq1:seq1, (mi1[k].dir==-1)?rqual1:qual1, mi1[k].cigar);
+			}
+			for (k=0; k<size2; k++)
+			{
+				mi2[k].score = calculateScore(mi2[k].loc, (mi2[k].dir==-1)?rseq2:seq2, (mi2[k].dir==-1)?rqual2:qual2, mi2[k].cigar);
+			}
+		}
+		if (pairedEndDiscordantMode)
+		{
+			for (j=0; j<size1; j++)
+			{
+				for(k = 0; k < size2; k++)
+				{
+					if(
+							(mi2[k].loc-mi1[j].loc >= minPairEndedDiscordantDistance &&
+							 mi2[k].loc-mi1[j].loc <= maxPairEndedDiscordantDistance &&
+							 mi1[j].dir > 0 && mi2[k].dir < 0 )
+							||
+							(mi1[j].loc-mi2[k].loc >= minPairEndedDiscordantDistance &&
+							 mi1[j].loc-mi2[k].loc <= maxPairEndedDiscordantDistance &&
+							 mi1[j].dir < 0 && mi2[k].dir > 0)
+					  )
+					{
+						//POSSIBLE CONCORDANT
+						if(_msf_readHasConcordantMapping[i] == 0)
+						{
+							setPairFullMappingInfo(i, mi1[j], mi2[k]);
+							_msf_readHasConcordantMapping[i] = 1;
+							_msf_seqList[i*2].hits[0] = 1;
+							_msf_seqList[i*2+1].hits[0] = 1;
+						}
+						else
+						{
+							if(bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err >= mi1[j].err + mi2[k].err)
+							{
+								if( bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err ==
+										mi1[j].err + mi2[k].err &&
+										findNearest(abs(bestHitMappingInfo[i*2+1].loc - bestHitMappingInfo[i*2].loc),
+											abs(mi2[k].loc - mi1[j].loc),
+											meanDistanceMapping
+											) == 0 )
+								{
+									continue;
+								}
+								setPairFullMappingInfo(i, mi1[j], mi2[k]);
+							}
+						}
+					}
+					//DISCORDANT TO TEMP FILE FOR POST PROCESSIING
+					else if(_msf_readHasConcordantMapping[i] == 0 &&
+						   	_msf_seqHits[i*2] != 0 &&
+							_msf_seqHits[i*2+1] != 0)
+					{
+						int tmp;
+						int rNo = i;
+						int loc = mi1[j].loc*mi1[j].dir;
+						int err = mi1[j].err;
+						float sc = mi1[j].score;
+						char l = strlen(_msf_refGenName);
+						if(_msf_discordantMapping[i*2] < maxDiscordantOutput)
+						{
+							tmp = fwrite(&rNo, sizeof(int), 1, out);
+							tmp = fwrite(&l, sizeof(char), 1, out);
+							tmp = fwrite(_msf_refGenName, sizeof(char), l, out);
+							tmp = fwrite(&loc, sizeof(int), 1, out);
+							tmp = fwrite(&err, sizeof(int), 1, out);
+							tmp = fwrite(&sc, sizeof(float), 1, out);
+							tmp = fwrite (&(mi1[j].cigarSize), sizeof(int), 1, out);
+							tmp = fwrite ((mi1[j].cigar), sizeof(char), mi1[j].cigarSize, out);
+							tmp = fwrite (&(mi1[j].mdSize), sizeof(int), 1, out);
+							tmp = fwrite ((mi1[j].md), sizeof(char), mi1[j].mdSize, out);
+							loc = mi2[k].loc*mi2[k].dir;
+							err = mi2[k].err;
+							sc = mi2[k].score;
+							tmp = fwrite(&loc, sizeof(int), 1, out);
+							tmp = fwrite(&err, sizeof(int), 1, out);
+							tmp = fwrite(&sc, sizeof(float), 1, out);
+							tmp = fwrite (&(mi2[k].cigarSize), sizeof(int), 1, out);
+							tmp = fwrite ((mi2[k].cigar), sizeof(char), mi2[k].cigarSize, out);
+							tmp = fwrite (&(mi2[k].mdSize), sizeof(int), 1, out);
+							tmp = fwrite ((mi2[k].md), sizeof(char), mi2[k].mdSize, out);
+							_msf_discordantMapping[i*2]++;
+						}
+						//SET THE BEST DISCORDANT
+						//BEGIN {Farhad Hormozdiari}
+						if( bestHitMappingInfo[i*2].loc == -1 &&
+								bestHitMappingInfo[i*2+1].loc == -1 &&
+								_msf_readHasConcordantMapping[i] == 0)
+						{
+							setPairFullMappingInfo(i, mi1[j], mi2[k]);
+							_msf_seqList[i*2].hits[0] = 1;
+							_msf_seqList[i*2+1].hits[0] = 1;
+						}
+						else if( bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err >= mi1[j].err + mi2[k].err
+								&&  _msf_readHasConcordantMapping[i] == 0)
+						{
+							if(bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err == mi1[j].err + mi2[k].err &&
+									findNearest( abs(bestHitMappingInfo[i*2+1].loc - bestHitMappingInfo[i*2].loc),
+										abs(mi1[j].loc - mi2[k].loc),
+										meanDistanceMapping
+										) == 0
+							  )
+							{
+								continue;
+							}
+							setPairFullMappingInfo(i, mi1[j], mi2[k]);
+						}
+						//END {Farhad Hormozdiari}
+					}
+				}
+			}
+		}
+		else
+		{
+			for (j=0; j<size1; j++)
+			{
+				for(k = 0; k < size2; k++)
+				{
+				  if((mi2[k].loc-mi1[j].loc >= minPairEndedDistance &&
+				      mi2[k].loc-mi1[j].loc <= maxPairEndedDistance &&
+				      mi1[j].dir > 0 && mi2[k].dir < 0)
+				     ||
+				     (mi1[j].loc-mi2[k].loc >= minPairEndedDistance &&
+				      mi1[j].loc-mi2[k].loc <= maxPairEndedDistance &&
+				      mi1[j].dir < 0 && mi2[k].dir > 0)
+					   )
+					{
+						char *seq;
+						char *qual;
+						char d1;
+						char d2;
+						int isize;
+						int proper=0;
+						// ISIZE CALCULATION
+						// The distance between outer edges
+						isize = abs(mi1[j].loc - mi2[k].loc)+SEQ_LENGTH-2;
+						if (mi1[j].loc - mi2[k].loc > 0)
+						{
+							isize *= -1;
+						}
+						d1 = (mi1[j].dir == -1)?1:0;
+						d2 = (mi2[k].dir == -1)?1:0;
+						//SET THE READ HAS CONCORDANT MAPPING
+						_msf_readHasConcordantMapping[i] = 1;
+						if ( d1 )
+						{
+							seq = rseq1;
+							qual = rqual1;
+						}
+						else
+						{
+							seq = seq1;
+							qual = qual1;
+						}
+						if ((mi1[j].loc < mi2[k].loc && !d1 && d2) ||
+						    (mi1[j].loc > mi2[k].loc && d1 && !d2) )
+						{
+							proper = 2;
+						}
+						else
+						{
+							proper = 0;
+						}
+						_msf_output.POS			= mi1[j].loc;
+						_msf_output.MPOS		= mi2[k].loc;
+						_msf_output.FLAG		= 1+proper+16*d1+32*d2+64;
+						_msf_output.ISIZE		= isize;
+						_msf_output.SEQ			= seq,
+						_msf_output.QUAL		= qual;
+						_msf_output.QNAME		= _msf_seqList[i*2].name;
+						_msf_output.RNAME		= _msf_refGenName;
+						_msf_output.MAPQ		= 255;
+						_msf_output.CIGAR		= cigar;
+						_msf_output.MRNAME		= "=";
+						_msf_output.optSize	= 2;
+						_msf_output.optFields	= _msf_optionalFields;
+						_msf_optionalFields[0].tag = "NM";
+						_msf_optionalFields[0].type = 'i';
+						_msf_optionalFields[0].iVal = mi1[j].err;
+						_msf_optionalFields[1].tag = "MD";
+						_msf_optionalFields[1].type = 'Z';
+						_msf_optionalFields[1].sVal = mi1[j].md;
+						if(!bestMode)
+							output(_msf_output);
+						if ( d2 )
+						{
+							seq = rseq2;
+							qual = rqual2;
+						}
+						else
+						{
+							seq = seq2;
+							qual = qual2;
+						}
+						_msf_output.POS			= mi2[k].loc;
+						_msf_output.MPOS		= mi1[j].loc;
+						_msf_output.FLAG		= 1+proper+16*d2+32*d1+128;
+						_msf_output.ISIZE		= -isize;
+						_msf_output.SEQ			= seq,
+						_msf_output.QUAL		= qual;
+						_msf_output.QNAME		= _msf_seqList[i*2].name;
+						_msf_output.RNAME		= _msf_refGenName;
+						_msf_output.MAPQ		= 255;
+						_msf_output.CIGAR		= cigar;
+						_msf_output.MRNAME		= "=";
+						_msf_output.optSize	= 2;
+						_msf_output.optFields	= _msf_optionalFields;
+						_msf_optionalFields[0].tag = "NM";
+						_msf_optionalFields[0].type = 'i';
+						_msf_optionalFields[0].iVal = mi2[k].err;;
+						_msf_optionalFields[1].tag = "MD";
+						_msf_optionalFields[1].type = 'Z';
+						_msf_optionalFields[1].sVal = mi2[k].md;
+						if(!bestMode)
+							output(_msf_output);
+						//SET THE BEST CONCORDANT
+						//BEGIN {Farhad Hormozdiari}
+							if(bestHitMappingInfo[i*2].loc == -1 && bestHitMappingInfo[i*2+1].loc == -1)
+							{
+								setPairFullMappingInfo(i, mi1[j], mi2[k]);
+							}
+							else
+							{
+								if(bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err >= mi1[j].err + mi2[k].err)
+								{
+									if( bestHitMappingInfo[i*2].err + bestHitMappingInfo[i*2+1].err == mi1[j].err + mi2[k].err &&
+									    findNearest(abs(bestHitMappingInfo[i*2+1].loc - bestHitMappingInfo[i*2].loc),
+										      abs(mi2[k].loc - mi1[j].loc),
+										      meanDistanceMapping
+									    ) == 0 )
+									{
+										continue;
+									}
+									setPairFullMappingInfo(i, mi1[j], mi2[k]);
+								}
+							}
+					//END   {Farhad Hormozdiari}
+					}
+				}
+			}
+		}
+	}
+	freeMem(rqual1, 0);
+	freeMem(rqual2, 0);
+	if (pairedEndDiscordantMode)
+	{
+		fclose(out);
+		fclose(out1);
+	}
+	for (i=0; i<_msf_openFiles; i++)
+	{
+		fclose(in1[i]);
+		fclose(in2[i]);
+		unlink(fname1[i]);
+		unlink(fname2[i]);
+	}
+	freeMem(mi1, sizeof(FullMappingInfo)*_msf_maxLSize);
+	freeMem(mi2, sizeof(FullMappingInfo)*_msf_maxRSize);
+	_msf_openFiles = 0;
+}
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+float str2int(char *str, int index1, int index2)
+{
+	char tmp[200];
+	strncpy(tmp, &str[index1], index2-index1);
+	tmp[index2-index1] = '\0';
+	return atol(tmp);
+}
+float calculateScore(int index, char *seq, char *qual,char *md)
+{
+	int i;
+	int j;
+	char *ref;
+	char *ver;
+	ref = _msf_refGen + index-1;
+	ver = seq;
+	float score = 1;
+	char tmp[200];
+	int value = 0;
+	int end = 0;
+	int index1 = 0;
+	int index2 = 0;
+	i=0;
+	while(1)
+	{
+		if(i>=strlen(md))
+			break;
+		index1 = i;
+		while(md[i] >='0' && md[i]<='9')
+		{
+			i++;
+		}
+		index2 = i;
+		value = str2int(md, index1,index2);
+		if(md[i]=='M')
+		{
+			for(j=0;j<value;j++)
+			{
+				tmp[end]='M';
+				end++;
+			}
+		}
+		else if(md[i]=='I')
+		{
+			for(j=0;j<value;j++)
+			{
+				tmp[end]='I';
+				end++;
+			}
+		}
+		else if(md[i] == 'D')
+		{
+			for(j=0;j<value;j++)
+			{
+				tmp[end]='D';
+				end++;
+			}
+		}
+		i++;
+	}
+	tmp[end] = '\0';
+	j =  0;
+	for (i = 0; i < end; i++)
+	{
+		if(tmp[i] == 'M')
+		{
+			if (*ref != *ver)
+			{
+				score *= 0.001 + 1/pow( 10, ((qual[j]-33)/10.0) );
+			}
+			ref++;
+			ver++;
+			j++;
+		}
+		else if(tmp[i] == 'I')
+		{
+			ver++;
+			j++;
+		}
+		else if(tmp[i] == 'D')
+		{
+			ref++;
+		}
+	}
+	return score;
+}
+int matoi(char *str, int start, int end)
+{
+	int i = 0;
+	char tmp[200];
+	for(i=0;i < end-start; i++)
+		tmp[i] = str[start+i];
+	tmp[i]='\0';
+	return atoi(tmp);
+}
+void convertCigarToMatrix(char *cigar, int cigar_size, char * matrix)
+{
+	int i = 0;
+	int j = 0;
+	int start = 0;
+	int size = 0;
+matrix[0] = '\0';
+	while(i < cigar_size)
+	{
+		if(cigar[i] >= '0' && cigar[i] <= '9')
+		{
+			start = i;
+			while(cigar[i] >= '0' && cigar[i] <= '9' && i < cigar_size)
+				i++;
+			int value = matoi(cigar, start, i);
+			for(j = 0;  j < value; j++)
+			{
+				if(cigar[i] == 'M')
+					matrix[size] = 'M';
+				else if(cigar[i] == 'D')
+					matrix[size] ='D';
+				else if(cigar[i] == 'I')
+					matrix[size] = 'I';
+				size++;
+			}
+		}
+		i++;
+	}
+	matrix[size] = '\0';
+}
+void convertMDToMatrix(char *md, int md_size, char * matrix)
+{
+	int i = 0;
+	int j = 0;
+	int start = 0;
+	int size = 0;
+matrix[0] = '\0';
+	while(i < md_size)
+	{
+		if(md[i] >= '0' && md[i] <= '9')
+		{
+			start = i;
+			while(md[i] >= '0' && md[i] <= '9' && i < md_size)
+				i++;
+			int value = matoi(md, start, i);
+			for(j = 0;  j < value; j++)
+			{
+				matrix[size] = 'M';
+				size++;
+			}
+			i--;
+		}
+		else if(md[i] == '^')
+		{
+			matrix[size] = 'D';
+			size++;
+		}
+		else
+		{
+			 matrix[size] = md[i];
+size++;
+		}
+		//size++;
+		i++;
+	}
+	matrix[size] = '\0';
+}
+void convertMDCigarToMatrix(char *cigar, int cigar_size, char *md, int md_size, char *matrix)
+{
+	int i = 0;
+	int j = 0;
+	int size = 0;
+	char tmp1[200];
+	char tmp2[200];
+	convertMDToMatrix(md,md_size, tmp2);
+	convertCigarToMatrix(cigar, cigar_size,tmp1);
+	while(i < strlen(tmp1))
+	{
+		if(tmp1[i]=='M')
+		{
+			if(j < strlen(tmp2))
+			{
+				if(tmp2[j]=='M')
+				{
+					matrix[size]='M';
+					size++;
+				}
+				if(tmp2[j]!='M')
+				{
+					matrix[size]=tmp2[j];
+					size++;
+				}
+			}
+			else
+			{
+				matrix[size]='M';
+				size++;
+			}
+		}
+		else if(tmp1[i] == 'D')
+		{
+			matrix[size]='D';
+			size++;
+			j++;
+			matrix[size]=tmp2[j];
+			size++;
+		}
+		else if(tmp1[i] == 'I')
+		{
+			matrix[size]='I';
+			size++;
+		}
+		i++;
+		if(j < strlen(tmp2))
+			j++;
+	}
+	if(strlen(tmp1))
+	matrix[size] = '\0';
+}
+void convertInsertion(char * in_matrix, char * seq, char *out_matrix)
+{
+	int i = 0;
+	int j = 0;
+	int size = 0;
+	while( i  < strlen(in_matrix))
+	{
+		if(in_matrix[i] == 'M')
+		{
+			out_matrix[size] = 'M';
+			size++;
+			j++;
+		}
+		else if(in_matrix[i] == 'D')
+		{
+			out_matrix[size] = 'D';
+			size++;
+			i++;
+			j++;
+			out_matrix[size] = seq[j];
+			j++;
+			size++;
+		}
+		else if(in_matrix[i] == 'I')
+		{
+			out_matrix[size] = 'I';
+			size++;
+			out_matrix[size] = seq[j];
+			size++;
+			j++;
+		}
+		else
+		{
+			out_matrix[size] = in_matrix[i];
+			size++;
+			j++;
+		}
+		i++;
+	}
+	out_matrix[size] = '\0';
+}
+/**********************************************/
+void outputPairedEndDiscPP()
+{
+	char tmp_matrix1[200];
+	char tmp_matrix2[200];
+	char matrix1[200];
+	char matrix2[200];
+	char cigar1[200];
+	char editString1[200];
+	char cigar2[200];
+	char editString2[200];
+	char seq1[SEQ_LENGTH+1];
+	char qual1[SEQ_LENGTH+1];
+	char seq2[SEQ_LENGTH+1];
+	char qual2[SEQ_LENGTH+1];
+	char genName[SEQ_LENGTH];
+	char fname1[FILE_NAME_LENGTH];
+	char fname2[FILE_NAME_LENGTH];
+	char l;
+	int l_size;
+	int loc1, loc2;
+	int err1, err2;
+	char dir1, dir2;
+	float sc1, sc2, lsc=0;
+	int flag = 0;
+	int rNo,lrNo = -1;
+	int tmp;
+	FILE *in, *out;
+	sprintf(fname1, "%s__%s__disc", mappingOutputPath, mappingOutput);
+	sprintf(fname2, "%s%s_DIVET.vh", mappingOutputPath, mappingOutput);
+	in   = fileOpen(fname1, "r");
+	out  = fileOpen(fname2, "w");
+	if (in != NULL)
+	{
+		flag = fread(&rNo, sizeof(int), 1, in);
+	}
+	else
+	{
+		flag  = 0;
+	}
+	seq1[SEQ_LENGTH]   = '\0';
+	qual1[SEQ_LENGTH]  = '\0';
+	seq2[SEQ_LENGTH]   = '\0';
+	qual2[SEQ_LENGTH]  = '\0';
+	while (flag)
+	{
+		tmp = fread(&l, sizeof(char), 1, in);
+		tmp = fread(genName, sizeof(char), l, in);
+		genName[(int)l]='\0';
+		tmp = fread(&loc1, sizeof(int), 1, in);
+		tmp = fread(&err1, sizeof(int), 1, in);
+		tmp = fread(&sc1, sizeof(float), 1, in);
+		//tmp = fwrite (&(mi2[k].cigarSize), sizeof(int), 1, out);
+		tmp = fread(&l_size, sizeof(int), 1, in);
+		tmp = fread(cigar1, sizeof(char), l_size, in);
+		cigar1[(int)l_size]='\0';
+		//tmp = fwrite ((mi2[k].cigar), sizeof(char), mi2[k].cigarSize, out);
+		//tmp = fwrite (&(mi2[k].mdSize), sizeof(int), 1, out);
+		tmp = fread(&l_size, sizeof(int), 1, in);
+		tmp = fread(editString1, sizeof(char), l_size, in);
+		editString1[(int)l_size]='\0';
+		//tmp = fwrite ((mi2[k].md), sizeof(char), mi2[k].mdSize, out);
+		tmp = fread(&loc2, sizeof(int), 1, in);
+		tmp = fread(&err2, sizeof(int), 1, in);
+		tmp = fread(&sc2, sizeof(float), 1, in);
+		tmp = fread(&l_size, sizeof(int), 1, in);
+		tmp = fread(cigar2, sizeof(char), l_size, in);
+		cigar2[(int)l_size]='\0';
+		tmp = fread(&l_size, sizeof(int), 1, in);
+		tmp = fread(editString2, sizeof(char), l_size, in);
+		editString2[(int)l_size]='\0';
+		convertMDCigarToMatrix(cigar1, strlen(cigar1), editString1, strlen(editString1), tmp_matrix1);
+		convertMDCigarToMatrix(cigar2, strlen(cigar2), editString2, strlen(editString2), tmp_matrix2);
+		if(_msf_readHasConcordantMapping[rNo] == 0)
+		{
+			dir1 = dir2 = 'F';
+			strncpy(seq1, _msf_seqList[rNo*2].seq, SEQ_LENGTH);
+			strncpy(seq2, _msf_seqList[rNo*2+1].seq, SEQ_LENGTH);
+			if (loc1 < 0)
+			{
+				dir1 = 'R';
+				loc1 = -loc1;
+				strncpy(seq1, _msf_seqList[rNo*2].rseq, SEQ_LENGTH);
+			}
+			if (loc2 < 0)
+			{
+				dir2 = 'R';
+				loc2 = -loc2;
+				strncpy(seq2, _msf_seqList[rNo*2+1].rseq, SEQ_LENGTH);
+			}
+			convertInsertion(tmp_matrix1, seq1, matrix1);
+			convertInsertion(tmp_matrix2, seq2, matrix2);
+			if (rNo != lrNo)
+			{
+				int j;
+				for (j=0; j<SEQ_LENGTH; j++)
+				{
+					lsc += _msf_seqList[rNo*2].qual[j]+_msf_seqList[rNo*2+1].qual[j];
+				}
+				lsc /= 2*SEQ_LENGTH;
+				lsc -= 33;
+				lrNo = rNo;
+			}
+			char event = '\0';
+			if ( dir1 == dir2 )
+			{
+				event = 'V';
+			}
+			else
+			{
+				if (loc1 < loc2)
+				{
+					if (dir1 == 'R' && dir2 == 'F')
+					{
+						event = 'E';
+					}
+					else if ( loc2 - loc1 >= maxPairEndedDiscordantDistance )
+					{
+						event = 'D';
+					}
+					else
+					{
+						event = 'I';
+					}
+				}
+				else if (loc2 < loc1)
+				{
+					if (dir2 == 'R' && dir1 == 'F')
+					{
+						event = 'E';
+					}
+					else if ( loc1 - loc2 >= maxPairEndedDiscordantDistance )
+					{
+						event = 'D';
+					}
+					else
+					{
+						event = 'I';
+					}
+				}
+			}
+			_msf_seqList[rNo*2].hits[0] = 2;
+			if(event != 'E')
+				fprintf(out, "%s\t%s\t%d\t%d\t%c\t%d\t%d\t%c\t%c\t%d\t%0.0f\t%e\n",
+					_msf_seqList[rNo*2].name, genName, loc1, (loc1+SEQ_LENGTH-1), dir1,
+					loc2, (loc2+SEQ_LENGTH-1), dir2, event, (err1+err2), lsc, sc1*sc2);
+		}
+		flag = fread(&rNo, sizeof(int), 1, in);
+	}
+	fclose(in);
+	fclose(out);
+	unlink(fname1);
+}
+void finalizeOEAReads(char *fileName)
+{
+	FILE *fp_out1;
+	FILE * in;
+	char genName[SEQ_LENGTH];
+	char fname1[FILE_NAME_LENGTH];
+	char fname2[FILE_NAME_LENGTH];
+	char l=0;
+	int loc1=0;
+	int err1;
+	char d;
+	float sc1=0;
+	int flag = 0;
+	int rNo=-1;
+	int tmp=0;
+	int cigarSize = 0;
+	int mdSize = 0;
+	char cigar[SEQ_LENGTH+1];
+	char md[SEQ_LENGTH+1];
+	char *seq1, *seq2, *qual1, *qual2;
+	char *rqual1, *rqual2;
+	seq1=NULL; seq2=NULL; qual1=NULL; qual2=NULL;
+	rqual1 = getMem(200*sizeof(char));
+	rqual2 = getMem(200*sizeof(char));
+	rqual1[0] = '\0';
+	rqual2[0] = '\0';
+	sprintf(fname1, "%s%s_OEA", mappingOutputPath, mappingOutput);
+	fp_out1 = fileOpen(fname1, "w");
+	in = NULL;
+	if (pairedEndDiscordantMode){
+	  sprintf(fname2, "%s__%s__oea", mappingOutputPath, mappingOutput);
+	  in   = fileOpen(fname2, "r");
+	}
+	if (in != NULL)
+	{
+		flag = fread(&rNo, sizeof(int), 1, in);
+	}
+	else
+	{
+		flag  = 0;
+	}
+	while (flag)
+	{
+		cigar[0] = '\0';
+		md[0] = '\0';
+		tmp = fread(&l, sizeof(char), 1, in);
+		tmp = fread(genName, sizeof(char), l, in);
+		genName[(int)l]='\0';
+		tmp = fread(&loc1, sizeof(int), 1, in);
+		tmp = fread(&err1, sizeof(int), 1, in);
+		tmp = fread(&sc1, sizeof(float), 1, in);
+		tmp = fread (&cigarSize, sizeof(int), 1, in);
+		tmp = fread (cigar, sizeof(char), cigarSize, in);
+		cigar[cigarSize] = '\0';
+		tmp = fread (&mdSize, sizeof(int), 1, in);
+		tmp = fread (md, sizeof(char), mdSize, in);
+		md[mdSize] = '\0';
+		d = 1;
+		if(loc1 < 0)
+		{
+			d = -1;
+			loc1 *= -1;
+			seq1 = _msf_seqList[rNo].rseq;
+			reverse(_msf_seqList[rNo].qual, rqual1, SEQ_LENGTH);
+			rqual1[SEQ_LENGTH] = '\0';
+		}
+		else
+		{
+			seq1 = _msf_seqList[rNo].seq;
+			qual1 = _msf_seqList[rNo].qual;
+		}
+		if(rNo % 2 == 0)
+		{
+			seq2 = _msf_seqList[rNo+1].seq;
+			qual2 = _msf_seqList[rNo+1].qual;
+		}
+		else
+		{
+			seq2 = _msf_seqList[rNo-1].seq;
+			qual2 = _msf_seqList[rNo-1].qual;
+		}
+		if(_msf_seqHits[rNo] != 0 && _msf_seqHits[(rNo%2==0)?rNo+1:rNo-1] == 0)
+		{
+			_msf_output.POS                 = loc1;
+			_msf_output.MPOS                = 0;
+			_msf_output.FLAG                = (rNo % 2 ==0)? 1+4+32*d+128  : 1+8+16*d+64 ;
+			_msf_output.ISIZE               = 0;
+			_msf_output.SEQ                 = seq1;
+			_msf_output.QUAL                = qual1;
+			_msf_output.QNAME               = _msf_seqList[rNo].name;
+			_msf_output.RNAME               = genName;
+			_msf_output.MAPQ                = 255;
+			_msf_output.CIGAR               = cigar;
+			_msf_output.MRNAME              = "=";
+			_msf_output.optSize     = 4;
+			_msf_output.optFields   = _msf_optionalFields;
+			_msf_optionalFields[0].tag = "NM";
+			_msf_optionalFields[0].type = 'i';
+			_msf_optionalFields[0].iVal = err1;
+			_msf_optionalFields[1].tag = "MD";
+			_msf_optionalFields[1].type = 'Z';
+			_msf_optionalFields[1].sVal = md;
+			//for the OEA reads
+			_msf_optionalFields[2].tag = "NS";
+			_msf_optionalFields[2].type = 'Z';
+			_msf_optionalFields[2].sVal = seq2;
+			_msf_optionalFields[3].tag = "NQ";
+			_msf_optionalFields[3].type = 'Z';
+			_msf_optionalFields[3].sVal = qual2;
+			outputSAM(fp_out1, _msf_output);
+			_msf_seqList[rNo].hits[0] = -1;
+			_msf_seqList[(rNo%2==0)?rNo+1:rNo-1].hits[0] = -1;
+		}
+		flag = fread(&rNo, sizeof(int), 1, in);
+	}
+	freeMem(rqual1, 0);
+	freeMem(rqual2, 0);
+	unlink(fname2);
+	fclose(fp_out1);
+}
+/*
+void outputOEA(char *fileName1, FILE * fp_out, int readSegment)
+{
+	int i = 0;
+	int j = 0;
+	char *index;
+	int size1 = 0;
+	FILE *fp1;
+	char geneFileName1[FILE_NAME_LENGTH];
+	char matrix[200];
+	char cigar[MAX_CIGAR_SIZE];
+	char editString[200];
+	FullMappingInfoLink *miL = getMem(_msf_seqListSize * sizeof(FullMappingInfoLink));
+	if(fileName1 != NULL)
+	{
+		fp1 = fileOpen(fileName1, "r");
+		index = strstr(fileName1, "__");
+		strncpy(geneFileName1, index + 2 * sizeof(char), strstr(index + 2, "__") - index - 2);
+		geneFileName1[strstr(index + 2, "__") - index - 2] = '\0';
+		for(i = 0; i < _msf_seqListSize / 2; i++)
+		{
+			fread(&size1, sizeof(int), 1, fp1);
+			miL[i].mi = getMem(size1 * sizeof(FullMappingInfo) );
+			miL[i].size = size1;
+			for(j = 0; j < size1; j++)
+			{
+				fread(&(miL[i].mi[j].loc), sizeof(int), 1, fp1);
+				fread (&(miL[i].mi[j].err), sizeof(int), 1, fp1);
+				fread (&(miL[i].mi[j].cigarSize), sizeof(int), 1, fp1);
+				fread ((miL[i].mi[j].cigar), sizeof(char), miL[i].mi[j].cigarSize+1, fp1);
+				fread (&(miL[i].mi[j].mdSize), sizeof(int), 1, fp1);
+				fread ((miL[i].mi[j].md), sizeof(char), miL[i].mi[j].mdSize+1, fp1);
+				miL[i].mi[j].dir = 1;
+				if(miL[i].mi[j].loc < 1)
+				{
+					miL[i].mi[j].loc *= -1;
+					miL[i].mi[j].dir = -1;
+				}
+			}
+			int tmpSize = (readSegment==0) ? _msf_seqHits[i*2+1]  : _msf_seqHits[i*2];
+			if(_msf_seqHits[i*2+readSegment] == 0 && size1 != 0 && _msf_oeaMapping[i*2+(readSegment == 0 ? 1: 0)] <= maxOEAOutput)
+			{
+				int d1 = 0;
+				char *seq, *qual;
+				char *seq1, *seq2, *rseq1, *rseq2, *qual1, *qual2;
+				char rqual1[SEQ_LENGTH+1], rqual2[SEQ_LENGTH+1];
+				rqual1[SEQ_LENGTH] = rqual2[SEQ_LENGTH] = '\0';
+				seq1 = _msf_seqList[i*2].seq;
+				rseq1 = _msf_seqList[i*2].rseq;
+				qual1 = _msf_seqList[i*2].qual;
+				reverse(_msf_seqList[i*2].qual, rqual1, SEQ_LENGTH);
+				seq2 = _msf_seqList[i*2+1].seq;
+				rseq2 = _msf_seqList[i*2+1].rseq;
+				qual2 = _msf_seqList[i*2+1].qual;
+				reverse(_msf_seqList[i*2+1].qual, rqual2, SEQ_LENGTH);
+				for(j = 0; j < size1 && _msf_oeaMapping[i*2+(readSegment == 0 ? 1: 0)] <= maxOEAOutput; j++)
+				{
+					d1 = (miL[i].mi[j].dir == -1)?1:0;
+					if(readSegment == 0)
+					{
+						if ( d1 )
+						{
+							seq = rseq2;
+							qual = rqual2;
+						}
+						else
+						{
+							seq = seq2;
+							qual = qual2;
+						}
+					}
+					else
+					{
+						if ( d1 )
+						{
+							seq = rseq1;
+							qual = rqual1;
+						}
+						else
+						{
+							seq = seq1;
+							qual = qual1;
+						}
+					}
+					_msf_oeaMapping[i*2+(readSegment == 0 ? 1: 0)]++;
+					_msf_output.POS                 = (readSegment==1)?miL[i].mi[j].loc:0;
+					_msf_output.MPOS                = (readSegment==0)?miL[i].mi[j].loc:0;
+					_msf_output.FLAG                = (readSegment==0)? 1+4+32*d1+128  : 1+8+16*d1+64 ;
+					_msf_output.ISIZE               = 0;
+					_msf_output.SEQ                 = seq,
+					_msf_output.QUAL                = qual;
+					_msf_output.QNAME               = _msf_seqList[i*2+(readSegment==0?1:0)].name;
+					_msf_output.RNAME               = geneFileName1;
+					_msf_output.MAPQ                = 255;
+					_msf_output.CIGAR               = miL[i].mi[j].cigar;
+					_msf_output.MRNAME              = "=";
+					//_msf_output.NSEQ                = (readSegment == 0)?seq1:seq2;
+					//_msf_output.NQUAL               = (readSegment == 0)?qual1:qual2;
+					_msf_output.optSize     = 4;
+					_msf_output.optFields   = _msf_optionalFields;
+					_msf_optionalFields[0].tag = "NM";
+					_msf_optionalFields[0].type = 'i';
+					_msf_optionalFields[0].iVal = miL[i].mi[j].err;
+					_msf_optionalFields[1].tag = "MD";
+					_msf_optionalFields[1].type = 'Z';
+					_msf_optionalFields[1].sVal = miL[i].mi[j].md;
+					//for the OEA reads
+					_msf_optionalFields[2].tag = "NS";
+					_msf_optionalFields[2].type = 'Z';
+					_msf_optionalFields[2].sVal = (readSegment == 0)?seq1:seq2;
+					_msf_optionalFields[3].tag = "NQ";
+					_msf_optionalFields[3].type = 'Z';
+					_msf_optionalFields[3].sVal = (readSegment == 0)?qual1:qual2;
+					outputSAM(fp_out, _msf_output);
+				}
+			}
+		}
+	}
+	for(i = 0; i < _msf_seqListSize / 2; i++)
+	{
+		freeMem(miL[i].mi, miL[i].size * sizeof(FullMappingInfo));
+	}
+	freeMem(miL, _msf_seqListSize * sizeof(FullMappingInfoLink));
+	fclose(fp1);
+}
+void finalizeOEAReads(char *fileName)
+{
+	int i = 0;
+	int k = 0;
+	FILE *fp_out1;
+	char fname1[200];
+	_msf_oeaMapping = getMem(_msf_seqListSize * sizeof(int));
+	for(i = 0; i < _msf_seqListSize; i++)
+	{
+		_msf_oeaMapping[i] = 0;
+	}
+	sprintf(fname1, "%s%s_OEA", mappingOutputPath, mappingOutput);
+	fp_out1 = fileOpen(fname1, "w");
+	for(i = 0; i < _msf_maxFile; i++)
+	{
+		for(k = 0; k < _msf_fileCount[i]; k++)
+		{
+			outputOEA(_msf_fileName[i][k][1], fp_out1, 0);
+		}// for k
+	} //for i
+	for(i = 0; i < _msf_maxFile; i++)
+	{
+		for(k = 0; k < _msf_fileCount[i]; k++)
+		{
+			outputOEA(_msf_fileName[i][k][0], fp_out1, 1);
+		}// for k
+	} //for i
+	fclose(fp_out1);
+}
+*/
+void outputTransChromosal(char *fileName1, char *fileName2, FILE * fp_out)
+{
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	char *index;
+	int size1 = 0;
+	int size2 = 0;
+	FILE *fp1 = NULL;
+	FILE *fp2 = NULL;
+	char geneFileName1[FILE_NAME_LENGTH];
+	char geneFileName2[FILE_NAME_LENGTH];
+	FullMappingInfoLink *miL = getMem(_msf_seqListSize * sizeof(FullMappingInfoLink));
+	FullMappingInfoLink *miR = getMem(_msf_seqListSize * sizeof(FullMappingInfoLink));
+	if(fileName1 != NULL && fileName2 != NULL)
+	{
+		fp1 = fileOpen(fileName1, "r");
+		fp2 = fileOpen(fileName2, "r");
+		index = strstr(fileName1, "__");
+		strncpy(geneFileName1, index + 2 * sizeof(char), strstr(index + 2, "__") - index - 2);
+		geneFileName1[strstr(index + 2, "__") - index - 2] = '\0';
+index = strstr(fileName2, "__");
+strncpy(geneFileName2, index + 2 * sizeof(char), strstr(index + 2, "__") - index - 2);
+		geneFileName2[strstr(index + 2, "__") - index - 2] = '\0';
+		for(i = 0; i < _msf_seqListSize / 2; i++)
+		{
+			fread(&size1, sizeof(int), 1, fp1);
+			fread(&size2, sizeof(int), 1, fp2);
+			miL[i].mi = getMem(size1 * sizeof(FullMappingInfo) );
+			miR[i].mi = getMem(size2 * sizeof(FullMappingInfo) );
+			miL[i].size = size1;
+			miR[i].size = size2;
+			for(j = 0; j < size1; j++)
+			{
+				fread(&(miL[i].mi[j].loc), sizeof(int), 1, fp1);
+				fread (&(miL[i].mi[j].err), sizeof(int), 1, fp1);
+fread (&(miL[i].mi[j].cigarSize), sizeof(int), 1, fp1);
+fread ((miL[i].mi[j].cigar), sizeof(char), miL[i].mi[j].cigarSize+1, fp1);
+fread (&(miL[i].mi[j].mdSize), sizeof(int), 1, fp1);
+fread ((miL[i].mi[j].md), sizeof(char), miL[i].mi[j].mdSize+1, fp1);
+				miL[i].mi[j].dir = 1;
+				if(miL[i].mi[j].loc < 1)
+				{
+					miL[i].mi[j].loc *= -1;
+					miL[i].mi[j].dir = -1;
+				}
+			}
+			for(k = 0; k < size2; k++)
+			{
+fread(&(miR[i].mi[k].loc), sizeof(int), 1, fp2);
+				fread (&(miR[i].mi[k].err), sizeof(int), 1, fp2);
+fread (&(miR[i].mi[k].cigarSize), sizeof(int), 1, fp2);
+fread ((miR[i].mi[k].cigar), sizeof(char), miR[i].mi[k].cigarSize+1, fp2);
+fread (&(miR[i].mi[k].mdSize), sizeof(int), 1, fp2);
+fread ((miR[i].mi[k].md), sizeof(char), miR[i].mi[k].mdSize+1, fp2);
+	 			miR[i].mi[k].dir = 1;
+if(miR[i].mi[k].loc < 1)
+{
+miR[i].mi[k].loc *= -1;
+miR[i].mi[k].dir = -1;
+}
+			}
+			if(_msf_readHasConcordantMapping[i] == 0 && size1 != 0 && size2 != 0 && (size1 * size2 < MAX_TRANS_CHROMOSAL_OUTPUT))
+			{
+				int d1 = 0;
+				int d2 = 0;
+				char *seq, *qual;
+				char *seq1, *seq2, *rseq1, *rseq2, *qual1, *qual2;
+		                char rqual1[SEQ_LENGTH+1], rqual2[SEQ_LENGTH+1];
+		                rqual1[SEQ_LENGTH] = rqual2[SEQ_LENGTH] = '\0';
+			        seq1 = _msf_seqList[i*2].seq;
+		                rseq1 = _msf_seqList[i*2].rseq;
+			        qual1 = _msf_seqList[i*2].qual;
+		                reverse(_msf_seqList[i*2].qual, rqual1, SEQ_LENGTH);
+		        seq2 = _msf_seqList[i*2+1].seq;
+		                rseq2 = _msf_seqList[i*2+1].rseq;
+		qual2 = _msf_seqList[i*2+1].qual;
+		                reverse(_msf_seqList[i*2+1].qual, rqual2, SEQ_LENGTH);
+				for(j = 0; j < size1; j++)
+				{
+					d1 = (miL[i].mi[j].dir == -1)?1:0;
+if ( d1 )
+					{
+						seq = rseq1;
+						qual = rqual1;
+					}
+					else
+					{
+						seq = seq1;
+						qual = qual1;
+					}
+					for(k = 0; k < size2; k++)
+					{
+	                        d2 = (miR[i].mi[k].dir == -1)?1:0;
+_msf_output.POS                 = miL[i].mi[j].loc;
+_msf_output.MPOS                = miR[i].mi[k].loc;
+_msf_output.FLAG                = 0;
+_msf_output.ISIZE               = 0;
+_msf_output.SEQ                 = seq,
+_msf_output.QUAL                = qual;
+_msf_output.QNAME               = _msf_seqList[i*2].name;
+_msf_output.RNAME               = geneFileName1;
+_msf_output.MAPQ                = 255;
+_msf_output.CIGAR               = miL[i].mi[j].cigar;
+_msf_output.MRNAME              = "=";
+_msf_output.optSize     = 2;
+_msf_output.optFields   = _msf_optionalFields;
+_msf_optionalFields[0].tag = "NM";
+_msf_optionalFields[0].type = 'i';
+_msf_optionalFields[0].iVal = miL[i].mi[j].err;
+_msf_optionalFields[1].tag = "MD";
+_msf_optionalFields[1].type = 'Z';
+_msf_optionalFields[1].sVal = miL[i].mi[j].md;
+						if ( d2 )
+{
+seq = rseq2;
+qual = rqual2;
+}
+else
+{
+seq = seq2;
+qual = qual2;
+}
+outputSAM(fp_out, _msf_output);
+					 	_msf_output.POS                 = miR[i].mi[k].loc;
+_msf_output.MPOS                = miL[i].mi[j].loc;
+_msf_output.FLAG                = 0;
+_msf_output.ISIZE               = 0;
+_msf_output.SEQ                 = seq,
+_msf_output.QUAL                = qual;
+_msf_output.QNAME               = _msf_seqList[i*2+1].name;
+_msf_output.RNAME               = geneFileName2;
+_msf_output.MAPQ                = 255;
+_msf_output.CIGAR               = miR[i].mi[k].cigar;
+_msf_output.MRNAME              = "=";
+_msf_output.optSize     = 2;
+_msf_output.optFields   = _msf_optionalFields;
+_msf_optionalFields[0].tag = "NM";
+_msf_optionalFields[0].type = 'i';
+_msf_optionalFields[0].iVal = miR[i].mi[k].err;
+_msf_optionalFields[1].tag = "MD";
+_msf_optionalFields[1].type = 'Z';
+_msf_optionalFields[1].sVal = miR[i].mi[k].md;
+outputSAM(fp_out, _msf_output);
+					}
+				}
+			}
+		}
+	}
+	for(i = 0; i < _msf_seqListSize / 2; i++)
+	{
+		freeMem(miL[i].mi, miL[i].size * sizeof(FullMappingInfo));
+		freeMem(miR[i].mi, miR[i].size * sizeof(FullMappingInfo));
+	}
+	freeMem(miL, _msf_seqListSize * sizeof(FullMappingInfoLink));
+freeMem(miR, _msf_seqListSize * sizeof(FullMappingInfoLink));
+	fclose(fp1);
+	fclose(fp2);
+}
+/*
+	if flag is 1 it will output all the possible trans chromsal mapping
+	otherwise only tmp file will be delete
+*/
+void outputAllTransChromosal(int flag)
+{
+	int i = 0;
+	int j = 0;
+	int k = 0;
+	int l = 0;
+	FILE *fp_out = NULL;
+	char fname1[200];
+	if(flag)
+	{
+	        fp_out = fileOpen(fname1, "w");
+	        sprintf(fname1, "%s%s_TRANSCHROMOSOMAL", mappingOutputPath, mappingOutput);
+	//	for(i = 0; i < _msf_maxFile; i++)
+	//	{
+			i  = 0;
+			for(j = i+1; j < _msf_maxFile; j++)
+			{
+				if(i != j)
+				{
+					for(k = 0; k < _msf_fileCount[i]; k++)
+					{
+						for(l = 0; l < _msf_fileCount[j]; l++)
+						{
+							outputTransChromosal(_msf_fileName[i][k][0], _msf_fileName[j][l][1], fp_out);
+						}// for l
+					}// for k
+				}// if
+			}// for j
+	//	} //for i
+	}
+	for(i = 0; i < _msf_maxFile; i++)
+	{
+		for(j = 0; j < _msf_fileCount[i]; j++)
+		{
+			unlink(_msf_fileName[i][j][0]);
+			unlink(_msf_fileName[i][j][1]);
+		}
+	}
+	if(flag)
+	  fclose(fp_out);
+}

Mercurial > repos > calkan > mrfast

comparison mrfast-2.1.0.4/MrFAST.c @ 0:7b3dc85dc7fd