Makefile

 CFLAGS = -Wall -O3 -std=c99
-LIBS = align.o alnfrags.o ankers.o assembly.o chain.o compdna.o compkmers.o compress.o decon.o ef.o filebuff.o frags.o hashmap.o hashmapindex.o hashmapkma.o hashmapkmers.o hashtable.o index.o kma.o kmapipe.o kmeranker.o kmers.o kmmap.o loadupdate.o makeindex.o mt1.o nw.o pherror.o printconsensus.o qseqs.o qualcheck.o runinput.o runkma.o sam.o savekmers.o seq2fasta.o seqmenttree.o seqparse.o shm.o sparse.o spltdb.o stdnuc.o stdstat.o tmp.o update.o updateindex.o updatescores.o valueshash.o vcf.o
+LIBS = align.o alnfrags.o ankers.o assembly.o chain.o compdna.o compkmers.o compress.o decon.o ef.o filebuff.o frags.o hashmap.o hashmapindex.o hashmapkma.o hashmapkmers.o hashtable.o index.o kma.o kmapipe.o kmeranker.o kmerlink.o kmers.o kmmap.o loadupdate.o makeindex.o mt1.o nw.o pherror.o printconsensus.o qseqs.o qualcheck.o runinput.o runkma.o sam.o savekmers.o seq2fasta.o seqmenttree.o seqparse.o shm.o sparse.o spltdb.o stdnuc.o stdstat.o tmp.o update.o updateindex.o updatescores.o valueshash.o vcf.o
 PROGS = kma kma_index kma_shm kma_update
 
 .c .o:
@@ -46,6 +46,7 @@ index.o: index.h compress.h decon.h hashmap.h hashmapkma.h loadupdate.h makeinde
 kma.o: kma.h ankers.h assembly.h chain.h hashmapkma.h kmers.h mt1.h penalties.h pherror.h qseqs.h runinput.h runkma.h savekmers.h sparse.h spltdb.h tmp.h version.h
 kmapipe.o: kmapipe.h pherror.h
 kmeranker.o: kmeranker.h penalties.h
+kmerlink.o: kmerlink.h kmeranker.h
 kmers.o: kmers.h ankers.h compdna.h hashmapkma.h kmapipe.h pherror.h qseqs.h savekmers.h spltdb.h
 kmmap.o: kmmap.h hashmapkma.h
 loadupdate.o: loadupdate.h pherror.h hashmap.h hashmapkma.h updateindex.h

align.c

@@ -210,7 +210,6 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 	start = chainSeedsPtr(points, q_len, t_len, kmersize, &aligned->mapQ);
 	score = points->score[start];
 	
-	//fprintf(stderr, "%d\t%d\t%d\t%d\n", t_len, mem_count, aligned->mapQ, score);
 	if(aligned->mapQ < mq || score < kmersize) {
 		Stat.score = 0;
 		Stat.len = 1;
@@ -250,6 +249,8 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 				NWstat = NW(template_index->seq, qseq, -1 - (t_s == 0), t_s, t_e, q_s, q_e, Frag_align, matrices);
 			} else {
 				NWstat = NW_band(template_index->seq, qseq, -1 - (t_s == 0), t_s, t_e, q_s, q_e, Frag_align, band, matrices);
+				/* here */
+				//NWstat = NW(template_index->seq, qseq, -1 - (t_s == 0), t_s, t_e, q_s, q_e, Frag_align, matrices);
 			}
 			
 			/* trim leading gaps */
@@ -343,7 +344,7 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 					NWstat = NW(template_index->seq, qseq, 0, t_s, t_e, q_s, q_e, Frag_align, matrices);
 				} else {
 					NWstat = NW_band(template_index->seq, qseq, 0, t_s, t_e, q_s, q_e, Frag_align, band, matrices);
-					//NWstat = NW(template_index->seq, qseq, 0, t_s, t_e, q_s, q_e, Frag_align);
+					//NWstat = NW(template_index->seq, qseq, 0, t_s, t_e, q_s, q_e, Frag_align, matrices);
 				}
 				
 				memcpy(aligned->t + Stat.len, Frag_align->t, NWstat.len);
@@ -377,7 +378,7 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 			NWstat = NW(template_index->seq, qseq, 1 + (t_e == t_len), t_s, t_e, q_s, q_e, Frag_align, matrices);
 		} else {
 			NWstat = NW_band(template_index->seq, qseq, 1 + (t_e == t_len), t_s, t_e, q_s, q_e, Frag_align, band, matrices);
-			//NWstat = NW(template_index->seq, qseq, 1 + (t_e == t_len), t_s, t_e, q_s, q_e, Frag_align);
+			//NWstat = NW(template_index->seq, qseq, 1 + (t_e == t_len), t_s, t_e, q_s, q_e, Frag_align, matrices);
 		}
 		/* trim trailing gaps */
 		if(t_e == t_len) {
@@ -390,12 +391,12 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 				--bias;
 			}
 			++bias;
-			/*
+			
+			
 			if(bias != NWstat.len) {
-				NWstat.score -= (W1 + (NWstat.len - bias) * U);
+				//NWstat.score -= (W1 + (NWstat.len - bias) * U);
 				NWstat.len = bias;
 			}
-			*/
 		}
 		
 		memcpy(aligned->t + Stat.len, Frag_align->t, NWstat.len);

assembly.c

@@ -274,7 +274,6 @@ void * assemble_KMA_threaded(void *arg) {
 	
 	static volatile int excludeIn[1] = {0}, excludeOut[1] = {0}, excludeMatrix[1] = {0}, mainTemplate = -2, thread_wait = 0;
 	static char *template_name;
-	static HashMap_index *template_index;
 	Assemble_thread *thread = arg;
 	int i, j, t_len, aln_len, start, end, bias, myBias, gaps, pos, spin, sam;
 	int read_score, depthUpdate, bestBaseScore, bestScore, template, asm_len;
@@ -295,6 +294,7 @@ void * assemble_KMA_threaded(void *arg) {
 	AssemInfo *matrix;
 	AlnPoints *points;
 	NWmat *NWmatrices;
+	HashMap_index *template_index;
 	
 	/* get input */
 	template = thread->template;
@@ -317,6 +317,7 @@ void * assemble_KMA_threaded(void *arg) {
 	sam = thread->sam;
 	spin = thread->spin;
 	thread_num = thread->thread_num;
+	template_index = thread->template_index;
 	
 	if(template != -2) {
 		/* all assemblies done, 
@@ -331,7 +332,6 @@ void * assemble_KMA_threaded(void *arg) {
 		/* Allocate assembly arrays */
 		lock(excludeMatrix);
 		template_name = thread->template_name;
-		template_index = thread->template_index;
 		t_len = template_index->len;
 		matrix->len = t_len;
 		if(matrix->size < (t_len << 1)) {
@@ -552,6 +552,7 @@ void * assemble_KMA_threaded(void *arg) {
 								header->seq[header->len - 1] = 0;
 								samwrite(qseq, header, 0, template_name, aligned, stats);
 							}
+							
 						} else if(sam) {
 							stats[1] = read_score;
 							stats[2] = start;
@@ -617,7 +618,6 @@ void * assemble_KMA_threaded(void *arg) {
 		
 		return NULL;
 	}
-	
 	/* diff */
 	/* pre on dense */
 	/* Pepare and make alignment on consensus */
@@ -697,6 +697,7 @@ void * assemble_KMA_threaded(void *arg) {
 		pos = assembly[pos].next;
 	}
 	
+
 	/* Trim alignment on consensus */
 	coverScore = 0;
 	bias = 0;

compress.c

@@ -20,6 +20,7 @@
 #include <limits.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <sys/param.h>
 #include "compress.h"
 #include "hashmap.h"
 #include "hashmapkma.h"
@@ -44,6 +45,7 @@ static void mmapinit(HashMapKMA *finalDB, long unsigned size, FILE *out) {
 	if(finalDB->exist == MAP_FAILED) {
 		ERROR();
 	}
+	posix_madvise(finalDB->exist, size, POSIX_MADV_SEQUENTIAL);
 	*finalDB->exist++ = finalDB->DB_size;
 	*finalDB->exist++ = finalDB->kmersize;
 	*finalDB->exist++ = finalDB->prefix_len;
@@ -72,8 +74,8 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 	void *data;
 	long unsigned i, j, check, size, v_size;
 	long unsigned index, t_index, v_index, new_index, null_index;
-	unsigned swap, *values;
-	short unsigned *values_s;
+	unsigned swap, *values, *finalV;
+	short unsigned *values_s, *finalV_s;
 	HashMapKMA *finalDB;
 	ValuesHash *shmValues;
 	ValuesTable *node, *next, *table;
@@ -227,6 +229,7 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 						if(finalDB->exist == MAP_FAILED) {
 							ERROR();
 						}
+						posix_madvise(finalDB->exist, size, POSIX_MADV_SEQUENTIAL);
 						finalDB->exist_l = (long unsigned *) (finalDB->exist + 3);
 						finalDB->exist_l += 5;
 						
@@ -470,6 +473,17 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 	/* move values */
 	if(finalDB->DB_size < USHRT_MAX) {
 		for(node = table; node != 0; node = next) {
+			next = node->next;
+			values_s = (short unsigned *)(node->values);
+			finalV_s = finalDB->values_s + node->v_index - 1;
+			i = 2 + *values_s--;
+			while(--i) {
+				*++finalV_s = *++values_s;
+			}
+			free(node->values);
+			free(node);
+			
+			/*
 			next = node->next;
 			values_s = (short unsigned *)(node->values);
 			for(i = node->v_index, j = 0; j <= *values_s; ++i, ++j) {
@@ -477,9 +491,21 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 			}
 			free(values_s);
 			free(node);
+			*/
 		}
 	} else {
 		for(node = table; node != 0; node = next) {
+			next = node->next;
+			values = node->values;
+			finalV = finalDB->values + node->v_index - 1;
+			i = 2 + *values--;
+			while(--i) {
+				*++finalV = *++values;
+			}
+			free(node->values);
+			free(node);
+			
+			/*
 			next = node->next;
 			values = node->values;
 			for(i = node->v_index, j = 0; j <= *values; ++i, ++j) {
@@ -487,6 +513,7 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 			}
 			free(values);
 			free(node);
+			*/
 		}
 	}
 	
@@ -537,8 +564,8 @@ HashMapKMA * compressKMA_DB(HashMap *templates, FILE *out) {
 HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 	
 	long unsigned i, j, v_index, new_index, null_index, size, v_size;
-	unsigned check, swap, *values;
-	short unsigned *values_s;
+	unsigned check, swap, *values, *finalV;
+	short unsigned *values_s, *finalV_s;
 	HashMapKMA *finalDB;
 	ValuesHash *shmValues;
 	ValuesTable *node, *next, *table;
@@ -653,6 +680,7 @@ HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 			if(finalDB->exist == MAP_FAILED) {
 				ERROR();
 			}
+			posix_madvise(finalDB->exist, size, POSIX_MADV_SEQUENTIAL);
 			finalDB->exist_l = (long unsigned *) (finalDB->exist + 3);
 			finalDB->exist_l += 5;
 		} else {
@@ -810,6 +838,7 @@ HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 		if(finalDB->values == MAP_FAILED) {
 			ERROR();
 		}
+		posix_madvise(finalDB->values, size, POSIX_MADV_SEQUENTIAL);
 		finalDB->values = ((void *) finalDB->values) + v_size;
 		if(finalDB->DB_size < USHRT_MAX) {
 			finalDB->values_s = (short unsigned *) finalDB->values;
@@ -819,6 +848,17 @@ HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 	/* move values */
 	if(finalDB->DB_size < USHRT_MAX) {
 		for(node = table; node != 0; node = next) {
+			next = node->next;
+			values_s = (short unsigned *)(node->values);
+			finalV_s = finalDB->values_s + node->v_index - 1;
+			i = 2 + *values_s--;
+			while(--i) {
+				*++finalV_s = *++values_s;
+			}
+			free(node->values);
+			free(node);
+			
+			/*
 			next = node->next;
 			values_s = (short unsigned *)(node->values);
 			for(i = node->v_index, j = 0; j <= *values_s; ++i, ++j) {
@@ -826,9 +866,21 @@ HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 			}
 			free(values_s);
 			free(node);
+			*/
 		}
 	} else {
 		for(node = table; node != 0; node = next) {
+			next = node->next;
+			values = node->values;
+			finalV = finalDB->values + node->v_index - 1;
+			i = 2 + *values--;
+			while(--i) {
+				*++finalV = *++values;
+			}
+			free(node->values);
+			free(node);
+			
+			/*
 			next = node->next;
 			values = node->values;
 			for(i = node->v_index, j = 0; j <= *values; ++i, ++j) {
@@ -836,6 +888,7 @@ HashMapKMA * compressKMA_megaDB(HashMap *templates, FILE *out) {
 			}
 			free(values);
 			free(node);
+			*/
 		}
 	}
 	/* dump final DB */

debian/changelog

+kma (1.2.12-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream version
+
+ -- Steffen Moeller <moeller@debian.org>  Mon, 23 Sep 2019 09:58:31 +0200
+
 kma (1.2.11-1) unstable; urgency=medium
 
   * Team upload.

debian/patches/hardening.patch

 Author: Andreas Tille <tille@debian.org>
 Last-Update: Thu, 07 Jun 2018 21:55:00 +0200
 Description: Propagate hardening options
+ Brought upstream as https://bitbucket.org/genomicepidemiology/kma/pull-requests/5
 
 Index: kma/Makefile
 ===================================================================
@@ -9,7 +10,7 @@ Index: kma/Makefile
 @@ -1,4 +1,4 @@
 -CFLAGS = -Wall -O3 -std=c99
 +CFLAGS += -Wall -O3 -std=c99
- LIBS = align.o alnfrags.o ankers.o assembly.o chain.o compdna.o compkmers.o compress.o decon.o ef.o filebuff.o frags.o hashmap.o hashmapindex.o hashmapkma.o hashmapkmers.o hashtable.o index.o kma.o kmapipe.o kmeranker.o kmers.o kmmap.o loadupdate.o makeindex.o mt1.o nw.o pherror.o printconsensus.o qseqs.o qualcheck.o runinput.o runkma.o sam.o savekmers.o seq2fasta.o seqmenttree.o seqparse.o shm.o sparse.o spltdb.o stdnuc.o stdstat.o tmp.o update.o updateindex.o updatescores.o valueshash.o vcf.o
+ LIBS = align.o alnfrags.o ankers.o assembly.o chain.o compdna.o compkmers.o compress.o decon.o ef.o filebuff.o frags.o hashmap.o hashmapindex.o hashmapkma.o hashmapkmers.o hashtable.o index.o kma.o kmapipe.o kmeranker.o kmerlink.o kmers.o kmmap.o loadupdate.o makeindex.o mt1.o nw.o pherror.o printconsensus.o qseqs.o qualcheck.o runinput.o runkma.o sam.o savekmers.o seq2fasta.o seqmenttree.o seqparse.o shm.o sparse.o spltdb.o stdnuc.o stdstat.o tmp.o update.o updateindex.o updatescores.o valueshash.o vcf.o
  PROGS = kma kma_index kma_shm kma_update
  
 @@ -8,16 +8,16 @@ PROGS = kma kma_index kma_shm kma_update

debian/upstream/metadata

 Reference:
  Author: Philip T. L. C. Clausen and Frank M. Aarestrup and Ole Lund
- Title: Rapid and precise alignment of raw reads against redundant databases with KMA
+ Title: >
+  Rapid and precise alignment of raw reads against redundant databases
+  with KMA
  Journal: BMC Bioinformatics
  Year: 2018
  Volume: 19
@@ -15,3 +17,5 @@ Registry:
    Entry: OMICS_31606
  - Name: bio.tools
    Entry: NA
+ - Name: conda:bioconda
+   Entry: kma

hashmapindex.c

@@ -37,17 +37,26 @@ typedef int key_t;
 void (*destroyPtr)(HashMap_index *) = &alignClean;
 HashMap_index * (*alignLoadPtr)(HashMap_index *, int, int, int, int, long unsigned, long unsigned) = &alignLoad_fly;
 
-void hashMap_index_initialize(HashMap_index *dest, int len, int kmerindex) {
+int hashMap_index_initialize(HashMap_index *dest, int len, int kmerindex) {
 	
 	dest->len = len;
-	dest->size = len << 1;
 	dest->kmerindex = kmerindex;
-	
+	if(dest->size < (len << 1)) {
+		if(dest->size) {
+			free(dest->index);
+			free(dest->seq);
+		}
+		dest->size = len << 1;
 		dest->index = calloc(dest->size, sizeof(int));
 		dest->seq = malloc(((len >> 5) + 1) * sizeof(long unsigned));
 		if(!dest->index || !dest->seq) {
 			ERROR();
 		}
+		
+		return 0;
+	}
+	
+	return 1;
 }
 
 void hashMap_index_set(HashMap_index *dest) {
@@ -248,8 +257,13 @@ HashMap_index * hashMap_index_load(HashMap_index *src, int seq, int index, int l
 	
 	if(src == 0) {
 		src = smalloc(sizeof(HashMap_index));
+		src->size = 0;
+	}
+	if(hashMap_index_initialize(src, len, kmersize)) {
+		src->size = len << 1;
+		memset(src->index, 0, src->size * sizeof(int));
 	}
-	hashMap_index_initialize(src, len, kmersize);
+	
 	read(seq, src->seq, ((src->len >> 5) + 1) * sizeof(long unsigned));
 	read(index, src->index, src->size * sizeof(int));
 	
@@ -267,8 +281,13 @@ HashMap_index * hashMap_index_build(HashMap_index *src, int seq, int len, int km
 	
 	if(src == 0) {
 		src = smalloc(sizeof(HashMap_index));
+		src->size = 0;
+	}
+	
+	if(hashMap_index_initialize(src, len, kmersize)) {
+		memset(src->index, 0, src->size * sizeof(int));
 	}
-	hashMap_index_initialize(src, len, kmersize);
+	
 	shifter = sizeof(long unsigned) * sizeof(long unsigned) - (src->kmerindex << 1);
 	
 	read(seq, src->seq, ((src->len >> 5) + 1) * sizeof(long unsigned));
@@ -375,10 +394,12 @@ HashMap_index * alignLoad_shm_initial(char *templatefilename, int file_len, int
 }
 
 void alignClean(HashMap_index *template_index) {
-	
+	/* Overwrite later */
+	/*
 	if(template_index) {
 		hashMap_index_destroy(template_index);
 	}
+	*/
 }
 
 void alignClean_shm(HashMap_index *template_index) {

hashmapindex.h

@@ -35,7 +35,7 @@ struct hashMap_index {
 extern void (*destroyPtr)(HashMap_index *);
 extern HashMap_index * (*alignLoadPtr)(HashMap_index *, int, int, int, int, long unsigned, long unsigned);
 
-void hashMap_index_initialize(HashMap_index *dest, int len, int kmerindex);
+int hashMap_index_initialize(HashMap_index *dest, int len, int kmerindex);
 void hashMap_index_set(HashMap_index *dest);
 void hashMap_index_destroy(HashMap_index *dest);
 int hashMap_index_get(const HashMap_index *dest, long unsigned key, unsigned shifter);

kma.c

@@ -173,7 +173,7 @@ int kma_main(int argc, char *argv[]) {
 	static unsigned nc, nf, shm, exhaustive;
 	static char *outputfilename, *templatefilename, **templatefilenames;
 	static char **inputfiles, **inputfiles_PE, **inputfiles_INT, ss;
-	static double ID_t, scoreT, evalue;
+	static double ID_t, scoreT, coverT, evalue;
 	static Penalties *rewards;
 	int i, j, args, exe_len, status, size, escape, tmp, step1, step2;
 	unsigned totFrags;
@@ -221,6 +221,7 @@ int kma_main(int argc, char *argv[]) {
 		mq = 0;
 		bcd = 1;
 		scoreT = 0.5;
+		coverT = 0.5;
 		ID_t = 1.0;
 		one2one = 0;
 		ss = 'q';
@@ -264,7 +265,6 @@ int kma_main(int argc, char *argv[]) {
 		deConPrintPtr = printPtr;
 		*/
 		
-		
 		/* PARSE COMMAND LINE OPTIONS */
 		args = 1;
 		while(args < argc) {
@@ -1137,7 +1137,7 @@ int kma_main(int argc, char *argv[]) {
 	} else if(step2) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);
 		strcpy(myTemplatefilename, templatefilename);
-		status = save_kmers_batch(myTemplatefilename, "-s1", shm, thread_num, exhaustive, rewards, ioStream, sam);
+		status = save_kmers_batch(myTemplatefilename, "-s1", shm, thread_num, exhaustive, rewards, ioStream, sam, scoreT, coverT);
 		free(myTemplatefilename);
 	} else if(sparse_run) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);

kmapipe.c

@@ -226,13 +226,14 @@ FILE * kmaPipeFork(const char *cmd, const char *type, FILE *ioStream, int *statu
 		for (last = 0, src = pidlist; src->fp != ioStream; last = src, src = src->next) {
 			if(!src) {
 				*status = 1;
+				unlock(lock);
 				return 0;
 			}
 		}
 		unlock(lock);
 		
 		/* close stream and get exit status */
-		*status = 1;
+		*status = 0;
 		#ifndef _WIN32
 		while ((pid = waitpid(src->pid, status, 0)) == -1 && errno == EINTR) {
 			usleep(100);
@@ -240,7 +241,6 @@ FILE * kmaPipeFork(const char *cmd, const char *type, FILE *ioStream, int *statu
 		#else
 		WaitForSingleObject(src->pid, INFINITE);
 		#endif
-		*status = 0;
 		fclose(ioStream);
 		
 		/* Remove the entry from the linked list. */

kmeranker.c

@@ -18,14 +18,15 @@
 */
 #define _XOPEN_SOURCE 600
 #include <stdlib.h>
+#include <math.h>
 #include "kmeranker.h"
 #include "penalties.h"
 
-KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int *bests, int *Score, int *extendScore, char *include) {
+KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int kmersize, int *bests, int *Score, int *extendScore, char *include) {
 	
 	/* get set of best templates and silences the chain, except for the initial anker */
-	int template, score, bestScore, Wl, W1, U, M, MM, gaps;
-	unsigned i, j, start, end, SU;
+	int i, j, template, score, bestScore, Wl, W1, U, M, MM, gaps;
+	int start, end, pos, SU;
 	unsigned *values;
 	short unsigned *values_s;
 	KmerAnker *node, *prev;
@@ -35,8 +36,9 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 		return 0;
 	} else if((SU = *include)) {
 		values_s = (short unsigned *) src->values;
-		i = *values_s + 1;
-		*bests += i;
+		i = *values_s;
+		*bests = i;
+		++i;
 		values_s += i;
 		bests += i;
 		while(--i) {
@@ -44,23 +46,27 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 		}
 	} else {
 		values = src->values;
-		i = *values + 1;
-		*bests += i;
+		i = *values;
+		*bests = i;
+		++i;
 		values += i;
 		bests += i;
 		while(--i) {
 			include[(*--bests = *--values)] ^= 1;
 		}
 	}
+	--bests;
 	M = rewards->M;
 	MM = rewards->MM;
 	U = rewards->U;
 	W1 = rewards->W1;
 	Wl = rewards->Wl;
 	bestScore = src->score;
+	values_s = 0;
+	values = 0;
 	
 	/* get chaining scores */
-	for(node = src; node != 0; node = node->next) {
+	for(node = src; node; node = node->next) {
 		if(SU) {
 			values_s = (short unsigned *) node->values;
 			i = *values_s + 1;
@@ -73,41 +79,44 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 		
 		start = node->start;
 		end = node->end;
-		
 		while(--i) {
 			template = SU ? *--values_s : *--values;
-			if(!include[template]) {
+			if(include[template]) {
 				score = Score[template];
-				gaps = start - extendScore[template];
-				if(gaps == 0) {
-					/* continued match */
-					score += M;
-				} else if(gaps < 0) {
-					/* deletion */
-					score += (W1 + abs(gaps + 1) * U);
-				} else if(gaps == 1) {
-					/* one unmatched between */
-					score += MM;
+				pos = extendScore[template];
+				gaps = pos - end;
+				
+				/* extend chain */
+				if(pos == 0) {
+					if(node->cStart) {
+						score = W1 + (node->cStart - 1) * U;
+						score = node->weight + (Wl < score ? score : Wl);
 					} else {
-					/* unmatched bases between */
-					if((gaps = (gaps - 2) * M) < 0) {
-						score += gaps;
-					}
-					score += (MM << 1);
+						score = node->weight;
 					}
-				
-				/* check local chain is needed */
-				if(score < Wl) {
-					score = Wl + node->weight;
-					/* mark as terminating */
-					extendScore[template] = 0;
 				} else {
-					score += node->weight;
-					/* mark descendant */
-					extendScore[template] = end;
+					if(gaps < 0) {
+						if(gaps == -kmersize) {
+							score += node->weight - (kmersize - 1) * M;
+						} else {
+							score += (W1 + (-gaps - 1) * U) + node->weight + gaps * M;
+						}
+					} else if(gaps == 0) {
+						score += node->weight + W1;
+					} else if(gaps <= 2) {
+						score += node->weight + gaps * MM;
+					} else if((MM * 2 + (gaps - 2) * M) < 0) {
+						score += node->weight + (MM * 2 + (gaps - 2) * M);
+					} else {
+						score += node->weight - log(gaps * M);
 					}
 					
-				/* update Scores */
+					/* verify extension */
+					if(score < 0 || score == bestScore) {
+						include[template] = 0;
+					}
+				}
+				extendScore[template] = start;
 				Score[template] = score;
 			}
 		}
@@ -124,6 +133,7 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 		if(Score[(template = bests[i])] == bestScore) {
 			bests[++j] = template;
 		}
+		
 		/* clear Score arrays */
 		Score[template] = 0;
 		extendScore[template] = 0;
@@ -134,18 +144,22 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 	return prev;
 }
 
-KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize) {
+KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize, double mrs) {
 	
 	KmerAnker *node, *prev;
 	
-	while(V_score->score < kmersize) {
-		++V_score;
+	if(!V_score || !V_score->score) {
+		return 0;
 	}
+	while((V_score->score < kmersize || V_score->score < mrs * (V_score->end - V_score->cStart)) && (V_score = V_score->descend));
 	
+	if(!V_score) {
+		return 0;
+	}
 	prev = V_score;
 	node = V_score->descend;
 	while(node) {
-		if(kmersize <= node->score) {
+		if(kmersize <= node->score && mrs * (node->end - node->cStart) <= node->score) {
 			prev->descend = node;
 			prev = node;
 		}
@@ -156,14 +170,6 @@ KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize) {
 	return V_score;
 }
 
-unsigned getStartAnker(KmerAnker *src) {
-	
-	while(src->next) {
-		src = src->next;
-	}
-	return src->start;
-}
-
 KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties) {
 	
 	KmerAnker *node, *prev, *best;
@@ -174,18 +180,19 @@ KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties) {
 		prev = prev->descend;
 	}
 	*src = prev;
-	best = prev;
+	if(!(best = prev)) {
+		return 0;
+	}
 	node = prev->descend;
-	
 	while(node) {
 		if(best->score < node->score) {
 			best = node;
 			*ties = 0;
 		} else if(best->score == node->score) {
-			if((node->end - node->start) < (best->end - best->start)) {
+			if((node->end - node->start) < (best->end - best->cStart)) {
 				best = node;
 				*ties = 0;
-			} else if((node->end - node->start) == (best->end - best->start)) {
+			} else if((node->end - node->start) == (best->end - best->cStart)) {
 				best = node;
 				++*ties;
 			}
@@ -200,11 +207,15 @@ KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties) {
 	return best;
 }
 
-KmerAnker * getTieAnker(KmerAnker *stop, KmerAnker *src, int score, int len) {
+KmerAnker * getTieAnker(int stop, KmerAnker *src, int score, int len) {
+	
+	if(!src || src->start == stop) {
+		return 0;
+	}
 	
 	/* search downwards */
-	while(--src != stop) {
-		if(src->score == score && (src->end - src->start) == len) {
+	while((--src)->start != stop) {
+		if(src->score == score && (src->end - src->cStart) == len) {
 			return src;
 		}
 	}

kmeranker.h

@@ -27,15 +27,15 @@ struct kmerAnker {
 	int weight;
 	unsigned start;
 	unsigned end;
-	unsigned flag;
+	unsigned cStart;
+	unsigned ties;
 	unsigned *values;
 	struct kmerAnker *descend; /* descending anker */
-	struct kmerAnker *next; /* next anker in chain */
+	struct kmerAnker *next; /* next ankers in chain */
 };
 #endif
 
-KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int *bests, int *Score, int *extendScore, char *include);
-KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize);
-unsigned getStartAnker(KmerAnker *src);
+KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int kmersize, int *bests, int *Score, int *extendScore, char *include);
+KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize, double mrs);
 KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties);
-KmerAnker * getTieAnker(KmerAnker *stop, KmerAnker *src, int score, int len);
+KmerAnker * getTieAnker(int stop, KmerAnker *src, int score, int len);

kmerlink.c

+/* Philip T.L.C. Clausen Jan 2017 plan@dtu.dk */
+
+/*
+ * Copyright (c) 2017, Philip Clausen, Technical University of Denmark
+ * All rights reserved.
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * 
+ *		http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+*/
+#define _XOPEN_SOURCE 600
+#include <stdlib.h>
+#include "kmerlink.h"
+#include "pherror.h"
+
+KmerLink * initKmerLink(unsigned size) {
+	
+	KmerLink *dest;
+	
+	dest = smalloc(sizeof(KmerLink));
+	dest->size = size;
+	dest->n = 0;
+	dest->list = smalloc(size * sizeof(KmerAnker));
+	
+	return dest;
+}
+
+void reallocKmerLink(KmerLink *src, unsigned newSize) {
+	
+	src->size = newSize;
+	src->list = realloc(src->list, newSize * sizeof(KmerAnker));
+	if(!src->list) {
+		ERROR();
+	}
+}
+
+KmerAnker * pushKmerLink(KmerLink *src, KmerAnker *node) {
+	
+	KmerAnker *dest;
+	
+	if(src->size == ++src->n) {
+		reallocKmerLink(src, src->size << 1);
+	}
+	dest = src->list + (src->n - 1);
+	*dest = *node;
+	
+	return dest;
+}
+
+KmerAnker * popKmerLink(KmerLink *src, int n) {
+	
+	if(n < src->n) {
+		src->n -= n;
+	} else {
+		src->n = 0;
+	}
+	
+	return src->list + src->n;
+}
+
+void destroyKmerLink(KmerLink *src) {
+	
+	free(src->list);
+	free(src);
+}

kmerlink.h

+/* Philip T.L.C. Clausen Jan 2017 plan@dtu.dk */
+
+/*
+ * Copyright (c) 2017, Philip Clausen, Technical University of Denmark
+ * All rights reserved.
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * 
+ *		http://www.apache.org/licenses/LICENSE-2.0
+ * 
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+*/
+#define _XOPEN_SOURCE 600
+#include "kmeranker.h"
+
+#ifndef KMERLINK
+#define KMERLINK 1;
+typedef struct kmerLink KmerLink;
+struct kmerLink {
+	unsigned size;
+	unsigned n;
+	KmerAnker *list;
+};
+#endif
+
+KmerLink * initKmerLink(unsigned size);
+void reallocKmerLink(KmerLink *src, unsigned newSize);
+KmerAnker * pushKmerLink(KmerLink *src, KmerAnker *node);
+KmerAnker * popKmerLink(KmerLink *src, int n);
+void destroyKmerLink(KmerLink *src);

kmers.c

@@ -47,7 +47,7 @@ typedef int key_t;
 #define shmctl(shmid, cmd, buf) fprintf(stderr, "sysV not available on Windows.\n")
 #endif
 
-int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam) {
+int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam, double mrs, double coverT) {
 	
 	int i, file_len, shmid, deCon, *bestTemplates, *template_lengths;
 	FILE *inputfile, *templatefile;
@@ -116,6 +116,7 @@ int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int th
 	if(printPtr == &print_ankers_spltDB || printPtr == &print_ankers_Sparse_spltDB) {
 		printPtr(0, 0, thread_num, 0, 0, 0);
 	}
+	template_lengths = 0;
 	if(kmerScan == &save_kmers_HMM) {
 		/* load lengths */
 		strcat(templatefilename, ".length.b");
@@ -138,8 +139,8 @@ int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int th
 		templatefilename[file_len] = 0;
 		fclose(templatefile);
 		save_kmers_HMM(templates, 0, &(int){thread_num}, template_lengths, 0, 0, 0, 0, 0, 0, 0, 0, 0);
-	} else {
-		template_lengths = 0;
+	} else if(kmerScan == &save_kmers_chain || kmerScan == &save_kmers_sparse_chain) {
+		kmerScan(0, 0, &(int){thread_num}, (int *)(&mrs), (int *)(&coverT), 0, 0, 0, 0, 0, 0, 0, 0);
 	}
 	
 	t1 = clock();
@@ -234,6 +235,8 @@ int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int th
 	}
 	if(kmerScan == &save_kmers_HMM && (shm & 4) == 0) {
 		free(template_lengths);
+	} else if(kmerScan == &save_kmers_chain || kmerScan == &save_kmers_sparse_chain) {
+		kmerScan(0, 0, &(int){thread_num}, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 	}
 	for(thread = threads; thread; thread = threads) {
 		threads = thread->next;

kmers.h

@@ -19,4 +19,4 @@
 #define _XOPEN_SOURCE 600
 #include "penalties.h"
 
-int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam);
+int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam, double mrs, double coverT);

nw.c

@@ -414,6 +414,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
 	
 	/* Perform banded NW */
 	pos[0] = 0;
+	pos[1] = 0;
 	en = 0;
 	c_pos = (t_len + q_len) >> 1;
 	for(m = t_len - 1, nuc_pos = t_e - 1; m >= 0; --m, --nuc_pos, --c_pos) {
@@ -532,6 +533,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
 		if(k < 0 && Stat.score <= D_ptr[n]) {
 			Stat.score = D_ptr[n];
 			pos[0] = m;
+			pos[1] = n;
 		}
 		
 		tmp = D_ptr;
@@ -546,12 +548,15 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
 	
 	/* get start position of alignment */
 	if(k < 0) {
+		q_pos = pos[0];
+		if(pos[0] == 0) {
 			pos[1] = en;
-		q_pos = 0;
+		}
 		if(k == -2) {
 			for(n = en; n < bq_len; ++n) {
 				if(D_prev[n] > Stat.score) {
 					Stat.score = D_prev[n];
+					q_pos = 0;
 					pos[0] = 0;
 					pos[1] = (aligned->start = n);
 				}
@@ -563,7 +568,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
 		pos[0] = 0;
 		pos[1] = en;
 	}
-	
+	aligned->start = q_pos;
 	
 	/* back tracking */
 	m = pos[0];
@@ -960,6 +965,7 @@ AlnScore NW_band_score(const long unsigned *template, const unsigned char *query
 	
 	/* Perform banded NW */
 	pos[0] = 0;
+	pos[1] = 0;
 	en = 0;
 	c_pos = (t_len + q_len) >> 1;
 	for(m = t_len - 1, nuc_pos = t_e - 1; m >= 0; --m, --nuc_pos, --c_pos) {
@@ -1078,6 +1084,7 @@ AlnScore NW_band_score(const long unsigned *template, const unsigned char *query
 		if(k < 0 && Stat.score <= D_ptr[n]) {
 			Stat.score = D_ptr[n];
 			pos[0] = m;
+			pos[1] = n;
 		}
 		
 		tmp = D_ptr;
@@ -1092,12 +1099,15 @@ AlnScore NW_band_score(const long unsigned *template, const unsigned char *query
 	
 	/* get start position of alignment */
 	if(k < 0) {
+		q_pos = pos[0];
+		if(pos[0] == 0) {
 			pos[1] = en;
-		q_pos = 0;
+		}
 		if(k == -2) {
 			for(n = en; n < bq_len; ++n) {
 				if(D_prev[n] > Stat.score) {
 					Stat.score = D_prev[n];
+					q_pos = 0;
 					pos[0] = 0;
 					pos[1] = n;
 				}

qseqs.c

@@ -42,8 +42,8 @@ void insertKmerBound(Qseqs *header, int start, int end) {
 	
 	int *seq;
 	
-	if((header->len + 3 * sizeof(int)) < header->size) {
-		header->size = (header->len + 3 * sizeof(int)) << 1;
+	if((header->len + 2 * sizeof(int)) < header->size) {
+		header->size = (header->len + 2 * sizeof(int)) << 1;
 		if(!(header->seq = realloc(header->seq, header->size))) {
 			ERROR();
 		}
@@ -52,7 +52,8 @@ void insertKmerBound(Qseqs *header, int start, int end) {
 	seq = (int *) (header->seq + header->len + 1);
 	*seq = start;
 	*++seq = end;
-	*++seq = 0;
-	header->len += (2 * sizeof(int) + 1);
+	/* here */
+	/* uncomment */
+	//header->len += (2 * sizeof(int));
 	
 }

runkma.c

@@ -215,6 +215,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 		index_in_no = fileno(index_in);
 		read(index_in_no, &kmersize, sizeof(int));
 	}
+	
 	templatefilename[file_len] = 0;
 	
 	/* allocate stuff */
@@ -1311,6 +1312,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 	AlnPoints *points;
 	NWmat *NWmatrices;
 	Assemble_thread *threads, *thread;
+	HashMap_index *template_index;
 	
 	/* get lengths and names */
 	file_len = strlen(templatefilename);
@@ -2095,6 +2097,16 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 			matrix->size = template_lengths[i];
 		}
 	}
+	
+	template_index = smalloc(sizeof(HashMap_index));
+	template_index->size = 0;
+	if(alignLoadPtr != alignLoad_fly_shm) {
+		hashMap_index_initialize(template_index, matrix->size, kmersize);
+	} else {
+		template_index->seq = 0;
+		template_index->index = 0;
+	}
+	
 	if(assembly_KMA_Ptr == &assemble_KMA_threaded) {
 		matrix->size <<= 1;
 	} else {
@@ -2153,7 +2165,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 		thread->points = seedPoint_init(delta, rewards);
 		thread->points->len = 0;
 		thread->spin = (sparse < 0) ? 10 : 100;
-		 
+		thread->template_index = template_index;
 		thread->next = threads;
 		threads = thread;
 		
@@ -2211,7 +2223,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 	thread->header = header;
 	thread->points = points;
 	thread->points->len = 0;
-	thread->template_index = 0;
+	thread->template_index = template_index;
 	thread->next = 0;
 	thread->spin = (sparse < 0) ? 10 : 100;
 	
@@ -2228,7 +2240,6 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 		fprintf(stderr, "# Progress:\t%3d%%\r", 0);
 		fflush(stderr);
 	}
-	
 	for(template = 1; template < DB_size; ++template) {
 		if(w_scores[template] > 0) {
 			if(progress) {
@@ -2269,7 +2280,6 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 				//status |= assemblyPtr(aligned_assem, template, template_fragments, fileCount, frag_out, aligned, gap_align, qseq, header, matrix, points, NWmatrices);
 				thread->template = template;
 				assembly_KMA_Ptr(thread);
-				
 				/* Depth, ID and coverage */
 				if(aligned_assem->cover > 0) {
 					coverScore = aligned_assem->cover;
@@ -2333,7 +2343,6 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 			seq_seeker += ((template_lengths[template] >> 5) + 1);
 		}
 	}
-	
 	if(progress) {
 		fprintf(stderr, "\n");
 	}