# HG changeset patch
# User calkan
# Date 1329838768 18000
# Node ID ec628ba338784b868b166685fbfaac81b739225a

Uploaded source code for mrsFAST

diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/CommandLineParser.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/CommandLineParser.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,373 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-01-29
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <getopt.h>
+#include <string.h>
+#include <ctype.h>
+#include "Common.h"
+#include "CommandLineParser.h"
+
+int						uniqueMode=1;
+int						indexingMode;
+int						searchingMode;
+int						bisulfiteMode;
+int						pairedEndMode;
+int						pairedEndDiscordantMode;
+int						pairedEndProfilingMode;
+int						seqCompressed;
+int						outCompressed;
+int						cropSize = 0;
+int						progressRep = 0;
+int						minPairEndedDistance=-1;
+int						maxPairEndedDistance=-1;
+int						minPairEndedDiscordantDistance=-1;
+int						maxPairEndedDiscordantDistance=-1;
+char					*seqFile1;
+char					*seqFile2;
+char					*mappingOutput = "output";
+char					*mappingOutputPath = "";
+char					*unmappedOutput = "";
+char					fileName[1000][2][FILE_NAME_LENGTH];
+int						fileCnt;
+unsigned char			errThreshold=2;
+unsigned char			maxHits=0;
+unsigned char			WINDOW_SIZE = 12;
+unsigned int			CONTIG_SIZE;
+unsigned int			CONTIG_MAX_SIZE;
+
+void printHelp();
+
+int parseCommandLine (int argc, char *argv[])
+{
+
+	int o;
+	int index;
+	char *fastaFile = NULL;
+	char *fastaOutputFile = NULL;
+	char *indexFile = NULL;
+	char *batchFile = NULL ;
+	int  batchMode = 0;
+	static struct option longOptions[] = 
+	{
+
+		//		{"bs",				no_argument,		&bisulfiteMode,		1},
+		{"pe",				no_argument,		&pairedEndMode,		1},
+		{"discordant-vh",	no_argument,		&pairedEndDiscordantMode,	1},
+		{"profile",			no_argument, 		&pairedEndProfilingMode,	1},
+		{"seqcomp",			no_argument,		&seqCompressed,		1},
+		{"outcomp",			no_argument,		&outCompressed,		1},
+		{"progress",		no_argument,		&progressRep,		1},
+		{"index",			required_argument,	0, 					'i'},
+		{"search",			required_argument,	0,					's'},
+		{"help",			no_argument,		0,					'h'},
+		{"version",			no_argument,		0,					'v'},
+		{"seq",				required_argument,	0,					'x'},
+		{"seq1",			required_argument,	0,					'x'},
+		{"seq2",			required_argument,	0,					'y'},
+		{"ws",				required_argument,  0,					'w'},
+		{"min",				required_argument,  0,					'l'},
+		{"max",				required_argument,  0,					'm'},
+		{"crop",			required_argument,  0,					'c'}
+
+	};
+
+
+
+	while ( (o = getopt_long ( argc, argv, "f:i:u:o:s:e:n:bhv", longOptions, &index))!= -1 )
+	{
+		switch (o)
+		{
+			case 'i':
+				indexingMode = 1;
+				fastaFile = optarg;
+				break;
+			case 's':
+				searchingMode = 1;
+				fastaFile = optarg;
+				break;
+			case 'b':
+				batchMode = 1;
+				break;
+			case 'c': 
+				cropSize = atoi(optarg);
+				break;
+			case 'w':
+				WINDOW_SIZE = atoi(optarg);
+				break;
+			case 'x':
+				seqFile1 = optarg;
+				break;
+			case 'y':
+				seqFile2 = optarg;
+				break;
+			case 'u':
+				unmappedOutput = optarg;
+				break;
+			case 'o':
+				mappingOutput = getMem(FILE_NAME_LENGTH);
+				mappingOutputPath = getMem(FILE_NAME_LENGTH);
+				stripPath (optarg, &mappingOutputPath, &mappingOutput);
+				break;
+			case 'n':
+				maxHits = atoi(optarg);
+				break;
+			case 'e':
+				errThreshold = atoi(optarg);
+				break;
+			case 'l':
+				minPairEndedDistance = atoi(optarg);
+				break;
+			case 'm':
+				maxPairEndedDistance = atoi(optarg);
+				break;					
+			case 'h':
+				printHelp();
+				return 0;
+				break;
+			case 'v':
+				fprintf(stdout, "%s.%s\n", versionNumber, versionNumberF);
+				return 0;
+				break;
+		}
+
+	}
+
+	if (indexingMode + searchingMode != 1)
+	{
+		fprintf(stdout, "ERROR: Indexing / Searching mode should be selected\n");
+		return 0;
+	}
+
+	if (WINDOW_SIZE > 14 || WINDOW_SIZE < 8)
+	{
+		fprintf(stdout, "ERROR: Window size should be in [8..14]\n");
+		return 0;
+	}
+
+
+	if ( indexingMode )
+	{
+		CONTIG_SIZE		= 15000000;
+		CONTIG_MAX_SIZE	= 40000000;
+
+		if (batchMode)
+		{
+			batchFile = fastaFile;
+			fastaFile = NULL;
+		}
+
+		if (batchFile == NULL && fastaFile == NULL)
+		{
+			fprintf(stdout, "ERROR: Reference(s) should be indicated for indexing\n");
+			return 0;
+		}
+
+		if (pairedEndDiscordantMode)
+		{
+			fprintf(stdout, "ERROR: --discordant cannot be used in indexing mode. \n");
+			return 0;
+		}
+
+	}
+
+
+	if ( searchingMode )
+	{
+		CONTIG_SIZE		= 300000000;
+		CONTIG_MAX_SIZE	= 300000000;
+
+
+		if (batchMode)
+		{
+			batchFile = fastaFile;
+			fastaFile = NULL;
+		}
+
+		if (batchFile == NULL && fastaFile == NULL)
+		{
+			fprintf(stdout, "ERROR: Index File(s) should be indiciated for searching\n");
+			return 0;
+		}
+
+		if (seqFile1 == NULL && seqFile2 == NULL)
+		{
+			fprintf(stdout, "ERROR: Please indicate a sequence file for searching.\n");
+			return 0;
+		}
+
+
+		if (!pairedEndMode && !bisulfiteMode && seqFile2 != NULL)
+		{
+			fprintf(stdout, "ERROR: Second File can be indicated in pairedend/bisulfite mode\n");
+			return 0;
+		}
+
+		if (pairedEndMode && (minPairEndedDistance <0 || maxPairEndedDistance < 0 || minPairEndedDistance > maxPairEndedDistance))
+		{
+			fprintf(stdout, "ERROR: Please enter a valid range for pairedend sequences.\n");
+			return 0;
+		}
+
+		if (pairedEndMode && seqFile1 == NULL)
+		{
+			fprintf(stdout, "ERROR: Please indicate the first file for pairedend search.\n");
+			return 0;
+		}
+
+		if (!pairedEndMode && pairedEndDiscordantMode)
+		{
+			fprintf(stdout, "ERROR: --discordant should be used with --pe");
+			return 0;
+		}
+
+		if (!pairedEndMode && pairedEndProfilingMode)
+		{
+			fprintf(stdout, "ERROR: --profile should be used with --pe");
+			return 0;
+		}
+	}
+
+	int i = 0;
+
+
+	if (batchMode)
+	{
+		FILE *fp = fileOpen(batchFile, "r");
+
+		if (fp == NULL)
+			return 0;
+
+		fileCnt  = 0;
+
+		while ( fgets(fileName[fileCnt][0], FILE_NAME_LENGTH, fp))
+		{
+			for (i = strlen(fileName[fileCnt][0])-1; i>=0; i--)
+				if ( !isspace(fileName[fileCnt][0][i]))
+					break;
+			fileName[fileCnt][0][i+1] = '\0';
+
+			if (strcmp(fileName[fileCnt][0], "") != 0)
+			{
+				if (bisulfiteMode)
+					sprintf(fileName[fileCnt][1], "%s.bsindex", fileName[fileCnt][0]); 
+				else
+					sprintf(fileName[fileCnt][1], "%s.index", fileName[fileCnt][0]); 
+				fileCnt++;
+			}
+		}
+	}
+	else
+	{
+		sprintf(fileName[fileCnt][0], "%s", fastaFile);
+		if (bisulfiteMode)
+			sprintf(fileName[fileCnt][1], "%s.bsindex", fileName[fileCnt][0]); 
+		else
+			sprintf(fileName[fileCnt][1], "%s.index", fileName[fileCnt][0]); 
+		fileCnt++;
+	}
+
+
+	if (pairedEndProfilingMode)
+	{
+
+		minPairEndedDistance = 0;
+		maxPairEndedDistance = 300000000;
+
+	}
+
+	if (pairedEndDiscordantMode)
+	{
+		minPairEndedDiscordantDistance = minPairEndedDistance;
+		maxPairEndedDiscordantDistance = maxPairEndedDistance;
+
+		minPairEndedDistance = 0;
+		maxPairEndedDistance = 300000000;
+	}
+
+	return 1;
+}
+
+
+void printHelp()
+{
+	char *errorType;
+	if (mrFAST)
+	{
+		fprintf(stdout,"mrFAST : Micro-Read Fast Alignment Search Tool.\n\n");
+		fprintf(stdout,"Usage: mrFAST [options]\n\n");
+		errorType="edit distance";
+	}
+	else
+	{
+		fprintf(stdout,"mrsFAST : Micro-Read Substitutions (only) Fast Alignment Search Tool.\n\n");
+		fprintf(stdout,"mrsFAST is a cache oblivious read mapping tool. mrsFAST capable of mapping\n");
+		fprintf(stdout,"single and paired end reads to the reference genome. Bisulfite treated \n");
+		fprintf(stdout,"sequences are not supported in this version. By default mrsFAST reports  \n");
+		fprintf(stdout,"the output in SAM format.\n\n");
+		fprintf(stdout,"Usage: mrsFAST [options]\n\n");
+		errorType="hamming distance";
+	}
+
+	fprintf(stdout,"General Options:\n");
+	fprintf(stdout," -v|--version\t\tCurrent Version.\n");
+	fprintf(stdout," -h\t\t\tShows the help file.\n");
+	fprintf(stdout,"\n\n");
+
+	fprintf(stdout,"Indexing Options:\n");
+	fprintf(stdout," --index [file]\t\tGenerate an index from the specified fasta file. \n");
+	fprintf(stdout," -b\t\t\tIndicates the indexing will be done in batch mode.\n\t\t\tThe file specified in --index should contain the \n\t\t\tlist of fasta files.\n");
+	fprintf(stdout," -ws [int]\t\tSet window size for indexing (default:12-min:8 max:14).\n");
+	//	fprintf(stdout," -bs \t\t\tBisulfite mode.");
+	fprintf(stdout,"\n\n");
+
+	fprintf(stdout,"Searching Options:\n");
+	fprintf(stdout," --search [file]\tSearch the specified genome. Index file should be \n\t\t\tin same directory as the fasta file.\n");
+	fprintf(stdout," -b\t\t\tIndicates the mapping will be done in batch mode. \n\t\t\tThe file specified in --search should contain the \n\t\t\tlist of fasta files.\n");
+	fprintf(stdout," --pe \t\t\tSearch will be done in Pairedend mode.\n");
+	//	fprintf(stdout," --bs \t\t\tSearch will be done in Bisulfite mode.\n");
+	fprintf(stdout," --seq [file]\t\tInput sequences in fasta/fastq format [file]. If \n\t\t\tpairend reads are interleaved, use this option.\n");
+	fprintf(stdout," --seq1 [file]\t\tInput sequences in fasta/fastq format [file] (First \n\t\t\tfile). Use this option to indicate the first file of \n\t\t\tpair-end reads. You can use this option alone in \n\t\t\tbisulfite mode. \n");
+	fprintf(stdout," --seq2 [file]\t\tInput sequences in fasta/fastq format [file] (Second \n\t\t\tfile). Use this option to indicate the second file of \n\t\t\tpair-end reads. You can use this option alone in \n\t\t\tbisulfite mode. \n");
+	fprintf(stdout," -o [file]\t\tOutput of the mapped sequences. The default is output.\n");
+	fprintf(stdout," --seqcomp \t\tIndicates that the input sequences are compressed(gz).\n");
+	fprintf(stdout," --outcomp \t\tIndicates that output file should be compressed(gz).\n");
+	//	fprintf(stdout," -u [file]\t\tSave unmapped sequences to the [file] in fasta format.\n");
+	fprintf(stdout," -n [int]\t\tMaximum number of locations reported for a sequence \n\t\t\t(default 0, all mappings). \n");
+	fprintf(stdout," -e [int]\t\t%s (default 2).\n", errorType);
+	fprintf(stdout," --min [int]\t\tMin inferred distance allowed between two pairend sequences.\n");
+	fprintf(stdout," --max [int]\t\tMax inferred distance allowed between two pairend sequences.\n");
+}
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/CommandLineParser.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/CommandLineParser.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,42 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __COMMAND_LINE_PARSER__
+#define __COMMAND_LINE_PARSER__
+
+int parseCommandLine (int argc, char *argv[]);
+
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Common.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Common.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,168 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-02-01
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <zlib.h>
+#include <string.h>
+#include "Common.h"
+
+
+unsigned short 			SEQ_LENGTH = 0;
+unsigned short 			QUAL_LENGTH = 0;
+long long				memUsage = 0;
+/**********************************************/
+FILE *fileOpen(char *fileName, char *mode)
+{
+	FILE *fp;
+	fp = fopen (fileName, mode);
+	if (fp == NULL)
+	{
+		fprintf(stdout, "Error: Cannot Open the file %s\n", fileName);
+		fflush(stdout);
+		exit(0);
+	}
+	return fp;
+}
+/**********************************************/
+gzFile fileOpenGZ(char *fileName, char *mode)
+{
+	gzFile gzfp;
+	gzfp = gzopen (fileName, mode);
+	if (gzfp == Z_NULL)
+	{
+		fprintf(stdout, "Error: Cannot Open the file %s\n", fileName);
+		fflush(stdout);
+		exit(0);
+	}
+	return gzfp;
+}
+/**********************************************/
+double getTime(void)
+{
+	struct timeval t;
+	gettimeofday(&t, NULL);
+	return t.tv_sec+t.tv_usec/1000000.0;
+}
+
+/**********************************************/
+char reverseCompleteChar(char c)
+{
+	char ret;
+	switch (c)
+	{
+		case 'A': 
+					ret = 'T';
+					break;
+		case 'T':
+					ret = 'A';
+					break;
+		case 'C':	
+					ret = 'G';
+					break;
+		case 'G':
+					ret = 'C';
+					break;
+		default:
+					ret = 'N';
+					break;
+	}
+	return ret;
+}
+/**********************************************/
+void reverseComplete (char *seq, char *rcSeq , int length)
+{
+	int i;
+	for (i=0; i<length; i++)
+	{
+		rcSeq[i]=reverseCompleteChar (seq[length-1-i]) ;
+	}
+}
+/**********************************************/
+void * getMem(size_t size)
+{
+	memUsage+=size;
+	return malloc(size);
+}
+/**********************************************/
+void freeMem(void *ptr, size_t size)
+{
+	memUsage-=size;
+	free(ptr);
+}
+/**********************************************/
+double getMemUsage()
+{
+	return memUsage/1048576.0;
+}
+/**********************************************/
+void reverse (char *seq, char *rcSeq , int length)
+{
+	int i;
+	for (i=0; i<length; i++)
+	{
+		rcSeq[i]=seq[length-1-i];
+	}
+}
+/**********************************************/
+void stripPath(char *full, char **path, char **fileName)
+{
+	int i;
+	int pos = -1;
+
+	for (i=strlen(full)-1; i>=0; i--)
+	{
+		if (full[i]=='/')
+		{
+			pos = i;
+			break;
+		}
+
+	}
+
+	if (pos != -1)
+	{
+		sprintf(*fileName, "%s%c", (full+pos+1), '\0');
+		full[pos+1]='\0';
+		sprintf(*path,"%s%c", full, '\0');
+	}
+	else
+	{
+		sprintf(*fileName, "%s%c", full, '\0');
+		sprintf(*path,"%c", '\0');
+	}
+}
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Common.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Common.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,98 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __COMMON__
+#define __COMMON__
+
+#include <zlib.h>
+
+#define SEQ_MAX_LENGTH		200			// Seq Max Length
+#define CONTIG_OVERLAP		200 		// No. of characters overlapped between contings
+#define CONTIG_NAME_SIZE	200			// Contig name max size
+#define FILE_NAME_LENGTH	400			// Filename Max Length
+#define DISCORDANT_CUT_OFF  800
+
+extern unsigned int		CONTIG_SIZE;
+extern unsigned int		CONTIG_MAX_SIZE;
+
+
+extern unsigned char	WINDOW_SIZE				;		// WINDOW SIZE for indexing/searching
+extern unsigned short	SEQ_LENGTH;						// Sequence(read) length
+extern unsigned short	QUAL_LENGTH;
+
+extern char				*versionNumber;
+extern char				*versionNumberF;
+extern unsigned char	mrFAST;
+
+
+extern int				uniqueMode;
+extern int				indexingMode;
+extern int				searchingMode;
+extern int				bisulfiteMode;
+extern int				pairedEndMode;
+extern int				pairedEndDiscordantMode;
+extern int				pairedEndProfilingMode;
+extern int				seqCompressed;
+extern int				outCompressed;
+extern int				cropSize;
+extern int				progressRep;
+extern char 			*seqFile1;
+extern char				*seqFile2;
+extern char				*seqUnmapped;
+extern char				*mappingOutput;
+extern char 			*mappingOutputPath;
+extern char				*unmappedOutput;
+extern unsigned char	seqFastq;
+extern unsigned char	errThreshold;
+extern unsigned char	maxHits;	
+extern int				minPairEndedDiscordantDistance;
+extern int				maxPairEndedDiscordantDistance;
+extern int				minPairEndedDistance;
+extern int				maxPairEndedDistance;
+extern char				fileName[1000][2][FILE_NAME_LENGTH];
+extern int				fileCnt;
+extern long long		memUsage;
+
+FILE	* fileOpen(char *fileName, char *mode);
+gzFile	fileOpenGZ(char *fileName, char *mode);
+double	getTime(void);
+void	reverseComplete (char *seq, char *rcSeq , int length);
+void	* getMem(size_t size);
+void	freeMem(void * ptr, size_t size);
+double	getMemUsage();
+void 	reverse (char *seq, char *rcSeq , int length);
+void 	stripPath(char *full, char **path, char **fileName);
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/HashTable.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/HashTable.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,428 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "Common.h"
+#include "RefGenome.h"
+#include "HashTable.h"
+/**********************************************/
+FILE		*_ih_fp					= NULL;
+IHashTable	*_ih_hashTable			= NULL;
+int 		_ih_maxHashTableSize	= 0;
+char		*_ih_refGen				= NULL;
+char		*_ih_refGenName			= NULL;
+long long	_ih_memUsage			= 0;
+int			_ih_refGenOff			= 0;
+/**********************************************/
+
+int hashVal(char *seq)
+{
+	int i=0;
+	int val=0, numericVal=0;
+
+	while(i<WINDOW_SIZE)
+	{
+		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;
+}
+/**********************************************/
+void freeIHashTableContent(IHashTable *hashTable, unsigned int maxSize)
+{
+	int i=0;
+	for (i=0; i<maxSize; i++)
+	{
+		if (hashTable[i].locs != NULL)
+		{
+			freeMem(hashTable[i].locs, (hashTable[i].locs[0]+1)*(sizeof(unsigned int)));
+			hashTable[i].locs = NULL;
+		}
+	}
+}
+/**********************************************/
+void initSavingIHashTable(char *fileName)
+{
+	int tmp;
+	_ih_fp = fileOpen(fileName, "w");
+	unsigned char bIndex = 0;				// Bisulfite Index
+	
+	// First Two bytes are indicating the type of the index & window size
+	tmp = fwrite(&bIndex, sizeof(bIndex), 1, _ih_fp);
+	tmp = fwrite(&WINDOW_SIZE, sizeof(WINDOW_SIZE), 1, _ih_fp);
+	
+}
+/**********************************************/
+void finalizeSavingIHashTable()
+{
+	fclose(_ih_fp);
+}
+/**********************************************/
+void saveIHashTable(IHashTable *hashTable,  unsigned int size, unsigned int maxSize, char *refGen, char *refGenName, int refGenOffset)
+{
+	int tmp;
+	
+	// Every Chunk starts with a byte indicating whether it has extra info;
+	unsigned char extraInfo = 0;
+	tmp = fwrite (&extraInfo, sizeof(extraInfo), 1, _ih_fp);
+
+	short len = strlen(refGenName);
+	tmp = fwrite(&len, sizeof(len), 1, _ih_fp);
+	tmp = fwrite(refGenName, sizeof(char), len, _ih_fp);
+
+	tmp = fwrite(&refGenOffset, sizeof(refGenOffset), 1, _ih_fp);
+	
+	unsigned int refGenLength = strlen(refGen);
+	tmp = fwrite(&refGenLength, sizeof(refGenLength), 1, _ih_fp);
+	tmp = fwrite(refGen, sizeof(char), refGenLength, _ih_fp);
+
+	tmp = fwrite(&size, sizeof(size), 1, _ih_fp);
+
+	int i=0,j=0;
+	unsigned char cnt=0;
+	for (i=0; i<maxSize; i++)
+	{
+		if (hashTable[i].locs != NULL)
+		{
+			tmp = fwrite(&i, sizeof(i), 1, _ih_fp);
+
+			if (hashTable[i].locs[0] < 250)
+			{
+				cnt = hashTable[i].locs[0];
+				tmp = fwrite(&cnt, sizeof(cnt), 1, _ih_fp);
+			}
+			else
+			{
+				cnt =0;
+				tmp = fwrite (&cnt, sizeof(cnt), 1, _ih_fp);
+				tmp = fwrite (&(hashTable[i].locs[0]), sizeof(hashTable[i].locs[0]), 1, _ih_fp);
+			}
+
+			for (j=1; j<=hashTable[i].locs[0]; j++)
+				tmp = fwrite(&(hashTable[i].locs[j]), sizeof(hashTable[i].locs[j]), 1, _ih_fp);
+		}
+	}
+}
+/**********************************************/
+unsigned int addIHashTableLocation(IHashTable *hashTable, int hv, int location)
+{
+	unsigned int	sizeInc				= 0;
+
+	if (hashTable[hv].locs == NULL)
+	{
+		sizeInc = 1;
+		hashTable[hv].locs = getMem (sizeof(unsigned int)*2);
+		hashTable[hv].locs[0]=1;
+		hashTable[hv].locs[1]=location;
+	}
+	else
+	{
+		int size = hashTable[hv].locs[0];
+		int i;
+		unsigned int *tmp = getMem( (size + 2) * sizeof(unsigned int) );
+
+		for (i = 0; i <= size; i++)
+		{
+			tmp[i] = hashTable[hv].locs[i];
+		}
+		size++;
+		tmp[0] = size;
+		tmp[size] = location;
+		
+		freeMem(hashTable[hv].locs, (hashTable[hv].locs[0]*(sizeof(unsigned int))));
+		hashTable[hv].locs = tmp;
+	}
+	return sizeInc;
+}
+/**********************************************/
+void generateIHashTable(char *fileName, char *indexName)
+{
+	double          startTime           = getTime();
+	unsigned int	hashTableSize		= 0;
+	unsigned int 	hashTableMaxSize	= pow(4, WINDOW_SIZE);
+	IHashTable		*hashTable			= getMem(sizeof(IHashTable)*hashTableMaxSize);
+	char 			*refGenName;
+	char			*refGen;
+	int				refGenOff			= 0;
+	int i, hv, l, flag;
+
+
+	for ( i = 0; i < hashTableMaxSize; i++)
+	{
+		hashTable[i].locs = NULL;
+	}
+
+	//Loading Fasta File
+	if (!initLoadingRefGenome(fileName))
+		return;		
+	initSavingIHashTable(indexName);
+	
+	fprintf(stdout, "Generating Index from %s", fileName);
+	fflush(stdout);
+
+	char *prev = getMem (CONTIG_NAME_SIZE);
+	prev[0]='\0';
+
+	do
+	{
+		flag = 	 loadRefGenome (&refGen, &refGenName, &refGenOff);	
+
+		if ( strcmp(prev, refGenName) != 0)
+		{
+			fprintf(stdout, "\n - %s ", refGenName);
+			fflush(stdout);
+			sprintf(prev, "%s", refGenName);
+		}
+		else
+		{
+			fprintf(stdout, ".");
+			fflush(stdout);
+		}
+		
+		l = strlen(refGen) - WINDOW_SIZE;
+
+		for (i=0; i < l; i++)
+		{
+			hv = hashVal(refGen + i);
+			if (hv != -1)
+			{
+				hashTableSize += addIHashTableLocation (hashTable, hv, i+1);
+			}
+		}
+
+		saveIHashTable(hashTable, hashTableSize, hashTableMaxSize, refGen, refGenName, refGenOff);
+		freeIHashTableContent(hashTable, hashTableMaxSize);
+		hashTableSize = 0;
+	} while (flag);
+
+	freeMem(prev, CONTIG_NAME_SIZE);
+	freeMem(hashTable, sizeof(IHashTable)*hashTableMaxSize);
+
+	finalizeLoadingRefGenome();
+	finalizeSavingIHashTable();
+
+	fprintf(stdout, "\nDONE in %0.2fs!\n", (getTime()-startTime));
+}
+
+/**********************************************/
+void finalizeLoadingIHashTable()
+{
+	freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
+	freeMem(_ih_hashTable, sizeof(IHashTable)* _ih_maxHashTableSize);
+	freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
+	freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
+	fclose(_ih_fp);
+}
+
+/**********************************************/
+int  loadIHashTable(double *loadTime)
+{
+	double startTime = getTime();
+	unsigned char extraInfo = 0;
+	short len;
+	unsigned int refGenLength;
+	unsigned int hashTableSize;
+	unsigned int tmpSize;
+	int tmp;
+	int i=0,j=0;
+
+	if ( fread(&extraInfo, sizeof(extraInfo), 1, _ih_fp) != sizeof(extraInfo) )
+		return 0;
+
+	freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
+	freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
+	freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
+
+	// Reading Chr Name
+	tmp = fread(&len, sizeof(len), 1, _ih_fp);
+	_ih_refGenName = getMem(sizeof(char)* (len+1));
+	tmp = fread(_ih_refGenName, sizeof(char), len, _ih_fp);
+	_ih_refGenName [len] ='\0';
+
+	tmp = fread(&_ih_refGenOff, sizeof (_ih_refGenOff), 1, _ih_fp);
+
+	// Reading Size and Content of Ref Genome
+	tmp = fread(&refGenLength, sizeof(refGenLength), 1, _ih_fp);
+	_ih_refGen = getMem(sizeof(char)*(refGenLength+1));
+	tmp = fread(_ih_refGen, sizeof(char), refGenLength, _ih_fp);
+	_ih_refGen[refGenLength]='\0';
+
+	//Reading Hashtable Size and Content
+	tmp = fread(&hashTableSize, sizeof(hashTableSize), 1, _ih_fp);
+	
+	unsigned int hv;
+	unsigned char cnt=0;
+	for (i=0; i<hashTableSize; i++)
+	{
+		tmp = fread(&hv, sizeof(hv), 1, _ih_fp);
+		tmp = fread(&cnt, sizeof(cnt), 1, _ih_fp);
+
+		if (cnt>0)
+		{
+			tmpSize = cnt;
+		}
+		else
+		{
+			tmp = fread(&tmpSize, sizeof(tmpSize), 1, _ih_fp);
+		}
+	
+		_ih_hashTable[hv].locs = getMem( sizeof(unsigned int)* (tmpSize+1) );
+		_ih_hashTable[hv].locs[0] = tmpSize;
+
+		for (j=1; j<=tmpSize; j++)
+			tmp = fread(&(_ih_hashTable[hv].locs[j]), sizeof(_ih_hashTable[hv].locs[j]), 1, _ih_fp);
+	}
+	*loadTime = getTime()-startTime;
+
+	return 1;
+}
+/**********************************************/
+unsigned int *getIHashTableCandidates(int hv)
+{
+	if ( hv != -1 )
+		return _ih_hashTable[hv].locs;
+	else 
+		return NULL;
+}
+/**********************************************/
+/**********************************************/
+/**********************************************/
+void configHashTable()
+{
+	if (WINDOW_SIZE <= 14)
+	{
+		generateHashTable = &generateIHashTable;
+		loadHashTable = &loadIHashTable;
+		finalizeLoadingHashTable = &finalizeLoadingIHashTable;
+		getCandidates = &getIHashTableCandidates;
+	}
+	else
+	{
+
+
+	}
+}
+/**********************************************/
+int initLoadingHashTable(char *fileName)
+{
+	int i;	
+	unsigned char bsIndex;
+	int tmp; 
+
+	_ih_fp = fileOpen(fileName, "r");	
+
+	if (_ih_fp == NULL)
+		return 0;
+
+	tmp = fread(&bsIndex, sizeof(bsIndex), 1, _ih_fp);
+	if (bsIndex)
+	{
+		fprintf(stdout, "Error: Wrong Type of Index indicated");
+		return 0;
+	}
+	
+	tmp = fread(&WINDOW_SIZE, sizeof(WINDOW_SIZE), 1, _ih_fp);
+
+	configHashTable();
+
+	if (_ih_maxHashTableSize != pow(4, WINDOW_SIZE))
+	{
+
+		if (_ih_hashTable != NULL)
+		{
+			freeIHashTableContent(_ih_hashTable, _ih_maxHashTableSize);
+			freeMem(_ih_hashTable, sizeof(IHashTable)* _ih_maxHashTableSize);
+			freeMem(_ih_refGen, strlen(_ih_refGen)+1) ;
+			freeMem(_ih_refGenName, strlen(_ih_refGenName)+1);
+		}
+
+		_ih_maxHashTableSize = pow(4, WINDOW_SIZE);
+
+		_ih_hashTable = getMem (sizeof(IHashTable) * _ih_maxHashTableSize);
+		for (i=0; i<_ih_maxHashTableSize; i++)
+			_ih_hashTable[i].locs = NULL;
+		_ih_refGen = getMem(1);
+		_ih_refGen[0]='\0';
+		_ih_refGenName = getMem(1);
+		_ih_refGenName[0] = '\0';
+	}
+
+	return 1;
+}
+/**********************************************/
+char *getRefGenome()
+{
+	return _ih_refGen;
+
+}
+/**********************************************/
+char *getRefGenomeName()
+{
+	return _ih_refGenName;
+
+}
+/**********************************************/
+int getRefGenomeOffset()
+{
+	return _ih_refGenOff;
+
+}
+/**********************************************/
+HashTable *getHashTable()
+{
+	return NULL;
+}
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/HashTable.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/HashTable.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,64 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __HASH_TABLE__
+#define __HASH_TABLE__
+
+typedef struct HashTable
+{
+	long long hv;
+	int *locs; 
+} HashTable;
+
+typedef struct 
+{
+	unsigned int *locs;
+} IHashTable;
+
+int				hashVal(char *seq);
+void			configHashTable();
+char			*getRefGenome();
+char			*getRefGenomeName();
+int				getRefGenomeOffset();
+int				initLoadingHashTable(char *fileName);
+HashTable		*getHashTable();
+
+void 			(*generateHashTable)(char *fileName, char *indexName);
+int				(*loadHashTable)(double *loadTime);
+void			(*finalizeLoadingHashTable)();
+unsigned int	*(*getCandidates)(int hv);
+
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Makefile
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Makefile	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,10 @@
+ALL: mrsfast
+
+LDFLAGS=-static
+LIBS=-lz -lm
+CFLAGS=-O3
+CC=gcc
+
+mrsfast: baseFAST.o MrsFAST.o Common.o CommandLineParser.o RefGenome.o HashTable.o Reads.o Output.o
+	$(CC) $^ -o $@ ${LDFLAGS} ${LIBS}
+	rm *.o
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/MrsFAST.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/MrsFAST.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,1757 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-01-29
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "Common.h"
+#include "Reads.h"
+#include "HashTable.h"
+#include "Output.h"
+#include "MrsFAST.h"
+#include "RefGenome.h"
+
+float calculateScore(int index, char *seq, char *qual, int *err);
+unsigned char		mrFAST = 0;
+char				*versionNumberF="0.2";
+
+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;
+
+ReadIndexTable		*_msf_rIndex = NULL;
+int					_msf_rIndexSize;
+int					_msf_rIndexMax;
+
+SAM					_msf_output;
+
+OPT_FIELDS			*_msf_optionalFields;
+
+char				*_msf_op;
+
+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;
+/**********************************************/
+int compare (const void *a, const void *b)
+{
+	return ((Pair *)a)->hv - ((Pair *)b)->hv;
+}
+/**********************************************/
+void preProcessReads()
+{
+	int i=0;
+	int j=0;
+	int pos = 0;
+
+	_msf_rIndexMax = -1;
+
+	int tmpSize = _msf_seqListSize*_msf_samplingLocsSize*2;
+	Pair *tmp = getMem(sizeof(Pair)*tmpSize);
+
+	for (i=0; i<_msf_seqListSize; i++)
+	{
+		for (j=0; j< _msf_samplingLocsSize; j++)
+		{
+
+			tmp[pos].hv = hashVal(_msf_seqList[i].seq+_msf_samplingLocs[j]);
+			tmp[pos].seqInfo = pos;
+			pos++;
+		}
+		for (j=0; j<_msf_samplingLocsSize; j++)
+		{
+			tmp[pos].hv = hashVal(_msf_seqList[i].rseq+_msf_samplingLocs[j]);
+			tmp[pos].seqInfo = pos;
+			pos++;
+		}
+	}
+
+	qsort(tmp, tmpSize, sizeof(Pair), compare);
+
+
+	int uniq = 0;
+	int prev = -2;
+	int beg = -1;
+	int end = -1;
+
+	for (i=0; i<tmpSize; i++)
+	{
+		if (prev != tmp[i].hv)
+		{
+			uniq ++;
+			prev = tmp[i].hv;
+		}
+	}
+
+	_msf_rIndexSize = uniq;
+	_msf_rIndex = getMem(sizeof(ReadIndexTable)*_msf_rIndexSize);
+	prev = -2;
+
+	j=0;
+	beg =0;
+	while (beg < tmpSize)
+	{
+		end = beg;
+		while (end+1<tmpSize && tmp[end+1].hv==tmp[beg].hv)
+			end++;
+
+		_msf_rIndex[j].hv = tmp[beg].hv;
+		_msf_rIndex[j].seqInfo = getMem(sizeof(int)*(end-beg+2));
+		_msf_rIndex[j].seqInfo[0] = end-beg+1;
+		if ((end-beg+1) > _msf_rIndexMax)
+			_msf_rIndexMax = end-beg+1;
+
+		for (i=1; i<=_msf_rIndex[j].seqInfo[0]; i++)
+		{
+			_msf_rIndex[j].seqInfo[i]=tmp[beg+i-1].seqInfo;
+		}
+		j++;
+		beg = end+1;
+	}
+	freeMem(tmp, sizeof(Pair)*tmpSize);
+}
+/**********************************************/
+void initFAST(Read *seqList, int seqListSize, int *samplingLocs, int samplingLocsSize, char *genFileName)
+{
+	if (_msf_optionalFields == NULL)
+	{
+		_msf_op = getMem(SEQ_LENGTH);
+		if (pairedEndMode)
+		{
+			_msf_optionalFields = getMem(4*sizeof(OPT_FIELDS));
+		}
+		else
+		{
+			_msf_optionalFields = getMem(2*sizeof(OPT_FIELDS));
+		}
+
+		int i;
+		for (i=0; i<200;i++)
+		{
+			sprintf(_msf_numbers[i],"%d%c",i, '\0');
+		}
+		sprintf(_msf_cigar, "%dM", SEQ_LENGTH);
+	}
+
+	if (_msf_samplingLocsEnds == NULL)
+	{
+		int i;
+		_msf_samplingLocs = samplingLocs;
+		_msf_samplingLocsSize = samplingLocsSize;
+
+		_msf_samplingLocsEnds = malloc(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();
+	}
+	if (_msf_refGenName == NULL)
+	{
+		_msf_refGenName = getMem(SEQ_LENGTH);
+	}
+	_msf_refGen =  getRefGenome();
+	_msf_refGenLength = strlen(_msf_refGen);
+	_msf_refGenOffset = getRefGenomeOffset();
+	sprintf(_msf_refGenName,"%s%c", getRefGenomeName(), '\0');
+
+	if (pairedEndMode && _msf_seqHits == NULL)
+	{
+		_msf_mappingInfo  = getMem(seqListSize * sizeof (MappingInfo));
+
+		int i=0;
+		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/2) * sizeof(int));
+
+
+		for (i=0; i<_msf_seqListSize/2; i++)
+		{
+			_msf_seqHits[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/2) * sizeof(int));
+	freeMem(_msf_refGenName, 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);
+}
+
+/**********************************************/
+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;
+
+	int matchCnt = 0;
+	int pp = 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;
+}
+
+/**********************************************/
+int calculateMD(int index, char *seq, int err, char **opSeq)
+{
+	int i;
+	char *ref;
+	char *ver;
+	short matchCnt = 0;
+	char *op = *opSeq;
+	int pp = 0;
+
+	ref = _msf_refGen + index-1;
+	ver = seq;
+
+	if (err>0 || err == -1 )
+	{
+
+		err = 0;
+		for (i=0; i < SEQ_LENGTH; i++)
+		{
+			if (* ref != *ver)
+			{
+				err++;
+				if (matchCnt)
+				{
+					if (matchCnt < 10)
+					{
+						op[pp++]=_msf_numbers[matchCnt][0];
+					}
+					else if (matchCnt < 100)
+					{
+						op[pp++]=_msf_numbers[matchCnt][0];
+						op[pp++]=_msf_numbers[matchCnt][1];
+					}
+					else
+					{
+						op[pp++]=_msf_numbers[matchCnt][0];
+						op[pp++]=_msf_numbers[matchCnt][1];
+						op[pp++]=_msf_numbers[matchCnt][2];
+					}
+
+					matchCnt = 0;
+				}
+				op[pp++]=*ref;
+			}
+			else
+			{
+				matchCnt++;
+			}
+			ref++;
+			ver++;
+		}
+
+	}
+	if (err == 0)
+	{
+		matchCnt = SEQ_LENGTH;
+	}
+
+	if (matchCnt>0)
+	{
+		if (matchCnt < 10)
+		{
+			op[pp++]=_msf_numbers[matchCnt][0];
+		}
+		else if (matchCnt < 100)
+		{
+			op[pp++]=_msf_numbers[matchCnt][0];
+			op[pp++]=_msf_numbers[matchCnt][1];
+		}
+		else
+		{
+			op[pp++]=_msf_numbers[matchCnt][0];
+			op[pp++]=_msf_numbers[matchCnt][1];
+			op[pp++]=_msf_numbers[matchCnt][2];
+		}
+	}
+	op[pp]='\0';
+
+	return err;
+}
+
+/**********************************************/
+void mapSingleEndSeqListBal(unsigned int *l1, int s1, unsigned int *l2, int s2, int dir)
+{
+
+	if (s1 == 0 || s2 == 0)
+	{
+		return;
+	}
+	else if (s1 == s2 && s1 <= 50)
+	{
+
+		int j = 0;
+		int z = 0;
+		int *locs;
+		int *seqInfo;
+		char *_tmpSeq, *_tmpQual;
+		char rqual[QUAL_LENGTH+1];
+		rqual[QUAL_LENGTH]='\0';
+
+		if (dir > 0)
+		{
+			locs		= (int *) l1;
+			seqInfo		= (int *) l2;
+		}
+		else
+		{
+			locs		= (int *) l2;
+			seqInfo		= (int *) l1;
+		}
+
+
+		for (j=0; j<s2; j++)
+		{
+			int re = _msf_samplingLocsSize * 2;
+			int r = seqInfo[j]/re;
+			if (maxHits!=0 && _msf_seqList[r].hits[0] == maxHits)
+			{
+				continue;
+			}
+
+			int x = seqInfo[j] % re;
+			int o = x % _msf_samplingLocsSize;
+			char d = (x/_msf_samplingLocsSize)?1:0;
+
+
+			if (d)
+			{
+				reverse(_msf_seqList[r].qual, rqual, QUAL_LENGTH);
+				_tmpQual = rqual;
+				_tmpSeq = _msf_seqList[r].rseq;
+			}
+			else
+			{
+				_tmpQual = _msf_seqList[r].qual;
+				_tmpSeq = _msf_seqList[r].seq;
+			}
+
+
+			for (z=0; z<s1; z++)
+			{
+				int genLoc = locs[z]-_msf_samplingLocs[o];
+
+
+				if ( genLoc < _msf_refGenBeg || genLoc > _msf_refGenEnd )
+					continue;
+
+				int err = -1;
+
+
+
+				err = verifySingleEnd(genLoc, _tmpSeq, o);
+
+
+
+				if (err != -1)
+				{
+					calculateMD(genLoc, _tmpSeq, err, &_msf_op);
+					mappingCnt++;
+					_msf_seqList[r].hits[0]++;
+
+					_msf_output.QNAME		= _msf_seqList[r].name;
+					_msf_output.FLAG		= 16 * d;
+					_msf_output.RNAME		= _msf_refGenName;
+					_msf_output.POS			= genLoc + _msf_refGenOffset;
+					_msf_output.MAPQ		= 255;
+					_msf_output.CIGAR		= _msf_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 = _msf_op;
+
+					output(_msf_output);
+
+
+					if (_msf_seqList[r].hits[0] == 1)
+					{
+						mappedSeqCnt++;
+					}
+
+					if ( maxHits == 0 )
+					{
+						_msf_seqList[r].hits[0] = 2;
+					}
+
+
+					if ( maxHits!=0 && _msf_seqList[r].hits[0] == maxHits)
+					{
+						completedSeqCnt++;
+						break;
+					}
+				}
+
+			}
+		}
+	}
+	else
+	{
+		int tmp1=s1/2, tmp2= s2/2;
+		if (tmp1 != 0)
+			mapSingleEndSeqListBal(l1, tmp1, l2+tmp2, s2-tmp2, dir);
+		mapSingleEndSeqListBal(l2+tmp2, s2-tmp2, l1+tmp1, s1-tmp1, -dir);
+		if (tmp2 !=0)
+			mapSingleEndSeqListBal(l1+tmp1, s1-tmp1, l2, tmp2, dir);
+		if (tmp1 + tmp2 != 0)
+			mapSingleEndSeqListBal(l2, tmp2, l1, tmp1, -dir);
+	}
+}
+
+
+/**********************************************/
+void mapSingleEndSeqListTOP(unsigned int *l1, int s1, unsigned int *l2, int s2)
+{
+	if (s1 < s2)
+	{
+		mapSingleEndSeqListBal(l1, s1, l2, s1,1);
+		mapSingleEndSeqListTOP(l1, s1, l2+s1, s2-s1);
+	}
+	else if (s1 > s2)
+	{
+		mapSingleEndSeqListBal(l1, s2, l2, s2,1);
+		mapSingleEndSeqListTOP(l1+s2, s1-s2, l2, s2);
+	}
+	else
+	{
+		mapSingleEndSeqListBal(l1, s1, l2, s2,1);
+	}
+}
+
+
+/**********************************************/
+void mapSingleEndSeqList(unsigned int *l1, int s1, unsigned int *l2, int s2)
+{
+	if ( s2/s1 <= 2)
+	{
+		int j = 0;
+		int z = 0;
+		int *locs = (int *) l1;
+		int *seqInfo = (int *) l2;
+		char *_tmpSeq, *_tmpQual;
+		char rqual[QUAL_LENGTH+1];
+		rqual[QUAL_LENGTH]='\0';
+
+		for (j=0; j<s2; j++)
+		{
+			int re = _msf_samplingLocsSize * 2;
+			int r = seqInfo[j]/re;
+			if (maxHits!=0 && _msf_seqList[r].hits[0] == maxHits)
+			{
+				continue;
+			}
+
+			int x = seqInfo[j] % re;
+			int o = x % _msf_samplingLocsSize;
+			char d = (x/_msf_samplingLocsSize)?1:0;
+
+
+			if (d)
+			{
+				reverse(_msf_seqList[r].qual, rqual, QUAL_LENGTH);
+				_tmpQual = rqual;
+				_tmpSeq = _msf_seqList[r].rseq;
+			}
+			else
+			{
+				_tmpQual = _msf_seqList[r].qual;
+				_tmpSeq = _msf_seqList[r].seq;
+			}
+
+
+			for (z=0; z<s1; z++)
+			{
+				int genLoc = locs[z]-_msf_samplingLocs[o];
+
+
+				if ( genLoc < _msf_refGenBeg || genLoc > _msf_refGenEnd )
+					continue;
+
+				int err = -1;
+
+
+
+				err = verifySingleEnd(genLoc, _tmpSeq, o);
+
+
+
+				if (err != -1)
+				{
+					calculateMD(genLoc, _tmpSeq, err, &_msf_op);
+					mappingCnt++;
+					_msf_seqList[r].hits[0]++;
+
+					_msf_output.QNAME		= _msf_seqList[r].name;
+					_msf_output.FLAG		= 16 * d;
+					_msf_output.RNAME		= _msf_refGenName;
+					_msf_output.POS			= genLoc + _msf_refGenOffset;
+					_msf_output.MAPQ		= 255;
+					_msf_output.CIGAR		= _msf_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 = _msf_op;
+
+					output(_msf_output);
+
+
+					if (_msf_seqList[r].hits[0] == 1)
+					{
+						mappedSeqCnt++;
+					}
+
+					if ( maxHits == 0 )
+					{
+						_msf_seqList[r].hits[0] = 2;
+					}
+
+
+					if ( maxHits!=0 && _msf_seqList[r].hits[0] == maxHits)
+					{
+						completedSeqCnt++;
+						break;
+					}
+				}
+
+			}
+		}
+	}
+	else if (s1 == 1)
+	{
+		//	fprintf(stderr, "1");
+		int tmp = s2/2;
+		mapSingleEndSeqList(l1, s1, l2, tmp);
+		mapSingleEndSeqList(l1, s1, l2+tmp, s2-tmp);
+	}
+	else if (s2 == 1)
+	{
+		//	fprintf(stderr, "2");
+		int tmp = s1/2;
+		mapSingleEndSeqList(l1, tmp, l2, s2);
+		mapSingleEndSeqList(l1+tmp, s1-tmp, l2, s2);
+	}
+	else
+	{
+		//	fprintf(stderr, "3");
+		int tmp1=s1/2, tmp2= s2/2;
+		mapSingleEndSeqList(l1, tmp1, l2, tmp2);
+		mapSingleEndSeqList(l1+tmp1, s1-tmp1, l2, tmp2);
+		mapSingleEndSeqList(l1+tmp1, s1-tmp1, l2+tmp2, s2-tmp2);
+		mapSingleEndSeqList(l1, tmp1, l2+tmp2, s2-tmp2);
+	}
+}
+/**********************************************/
+int	 mapSingleEndSeq()
+{
+	int i = 0;
+	unsigned int *locs = NULL;
+	unsigned int *seqInfo = NULL;
+	while ( i < _msf_rIndexSize )
+	{
+
+		locs = getCandidates (_msf_rIndex[i].hv);
+		if ( locs != NULL)
+		{
+			seqInfo  = _msf_rIndex[i].seqInfo;
+			mapSingleEndSeqListTOP (locs+1, locs[0], seqInfo+1, seqInfo[0]);
+		}
+		i++;
+	}
+	return 1;
+}
+
+
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+int compareOut (const void *a, const void *b)
+{
+	FullMappingInfo *aInfo = (FullMappingInfo *)a;
+	FullMappingInfo *bInfo = (FullMappingInfo *)b;
+	return aInfo->loc - bInfo->loc;
+}
+
+/**********************************************/
+void mapPairedEndSeqList(unsigned int *l1, int s1, unsigned int *l2, int s2)
+{
+	if ( s2/s1 <= 2)
+	{
+		int j = 0;
+		int z = 0;
+		int *locs = (int *) l1;
+		int *seqInfo = (int *) l2;
+		char *_tmpSeq, *_tmpQual;
+		char rqual[QUAL_LENGTH+1];
+		rqual[QUAL_LENGTH]='\0';
+
+		for (j=0; j<s2; j++)
+		{
+			int re = _msf_samplingLocsSize * 2;
+			int r = seqInfo[j]/re;
+
+			if (pairedEndDiscordantMode && (_msf_seqList[r].hits[0] == 1 || (_msf_seqHits[r/2] > DISCORDANT_CUT_OFF) ))
+			{
+				continue;
+			}
+
+			int x = seqInfo[j] % re;
+			int o = x % _msf_samplingLocsSize;
+			char d = (x/_msf_samplingLocsSize)?-1:1;
+
+
+			if (d==-1)
+			{
+				_tmpSeq = _msf_seqList[r].rseq;
+			}
+			else
+			{
+				_tmpSeq = _msf_seqList[r].seq;
+			}
+
+
+			for (z=0; z<s1; z++)
+			{
+				int genLoc = locs[z]-_msf_samplingLocs[o];
+
+
+				if ( genLoc < _msf_refGenBeg || genLoc > _msf_refGenEnd )
+					continue;
+
+				int err = -1;
+
+
+
+				err = verifySingleEnd(genLoc, _tmpSeq, o);
+
+
+				if (err != -1)
+				{
+					MappingLocations *parent = NULL;
+					MappingLocations *child = _msf_mappingInfo[r].next;
+
+					genLoc+= _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;
+						if (parent == NULL)
+							_msf_mappingInfo[r].next = tmp;
+						else
+							parent->next = tmp;
+					}
+					else
+					{
+						child->loc[_msf_mappingInfo[r].size % MAP_CHUNKS] = genLoc * d;
+					}
+
+
+
+					_msf_mappingInfo[r].size++;
+
+				}
+
+			}
+		}
+	}
+	else if (s1 == 1)
+	{
+		int tmp = s2/2;
+		mapPairedEndSeqList(l1, s1, l2, tmp);
+		mapPairedEndSeqList(l1, s1, l2+tmp, s2-tmp);
+	}
+	else if (s2 == 1)
+	{
+		int tmp = s1/2;
+		mapPairedEndSeqList(l1, tmp, l2, s2);
+		mapPairedEndSeqList(l1+tmp, s1-tmp, l2, s2);
+	}
+	else
+	{
+		int tmp1=s1/2, tmp2= s2/2;
+		mapPairedEndSeqList(l1, tmp1, l2, tmp2);
+		mapPairedEndSeqList(l1+tmp1, s1-tmp1, l2, tmp2);
+		mapPairedEndSeqList(l1+tmp1, s1-tmp1, l2+tmp2, s2-tmp2);
+		mapPairedEndSeqList(l1, tmp1, l2+tmp2, s2-tmp2);
+	}
+}
+
+/**********************************************/
+int	 mapPairedEndSeq()
+{
+	int i = 0;
+	unsigned int *locs = NULL;
+	unsigned int *seqInfo = NULL;
+	while ( i < _msf_rIndexSize )
+	{
+		locs = getCandidates (_msf_rIndex[i].hv);
+		if ( locs != NULL)
+		{
+			seqInfo  = _msf_rIndex[i].seqInfo;
+			mapPairedEndSeqList(locs+1, locs[0], seqInfo+1, seqInfo[0]);
+
+		}
+		i++;
+	}
+
+
+	char fname1[FILE_NAME_LENGTH];
+	char fname2[FILE_NAME_LENGTH];
+	MappingLocations *cur, *tmp;
+	int tmpOut;
+	int j;
+	int lmax=0, rmax=0;
+
+	sprintf(fname1, "%s__%s__%d__1",mappingOutputPath, mappingOutput, _msf_openFiles);
+	sprintf(fname2, "%s__%s__%d__2",mappingOutputPath, 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;
+				}
+				tmpOut = fwrite(&(cur->loc[j % MAP_CHUNKS]), sizeof(int), 1, out);
+			}
+			_msf_mappingInfo[i].size = 0;
+			//	_msf_mappingInfo[i].next = NULL;
+		}
+	}
+
+	_msf_maxLSize += lmax;
+	_msf_maxRSize += rmax;
+
+	fclose(out1);
+	fclose(out2);
+
+	//fprintf(stdout, "%d %d\n", _msf_maxLSize, _msf_maxRSize);
+
+}
+
+/**********************************************/
+void outputPairedEnd()
+{
+
+	char *curGen;
+	char *curGenName;
+	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, *out1, *out2;
+
+	char fname3[FILE_NAME_LENGTH];
+	char fname4[FILE_NAME_LENGTH];
+	char fname5[FILE_NAME_LENGTH];
+
+	if (pairedEndDiscordantMode)
+	{
+		sprintf(fname3, "%s__%s__disc", mappingOutputPath, mappingOutput);
+		sprintf(fname4, "%s__%s__oea1", mappingOutputPath, mappingOutput);
+		sprintf(fname5, "%s__%s__oea2", mappingOutputPath, mappingOutput);
+		out = fileOpen(fname3, "a");
+		out1 = fileOpen(fname4, "a");
+		out2 = fileOpen(fname5, "a");
+	}
+
+
+
+	int i;
+
+	FullMappingInfo *mi1 = getMem(sizeof(FullMappingInfo) * _msf_maxLSize);
+	FullMappingInfo *mi2 = getMem(sizeof(FullMappingInfo) * _msf_maxRSize);
+
+
+	for (i=0; i<_msf_openFiles; i++)
+	{
+		sprintf(fname1[i], "%s__%s__%d__1", mappingOutputPath, mappingOutput, i);
+		sprintf(fname2[i], "%s__%s__%d__2", mappingOutputPath, mappingOutput, i);
+		in1[i] = fileOpen(fname1[i], "r");
+		in2[i] = fileOpen(fname2[i], "r");
+	}
+
+
+	int size;
+	int j, k;
+	int size1, size2;
+
+	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]);
+					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]);
+
+					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;
+			}
+		}
+
+		//if (i == 6615)
+		//	fprintf(stdout, "%d: %s %d %d ",i, _msf_seqList[i*2].name, size1, size2);	
+
+		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;
+							size1=0;
+							size2=0;
+							break;
+						}
+					}
+					k++;
+				}
+			}
+
+			_msf_seqHits[i]+= size1*size2;
+			if (_msf_seqHits[i]> DISCORDANT_CUT_OFF)
+			{
+				_msf_seqList[i*2].hits[0]=1;
+				_msf_seqList[i*2+1].hits[0]=1;
+				size1=0;
+				size2=0;
+			}
+		}
+		//if (i == 6615)
+		//	fprintf(stdout, "%d %d\n", size1, size2);	
+
+		char *seq1, *seq2, *rseq1, *rseq2, *qual1, *qual2;
+		char rqual1[QUAL_LENGTH+1], rqual2[QUAL_LENGTH+1];
+		rqual1[QUAL_LENGTH] = rqual2[QUAL_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, QUAL_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, QUAL_LENGTH);
+
+
+		if (pairedEndDiscordantMode)
+		{
+			for (k=0; k<size1; k++)
+			{
+				int tm = -1;
+				mi1[k].score = calculateScore(mi1[k].loc, (mi1[k].dir==-1)?rseq1:seq1, (mi1[k].dir==-1)?rqual1:qual1, &tm);
+				mi1[k].err = tm;
+			}
+
+			for (k=0; k<size2; k++)
+			{
+				int tm = -1;
+				mi2[k].score = calculateScore(mi2[k].loc, (mi2[k].dir==-1)?rseq2:seq2, (mi2[k].dir==-1)?rqual2:qual2, &tm);
+				mi2[k].err = tm;
+			}
+
+		}
+		else
+		{
+			for (k=0; k<size1; k++)
+			{
+				mi1[k].err = calculateMD(mi1[k].loc, (mi1[k].dir==-1)?rseq1:seq1, -1, &_msf_op);
+				sprintf(mi1[k].md, "%s", _msf_op);
+			}
+
+			for (k=0; k<size2; k++)
+			{
+				mi2[k].err = calculateMD(mi2[k].loc, (mi2[k].dir==-1)?rseq2:seq2, -1, &_msf_op);
+				sprintf(mi2[k].md, "%s", _msf_op);
+			}
+		}
+		pos = 0;
+
+		for (j=0; j<size1; j++)
+		{
+			lm = mi1[j].loc - maxPairEndedDistance + 1;
+			ll = mi1[j].loc - minPairEndedDistance + 1;
+			rl = mi1[j].loc + minPairEndedDistance - 1;
+			rm = mi1[j].loc + maxPairEndedDistance - 1;
+
+			//fprintf(stdout, "%d %d %d %d %d\n",lm, ll,mi1[j].loc ,rl, rm); 
+
+			while (pos<size2 && mi2[pos].loc < lm)
+			{
+				pos++;
+			}
+
+			//fprintf(stdout, "POS: %d %d \n", pos, mi2[pos].loc);
+
+			k = pos;
+			while (k<size2 && mi2[k].loc <= rm)
+			{
+				if (mi2[k].loc <= ll || mi2[k].loc >= rl)
+				{
+					if (pairedEndDiscordantMode)
+					{
+						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);
+
+						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(char), 1, out);
+						tmp = fwrite(&sc, sizeof(float), 1, 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(char), 1, out);
+						tmp = fwrite(&sc, sizeof(float), 1, out);
+					} // end discordant
+					else
+					{ //start sampe
+						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-1;												
+						if (mi1[j].loc - mi2[k].loc > 0)
+						{
+							isize *= -1;
+						}
+
+						d1 = (mi1[j].dir == -1)?1:0;
+						d2 = (mi2[k].dir == -1)?1:0;
+
+						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		= _msf_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;
+
+
+						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		= _msf_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;
+
+						output(_msf_output);
+					} //end sampe
+				}
+				k++;
+			}
+		}
+	}
+
+	if (pairedEndDiscordantMode)
+	{
+		fclose(out);
+	}
+
+
+	freeMem(mi1, sizeof(FullMappingInfo)*_msf_maxLSize);
+	freeMem(mi2, sizeof(FullMappingInfo)*_msf_maxRSize);
+
+	for (i=0; i<_msf_openFiles; i++)
+	{
+		fclose(in1[i]);
+		fclose(in2[i]);
+		//fprintf(stdout, "%s %s \n", fname1[i], fname2[i]);
+		unlink(fname1[i]);
+		unlink(fname2[i]);
+	}
+	_msf_openFiles = 0;
+
+}
+
+/**********************************************/
+/**********************************************/
+/**********************************************/
+/**********************************************/
+float calculateScore(int index, char *seq, char *qual, int *err)
+{
+	int i;
+	char *ref;
+	char *ver;
+
+	ref = _msf_refGen + index-1;
+	ver = seq;
+	float score = 1;
+
+	if (*err > 0 || *err == -1)
+	{
+		*err = 0;
+
+		for (i=0; i < SEQ_LENGTH; i++)
+		{
+			if (*ref != *ver)
+			{
+				//fprintf(stdout, "%c %c %d", *ref, *ver, *err);
+				(*err)++;
+				score *= 0.001 + 1/pow( 10, ((qual[i]-33)/10.0) );
+			}
+			ref++;
+			ver++;
+		}
+
+	}
+	return score;
+}
+
+/**********************************************/
+void outputPairedEndDiscPP()
+{
+	char genName[SEQ_LENGTH];
+	char fname1[FILE_NAME_LENGTH];
+	char fname2[FILE_NAME_LENGTH];
+	char fname3[FILE_NAME_LENGTH];
+	char fname4[FILE_NAME_LENGTH];
+	char fname5[FILE_NAME_LENGTH];
+	char fname6[FILE_NAME_LENGTH];
+	char l;
+	int loc1, loc2;
+	char err1, err2;
+	char dir1, dir2;
+	float sc1, sc2, lsc=0;
+	int flag = 0;
+	int rNo,lrNo = -1;
+	int tmp;
+	FILE *in, *in1, *in2, *out, *out1, *out2;
+
+	sprintf(fname1, "%s__%s__disc", mappingOutputPath, mappingOutput);
+	sprintf(fname2, "%s__%s__oea1", mappingOutputPath, mappingOutput);
+	sprintf(fname3, "%s__%s__oea2", mappingOutputPath, mappingOutput);
+	sprintf(fname4, "%s%s_DIVET.vh", mappingOutputPath, mappingOutput);
+	sprintf(fname5, "%s%s_OEA1.vh", mappingOutputPath, mappingOutput);
+	sprintf(fname6, "%s%s_OEA2.vh", mappingOutputPath, mappingOutput);
+
+	in   = fileOpen(fname1, "r");
+	in1  = fileOpen(fname2, "r");
+	in2  = fileOpen(fname3, "r");
+	out  = fileOpen(fname4, "w");
+	out1 = fileOpen(fname5, "w");
+	out2 = fileOpen(fname6, "w");
+	if (in != NULL)
+	{
+		flag = fread(&rNo, sizeof(int), 1, in);
+	}
+	else
+	{
+		flag  = 0;
+	}
+
+
+	while (flag)
+	{
+
+		tmp = fread(&l, sizeof(char), 1, in);
+		tmp = fread(genName, sizeof(char), l, in);
+		genName[l]='\0';
+		tmp = fread(&loc1, sizeof(int), 1, in);
+		tmp = fread(&err1, sizeof(char), 1, in);
+		tmp = fread(&sc1, sizeof(float), 1, in);
+		tmp = fread(&loc2, sizeof(int), 1, in);
+		tmp = fread(&err2, sizeof(char), 1, in);
+		tmp = fread(&sc2, sizeof(float), 1, in);
+
+		//if (rNo ==6615)
+		//	fprintf(stdout, "%s %d: %d %0.20f %d %d %0.20f\n", genName, loc1, err1, sc1, loc2, err2, sc2);
+
+		if (_msf_seqList[rNo*2].hits[0] % 2 == 0 && _msf_seqHits[rNo] < DISCORDANT_CUT_OFF)
+		{
+			dir1 = dir2 = 'F';
+
+			if (loc1 < 0)
+			{
+				dir1 = 'R';
+				loc1 = -loc1;
+			}
+
+			if (loc2 < 0)
+			{
+				dir2 = 'R';
+				loc2 = -loc2;
+			}
+
+			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;
+			}
+
+			int inv = 0;
+			int eve = 0;
+			int dist = 0;
+			char event;
+
+			//fprintf(stdout, "%c %c ", dir1, dir2);
+
+			if ( dir1 == dir2 )
+			{
+				event = 'V';
+				//fprintf(stdout, "Inverstion \n");
+			}
+			else
+			{
+				if (loc1 < loc2)
+				{
+
+					//fprintf(stdout, "< %d ", loc2-loc1-SEQ_LENGTH);
+
+					if (dir1 == 'R' && dir2 == 'F')
+					{
+						event = 'E';
+
+						//fprintf(stdout, "Everted \n");
+					}
+					else if ( loc2 - loc1 >= maxPairEndedDiscordantDistance )
+					{
+						event = 'D';
+						//fprintf(stdout, "Deletion \n");
+					}
+					else
+					{
+						event = 'I';
+						//fprintf(stdout, "Insertion \n");
+					}
+				}
+				else if (loc2 < loc1)
+				{
+					//fprintf(stdout, "> %d ", loc1-loc2-SEQ_LENGTH);
+					if (dir2 == 'R' && dir1 == 'F')
+					{
+						event = 'E';
+						//fprintf(stdout, "Everted \n");
+					}
+					else if ( loc1 - loc2 >= maxPairEndedDiscordantDistance )
+					{
+						event = 'D';
+						//fprintf(stdout, "Deletion \n");
+					}
+					else
+					{
+						event = 'I';
+						//fprintf(stdout, "Insertion \n");
+					}
+				}
+			}
+			_msf_seqList[rNo*2].hits[0] = 2;
+			fprintf(out, "%s\t%s\t%d\t%d\t%c\t%s\t%d\t%d\t%c\t%c\t%d\t%0.0f\t%0.20f\n",
+					_msf_seqList[rNo*2].name, genName, loc1, (loc1+SEQ_LENGTH-1), dir1, genName, loc2, (loc2+SEQ_LENGTH-1), dir2, event, (err1+err2), lsc, sc1*sc2);
+		}
+		flag = fread(&rNo, sizeof(int), 1, in);
+
+	}
+
+	/*
+	   MappingInfoNode *lr[_msf_seqListSize/2];
+	   MappingInfoNode *rr[_msf_seqListSize/2];
+	   MappingInfoNode *cur, *tmpDel, *cur2;
+
+
+	   int ls[_msf_seqListSize/2];
+	   int rs[_msf_seqListSize/2];
+
+
+	   int i=0;
+
+	   for (i = 0; i<_msf_seqListSize/2; i++)
+	   {
+	   lr[i] = rr[i] = NULL;
+	   ls[i] = rs[i] = 0;
+	   }
+
+
+
+	   if (in1 != NULL)
+	   {
+	   flag = fread(&rNo, sizeof(int), 1, in1);
+	   }
+	   else
+	   {
+	   flag  = 0;
+	   }
+
+
+	   while (flag)
+	   {
+	   tmp = fread(&loc1, sizeof(int), 1, in1);
+	   tmp = fread(&err1, sizeof(char), 1, in1);
+	   tmp = fread(&sc1, sizeof(float), 1, in1);
+	   tmp = fread(&l, sizeof(char), 1, in1);
+	   tmp = fread(genName, sizeof(char), l, in1);
+	   genName[l]='\0';
+
+	   if (_msf_seqList[rNo*2].hits[0] == 0)
+	   {
+
+	   if ( ls[rNo] < DISCORDANT_CUT_OFF )
+	   {
+	   ls[rNo]++;
+
+	   cur = lr[rNo];
+
+	   if (cur !=NULL)
+	   {
+	   if (err1 == cur->err)
+	   {
+	   MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+	   nr->loc = loc1;
+	   nr->err = err1;
+	   nr->score = sc1;
+	   nr->next = lr[rNo];
+	   sprintf(nr->chr,"%s", genName);
+	   lr[rNo] = nr;
+	   }
+	   else if (err1 < cur->err)
+	   {
+	   MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+	   nr->loc = loc1;
+	   nr->err = err1;
+	   nr->score = sc1;
+	   sprintf(nr->chr,"%s", genName);
+	   nr->next = NULL;
+	   lr[rNo] = nr;
+	while (cur!=NULL)
+	{
+		tmpDel = cur;
+		cur = cur->next;
+		freeMem(tmpDel, sizeof(MappingInfoNode));
+	}
+}
+}
+else
+{
+
+	MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+	nr->loc = loc1;
+	nr->err = err1;
+	nr->score = sc1;
+	sprintf(nr->chr,"%s", genName);
+	nr->next = NULL;
+	lr[rNo] = nr;
+}
+
+if (ls[rNo] > DISCORDANT_CUT_OFF)
+{
+	cur = lr[rNo];
+	while (cur!=NULL)
+	{
+		tmpDel = cur;
+		cur = cur->next;
+		freeMem(tmpDel, sizeof(MappingInfoNode));
+	}
+}
+}
+
+}
+flag = fread(&rNo, sizeof(int), 1, in1);
+
+}
+
+
+if (in2 != NULL)
+{
+	flag = fread(&rNo, sizeof(int), 1, in2);
+}
+else
+{
+	flag  = 0;
+}
+
+
+while (flag)
+{
+	tmp = fread(&loc1, sizeof(int), 1, in2);
+	tmp = fread(&err1, sizeof(char), 1, in2);
+	tmp = fread(&sc1, sizeof(float), 1, in2);
+	tmp = fread(&l, sizeof(char), 1, in2);
+	tmp = fread(genName, sizeof(char), l, in2);
+	genName[l]='\0';
+
+	if (_msf_seqList[rNo*2].hits[0] == 0)
+	{
+
+		if ( rs[rNo] < DISCORDANT_CUT_OFF )
+		{
+			rs[rNo]++;
+
+			cur = rr[rNo];
+
+			if (cur !=NULL)
+			{
+				if (err1 == cur->err)
+				{
+					MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+					nr->loc = loc1;
+					nr->err = err1;
+					nr->score = sc1;
+					nr->next = rr[rNo];
+					sprintf(nr->chr,"%s", genName);
+					rr[rNo] = nr;
+				}
+				else if (err1 < cur->err)
+				{
+					MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+					nr->loc = loc1;
+					nr->err = err1;
+					nr->score = sc1;
+					sprintf(nr->chr,"%s", genName);
+					nr->next = NULL;
+					rr[rNo] = nr;
+					while (cur!=NULL)
+					{
+						tmpDel = cur;
+						cur = cur->next;
+						freeMem(tmpDel, sizeof(MappingInfoNode));
+					}
+				}
+			}
+			else
+			{
+
+				MappingInfoNode *nr = getMem(sizeof(MappingInfoNode));
+
+				nr->loc = loc1;
+				nr->err = err1;
+				nr->score = sc1;
+				sprintf(nr->chr,"%s", genName);
+				nr->next = NULL;
+				rr[rNo] = nr;
+			}
+
+			if (rs[rNo] > DISCORDANT_CUT_OFF)
+			{
+				cur = rr[rNo];
+				while (cur!=NULL)
+				{
+					tmpDel = cur;
+					cur = cur->next;
+					freeMem(tmpDel, sizeof(MappingInfoNode));
+				}
+			}
+		}
+	}
+	flag = fread(&rNo, sizeof(int), 1, in2);
+
+}
+
+
+for (i=0; i<_msf_seqListSize/2; i++)
+{
+	int j;
+	for (j=0; j<SEQ_LENGTH; j++)
+	{
+		lsc += _msf_seqList[i*2].qual[j]+_msf_seqList[i*2+1].qual[j];
+	}
+	lsc /= 2*SEQ_LENGTH;
+	lsc -= 33;
+	if (ls[i] * rs[i] < DISCORDANT_CUT_OFF && ls[i] & rs[i] > 0)
+	{
+		cur = lr[i];
+		while (cur != NULL)
+		{
+			cur2 = rr[i];
+			if (cur->loc < 0)
+			{
+				dir1 = 'R';
+				loc1 = -cur->loc;
+			}
+			else
+			{
+				dir1 = 'F';
+				loc1 = cur->loc;
+			}
+			while (cur2 != NULL)
+			{
+
+				if (cur2->loc < 0)
+				{
+					dir2 = 'R';
+					loc2 = -cur2->loc;
+				}
+				else
+				{
+					dir2 = 'F';
+					loc2 = cur2->loc;
+				}
+
+				fprintf(out, "%s\t%s\t%d\t%d\t%c\t%s\t%d\t%d\t%c\t%c\t%d\t%0.0f\t%0.20f\n",
+						_msf_seqList[i*2].name, cur->chr, loc1, (loc1+SEQ_LENGTH-1), dir1, cur2->chr, loc2, (loc2+SEQ_LENGTH-1), dir2, 'T', (cur->err+cur2->err), lsc, cur->score*cur2->score);
+				cur2 = cur2->next;
+			}
+			cur = cur->next;
+		}
+	}
+
+}*/
+
+
+fclose(in);
+fclose(in1);
+fclose(in2);
+fclose(out);
+fclose(out1);
+fclose(out2);
+
+unlink(fname1);
+unlink(fname2);
+unlink(fname3);
+unlink(fname5);
+unlink(fname6);
+}
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/MrsFAST.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/MrsFAST.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,103 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __MRS_FAST__
+#define __MRS_FAST__
+
+#include "Reads.h"
+
+#define MAP_CHUNKS 15
+
+// Pair is used to pre-processing and making the read index table
+typedef struct
+{
+	int hv;
+	int seqInfo;
+} Pair;
+
+typedef struct
+{
+	int hv;
+	unsigned int *seqInfo;
+} ReadIndexTable;
+
+
+typedef struct mn
+{
+	int loc;
+	char dir;
+	char err;
+	float score;
+	char md[40];
+	char chr[10];
+} FullMappingInfo;
+
+typedef struct lc
+{
+	int loc[MAP_CHUNKS];
+	struct lc *next;
+} MappingLocations;
+
+typedef struct inf
+{
+	int size;
+	MappingLocations *next;
+} MappingInfo;
+
+typedef struct 
+{
+	FILE * fp;
+	char name[400];
+} FILE_STRUCT;
+
+extern long long			verificationCnt;
+extern long long			mappingCnt;
+extern long long			mappedSeqCnt;
+extern long long			completedSeqCnt;
+
+void initFAST(	Read *seqList,
+				int seqListSize,
+				int *samplingLocs,
+				int samplingLocsSize, 
+				char *fileName);
+
+void finalizeFAST();
+
+int mapSingleEndSeq();
+int mapPaiedEndSeq();
+void outputPairedEnd();
+void outputPairedEndDiscPP();
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Output.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Output.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,175 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <zlib.h>
+#include <string.h>
+#include "Common.h"
+#include "Output.h"
+
+FILE			*_out_fp;
+gzFile			_out_gzfp;
+
+
+
+void finalizeGZOutput()
+{
+	gzclose(_out_gzfp);
+}
+
+void finalizeTXOutput()
+{
+	fclose(_out_fp);
+}
+
+
+void gzOutputQ(SAM map)
+{
+	gzprintf(_out_gzfp, "%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s", 
+			map.QNAME, 
+			map.FLAG,
+			map.RNAME, 
+			map.POS,
+			map.MAPQ,
+			map.CIGAR,
+			map.MRNAME,
+			map.MPOS,
+			map.ISIZE,
+			map.SEQ,
+			map.QUAL);
+	
+	int i;
+
+	for ( i = 0; i < map.optSize; i++)
+	{
+		switch (map.optFields[i].type)
+		{
+			case 'A':
+						gzprintf(_out_gzfp, "\t%s:%c:%c", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].cVal);
+						break;
+			case 'i':
+						gzprintf(_out_gzfp, "\t%s:%c:%d", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].iVal);
+						break;
+			case 'f':
+						gzprintf(_out_gzfp, "\t%s:%c:%f", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].fVal);
+						break;
+			case 'Z':
+			case 'H':
+						gzprintf(_out_gzfp, "\t%s:%c:%s", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].sVal);
+						break;
+		}
+	}
+	gzprintf(_out_gzfp, "\n");
+}
+
+void outputQ(SAM map)
+{
+
+	fprintf(_out_fp, "%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s", 
+			map.QNAME, 
+			map.FLAG,
+			map.RNAME, 
+			map.POS,
+			map.MAPQ,
+			map.CIGAR,
+			map.MRNAME,
+			map.MPOS,
+			map.ISIZE,
+			map.SEQ,
+			map.QUAL);
+
+	
+	int i;
+
+	for ( i = 0; i < map.optSize; i++)
+	{
+		switch (map.optFields[i].type)
+		{
+			case 'A':
+						fprintf(_out_fp, "\t%s:%c:%c", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].cVal);
+						break;
+			case 'i':
+						fprintf(_out_fp, "\t%s:%c:%d", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].iVal);
+						break;
+			case 'f':
+						fprintf(_out_fp, "\t%s:%c:%f", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].fVal);
+						break;
+			case 'Z':
+			case 'H':
+						fprintf(_out_fp, "\t%s:%c:%s", map.optFields[i].tag, map.optFields[i].type, map.optFields[i].sVal);
+						break;
+		}
+	}
+	
+	fprintf(_out_fp, "\n");
+}
+
+int initOutput ( char *fileName, int compressed)
+{
+	if (compressed)
+	{
+		char newFileName[strlen(mappingOutputPath)+strlen(fileName)+4];
+		sprintf(newFileName, "%s%s.gz", mappingOutputPath, fileName);
+		_out_gzfp = fileOpenGZ(newFileName, "w1f");
+		if (_out_gzfp == Z_NULL)
+		{
+			return 0;
+		}
+	
+		finalizeOutput = &finalizeGZOutput;
+
+		output = &gzOutputQ;
+	}
+	else
+	{
+	
+		char newFileName[strlen(mappingOutputPath)+strlen(fileName)];
+		sprintf(newFileName, "%s%s", mappingOutputPath, fileName);
+
+		_out_fp = fileOpen(newFileName, "w");
+		if (_out_fp == NULL)
+		{
+			return 0;
+		}
+
+		finalizeOutput = &finalizeTXOutput;
+		output = &outputQ;
+	}
+	return 1;
+}
+
+
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Output.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Output.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,76 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __OUTPUT__
+#define __OUTPUT__
+
+#define	FORWARD 0
+#define REVERSE 1
+
+typedef struct
+{
+	char		*tag;
+	char		type;
+	char		cVal;
+	int			iVal;
+	float		fVal;
+	char		*sVal;
+} OPT_FIELDS;
+
+typedef struct  
+{
+	char			*QNAME;
+	short			FLAG;
+	char			*RNAME;
+	int				POS;
+	unsigned char	MAPQ;
+	char			*CIGAR;
+	char			*MRNAME;
+	int				MPOS;
+	int				ISIZE;
+	char			*SEQ;
+	char			*QUAL;
+	
+	int			optSize;
+	OPT_FIELDS	*optFields;
+} SAM;
+
+int initOutput(char *fileName, int compressed);
+void (*finalizeOutput)();
+void (*output)(SAM map);
+
+
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Reads.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Reads.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,559 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <zlib.h>
+#include "Common.h"
+#include "Reads.h"
+
+
+
+FILE *_r_fp1;
+FILE *_r_fp2;
+gzFile _r_gzfp1;
+gzFile _r_gzfp2;
+Read *_r_seq;
+int _r_seqCnt;
+int *_r_samplingLocs;
+
+/**********************************************/
+char *(*readFirstSeq)(char *);
+char *(*readSecondSeq)(char *);
+/**********************************************/
+char *readFirstSeqTXT( char *seq )
+{
+	return fgets(seq, SEQ_MAX_LENGTH, _r_fp1);
+}
+
+/**********************************************/
+char *readSecondSeqTXT( char *seq )
+{
+	return fgets(seq, SEQ_MAX_LENGTH, _r_fp2);
+}
+/**********************************************/
+char *readFirstSeqGZ( char *seq )
+{
+	return gzgets(_r_gzfp1, seq, SEQ_MAX_LENGTH);
+}
+
+/**********************************************/
+char *readSecondSeqGZ( char *seq )
+{
+	return gzgets(_r_gzfp2, seq, SEQ_MAX_LENGTH);
+}
+/**********************************************/
+int readAllReads(char *fileName1,
+						char *fileName2,
+						int compressed,
+						unsigned char *fastq,
+						unsigned char pairedEnd,
+						Read **seqList,
+						unsigned int *seqListSize)
+{
+	double startTime=getTime();
+
+	char seq1[SEQ_MAX_LENGTH];
+	char rseq1[SEQ_MAX_LENGTH];
+	char name1[SEQ_MAX_LENGTH];
+	char qual1[SEQ_MAX_LENGTH];
+	char seq2[SEQ_MAX_LENGTH];
+	char rseq2[SEQ_MAX_LENGTH];
+	char name2[SEQ_MAX_LENGTH];
+	char qual2[SEQ_MAX_LENGTH];
+
+	char dummy[SEQ_MAX_LENGTH];
+	char ch;
+	int err1, err2;
+	int nCnt;
+	int discarded = 0;
+	int seqCnt = 0;
+	int maxCnt = 0;
+	int i;
+	Read *list = NULL;
+
+
+	if (!compressed)
+	{
+		_r_fp1 = fileOpen( fileName1, "r");
+
+		if (_r_fp1 == NULL)
+		{
+			return 0;
+		}
+
+		ch = fgetc(_r_fp1);
+
+		if ( pairedEnd && fileName2 != NULL )
+		{
+			_r_fp2 = fileOpen ( fileName2, "r" );
+			if (_r_fp2 == NULL)
+			{
+				return 0;
+			}
+		}
+		else
+		{
+			_r_fp2 = _r_fp1;
+		}
+
+		readFirstSeq = &readFirstSeqTXT;
+		readSecondSeq = &readSecondSeqTXT;
+	}
+	else
+	{
+
+		_r_gzfp1 = fileOpenGZ (fileName1, "r");
+
+		if (_r_gzfp1 == NULL)
+		{
+			return 0;
+		}
+
+		ch = gzgetc(_r_gzfp1);
+
+		if ( pairedEnd && fileName2 != NULL )
+		{
+			_r_fp2 = fileOpenGZ ( fileName2, "r" );
+			if (_r_fp2 == NULL)
+			{
+				return 0;
+			}
+		}
+		else
+		{
+			_r_fp2 = _r_fp1;
+		}
+
+		readFirstSeq = &readFirstSeqGZ;
+		readSecondSeq = &readSecondSeqGZ;
+	}
+
+	if (ch == '>')
+		*fastq = 0;
+	else
+		*fastq = 1;
+
+	// Counting the number of lines in the file
+	while (readFirstSeq(dummy)) maxCnt++;
+
+	if (!compressed)
+	{
+		rewind(_r_fp1);
+	}
+	else
+	{
+		gzrewind(_r_gzfp1);
+	}
+
+	// Calculating the Maximum # of sequences
+	if (*fastq)
+	{
+		maxCnt /= 4;
+	}
+	else
+	{
+		maxCnt /= 2;
+	}
+
+
+
+	if (pairedEnd && fileName2 != NULL )
+		maxCnt *= 2;
+
+	list = getMem(sizeof(Read)*maxCnt);
+
+	while( readFirstSeq(name1) )
+	{
+		err1 = 0;
+		err2 = 0;
+		readFirstSeq(seq1);
+		name1[strlen(name1)-1] = '\0';
+		for (i=0; i<strlen(name1);i++)
+		{
+			if (name1[i] == ' ')
+			{
+				name1[i] = '\0';
+				break;
+			}
+
+		}
+
+		if ( *fastq )
+		{
+			readFirstSeq(dummy);
+			readFirstSeq(qual1);
+			qual1[strlen(qual1)-1] = '\0';
+		}
+		else
+		{
+			sprintf(qual1, "*");
+		}
+
+
+		// Cropping
+		if (cropSize > 0)
+		{
+			seq1[cropSize] = '\0';
+			if ( *fastq )
+				qual1[cropSize] = '\0';
+		}
+
+
+		nCnt = 0;
+		for (i=0; i<strlen(seq1); i++)
+		{
+			seq1[i] = toupper (seq1[i]);
+			if (seq1[i] == 'N')
+			{
+				nCnt++;
+			}
+			else if (isspace(seq1[i]))
+			{
+
+				seq1[i] = '\0';
+				break;
+			}
+		}
+
+		if (nCnt > errThreshold)
+		{
+			err1 = 1;
+		}
+
+		// Reading the second seq of pair-ends
+		if (pairedEnd)
+		{
+			readSecondSeq(name2);
+			readSecondSeq(seq2);
+			name2[strlen(name2)-1] = '\0';
+			for (i=0; i<strlen(name2);i++)
+			{
+				if (name2[i] == ' ')
+				{
+					name2[i] = '\0';
+					break;
+				}
+
+			}
+
+			if ( *fastq )
+			{
+				readSecondSeq(dummy);
+				readSecondSeq(qual2);
+
+				qual2[strlen(qual2)-1] = '\0';
+			}
+			else
+			{
+				sprintf(qual2, "*");
+			}
+
+
+			// Cropping
+			if (cropSize > 0)
+			{
+				seq2[cropSize] = '\0';
+				if ( *fastq )
+					qual2[cropSize] = '\0';
+			}
+
+
+			nCnt = 0;
+			for (i=0; i<strlen(seq2); i++)
+			{
+				seq2[i] = toupper (seq2[i]);
+				if (seq2[i] == 'N')
+				{
+					nCnt++;
+
+				}
+				else if (isspace(seq2[i]))
+				{
+					seq2[i] = '\0';
+				}
+			}
+			if (nCnt > errThreshold)
+			{
+				err2 = 1;
+			}
+		}
+
+		if (!pairedEnd && !err1)
+		{
+
+			int _mtmp = strlen(seq1);
+			list[seqCnt].hits = getMem (1+3*_mtmp+3+strlen(name1)+1);
+			list[seqCnt].seq = list[seqCnt].hits + 1;
+			list[seqCnt].rseq = list[seqCnt].seq + _mtmp+1;
+			list[seqCnt].qual = list[seqCnt].rseq + _mtmp+1;
+			list[seqCnt].name = list[seqCnt].qual + _mtmp+1;
+
+
+			reverseComplete(seq1, rseq1, _mtmp);
+			rseq1[_mtmp] =  '\0';
+			int i;
+
+			list[seqCnt].hits[0] = 0;
+
+			for (i=0; i<=_mtmp; i++)
+			{
+				list[seqCnt].seq[i] = seq1[i];
+				list[seqCnt].rseq[i] = rseq1[i] ;
+				list[seqCnt].qual[i] = qual1[i];
+			}
+			sprintf(list[seqCnt].name,"%s%c", ((char*)name1)+1,'\0');
+
+			seqCnt++;
+
+		}
+		else if (pairedEnd && !err1 && !err2)
+		{
+			// Naming Conventions X/1, X/2 OR X
+			int tmplen = strlen(name1);
+			if (strcmp(name1, name2) != 0)
+			{
+				tmplen = strlen(name1)-2;
+			}
+		
+			if (strcmp(name1, "@IL11_266:2:1:922:509/1") == 0)
+			{
+				fprintf(stdout, "%d\n", seqCnt);
+			}
+			//first seq
+			int _mtmp = strlen(seq1);
+			list[seqCnt].hits = getMem (1+3*_mtmp+3+tmplen+1);
+			list[seqCnt].seq = list[seqCnt].hits + 1;
+			list[seqCnt].rseq = list[seqCnt].seq + _mtmp+1;
+			list[seqCnt].qual = list[seqCnt].rseq + _mtmp+1;
+			list[seqCnt].name = list[seqCnt].qual + _mtmp+1;
+
+			reverseComplete(seq1, rseq1, _mtmp);
+			rseq1[_mtmp] =  '\0';
+			int i;
+
+			list[seqCnt].hits[0] = 0;
+
+			for (i=0; i<=_mtmp; i++)
+			{
+				list[seqCnt].seq[i] = seq1[i];
+				list[seqCnt].rseq[i] = rseq1[i] ;
+				list[seqCnt].qual[i] = qual1[i];
+			}
+
+
+			name1[tmplen]='\0';
+			sprintf(list[seqCnt].name,"%s%c", ((char*)name1)+1,'\0');
+
+
+			seqCnt++;
+
+			//second seq
+			list[seqCnt].hits = getMem (1+3*_mtmp+3+tmplen+1);
+			list[seqCnt].seq = list[seqCnt].hits + 1;
+			list[seqCnt].rseq = list[seqCnt].seq + _mtmp+1;
+			list[seqCnt].qual = list[seqCnt].rseq + _mtmp+1;
+			list[seqCnt].name = list[seqCnt].qual + _mtmp+1;
+
+			reverseComplete(seq2, rseq2, _mtmp);
+			rseq2[_mtmp] =  '\0';
+
+			list[seqCnt].hits[0] = 0;
+
+			for (i=0; i<=_mtmp; i++)
+			{
+				list[seqCnt].seq[i] = seq2[i];
+				list[seqCnt].rseq[i] = rseq2[i] ;
+				list[seqCnt].qual[i] = qual2[i];
+			}
+
+
+			name2[tmplen]='\0';
+			sprintf(list[seqCnt].name,"%s%c", ((char*)name2)+1,'\0');
+
+
+			seqCnt++;
+
+		}
+		else
+		{
+			discarded++;
+		}
+	}
+
+	if (seqCnt > 0)
+	{
+		QUAL_LENGTH = SEQ_LENGTH = strlen(list[0].seq);
+		if (! *fastq)
+		{
+			QUAL_LENGTH = 1;
+		}
+		//fprintf(stderr, "%d %d\n", SEQ_LENGTH, QUAL_LENGTH);
+	}
+	else
+	{
+		fprintf(stdout, "ERR: No reads can be found for mapping\n");
+		return 0;
+	}
+
+
+	if (pairedEnd)
+	{
+//		seqCnt /= 2;
+	}
+
+
+	// Closing Files
+	if (!compressed)
+	{
+		fclose(_r_fp1);
+		if ( pairedEnd && fileName2 != NULL )
+		{
+			fclose(_r_fp2);
+		}
+	}
+	else
+	{
+		gzclose(_r_gzfp1);
+		if ( pairedEnd && fileName2 != NULL)
+		{
+			gzclose(_r_fp2);
+		}
+	}
+
+	*seqList = list;
+	*seqListSize = seqCnt;
+
+	_r_seq = list;
+	_r_seqCnt = seqCnt;
+
+	fprintf(stdout, "%d sequences are read in %0.2f. (%d discarded) [Mem:%0.2f M]\n", seqCnt, (getTime()-startTime), discarded, getMemUsage());
+	//totalLoadingTime+=getTime()-startTime;
+
+	return 1;
+}
+/**********************************************/
+void loadSamplingLocations(int **samplingLocs, int * samplingLocsSize)
+{
+	int i;
+	int samLocsSize = errThreshold + 1;
+	int *samLocs = getMem(sizeof(int)*samLocsSize);
+
+	for (i=0; i<samLocsSize; i++)
+	{
+		samLocs[i] = (SEQ_LENGTH / samLocsSize) *i;
+		if ( samLocs[i] + WINDOW_SIZE > SEQ_LENGTH)
+			samLocs[i] = SEQ_LENGTH - WINDOW_SIZE;
+	}
+
+	// Outputing the sampling locations
+
+/*	int j;
+ 	for (i=0; i<SEQ_LENGTH; i++)
+	{
+		fprintf(stdout, "-");
+	}
+	fprintf(stdout, "\n");
+
+	for ( i=0; i<samLocsSize; i++ )
+	{
+		for ( j=0; j<samLocs[i]; j++ )
+		{
+			fprintf(stdout," ");
+		}
+		for (j=0; j<WINDOW_SIZE; j++)
+		{
+			fprintf(stdout,"+");
+		}
+		fprintf(stdout, "\n");
+		fflush(stdout);
+	}
+	for ( i=0; i<SEQ_LENGTH; i++ )
+	{
+		fprintf(stdout, "-");
+	}
+	fprintf(stdout, "\n");*/
+	*samplingLocs = samLocs;
+	*samplingLocsSize = samLocsSize;
+	_r_samplingLocs = samLocs;
+}
+
+void finalizeReads(char *fileName)
+{
+	FILE *fp1=NULL;
+
+	if (fileName != NULL)
+	{
+		fp1 = fileOpen(fileName, "w");
+	}
+	if (pairedEndMode)
+		_r_seqCnt /=2;
+
+	int i=0;
+	for (i = 0; i < _r_seqCnt; i++)
+	{
+		if (pairedEndMode && _r_seq[2*i].hits[0] == 0 &&  strcmp(_r_seq[2*i].qual,"*")!=0)
+		{
+			fprintf(fp1,"@%s/1\n%s\n+\n%s\n@%s/2\n%s\n+\n%s\n", _r_seq[i*2].name, _r_seq[i*2].seq, _r_seq[i*2].qual, _r_seq[i*2].name, _r_seq[i*2+1].seq, _r_seq[i*2+1].qual);
+		}
+		else if (pairedEndMode && _r_seq[2*i].hits[0] == 0)
+		{
+			fprintf(fp1, ">%s/1\n%s\n>%s/2\n%s\n", _r_seq[i*2].name, _r_seq[i*2].seq, _r_seq[i*2].name, _r_seq[i*2+1].seq);
+		}
+		else if (_r_seq[i].hits[0] == 0 && strcmp(_r_seq[i].qual, "*")!=0)
+		{
+			fprintf(fp1,"@%s\n%s\n+\n%s\n", _r_seq[i].name, _r_seq[i].seq, _r_seq[i].qual);
+		}
+		else if (_r_seq[i].hits[0] == 0)
+		{
+			fprintf(fp1,">%s\n%s\n", _r_seq[i].name, _r_seq[i].seq);
+		}
+	}
+
+	fclose(fp1);
+	if (pairedEndMode)
+		_r_seqCnt *= 2;
+
+	for (i = 0; i < _r_seqCnt; i++)
+	{
+		freeMem(_r_seq[i].hits,0);
+	}
+
+
+	freeMem(_r_seq,0);
+	freeMem(_r_samplingLocs,0);
+}
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/Reads.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/Reads.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,52 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef __READ__
+#define __READ__
+
+typedef struct
+{
+	char *name;
+	char *seq;
+	char *rseq;
+	char *qual;
+	char *hits;
+} Read;
+
+int readAllReads(char *fileName1, char *fileName2, int compressed, unsigned char *fastq, unsigned char pe, Read **seqList, unsigned int *seqListSize);
+void loadSamplingLocations(int **samplingLocs, int *samplingLocsSize);
+void finalizeReads(char *fileName);
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/RefGenome.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/RefGenome.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,161 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include "Common.h"
+#include "RefGenome.h"
+
+FILE *_rg_fp;
+char *_rg_gen;
+char *_rg_name;
+int _rg_offset;
+int _rg_contGen; 
+
+/**********************************************/
+int initLoadingRefGenome(char *fileName)
+{
+	char ch;
+	_rg_fp = fileOpen (fileName, "r");
+	if (fscanf(_rg_fp, "%c", &ch))
+	{
+		if (ch == '>')
+		{
+			_rg_contGen = 0;
+			_rg_offset = 0;
+			_rg_gen = getMem(CONTIG_MAX_SIZE);
+			_rg_name = getMem(CONTIG_NAME_SIZE);
+			return 1;
+		}
+	}
+	return 0;
+}
+/**********************************************/
+void finalizeLoadingRefGenome()
+{
+	freeMem(_rg_gen, CONTIG_MAX_SIZE);
+	freeMem(_rg_name, CONTIG_NAME_SIZE); 
+	fclose(_rg_fp);
+}
+/**********************************************/
+int loadRefGenome(char **refGen, char **refGenName, int *refGenOff)
+{
+	char ch;
+	int i;
+	int returnVal = 0;
+	int actualSize=0;
+	int size;
+	char *tmp;
+	
+	// New Conting 
+	if (!_rg_contGen)
+	{
+		size = 0;
+		tmp = fgets(_rg_name, SEQ_MAX_LENGTH, _rg_fp);
+		int k;
+		for (k=0; k<strlen(_rg_name);k++)
+		{
+			if (_rg_name[k] == ' ')
+			{
+				_rg_name[k]='\0';
+				break;
+			}
+		}
+			
+	}
+	else
+	{
+		size=strlen(_rg_gen);
+		for( i = 0 ; i < CONTIG_OVERLAP ; i++ )
+		{
+			
+			_rg_gen[i] = _rg_gen[size-CONTIG_OVERLAP+i];
+			if (_rg_gen[i] != 'N')
+				actualSize++;
+		}
+		size = CONTIG_OVERLAP;
+	}
+	while( fscanf(_rg_fp, "%c", &ch) > 0 )
+	{
+		if (ch == '>')
+		{
+			_rg_contGen = 0;
+			returnVal = 1;
+			break;
+		}
+		else if (!isspace(ch))
+		{
+			ch = toupper(ch);
+			_rg_gen[size++] = ch;
+			if (ch != 'N')
+			{
+				actualSize++;
+			}
+			if (actualSize == CONTIG_SIZE || size == CONTIG_MAX_SIZE)
+			{
+				_rg_contGen = 1;
+				returnVal=1;
+				break;
+			}
+		}
+
+	}
+
+	_rg_gen[size] = '\0';
+	for (i=strlen(_rg_name)-1; i >= 0; i--)
+		if (!isspace(_rg_name[i]))
+			break;
+	_rg_name[i+1] = '\0';
+
+	*refGenOff = _rg_offset;
+	*refGenName = _rg_name;
+	*refGen = _rg_gen;
+
+	if (_rg_contGen == 1)
+	{
+		_rg_offset += size-CONTIG_OVERLAP;
+	}
+	else
+	{
+		_rg_offset = 0;
+	}
+	
+	
+	return returnVal;
+}
+/**********************************************/
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/RefGenome.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/RefGenome.h	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,44 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-12-08
+ */
+
+
+#ifndef _REF_GENOME_
+#define _REF_GENOME_
+
+int		initLoadingRefGenome(char *fileName);
+void	finalizeLoadingRefGenome();
+int		loadRefGenome(char **refGen, char **refGenName, int *refGenOff);
+
+#endif
diff -r 000000000000 -r ec628ba33878 mrsfast-2.3.0.2/baseFAST.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mrsfast-2.3.0.2/baseFAST.c	Tue Feb 21 10:39:28 2012 -0500
@@ -0,0 +1,398 @@
+/*
+ * Copyright (c) <2008 - 2009>, 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 name of the <ORGANIZATION> 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.
+ */
+
+/*
+ * Author         : Faraz Hach
+ * Email          : fhach AT cs DOT sfu
+ * Last Update    : 2009-02-01
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "Common.h"
+#include "CommandLineParser.h"
+#include "Reads.h"
+#include "Output.h"
+#include "HashTable.h"
+#include "MrsFAST.h"
+
+char 				*versionNumber = "2.3";			// Current Version
+unsigned char		seqFastq;
+
+int main(int argc, char *argv[])
+{
+	if (!parseCommandLine(argc, argv))
+		return 1;
+
+	configHashTable();
+	/****************************************************
+	 * INDEXING
+	 ***************************************************/
+	if (indexingMode)
+	{
+		int i;
+		/********************************
+		 * Regulard Mode
+		 ********************************/
+		if (!bisulfiteMode)
+		{
+			configHashTable();
+			for (i = 0; i < fileCnt; i++)
+			{
+				generateHashTable(fileName[i][0], fileName[i][1]);
+			}
+		}
+		else
+		/********************************
+		 * Bisulfite Mode
+		 ********************************/
+		{ // TODO
+		}
+	}
+	/****************************************************
+	 * SEARCHING
+	 ***************************************************/
+	else
+	{
+		Read *seqList;
+		unsigned int seqListSize;
+		int fc;
+		int samplingLocsSize;
+		int *samplingLocs;
+		double totalLoadingTime = 0;
+		double totalMappingTime = 0;
+		double startTime;
+		double loadingTime;
+		double mappingTime;
+		double lstartTime;
+		double ppTime;
+		double tmpTime;;
+		char *prevGen = getMem(CONTIG_NAME_SIZE);
+		prevGen[0]='\0';
+		char *curGen;
+		int	flag;
+		double maxMem=0;
+		char fname1[FILE_NAME_LENGTH];
+		char fname2[FILE_NAME_LENGTH];
+		char fname3[FILE_NAME_LENGTH];
+		char fname4[FILE_NAME_LENGTH];
+		char fname5[FILE_NAME_LENGTH];
+		// Loading Sequences & Sampling Locations
+		startTime = getTime();
+		if (bisulfiteMode && !pairedEndMode && seqFile1 == NULL)
+		{
+		    //TODO
+		}
+		else
+		{
+			if (!readAllReads(seqFile1, seqFile2, seqCompressed, &seqFastq, pairedEndMode, &seqList, &seqListSize))
+			{
+				return 1;
+			}
+		}
+
+		//loadSamplingLocations(&samplingLocs, &samplingLocsSize);
+		totalLoadingTime += getTime()-startTime;
+
+
+
+
+		if (pairedEndMode)
+		{
+			//Switching to Inferred Size 
+			minPairEndedDistance = minPairEndedDistance - SEQ_LENGTH + 2;
+			maxPairEndedDistance = maxPairEndedDistance - SEQ_LENGTH + 2;
+			if (pairedEndDiscordantMode)
+			{
+				maxPairEndedDiscordantDistance = maxPairEndedDiscordantDistance - SEQ_LENGTH + 2;
+				minPairEndedDiscordantDistance = minPairEndedDiscordantDistance - SEQ_LENGTH + 2;
+			}
+			
+			/* The size between the ends;
+			minPairEndedDistance = minPairEndedDistance + SEQ_LENGTH + 1;
+			maxPairEndedDistance = maxPairEndedDistance + SEQ_LENGTH + 1;
+			if (pairedEndDiscordantMode)
+			{
+				maxPairEndedDiscordantDistance = maxPairEndedDiscordantDistance + SEQ_LENGTH + 1;
+				minPairEndedDiscordantDistance = minPairEndedDiscordantDistance + SEQ_LENGTH + 1;
+			}*/
+			sprintf(fname1, "%s__%s__1", mappingOutputPath, mappingOutput);
+			sprintf(fname2, "%s__%s__2", mappingOutputPath, mappingOutput);
+			sprintf(fname3, "%s__%s__disc", mappingOutputPath, mappingOutput);
+			sprintf(fname4, "%s__%s__oea1", mappingOutputPath, mappingOutput);
+			sprintf(fname5, "%s__%s__oea2", mappingOutputPath, mappingOutput);
+			unlink(fname1);
+			unlink(fname2);
+			unlink(fname3);
+			unlink(fname4);
+			unlink(fname5);
+		}
+
+		// Preparing output
+		initOutput(mappingOutput, outCompressed);
+
+		fprintf(stdout, "-----------------------------------------------------------------------------------------------------------\n");
+		fprintf(stdout, "| %15s | %15s | %15s | %15s | %15s %15s |\n","Genome Name","Loading Time", "Mapping Time", "Memory Usage(M)","Total Mappings","Mapped reads");
+		fprintf(stdout, "-----------------------------------------------------------------------------------------------------------\n");
+
+		/********************************
+		 * Regular Mode
+		 ********************************/
+		if (!bisulfiteMode)
+		{
+			if (!pairedEndMode)
+			{
+				for (fc = 0; fc <fileCnt; fc++)
+				{
+					if (!initLoadingHashTable(fileName[fc][1]))
+					{
+						return 1;
+					}
+
+					loadSamplingLocations(&samplingLocs, &samplingLocsSize);
+
+					mappingTime = 0;
+					loadingTime = 0;
+					prevGen[0] = '\0';
+					flag = 1;
+					
+					do
+					{
+						flag = loadHashTable ( &tmpTime );  			// Reading a fragment
+						curGen = getRefGenomeName();
+
+						// First Time
+						if (flag && prevGen[0]== '\0')
+						{
+							sprintf(prevGen, "%s", curGen);
+						}
+
+						if ( !flag || strcmp(prevGen, curGen)!=0)
+						{
+
+							fprintf(stdout, "| %15s | %15.2f | %15.2f | %15.2f | %15lld %15lld |\n",
+									prevGen,loadingTime, mappingTime, maxMem, mappingCnt , mappedSeqCnt);
+							fflush(stdout);
+
+							totalMappingTime += mappingTime;
+							totalLoadingTime += loadingTime;
+
+							loadingTime = 0;
+							mappingTime = 0;
+							maxMem = 0;
+
+							if (!flag)
+							{
+								break;
+							}
+						}
+						else if (progressRep && mappingTime != 0)
+						{
+							fprintf(stdout, "| %15s | %15.2f | %15.2f | %15.2f | %15lld %15lld |\n",
+									prevGen,loadingTime, mappingTime, maxMem, mappingCnt , mappedSeqCnt);
+							fflush(stdout);
+						}
+
+						sprintf(prevGen, "%s", curGen);
+
+						loadingTime += tmpTime;
+						lstartTime = getTime();
+
+						initFAST(seqList, seqListSize, samplingLocs, samplingLocsSize, fileName[fc][0]);
+
+						mapSingleEndSeq();
+
+						mappingTime += getTime() - lstartTime;
+						if (maxMem < getMemUsage())
+						{
+							maxMem = getMemUsage();
+						}
+					} while (flag);
+
+				} // end for;
+				finalizeFAST();
+				finalizeLoadingHashTable();
+
+			}
+			// Pairedend Mapping Mode
+			else
+			{
+
+				for (fc = 0; fc <fileCnt; fc++)
+				{
+					if (!initLoadingHashTable(fileName[fc][1]))
+					{
+						return 1;
+					}
+					
+
+					loadSamplingLocations(&samplingLocs, &samplingLocsSize);
+					
+					mappingTime = 0;
+					loadingTime = 0;
+					prevGen[0] = '\0';
+					flag = 1;
+
+					do
+					{
+						flag = loadHashTable ( &tmpTime );  			// Reading a fragment
+						curGen = getRefGenomeName();
+
+						// First Time
+						if (flag && prevGen[0]== '\0')
+						{
+							sprintf(prevGen, "%s", curGen);
+						}
+
+						if ( !flag || strcmp(prevGen, curGen)!=0)
+						{
+
+							// DISCORDANT
+							lstartTime = getTime();					
+							outputPairedEnd();
+							mappingTime += getTime() - lstartTime;
+							//DISCORDANT			
+
+							fprintf(stdout, "| %15s | %15.2f | %15.2f | %15.2f | %15lld %15lld |\n",
+									prevGen,loadingTime, mappingTime, maxMem, mappingCnt , mappedSeqCnt);
+							fflush(stdout);
+
+							totalMappingTime += mappingTime;
+							totalLoadingTime += loadingTime;
+
+							loadingTime = 0;
+							mappingTime = 0;
+							maxMem = 0;
+
+							if (!flag)
+							{
+								break;
+							}
+						}
+						else if (progressRep && mappingTime != 0)
+						{
+							fprintf(stdout, "| %15s | %15.2f | %15.2f | %15.2f | %15lld %15lld |\n",
+									prevGen,loadingTime, mappingTime, maxMem, mappingCnt , mappedSeqCnt);
+							fflush(stdout);
+						}
+
+						sprintf(prevGen, "%s", curGen);
+
+						loadingTime += tmpTime;
+						lstartTime = getTime();
+						
+						initFAST(seqList, seqListSize, samplingLocs, samplingLocsSize, fileName[fc][0]);
+						
+						mapPairedEndSeq();
+							
+						mappingTime += getTime() - lstartTime;
+						if (maxMem < getMemUsage())
+						{
+							maxMem = getMemUsage();
+						}
+					} while (flag);
+
+				} // end for;
+
+				finalizeLoadingHashTable();
+				if (pairedEndDiscordantMode)
+				{
+					lstartTime = getTime();							
+					outputPairedEndDiscPP();
+					ppTime = getTime() - lstartTime;
+				}
+				finalizeFAST();
+			}
+		}
+		/********************************
+		 * Bisulfite Mode
+		 ********************************/
+		{
+			//TODO
+		}
+
+
+		finalizeOutput();
+
+		fprintf(stdout, "-----------------------------------------------------------------------------------------------------------\n");
+		fprintf(stdout, "%19s%16.2f%18.2f\n\n", "Total:",totalLoadingTime, totalMappingTime);
+		if (pairedEndDiscordantMode)
+			fprintf(stdout, "Post Processing Time: %18.2f \n", ppTime);
+		fprintf(stdout, "%-30s%10.2f\n","Total Time:", totalMappingTime+totalLoadingTime);
+		fprintf(stdout, "%-30s%10d\n","Total No. of Reads:", seqListSize);
+		fprintf(stdout, "%-30s%10lld\n","Total No. of Mappings:", mappingCnt);
+		fprintf(stdout, "%-30s%10.0f\n\n","Avg No. of locations verified:", ceil((float)verificationCnt/seqListSize));
+
+		int cof = (pairedEndMode)?2:1;
+
+		if (progressRep && maxHits != 0)
+		{
+			int frequency[maxHits+1];
+			int i;
+			for ( i=0 ; i <= maxHits; i++)
+			{
+				frequency[i] = 0;
+			}
+
+
+			for (fc = 0; fc < seqListSize; fc++)
+			{
+				frequency[*(seqList[fc*cof].hits)]++;
+			}
+			frequency[maxHits] = completedSeqCnt;
+			for ( i=0 ; i <= maxHits; i++)
+			{
+				fprintf(stdout, "%-30s%10d%10d%10.2f%%\n","Reads Mapped to ", i, frequency[i], 100*(float)frequency[i]/(float)seqListSize);
+			}
+		}
+
+
+		if (strcmp(unmappedOutput, "") == 0)
+		{
+			char fn[strlen(mappingOutputPath) + strlen(mappingOutput) + 6 ];
+			sprintf(fn, "%s%s.nohit", mappingOutputPath, mappingOutput );
+			finalizeReads(fn);
+		}
+		else
+		{
+			finalizeReads(unmappedOutput);
+		}
+
+
+
+		freeMem(prevGen, CONTIG_NAME_SIZE);
+	}
+
+
+
+	return 1;
+}