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 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
+CFLAGS ?= -Wall -O3
+CFLAGS += -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
 PROGS = kma kma_index kma_shm kma_update
 
 .c .o:
@@ -8,16 +9,16 @@ PROGS = kma kma_index kma_shm kma_update
 all: $(PROGS)
 
 kma: main.c libkma.a
-	$(CC) $(CFLAGS) -o $@ main.c libkma.a -lm -lpthread -lz
+	$(CC) $(CFLAGS) -o $@ main.c libkma.a -lm -lpthread -lz $(LDFLAGS)
 
 kma_index: kma_index.c libkma.a
-	$(CC) $(CFLAGS) -o $@ kma_index.c libkma.a -lm -lz
+	$(CC) $(CFLAGS) -o $@ kma_index.c libkma.a -lm -lz $(LDFLAGS)
 
 kma_shm: kma_shm.c libkma.a
-	$(CC) $(CFLAGS) -o $@ kma_shm.c libkma.a
+	$(CC) $(CFLAGS) -o $@ kma_shm.c libkma.a $(LDFLAGS)
 
 kma_update: kma_update.c libkma.a
-	$(CC) $(CFLAGS) -o $@ kma_update.c libkma.a
+	$(CC) $(CFLAGS) -o $@ kma_update.c libkma.a $(LDFLAGS)
 
 libkma.a: $(LIBS)
 	$(AR) -csru $@ $(LIBS)
@@ -46,7 +47,6 @@ 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

README.md

@@ -10,7 +10,7 @@ cd kma && make
 ```
 
 # Introduction #
-KMA is mapping a method designed to map raw reads directly against redundant databases, in an 
+KMA is a mapping method designed to map raw reads directly against redundant databases, in an 
 ultra-fast manner using seed and extend. KMA is particulary good at aligning 
 high quality reads against highly redundant databases, where unique matches often does 
 not exist. It works for long low quality reads as well, such as those from Nanopore. 

align.c

@@ -237,9 +237,11 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 		q_s = 0;
 		q_e = i;
 		if((q_e << 1) < t_e || (q_e + 64) < t_e) { // big leading template gap, cut down
-			t_s = t_e - MIN(64, (q_e << 1));
+			//t_s = t_e - MIN(64, (q_e << 1));
+			t_s = t_e - (q_e + (q_e < 64 ? q_e : 64));
 		} else if((t_e << 1) < q_e || (t_e + 64) < q_e) { // big leading query gap, cut down
-			q_s = q_e - MIN(64, (t_e << 1));
+			//q_s = q_e - MIN(64, (t_e << 1));
+			q_s = q_e - (t_e + (t_e < 64 ? t_e : 64));
 		}
 		
 		/* align */
@@ -249,7 +251,6 @@ 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);
 			}
 			
@@ -366,9 +367,13 @@ AlnScore KMA(const HashMap_index *template_index, const unsigned char *qseq, int
 	q_e = q_len;
 	t_e = t_len;
 	if(((q_len - q_s) << 1) < (t_len - t_s) || (q_len - q_s + 64) < (t_len - t_s)) { // big trailing template gap, cut down
-		t_e = t_s + MIN(64, ((q_len - q_s) << 1));
+		//t_e = t_s + MIN(64, ((q_len - q_s) << 1));
+		t_e = q_len - q_s;
+		t_e = t_s + (t_e + (t_e < 64 ? t_e : 64));
 	} else if(((t_len - t_s) << 1) < (q_len - q_s) || (t_len - t_s + 64) < (q_len - q_s)) { // big leading query gap, cut down
-		q_e = q_s + MIN(64, ((t_len - t_s) << 1));
+		//q_e = q_s + MIN(64, ((t_len - t_s) << 1));
+		q_e = t_len - t_s;
+		q_e = q_s + (q_e + (q_e < 64 ? q_e : 64));
 	}
 	
 	/* align trailing gap */
@@ -582,9 +587,11 @@ AlnScore KMA_score(const HashMap_index *template_index, const unsigned char *qse
 		q_s = 0;
 		q_e = i;
 		if((q_e << 1) < t_e || (q_e + 64) < t_e) { // big leading template gap, cut down
-			t_s = t_e - MIN(64, (q_e << 1));
+			//t_s = t_e - MIN(64, (q_e << 1));
+			t_s = t_e - (q_e + (q_e < 64 ? q_e : 64));
 		} else if((t_e << 1) < q_e || (t_e + 64) < q_e) { // big leading query gap, cut down
-			q_s = q_e - MIN(64, (t_e << 1));
+			//q_s = q_e - MIN(64, (t_e << 1));
+			q_s = q_e - (t_e + (t_e < 64 ? t_e : 64));
 		}
 		
 		/* align */
@@ -675,9 +682,13 @@ AlnScore KMA_score(const HashMap_index *template_index, const unsigned char *qse
 	q_e = q_len;
 	t_e = t_len;
 	if((t_len - t_s) > (q_len - q_s + 64) || (t_len - t_s) > ((q_len - q_s) << 1)) { // big trailing template gap, cut down
-		t_e = t_s + MIN(64, ((q_len - q_s) << 1));
+		//t_e = t_s + MIN(64, ((q_len - q_s) << 1));
+		t_e = q_len - q_s;
+		t_e = t_s + (t_e + (t_e < 64 ? t_e : 64));
 	} else if ((q_len - q_s) > (t_len - t_s + 64) || (q_len - q_s) > ((t_len - t_s) << 1)) { // big leading query gap, cut down
-		q_e = q_s + MIN(64, ((t_len - t_s) << 1));
+		//q_e = q_s + MIN(64, ((t_len - t_s) << 1));
+		q_e = t_len - t_s;
+		q_e = q_s + (q_e + (q_e < 64 ? q_e : 64));
 	}
 	
 	/* align trailing gap */

assembly.c

@@ -553,7 +553,7 @@ void * assemble_KMA_threaded(void *arg) {
 								samwrite(qseq, header, 0, template_name, aligned, stats);
 							}
 							
-						} else if(sam) {
+						} else if(sam && !(sam & 2096)) {
 							stats[1] = read_score;
 							stats[2] = start;
 							stats[3] = end;
@@ -567,7 +567,7 @@ void * assemble_KMA_threaded(void *arg) {
 								samwrite(qseq, header, 0, template_name, 0, stats);
 							}
 						}
-					} else if(sam) {
+					} else if(sam && !(sam & 2096)) {
 						stats[4] |= 4;
 						stats[1] = stats[4];
 						header->seq[header->len - 1] = 0;
@@ -599,6 +599,7 @@ void * assemble_KMA_threaded(void *arg) {
 				++file_i;
 			}
 		}
+		
 		lock(excludeIn);
 		--thread_wait;
 		unlock(excludeIn);
@@ -618,6 +619,7 @@ void * assemble_KMA_threaded(void *arg) {
 		
 		return NULL;
 	}
+	
 	/* diff */
 	/* pre on dense */
 	/* Pepare and make alignment on consensus */
@@ -697,7 +699,6 @@ void * assemble_KMA_threaded(void *arg) {
 		pos = assembly[pos].next;
 	}
 	
-
 	/* Trim alignment on consensus */
 	coverScore = 0;
 	bias = 0;
@@ -724,7 +725,6 @@ void * assemble_KMA_threaded(void *arg) {
 	aligned_assem->depthVar = depthVar;
 	aligned_assem->len = asm_len;
 	aligned_assem->aln_len = aln_len;
-	
 	return NULL;
 }
 
@@ -971,8 +971,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
 								header->seq[header->len - 1] = 0;
 								samwrite(qseq, header, 0, template_name, aligned, stats);
 							}
-						} else if(sam) {
-							/* here */
+						} else if(sam && !(sam & 2096)) {
 							stats[1] = read_score;
 							stats[2] = start;
 							stats[3] = end;
@@ -987,7 +986,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
 								samwrite(qseq, header, 0, template_name, 0, stats);
 							}
 						}
-					} else if(sam) {
+					} else if(sam && !(sam & 2096)) {
 						stats[4] |= 4;
 						stats[1] = stats[4];
 						header->seq[header->len - 1] = 0;
@@ -1115,10 +1114,29 @@ void * assemble_KMA_dense_threaded(void *arg) {
 	return NULL;
 }
 
-void skip_assemble_KMA(int template, int sam, int t_len, char *template_name, int file_count, FILE **files, Assem *aligned_assem, Qseqs *qseq, Qseqs *header) {
+void * skip_assemble_KMA(void *arg) {
+	
+	Assemble_thread *thread = arg;
+	int template, t_len, sam, nextTemplate, file_count, file_i, status;
+	int stats[5], buffer[8];
+	char *template_name;
+	FILE *file, **files;
+	Assem *aligned_assem;
+	Qseqs *qseq, *header;
 	
-	int nextTemplate, file_i, status, stats[5], buffer[8];
-	FILE *file;
+	if((template = thread->template) < 0) {
+		return NULL;
+	}
+	
+	/* get input */
+	sam = thread->sam;
+	t_len = thread->template_index->len;
+	template_name = thread->template_name;
+	file_count = thread->file_count;
+	files = thread->files;
+	aligned_assem = thread->aligned_assem;
+	qseq = thread->qseq;
+	header = thread->header;
 	
 	aligned_assem->cover = 0;
 	aligned_assem->depth = 0;
@@ -1188,4 +1206,8 @@ void skip_assemble_KMA(int template, int sam, int t_len, char *template_name, in
 			++file_i;
 		}
 	}
+	
+	aligned_assem->cover = 0; 
+	aligned_assem->aln_len = (1 - exp((-1.0) * aligned_assem->depth / t_len)) * t_len; // expected coverage from depth
+	return NULL;
 }

assembly.h

@@ -96,4 +96,4 @@ unsigned char nanoCaller(unsigned char bestNuc, unsigned char tNuc, int bestScor
 unsigned char refNanoCaller(unsigned char bestNuc, unsigned char tNuc, int bestScore, int depthUpdate, double evalue, Assembly *calls);
 void * assemble_KMA_threaded(void *arg);
 void * assemble_KMA_dense_threaded(void *arg);
-void skip_assemble_KMA(int template, int sam, int t_len, char *template_name, int file_count, FILE **files, Assem *aligned_assem, Qseqs *qseq, Qseqs *header);
+void * skip_assemble_KMA(void *arg);

compdna.c

@@ -106,12 +106,12 @@ int compDNAref(CompDNA *compressor, unsigned char *qseq, int seqlen) {
 		++bias;
 	}
 	seqlen -= bias;
+	
 	/* trim trailing N's */
-	--seqlen;
-	while(seq[seqlen] == 4) {
-		--seqlen;
-	}
+	if(seqlen) {
+		while(seq[--seqlen] == 4);
 		++seqlen;
+	}
 	
 	compressor->seqlen = seqlen;
 	if(seqlen & 31) {

debian/changelog

+kma (1.2.15-1) unstable; urgency=medium
+
+  * Team upload.
+
+  * New upstream version
+  * Standards-Version: 4.4.1
+  * Hardening-patch adopted upstream
+
+ -- Steffen Moeller <moeller@debian.org>  Fri, 22 Nov 2019 09:29:54 +0100
+
 kma (1.2.12-1) unstable; urgency=medium
 
   * Team upload.

debian/control

@@ -5,7 +5,7 @@ Section: science
 Priority: optional
 Build-Depends: debhelper-compat (= 12),
                zlib1g-dev
-Standards-Version: 4.4.0
+Standards-Version: 4.4.1
 Vcs-Browser: https://salsa.debian.org/med-team/kma
 Vcs-Git: https://salsa.debian.org/med-team/kam.git
 Homepage: https://bitbucket.org/genomicepidemiology/kma

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
-===================================================================
---- kma.orig/Makefile
-+++ 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 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
- all: $(PROGS)
- 
- kma: main.c libkma.a
--	$(CC) $(CFLAGS) -o $@ main.c libkma.a -lm -lpthread -lz
-+	$(CC) $(CFLAGS) -o $@ main.c libkma.a -lm -lpthread -lz $(LDFLAGS)
- 
- kma_index: kma_index.c libkma.a
--	$(CC) $(CFLAGS) -o $@ kma_index.c libkma.a -lm -lz
-+	$(CC) $(CFLAGS) -o $@ kma_index.c libkma.a -lm -lz $(LDFLAGS)
- 
- kma_shm: kma_shm.c libkma.a
--	$(CC) $(CFLAGS) -o $@ kma_shm.c libkma.a
-+	$(CC) $(CFLAGS) -o $@ kma_shm.c libkma.a $(LDFLAGS)
- 
- kma_update: kma_update.c libkma.a
--	$(CC) $(CFLAGS) -o $@ kma_update.c libkma.a
-+	$(CC) $(CFLAGS) -o $@ kma_update.c libkma.a $(LDFLAGS)
- 
- libkma.a: $(LIBS)
- 	$(AR) -csru $@ $(LIBS)

debian/patches/series

-hardening.patch

filebuff.c

@@ -24,6 +24,8 @@
 #include "pherror.h"
 #include "qseqs.h"
 
+int (*buffFileBuff)(FileBuff *) = &BuffgzFileBuff;
+
 int BuffgzFileBuff(FileBuff *dest) {
 	
 	int status;
@@ -65,6 +67,7 @@ int BuffgzFileBuff(FileBuff *dest) {
 		dest->bytes = 0;
 		dest->next = dest->buffer;
 		fprintf(stderr, "Gzip error %d\n", status);
+		errno |= status;
 	}
 	
 	return dest->bytes;
@@ -139,10 +142,13 @@ void gzcloseFileBuff(FileBuff *dest) {
 	int status;
 	if((status = inflateEnd(dest->strm)) != Z_OK) {
 		fprintf(stderr, "Gzip error %d\n", status);
+		errno |= status;
 	}
 	if(dest->z_err != Z_STREAM_END) {
 		fprintf(stderr, "Unexpected end of file\n");
+		errno |= status;
 	}
+	fclose(dest->file);
 	dest->file = 0;
 	dest->strm->avail_out = 0;
 }

filebuff.h

@@ -38,6 +38,8 @@ struct fileBuff {
 #define ENABLE_ZLIB_GZIP 32
 #endif
 
+/* pointer to load buffer from a regular or gz file stream */
+extern int (*buffFileBuff)(FileBuff *);
 int BuffgzFileBuff(FileBuff *dest);
 void init_gzFile(FileBuff *inputfile);
 FileBuff * setFileBuff(int buffSize);

frags.c

@@ -16,6 +16,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
 */
+#define _XOPEN_SOURCE 600
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>

hashmapindex.c

@@ -327,6 +327,14 @@ HashMap_index * alignLoad_fly_build_mem(HashMap_index *dest, int seq_in, int ind
 	return hashMap_index_build(dest, seq_in, len, kmersize);
 }
 
+HashMap_index * alignLoad_skip(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index) {
+	
+	dest->len = len;
+	dest->kmerindex = kmersize;
+	
+	return dest;
+}
+
 HashMap_index * alignLoad_fly_shm(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index) {
 	
 	static HashMap_index *src = 0;
@@ -350,9 +358,6 @@ HashMap_index * alignLoad_fly_shm(HashMap_index *dest, int seq_in, int index_in,
 
 HashMap_index * alignLoad_shm_initial(char *templatefilename, int file_len, int seq_in, int index_in, int kmersize) {
 	
-	/* here */
-	/* join all of them as one */
-	
 	key_t key;
 	int shmid;
 	HashMap_index *dest;

hashmapindex.h

@@ -51,6 +51,7 @@ HashMap_index * alignLoad_fly(HashMap_index *dest, int seq_in, int index_in, int
 HashMap_index * alignLoad_fly_mem(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index);
 HashMap_index * alignLoad_fly_build(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index);
 HashMap_index * alignLoad_fly_build_mem(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index);
+HashMap_index * alignLoad_skip(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index);
 HashMap_index * alignLoad_fly_shm(HashMap_index *dest, int seq_in, int index_in, int len, int kmersize, long unsigned seq_index, long unsigned index_index);
 HashMap_index * alignLoad_shm_initial(char *templatefilename, int file_len, int seq_in, int index_in, int kmersize);
 void alignClean(HashMap_index *template_index);

hashtable.c

@@ -51,7 +51,7 @@ int intpos_bin(const unsigned *str1, const int str2) {
 	return -1;
 }
 
-HashTable * collect_Kmers(const HashMapKMA *templates, int *Scores, int *Scores_tot, HashMap_kmers *foundKmers, Hit *hits) {
+HashTable * collect_Kmers(const HashMapKMA *templates, unsigned *Scores, unsigned *Scores_tot, HashMap_kmers *foundKmers, Hit *hits) {
 	
 	int template, SU;
 	unsigned i, j, *value;
@@ -118,7 +118,7 @@ HashTable * collect_Kmers(const HashMapKMA *templates, int *Scores, int *Scores_
 	return kmerList;
 }
 
-HashTable ** collect_Kmers_deCon(const HashMapKMA *templates, int *Scores, int *Scores_tot, HashMap_kmers *foundKmers, Hit *hits, int contamination) {
+HashTable ** collect_Kmers_deCon(const HashMapKMA *templates, unsigned *Scores, unsigned *Scores_tot, HashMap_kmers *foundKmers, Hit *hits, int contamination) {
 	
 	int template, SU;
 	unsigned i, j, n, *value;
@@ -219,7 +219,7 @@ HashTable ** collect_Kmers_deCon(const HashMapKMA *templates, int *Scores, int *
 	return Returner;
 }
 
-HashTable * withDraw_Kmers(int *Scores, int *Scores_tot, HashTable *kmerList, int template, Hit *hits) {
+HashTable * withDraw_Kmers(unsigned *Scores, unsigned *Scores_tot, HashTable *kmerList, int template, Hit *hits) {
 	
 	unsigned i;
 	HashTable *node, *prev;
@@ -267,7 +267,7 @@ HashTable * withDraw_Kmers(int *Scores, int *Scores_tot, HashTable *kmerList, in
 	return kmerList;
 }
 
-Hit withDraw_Contamination(int *Scores, int *Scores_tot, HashTable *kmerList, HashTable *deConTable, int template, Hit hits) {
+Hit withDraw_Contamination(unsigned *Scores, unsigned *Scores_tot, HashTable *kmerList, HashTable *deConTable, int template, Hit hits) {
 	
 	unsigned i, belong;
 	HashTable *node, *prev;

hashtable.h

@@ -31,14 +31,14 @@ struct hashTable {
 };
 
 struct hit {
-	unsigned n;
-	unsigned tot;
+	long unsigned n;
+	long unsigned tot;
 };
 #define HASHTABLE 1
 #endif
 
 int intpos_bin(const unsigned *str1, const int str2);
-HashTable * collect_Kmers(const HashMapKMA *templates, int *Scores, int *Scores_tot, HashMap_kmers *foundKmers, Hit *hits);
-HashTable ** collect_Kmers_deCon(const HashMapKMA *templates, int *Scores, int *Scores_tot, HashMap_kmers *foundKmers, Hit *hits, int contamination);
-HashTable * withDraw_Kmers(int *Scores, int *Scores_tot, HashTable *kmerList, int template, Hit *hits);
-Hit withDraw_Contamination(int *Scores, int *Scores_tot, HashTable *kmerList, HashTable *deConTable, int template, Hit hits);
+HashTable * collect_Kmers(const HashMapKMA *templates, unsigned *Scores, unsigned *Scores_tot, HashMap_kmers *foundKmers, Hit *hits);
+HashTable ** collect_Kmers_deCon(const HashMapKMA *templates, unsigned *Scores, unsigned *Scores_tot, HashMap_kmers *foundKmers, Hit *hits, int contamination);
+HashTable * withDraw_Kmers(unsigned *Scores, unsigned *Scores_tot, HashTable *kmerList, int template, Hit *hits);
+Hit withDraw_Contamination(unsigned *Scores, unsigned *Scores_tot, HashTable *kmerList, HashTable *deConTable, int template, Hit hits);

kma.c

@@ -119,11 +119,12 @@ static void helpMessage(int exeStatus) {
 	fprintf(helpOut, "#\t-ex_mode\tSearh kmers exhaustively\tFalse\n");
 	fprintf(helpOut, "#\t-ef\t\tPrint additional features\tFalse\n");
 	fprintf(helpOut, "#\t-vcf\t\tMake vcf file, 2 to apply FT\tFalse/0\n");
-	fprintf(helpOut, "#\t-sam\t\tOutput sam to stdout, 4 to \n#\t\t\tonly outputmapped reads\t\tFalse/0\n");
+	fprintf(helpOut, "#\t-sam\t\tOutput sam to stdout, 4 to \n#\t\t\tonly output mapped reads, \n#\t\t\t2096 for aligned\t\tFalse/0\n");
 	fprintf(helpOut, "#\t-nc\t\tNo consensus file\t\tFalse\n");
 	fprintf(helpOut, "#\t-nf\t\tNo frag file\t\tFalse\n");
 	fprintf(helpOut, "#\t-deCon\t\tRemove contamination\t\tFalse\n");
 	fprintf(helpOut, "#\t-dense\t\tDo not allow insertions\n#\t\t\tin assembly\t\t\tFalse\n");
+	fprintf(helpOut, "#\t-sasm\t\tSkip assembly\t\t\tFalse\n");
 	fprintf(helpOut, "#\t-ref_fsa\tConsensus sequnce will\n#\t\t\thave \"n\" instead of gaps\tFalse\n");
 	fprintf(helpOut, "#\t-matrix\t\tPrint assembly matrix\t\tFalse\n");
 	fprintf(helpOut, "#\t-a\t\tPrint all best mappings\t\tFalse\n");
@@ -156,6 +157,11 @@ static void helpMessage(int exeStatus) {
 	fprintf(helpOut, "#\t-gapopen\tPenalty for gap opening\t\t-3\n");
 	fprintf(helpOut, "#\t-gapextend\tPenalty for gap extension\t-1\n");
 	fprintf(helpOut, "#\t-per\t\tReward for pairing reads\t7\n");
+	/* here */
+	//fprintf(helpOut, "#\t-localopen\t\tPenalty for openning a local alignment\t-6\n");
+	fprintf(helpOut, "#\t-Npenalty\tPenalty matching N\t\t-1\n");
+	fprintf(helpOut, "#\t-transition\tPenalty for transition\t\t-2\n");
+	fprintf(helpOut, "#\t-transversion\tPenalty for transversion\t-2\n");
 	fprintf(helpOut, "#\t-cge\t\tSet CGE penalties and rewards\tFalse\n");
 	fprintf(helpOut, "#\t-t\t\tNumber of threads\t\t1\n");
 	fprintf(helpOut, "#\t-v\t\tVersion\n");
@@ -245,25 +251,6 @@ int kma_main(int argc, char *argv[]) {
 		Mt1 = 0;
 		inputfiles = 0;
 		deConPrintPtr = printPtr;
-		/*
-		assembly_KMA_Ptr = &assemble_KMA_threaded;
-		cmp = &cmp_or;
-		kmerScan = &save_kmers_HMM;
-		get_kmers_for_pair_ptr = &get_kmers_for_pair;
-		save_kmers_pair = &save_kmers_unionPair;
-		alnFragsPE = &alnFragsUnionPE;
-		printPairPtr = &printPair;
-		printPtr = &print_ankers;
-		printFsa_pair_ptr = &printFsa_pair;
-		ankerPtr = &ankerAndClean;
-		alignLoadPtr = &alignLoad_fly;
-		destroyPtr = &alignClean;
-		printFsa_ptr = &printFsa;
-		significantBase = &significantNuc; //-bc
-		baseCall = &baseCaller;
-		chainSeedsPtr = &chainSeeds;
-		deConPrintPtr = printPtr;
-		*/
 		
 		/* PARSE COMMAND LINE OPTIONS */
 		args = 1;
@@ -511,6 +498,8 @@ int kma_main(int argc, char *argv[]) {
 				}
 			} else if(strcmp(argv[args], "-dense") == 0) {
 				assembly_KMA_Ptr = &assemble_KMA_dense_threaded;
+			} else if(strcmp(argv[args], "-sasm") == 0) {
+				assembly_KMA_Ptr = &skip_assemble_KMA;
 			} else if(strcmp(argv[args], "-matrix") == 0) {
 				print_matrix = 1;
 			} else if(strcmp(argv[args], "-a") == 0) {
@@ -525,7 +514,6 @@ int kma_main(int argc, char *argv[]) {
 			} else if(strcmp(argv[args], "-ck") == 0) {
 				get_kmers_for_pair_ptr = &get_kmers_for_pair_count;
 			} else if(strcmp(argv[args], "-proxi") == 0) {
-				/* here */
 				if(++args < argc) {
 					support = strtod(argv[args], &exeBasic);
 					if(*exeBasic != 0 || support < 0 || 1 < support) {
@@ -664,7 +652,6 @@ int kma_main(int argc, char *argv[]) {
 					}
 				}
 			} else if(strcmp(argv[args], "-localopen") == 0) {
-				/* here */
 				/* add to help */
 				++args;
 				if(args < argc) {
@@ -676,7 +663,6 @@ int kma_main(int argc, char *argv[]) {
 					}
 				}
 			} else if(strcmp(argv[args], "-Npenalty") == 0) {
-				/* here */
 				/* add to help */
 				++args;
 				if(args < argc) {
@@ -698,7 +684,6 @@ int kma_main(int argc, char *argv[]) {
 					}
 				}
 			} else if(strcmp(argv[args], "-transition") == 0) {
-				/* here */
 				/* add to help */
 				++args;
 				if(args < argc) {
@@ -710,7 +695,6 @@ int kma_main(int argc, char *argv[]) {
 					}
 				}
 			} else if(strcmp(argv[args], "-transversion") == 0) {
-				/* here */
 				/* add to help */
 				++args;
 				if(args < argc) {

kmeranker.c

@@ -25,8 +25,8 @@
 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 i, j, template, score, bestScore, Wl, W1, U, M, MM, gaps;
-	int start, end, pos, SU;
+	int i, j, template, score, tmpScore, bestScore, Wl, W1, U, M, MM, gaps;
+	int start, end, pos, SU, nextAnker;
 	unsigned *values;
 	short unsigned *values_s;
 	KmerAnker *node, *prev;
@@ -64,9 +64,11 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 	bestScore = src->score;
 	values_s = 0;
 	values = 0;
+	nextAnker = 1;
+	prev = src;
 	
 	/* get chaining scores */
-	for(node = src; node; node = node->next) {
+	for(node = src; nextAnker; --node) {
 		if(SU) {
 			values_s = (short unsigned *) node->values;
 			i = *values_s + 1;
@@ -88,12 +90,7 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 				
 				/* extend chain */
 				if(pos == 0) {
-					if(node->cStart) {
-						score = W1 + (node->cStart - 1) * U;
-						score = node->weight + (Wl < score ? score : Wl);
-					} else {
 					score = node->weight;
-					}
 				} else {
 					if(gaps < 0) {
 						if(gaps == -kmersize) {
@@ -108,24 +105,33 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 					} else if((MM * 2 + (gaps - 2) * M) < 0) {
 						score += node->weight + (MM * 2 + (gaps - 2) * M);
 					} else {
-						score += node->weight - log(gaps * M);
+						score += node->weight - gaps * M;
+					}
+					/* mark as used */
+					node->score = 0;
 				}
-					
 				/* verify extension */
-					if(score < 0 || score == bestScore) {
+				if(score < 0) {
 					include[template] = 0;
+				} else if(bestScore <= score) {
+					/* test if anker is finalising the chain */
+					if(node->start) {
+						tmpScore = W1 + (node->start - 1) * U;
+						tmpScore = score + (Wl < tmpScore ? tmpScore : Wl);
+					} else {
+						tmpScore = score;
+					}
+					if(tmpScore == bestScore) {
+						score = bestScore;
+						nextAnker = 0;
+						prev = node;
 					}
 				}
 				extendScore[template] = start;
 				Score[template] = score;
 			}
 		}
-		
-		/* mark as used */
-		node->score = 0;
-		prev = node;
 	}
-	src->score = bestScore;
 	
 	/* get best templates */
 	j = 0;
@@ -144,14 +150,14 @@ KmerAnker * getBestChainTemplates(KmerAnker *src, const Penalties *rewards, int
 	return prev;
 }
 
-KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize, double mrs) {
+KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize) {
 	
 	KmerAnker *node, *prev;
 	
 	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));
+	while(V_score->score < kmersize && (V_score = V_score->descend));
 	
 	if(!V_score) {
 		return 0;
@@ -159,7 +165,7 @@ KmerAnker * pruneAnkers(KmerAnker *V_score, int kmersize, double mrs) {
 	prev = V_score;
 	node = V_score->descend;
 	while(node) {
-		if(kmersize <= node->score && mrs * (node->end - node->cStart) <= node->score) {
+		if(kmersize <= node->score) {
 			prev->descend = node;
 			prev = node;
 		}
@@ -189,14 +195,10 @@ KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties) {
 			best = node;
 			*ties = 0;
 		} else if(best->score == node->score) {
-			if((node->end - node->start) < (best->end - best->cStart)) {
-				best = node;
-				*ties = 0;
-			} else if((node->end - node->start) == (best->end - best->cStart)) {
 			best = node;
 			++*ties;
 		}
-		} else if(node->score != 0) {
+		if(node->score != 0) {
 			prev->descend = node;
 			prev = node;
 		}
@@ -207,15 +209,15 @@ KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties) {
 	return best;
 }
 
-KmerAnker * getTieAnker(int stop, KmerAnker *src, int score, int len) {
+KmerAnker * getTieAnker(int stop, KmerAnker *src, int score) {
 	
 	if(!src || src->start == stop) {
 		return 0;
 	}
 	
 	/* search downwards */
-	while((--src)->start != stop) {
-		if(src->score == score && (src->end - src->cStart) == len) {
+	while(stop < (--src)->start) {
+		if(src->score == score) {
 			return src;
 		}
 	}

kmeranker.h

@@ -27,15 +27,13 @@ struct kmerAnker {
 	int weight;
 	unsigned start;
 	unsigned end;
-	unsigned cStart;
 	unsigned ties;
 	unsigned *values;
 	struct kmerAnker *descend; /* descending anker */
-	struct kmerAnker *next; /* next ankers in chain */
 };
 #endif
 
 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 * pruneAnkers(KmerAnker *V_score, int kmersize);
 KmerAnker * getBestAnker(KmerAnker **src, unsigned *ties);
-KmerAnker * getTieAnker(int stop, KmerAnker *src, int score, int len);
+KmerAnker * getTieAnker(int stop, KmerAnker *src, int score);