[med-svn] [Git][med-team/kma][master] 4 commits: New upstream version 1.2.19

Steffen Möller gitlab at salsa.debian.org
Fri Dec 27 19:14:17 GMT 2019



Steffen Möller pushed to branch master at Debian Med / kma


Commits:
a28edffe by Steffen Moeller at 2019-12-27T19:12:34Z
New upstream version 1.2.19
- - - - -
94665d08 by Steffen Moeller at 2019-12-27T19:12:34Z
routine-update: New upstream version

- - - - -
4c7d79d3 by Steffen Moeller at 2019-12-27T19:12:35Z
Update upstream source from tag 'upstream/1.2.19'

Update to upstream version '1.2.19'
with Debian dir 95c2ef97f92ac5f5045cea7b50130379b01793ae
- - - - -
8e34dd7c by Steffen Moeller at 2019-12-27T19:12:51Z
routine-update: Ready to upload to unstable

- - - - -


23 changed files:

- alnfrags.c
- alnfrags.h
- assembly.c
- assembly.h
- debian/changelog
- ef.c
- ef.h
- kma.c
- kmers.c
- kmers.h
- mt1.c
- mt1.h
- runinput.c
- runinput.h
- runkma.c
- runkma.h
- savekmers.c
- spltdb.c
- spltdb.h
- stdstat.c
- stdstat.h
- vcf.c
- version.h


Changes:

=====================================
alnfrags.c
=====================================
@@ -35,9 +35,9 @@
 #include "threader.h"
 #include "updatescores.h"
 
-int (*alnFragsPE)(HashMap_index**, int*, int*, int, double, CompDNA*, CompDNA*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, int, long*, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*) = alnFragsUnionPE;
+int (*alnFragsPE)(HashMap_index**, int*, int*, int, double, int, CompDNA*, CompDNA*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, int, long*, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*) = alnFragsUnionPE;
 
-int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int rc_flag, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, int q_len, int kmersize, Qseqs *header, int *bestTemplates, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *best_read_score, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, int rc_flag, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, int q_len, int kmersize, Qseqs *header, int *bestTemplates, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *best_read_score, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
 	
 	int t_i, template, read_score, bestHits, aln_len;
 	int start, end, W1;
@@ -100,9 +100,10 @@ int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *tem
 		//read_score += (((start != 0) + (end != template_lengths[abs(template)])) * W1);
 		
 		/* Get normed score */
-		if(aln_len > 0) {
+		if(minlen <= aln_len) {
 			score = 1.0 * read_score / aln_len;
 		} else {
+			read_score = 0;
 			score = 0;
 		}
 		
@@ -142,7 +143,7 @@ int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *tem
 	return 1;
 }
 
-int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
 	
 	int t_i, template, read_score;
 	int compScore, bestHits, bestHits_r, aln_len, start, end, rc, W1;
@@ -198,13 +199,12 @@ int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int
 		alnStat = KMA_score(templates_index[template], qseq, qseq_comp->seqlen, qseq_comp, mq, scoreT, points, NWmatrices);
 		
 		/* get read score */
-		if(0 < alnStat.score) {
-			aln_len = alnStat.len;
-			
+		aln_len = alnStat.len;
+		read_score = alnStat.score;
+		if(minlen <= aln_len && 0 < read_score) {
 			start = alnStat.pos;
 			end = alnStat.pos + alnStat.len - alnStat.gaps;
 			
-			read_score = alnStat.score;
 			if(start == 0 && end == template_lengths[abs(template)]) {
 				read_score += abs(W1);
 			}
@@ -227,13 +227,12 @@ int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int
 		
 		alnStat = KMA_score(templates_index[template], qseq_r, qseq_r_comp->seqlen, qseq_r_comp, mq, scoreT, points, NWmatrices);
 		/* get read score */
-		if(0 < alnStat.score) {
-			aln_len = alnStat.len;
-			
+		aln_len = alnStat.len;
+		read_score = alnStat.score;
+		if(minlen <= aln_len && 0 < read_score) {
 			start = alnStat.pos;
 			end = alnStat.pos + alnStat.len - alnStat.gaps;
 			
-			read_score = alnStat.score;
 			if(start == 0 && end == template_lengths[abs(template)]) {
 				read_score += abs(W1);
 			}
@@ -421,7 +420,7 @@ int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int
 	return 3;
 }
 
-int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
 	
 	int t_i, template, read_score;
 	int compScore, bestHits, bestHits_r, aln_len, start, end, rc, W1, PE;
@@ -478,13 +477,12 @@ int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, i
 		alnStat = KMA_score(templates_index[template], qseq, qseq_comp->seqlen, qseq_comp, mq, scoreT, points, NWmatrices);
 		
 		/* get read score */
-		if(0 < alnStat.score) {
-			aln_len = alnStat.len;
-			
+		aln_len = alnStat.len;
+		read_score = alnStat.score;
+		if(minlen <= aln_len && 0 < read_score) {
 			start = alnStat.pos;
 			end = alnStat.pos + alnStat.len - alnStat.gaps;
 			
-			read_score = alnStat.score;
 			if(start == 0 && end == template_lengths[abs(template)]) {
 				read_score += abs(W1);
 			}
@@ -507,13 +505,12 @@ int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, i
 		
 		alnStat = KMA_score(templates_index[template], qseq_r, qseq_r_comp->seqlen, qseq_r_comp, mq, scoreT, points, NWmatrices);
 		/* get read score */
-		if(0 < alnStat.score) {
-			aln_len = alnStat.len;
-			
+		aln_len = alnStat.len;
+		read_score = alnStat.score;
+		if(minlen <= aln_len && 0 < read_score) {
 			start = alnStat.pos;
 			end = alnStat.pos + alnStat.len - alnStat.gaps;
 			
-			read_score = alnStat.score;
 			if(start == 0 && end == template_lengths[abs(template)]) {
 				read_score += abs(W1);
 			}
@@ -698,7 +695,7 @@ int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, i
 	return 3;
 }
 
-int alnFragsForcePE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsForcePE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
 	
 	int t_i, template, read_score, bestHits, aln_len, W1;
 	int start, end, rc;
@@ -751,11 +748,11 @@ int alnFragsForcePE(HashMap_index **templates_index, int *matched_templates, int
 		
 		/* align qseq */
 		alnStat = KMA_score(templates_index[template], qseq, qseq_comp->seqlen, qseq_comp, mq, scoreT, points, NWmatrices);
-		if(0 < alnStat.score) {
+		if(0 < alnStat.score && minlen <= alnStat.len) {
 			alnStat_r = KMA_score(templates_index[template], qseq_r, qseq_r_comp->seqlen, qseq_r_comp, mq, scoreT, points, NWmatrices);
 			
 			/* get read score */
-			if(0 < alnStat_r.score) {
+			if(0 < alnStat_r.score && minlen <= alnStat_r.len) {
 				aln_len = alnStat.len + alnStat_r.len;
 				
 				/* Handle negative insertsizes caused by trimming,
@@ -834,9 +831,10 @@ void * alnFrags_threaded(void * arg) {
 	
 	static volatile int excludeIn[1] = {0}, excludeOut[1] = {0}, excludeDB[1] = {0};
 	Aln_thread *thread = arg;
-	int rc_flag, read_score, delta, index_in, seq_in, kmersize, flag, flag_r;
-	int mq, sam, unmapped, best_read_score, stats[2], *matched_templates, *bestTemplates;
-	int *best_start_pos, *best_end_pos, *template_lengths, *bestTemplates_r;
+	int rc_flag, read_score, delta, index_in, seq_in, kmersize, minlen;
+	int flag, flag_r, mq, sam, unmapped, best_read_score, stats[2];
+	int *matched_templates, *bestTemplates, *bestTemplates_r;
+	int *template_lengths, *best_start_pos, *best_end_pos;
 	long *index_indexes, *seq_indexes;
 	long unsigned *alignment_scores, *uniq_alignment_scores;
 	double scoreT;
@@ -872,6 +870,7 @@ void * alnFrags_threaded(void * arg) {
 	points = thread->points;
 	NWmatrices = thread->NWmatrices;
 	kmersize = thread->kmersize;
+	minlen = thread->minlen;
 	mq = thread->mq;
 	sam = thread->sam;
 	scoreT = thread->scoreT;
@@ -910,9 +909,9 @@ void * alnFrags_threaded(void * arg) {
 		
 		if(kmersize <= qseq->len) {
 			if(read_score && kmersize <= qseq_r->len) { // PE
-				unmapped = alnFragsPE(templates_index, matched_templates, template_lengths, mq, scoreT, qseq_comp, qseq_r_comp, qseq->seq, qseq_r->seq, header, header_r, kmersize, bestTemplates, bestTemplates_r, alignment_scores, uniq_alignment_scores, best_start_pos, best_end_pos, &flag, &flag_r, &best_read_score, &read_score, seq_in, index_in, seq_indexes, index_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
+				unmapped = alnFragsPE(templates_index, matched_templates, template_lengths, mq, scoreT, minlen, qseq_comp, qseq_r_comp, qseq->seq, qseq_r->seq, header, header_r, kmersize, bestTemplates, bestTemplates_r, alignment_scores, uniq_alignment_scores, best_start_pos, best_end_pos, &flag, &flag_r, &best_read_score, &read_score, seq_in, index_in, seq_indexes, index_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
 			} else { // SE
-				unmapped = alnFragsSE(templates_index, matched_templates, template_lengths, mq, scoreT, rc_flag, qseq_comp, qseq_r_comp, qseq->seq, qseq_r->seq, qseq->len, kmersize, header, bestTemplates, alignment_scores, uniq_alignment_scores, best_start_pos, best_end_pos, &flag, &best_read_score, seq_in, index_in, seq_indexes, index_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
+				unmapped = alnFragsSE(templates_index, matched_templates, template_lengths, mq, scoreT, minlen, rc_flag, qseq_comp, qseq_r_comp, qseq->seq, qseq_r->seq, qseq->len, kmersize, header, bestTemplates, alignment_scores, uniq_alignment_scores, best_start_pos, best_end_pos, &flag, &best_read_score, seq_in, index_in, seq_indexes, index_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
 			}
 		} else {
 			unmapped = 0;


=====================================
alnfrags.h
=====================================
@@ -45,6 +45,7 @@ struct aln_thread {
 	int index_in;
 	int seq_in;
 	int kmersize;
+	int minlen;
 	int mq;
 	int sam;
 	double scoreT;
@@ -62,9 +63,9 @@ struct aln_thread {
 #define ALNTHREAD 1;
 #endif
 
-extern int (*alnFragsPE)(HashMap_index**, int*, int*, int, double, CompDNA*, CompDNA*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, int, long*, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*);
-int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int rc_flag, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, int q_len, int kmersize, Qseqs *header, int *bestTemplates, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *best_read_score, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsForcePE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+extern int (*alnFragsPE)(HashMap_index**, int*, int*, int, double, int, CompDNA*, CompDNA*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, int, long*, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*);
+int alnFragsSE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, int rc_flag, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, int q_len, int kmersize, Qseqs *header, int *bestTemplates, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *best_read_score, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsUnionPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsPenaltyPE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsForcePE(HashMap_index **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, unsigned char *qseq, unsigned char *qseq_r, Qseqs *header, Qseqs *header_r, int kmersize, int *bestTemplates, int *bestTemplates_r, long unsigned *alignment_scores, long unsigned *uniq_alignment_scores, int *best_start_pos, int *best_end_pos, int *flag, int *flag_r, int *best_read_score, int *best_read_score_r, int seq_in, int index_in, long *seq_indexes, long *index_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
 void * alnFrags_threaded(void * arg);


=====================================
assembly.c
=====================================
@@ -279,7 +279,7 @@ void * assemble_KMA_threaded(void *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;
 	int nextTemplate, file_i, file_count, delta, thread_num, mq, status, bcd;
-	int stats[5], buffer[8];
+	int minlen, stats[5], buffer[8];
 	unsigned coverScore;
 	long unsigned depth, depthVar;
 	const char bases[6] = "ACGTN-";
@@ -312,6 +312,7 @@ void * assemble_KMA_threaded(void *arg) {
 	NWmatrices = thread->NWmatrices;
 	delta = qseq->size;
 	mq = thread->mq;
+	minlen = thread->minlen;
 	scoreT = thread->scoreT;
 	evalue = thread->evalue;
 	bcd = thread->bcd;
@@ -357,6 +358,8 @@ void * assemble_KMA_threaded(void *arg) {
 		
 		/* start threads */
 		aligned_assem->score = 0;
+		aligned_assem->fragmentCountAln = 0;
+		aligned_assem->readCountAln = 0;
 		mainTemplate = template;
 		thread_wait = thread_num;
 		unlock(excludeMatrix);
@@ -446,9 +449,10 @@ void * assemble_KMA_threaded(void *arg) {
 						
 						/* Get normed score */
 						read_score = alnStat.score;
-						if(0 < aln_len) {
+						if(minlen <= aln_len) {
 							score = 1.0 * read_score / aln_len;
 						} else {
+							read_score = 0;
 							score = 0;
 						}
 						
@@ -463,6 +467,10 @@ void * assemble_KMA_threaded(void *arg) {
 							//lock(excludeMatrix);
 							lockTime(excludeMatrix, 10)
 							aligned_assem->score += read_score;
+							if(!(stats[4] & 2) || (stats[4] & 64)) {
+								++aligned_assem->fragmentCountAln;
+							}
+							++aligned_assem->readCountAln;
 							
 							/* diff */
 							i = 0;
@@ -738,7 +746,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
 	Assemble_thread *thread = arg;
 	int i, j, t_len, aln_len, start, end, file_i, file_count, template, spin;
 	int pos, read_score, bestScore, depthUpdate, bestBaseScore, nextTemplate;
-	int sam, thread_num, mq, status, bcd, stats[5], buffer[8];
+	int sam, thread_num, mq, status, bcd, minlen, stats[5], buffer[8];
 	unsigned coverScore, delta;
 	long unsigned depth, depthVar;
 	const char bases[6] = "ACGTN-";
@@ -770,6 +778,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
 	NWmatrices = thread->NWmatrices;
 	delta = qseq->size;
 	mq = thread->mq;
+	minlen = thread->minlen;
 	scoreT = thread->scoreT;
 	evalue = thread->evalue;
 	bcd = thread->bcd;
@@ -834,6 +843,8 @@ void * assemble_KMA_dense_threaded(void *arg) {
 		
 		/* start threads */
 		aligned_assem->score = 0;
+		aligned_assem->fragmentCountAln = 0;
+		aligned_assem->readCountAln = 0;
 		mainTemplate = template;
 		thread_wait = thread_num;
 		unlock(excludeOut);
@@ -927,11 +938,11 @@ void * assemble_KMA_dense_threaded(void *arg) {
 						
 						/* Get normed score */
 						read_score = alnStat.score;
-						if(0 < aln_len) {
+						if(minlen <= aln_len) {
 							score = 1.0 * read_score / aln_len;
 						} else {
-							score = 0;
 							read_score = 0;
+							score = 0;
 						}
 						
 						if(0 < read_score && scoreT <= score) {
@@ -946,6 +957,10 @@ void * assemble_KMA_dense_threaded(void *arg) {
 							//lock(excludeMatrix);
 							lockTime(excludeMatrix, 10)
 							aligned_assem->score += read_score;
+							if(!(stats[4] & 2) || (stats[4] & 64)) {
+								++aligned_assem->fragmentCountAln;
+							}
+							++aligned_assem->readCountAln;
 							
 							/* diff */
 							for(i = 0, pos = start; i < aln_len; ++i) {
@@ -1150,6 +1165,8 @@ void * skip_assemble_KMA(void *arg) {
 	aligned_assem->q[0] = 0;
 	aligned_assem->len = t_len;
 	aligned_assem->aln_len = t_len;
+	aligned_assem->fragmentCountAln = 0;
+	aligned_assem->readCountAln = 0;
 	
 	/* load reads of this template */
 	file_i = 0;
@@ -1213,5 +1230,6 @@ void * skip_assemble_KMA(void *arg) {
 	
 	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
=====================================
@@ -42,6 +42,8 @@ struct assem {
 	unsigned len;
 	unsigned aln_len;
 	unsigned size;
+	unsigned fragmentCountAln;
+	unsigned readCountAln;
 };
 
 struct assembly {
@@ -62,6 +64,7 @@ struct assemble_thread {
 	int file_count;
 	int spin;
 	int mq;
+	int minlen;
 	int bcd;
 	int sam;
 	int thread_num;


=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+kma (1.2.19-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream version
+
+ -- Steffen Moeller <moeller at debian.org>  Fri, 27 Dec 2019 20:12:38 +0100
+
 kma (1.2.16-1) unstable; urgency=medium
 
   * Team upload.


=====================================
ef.c
=====================================
@@ -41,10 +41,10 @@ void initExtendedFeatures(FILE *out, char *templatefilename, unsigned totFrags,
 	strftime(Date, sizeof(Date), "%Y-%m-%d", tm);
 	fprintf(out, "## date\t%s\n", Date);
 	fprintf(out, "## command\t%s\n", cmd);
-	fprintf(out, "# refSequence\treadCount\tfragmentCount\tmapScoreSum\trefCoveredPositions\trefConsensusSum\tbpTotal\tdepthVariance\tnucHighDepthVariance\tdepthMax\tsnpSum\tinsertSum\tdeletionSum\n");
+	fprintf(out, "# refSequence\treadCount\tfragmentCount\tmapScoreSum\trefCoveredPositions\trefConsensusSum\tbpTotal\tdepthVariance\tnucHighDepthVariance\tdepthMax\tsnpSum\tinsertSum\tdeletionSum\treadCountAln\tfragmentCountAln\n");
 }
 
-void getExtendedFeatures(char *template_name, AssemInfo *matrix, long unsigned *template_seq, int t_len, Assem *aligned_assem, unsigned fragmentCount, unsigned readCount, FILE *outfile) {
+void getExtendedFeatures(char *template_name, AssemInfo *matrix, long unsigned *template_seq, int t_len, Assem *aligned_assem, unsigned fragmentCount, unsigned readCount, unsigned fragmentCountAln, unsigned readCountAln, FILE *outfile) {
 	
 	unsigned pos, depthUpdate, maxDepth, nucHighVarSum;
 	long unsigned snpSum, insertSum, deletionSum;
@@ -94,10 +94,10 @@ void getExtendedFeatures(char *template_name, AssemInfo *matrix, long unsigned *
 			}
 		} while((pos = assembly[pos].next) != 0);
 		
-		fprintf(outfile, "%s\t%u\t%u\t%lu\t%u\t%u\t%lu\t%f\t%u\t%u\t%lu\t%lu\t%lu\n", template_name, readCount, fragmentCount, aligned_assem->score, aligned_assem->aln_len, aligned_assem->cover, aligned_assem->depth, (double) var, nucHighVarSum, maxDepth, snpSum, insertSum, deletionSum);
+		fprintf(outfile, "%s\t%u\t%u\t%lu\t%u\t%u\t%lu\t%f\t%u\t%u\t%lu\t%lu\t%lu\t%u\t%u\n", template_name, readCount, fragmentCount, aligned_assem->score, aligned_assem->aln_len, aligned_assem->cover, aligned_assem->depth, (double) var, nucHighVarSum, maxDepth, snpSum, insertSum, deletionSum, readCountAln, fragmentCountAln);
 	} else if(aligned_assem) {
-		fprintf(outfile, "%s\t%u\t%u\t%u\t%u\t%u\t%lu\t%f\t%u\t%u\t%u\t%u\t%u\n", template_name, readCount, fragmentCount, 0, 0, 0, aligned_assem->depth, 0.0, 0, 0, 0, 0, 0);
+		fprintf(outfile, "%s\t%u\t%u\t%u\t%u\t%u\t%lu\t%f\t%u\t%u\t%u\t%u\t%u\t%u\t%u\n", template_name, readCount, fragmentCount, 0, 0, 0, aligned_assem->depth, 0.0, 0, 0, 0, 0, 0, readCountAln, fragmentCountAln);
 	} else {
-		fprintf(outfile, "%s\t%u\t%u\t%u\t%u\t%u\t%u\t%f\t%u\t%u\t%u\t%u\t%u\n", template_name, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0, 0, 0);
+		fprintf(outfile, "%s\t%u\t%u\t%u\t%u\t%u\t%u\t%f\t%u\t%u\t%u\t%u\t%u\t%u\t%u\n", template_name, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0, 0, 0, 0, 0);
 	}
 }


=====================================
ef.h
=====================================
@@ -21,4 +21,4 @@
 #include "assembly.h"
 
 void initExtendedFeatures(FILE *out, char *templatefilename, unsigned totFrags, char *cmd);
-void getExtendedFeatures(char *template_name, AssemInfo *matrix, long unsigned *template_seq, int t_len, Assem *aligned_assem, unsigned fragmentCount, unsigned readCount, FILE *outfile);
+void getExtendedFeatures(char *template_name, AssemInfo *matrix, long unsigned *template_seq, int t_len, Assem *aligned_assem, unsigned fragmentCount, unsigned readCount, unsigned fragmentCountAln, unsigned readCountAln, FILE *outfile);


=====================================
kma.c
=====================================
@@ -112,6 +112,7 @@ static void helpMessage(int exeStatus) {
 	fprintf(helpOut, "#\t-ipe\t\tInput paired end file name(s)\n");
 	fprintf(helpOut, "#\t-int\t\tInput interleaved file name(s)\n");
 	fprintf(helpOut, "#\t-k\t\tKmersize\t\t\t%s\n", "DB defined");
+	fprintf(helpOut, "#\t-ml\t\tMinimum alignment length\t%d\n", 16);
 	fprintf(helpOut, "#\t-p\t\tp-value\t\t\t\t0.05\n");
 	fprintf(helpOut, "#\t-ConClave\tConClave version\t\t1\n");
 	fprintf(helpOut, "#\t-mem_mode\tUse kmers to choose best\n#\t\t\ttemplate, and save memory\tFalse\n");
@@ -176,7 +177,7 @@ int kma_main(int argc, char *argv[]) {
 	static int minPhred, fiveClip, sparse_run, mem_mode, Mt1, ConClave, bcd;
 	static int fileCounter, fileCounter_PE, fileCounter_INT, targetNum, vcf;
 	static int extendedFeatures, spltDB, mq, thread_num, kmersize, one2one;
-	static int ref_fsa, print_matrix, print_all, sam, Ts, Tv, **d;
+	static int ref_fsa, print_matrix, print_all, sam, Ts, Tv, minlen, **d;
 	static unsigned nc, nf, shm, exhaustive;
 	static char *outputfilename, *templatefilename, **templatefilenames;
 	static char **inputfiles, **inputfiles_PE, **inputfiles_INT, ss;
@@ -222,6 +223,7 @@ int kma_main(int argc, char *argv[]) {
 		print_all = 0;
 		ref_fsa = 0;
 		kmersize = 0;
+		minlen = 16;
 		evalue = 0.05;
 		exhaustive = 0;
 		shm = 0;
@@ -470,6 +472,15 @@ int kma_main(int argc, char *argv[]) {
 						kmersize = 32;
 					}
 				}
+			} else if(strcmp(argv[args], "-ml") == 0) {
+				++args;
+				if(args < argc) {
+					minlen = strtoul(argv[args], &exeBasic, 10);
+					if(*exeBasic != 0) {
+						fprintf(stderr, "# Invalid kmersize parsed\n");
+						exit(4);
+					}
+				}
 			} else if(strcmp(argv[args], "-mp") == 0) {
 				++args;
 				if(args < argc) {
@@ -1090,17 +1101,17 @@ int kma_main(int argc, char *argv[]) {
 			
 			/* SE */
 			if(fileCounter > 0) {
-				totFrags += run_input(inputfiles, fileCounter, minPhred, fiveClip, kmersize, to2Bit, ioStream);
+				totFrags += run_input(inputfiles, fileCounter, minPhred, fiveClip, minlen, to2Bit, ioStream);
 			}
 			
 			/* PE */
 			if(fileCounter_PE > 0) {
-				totFrags += run_input_PE(inputfiles_PE, fileCounter_PE, minPhred, fiveClip, kmersize, to2Bit, ioStream);
+				totFrags += run_input_PE(inputfiles_PE, fileCounter_PE, minPhred, fiveClip, minlen, to2Bit, ioStream);
 			}
 			
 			/* INT */
 			if(fileCounter_INT > 0) {
-				totFrags += run_input_INT(inputfiles_INT, fileCounter_INT, minPhred, fiveClip, kmersize, to2Bit, ioStream);
+				totFrags += run_input_INT(inputfiles_INT, fileCounter_INT, minPhred, fiveClip, minlen, to2Bit, ioStream);
 			}
 			
 			if(Mt1) {
@@ -1119,13 +1130,13 @@ int kma_main(int argc, char *argv[]) {
 	} else if(Mt1) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);
 		strcpy(myTemplatefilename, templatefilename);
-		runKMA_Mt1(myTemplatefilename, outputfilename, strjoin(argv, argc), kmersize, rewards, ID_t, mq, scoreT, evalue, bcd, Mt1, ref_fsa, print_matrix, vcf, sam, nc, nf, thread_num);
+		runKMA_Mt1(myTemplatefilename, outputfilename, strjoin(argv, argc), kmersize, minlen, rewards, ID_t, mq, scoreT, evalue, bcd, Mt1, ref_fsa, print_matrix, vcf, sam, nc, nf, thread_num);
 		free(myTemplatefilename);
 		fprintf(stderr, "# Closing files\n");
 	} else if(step2) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);
 		strcpy(myTemplatefilename, templatefilename);
-		status = save_kmers_batch(myTemplatefilename, "-s1", shm, thread_num, exhaustive, rewards, ioStream, sam, scoreT, coverT);
+		status = save_kmers_batch(myTemplatefilename, "-s1", shm, thread_num, exhaustive, rewards, ioStream, sam, minlen, scoreT, coverT);
 		free(myTemplatefilename);
 	} else if(sparse_run) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);
@@ -1138,11 +1149,11 @@ int kma_main(int argc, char *argv[]) {
 		myTemplatefilename = smalloc(strlen(templatefilename) + 64);
 		strcpy(myTemplatefilename, templatefilename);
 		if(spltDB == 0 && targetNum != 1) {
-			status = runKMA_spltDB(templatefilenames, targetNum, outputfilename, argc, argv, ConClave, kmersize, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
+			status = runKMA_spltDB(templatefilenames, targetNum, outputfilename, argc, argv, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
 		} else if(mem_mode) {
-			status = runKMA_MEM(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
+			status = runKMA_MEM(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
 		} else {
-			status = runKMA(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
+			status = runKMA(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, sam, nc, nf, shm, thread_num);
 		}
 		free(myTemplatefilename);
 		fprintf(stderr, "# Closing files\n");


=====================================
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, double mrs, double coverT) {
+int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam, int minlen, double mrs, double coverT) {
 	
 	int i, file_len, shmid, deCon, *bestTemplates, *template_lengths;
 	FILE *inputfile, *templatefile;
@@ -154,12 +154,12 @@ 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);
+		save_kmers_HMM(templates, 0, &(int){thread_num}, template_lengths, 0, 0, 0, 0, 0, 0, minlen, 0, 0);
 	}
 	/* here */
 	/*
 	else if(kmerScan == &save_kmers_chain || kmerScan == &save_kmers_sparse_chain) {
-		kmerScan(0, 0, &(int){thread_num}, (int *)(&coverT), (int *)(&mrs), 0, 0, 0, 0, 0, 0, 0, 0);
+		kmerScan(0, 0, &(int){thread_num}, (int *)(&coverT), (int *)(&mrs), 0, 0, 0, 0, 0, minlen, 0, 0);
 	}
 	*/
 	


=====================================
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, double mrs, double coverT);
+int save_kmers_batch(char *templatefilename, char *exePrev, unsigned shm, int thread_num, const int exhaustive, Penalties *rewards, FILE *out, int sam, int minlen, double mrs, double coverT);


=====================================
mt1.c
=====================================
@@ -79,7 +79,7 @@ void printFsa_pairMt1(Qseqs *header, Qseqs *qseq, Qseqs *header_r, Qseqs *qseq_r
 	
 }
 
-void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int kmersize, Penalties *rewards, double ID_t, int mq, double scoreT, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int sam, int nc, int nf, int thread_num) {
+void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int kmersize, int minlen, Penalties *rewards, double ID_t, int mq, double scoreT, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int sam, int nc, int nf, int thread_num) {
 	
 	int i, j, aln_len, t_len, coverScore, file_len, DB_size, delta;
 	int *template_lengths;
@@ -271,6 +271,7 @@ void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int
 		thread->num = i;
 		thread->thread_num = thread_num;
 		thread->mq = mq;
+		thread->minlen = minlen;
 		thread->scoreT = scoreT;
 		thread->evalue = evalue;
 		thread->bcd = bcd;
@@ -330,6 +331,7 @@ void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int
 	thread->num = 0;
 	thread->thread_num = thread_num;
 	thread->mq = mq;
+	thread->minlen = minlen;
 	thread->scoreT = scoreT;
 	thread->evalue = evalue;
 	thread->bcd = bcd;


=====================================
mt1.h
=====================================
@@ -24,4 +24,4 @@
 
 void printFsaMt1(Qseqs *header, Qseqs *qseq, CompDNA *compressor, FILE *out);
 void printFsa_pairMt1(Qseqs *header, Qseqs *qseq, Qseqs *header_r, Qseqs *qseq_r, CompDNA *compressor, FILE *out);
-void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int kmersize, Penalties *rewards, double ID_t, int mq, double scoreT, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int sam, int nc, int nf, int thread_num);
+void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int kmersize, int minlen, Penalties *rewards, double ID_t, int mq, double scoreT, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int sam, int nc, int nf, int thread_num);


=====================================
runinput.c
=====================================
@@ -28,7 +28,7 @@
 void (*printFsa_ptr)(Qseqs*, Qseqs*, CompDNA*, FILE*) = &printFsa;
 void (*printFsa_pair_ptr)(Qseqs*, Qseqs*, Qseqs*, Qseqs*, CompDNA*, FILE*) = &printFsa_pair;
 
-long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out) {
+long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out) {
 	
 	int fileCounter, phredCut, start, end;
 	unsigned FASTQ;
@@ -84,7 +84,7 @@ long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int five
 				*/
 				qseq->len = end - start;
 				/* print */
-				if(qseq->len > kmersize) {
+				if(qseq->len > minlen) {
 					/* dump seq */
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
@@ -106,7 +106,7 @@ long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int five
 					++start;
 				}
 				qseq->len = end - start;
-				if(qseq->len > kmersize) {
+				if(qseq->len > minlen) {
 					/* dump seq */
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
@@ -133,7 +133,7 @@ long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int five
 	return count;
 }
 
-long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out) {
+long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out) {
 	
 	int fileCounter, phredCut, start, start2, end;
 	unsigned FASTQ, FASTQ2;
@@ -230,19 +230,19 @@ long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int f
 				qseq2->len = end - start2;
 				
 				/* print */
-				if(qseq->len > kmersize && qseq2->len > kmersize) {
+				if(qseq->len > minlen && qseq2->len > minlen) {
 					qseq->seq += start;
 					qseq2->seq += start2;
 					printFsa_pair_ptr(header, qseq, header2, qseq2, compressor, out);
 					qseq->seq -= start;
 					qseq2->seq -= start2;
 					++count;
-				} else if(qseq->len > kmersize) {
+				} else if(qseq->len > minlen) {
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
 					qseq->seq -= start;
 					++count;
-				} else if(qseq2->len > kmersize) {
+				} else if(qseq2->len > minlen) {
 					qseq2->seq += start2;
 					printFsa_ptr(header2, qseq2, compressor, out);
 					qseq2->seq -= start2;
@@ -276,19 +276,19 @@ long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int f
 				qseq2->len = end - start2;
 				
 				/* print */
-				if(qseq->len > kmersize && qseq2->len > kmersize) {
+				if(qseq->len > minlen && qseq2->len > minlen) {
 					qseq->seq += start;
 					qseq2->seq += start2;
 					printFsa_pair_ptr(header, qseq, header2, qseq2, compressor, out);
 					qseq->seq -= start;
 					qseq2->seq -= start2;
 					++count;
-				} else if(qseq->len > kmersize) {
+				} else if(qseq->len > minlen) {
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
 					qseq->seq -= start;
 					++count;
-				} else if(qseq2->len > kmersize) {
+				} else if(qseq2->len > minlen) {
 					qseq2->seq += start2;
 					printFsa_ptr(header2, qseq2, compressor, out);
 					qseq2->seq -= start2;
@@ -326,7 +326,7 @@ long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int f
 	return count;
 }
 
-long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out) {
+long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out) {
 	
 	int fileCounter, phredCut, start, start2, end;
 	unsigned FASTQ;
@@ -410,19 +410,19 @@ long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int
 				qseq2->len = end - start2;
 				
 				/* print */
-				if(qseq->len > kmersize && qseq2->len > kmersize) {
+				if(qseq->len > minlen && qseq2->len > minlen) {
 					qseq->seq += start;
 					qseq2->seq += start2;
 					printFsa_pair_ptr(header, qseq, header2, qseq2, compressor, out);
 					qseq->seq -= start;
 					qseq2->seq -= start2;
 					++count;
-				} else if(qseq->len > kmersize) {
+				} else if(qseq->len > minlen) {
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
 					qseq->seq -= start;
 					++count;
-				} else if(qseq2->len > kmersize) {
+				} else if(qseq2->len > minlen) {
 					qseq2->seq += start2;
 					printFsa_ptr(header2, qseq2, compressor, out);
 					qseq2->seq -= start2;
@@ -457,19 +457,19 @@ long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int
 				qseq2->len = end - start2;
 				
 				/* print */
-				if(qseq->len > kmersize && qseq2->len > kmersize) {
+				if(qseq->len > minlen && qseq2->len > minlen) {
 					qseq->seq += start;
 					qseq2->seq += start2;
 					printFsa_pair_ptr(header, qseq, header2, qseq2, compressor, out);
 					qseq->seq -= start;
 					qseq2->seq -= start2;
 					++count;
-				} else if(qseq->len > kmersize) {
+				} else if(qseq->len > minlen) {
 					qseq->seq += start;
 					printFsa_ptr(header, qseq, compressor, out);
 					qseq->seq -= start;
 					++count;
-				} else if(qseq2->len > kmersize) {
+				} else if(qseq2->len > minlen) {
 					qseq2->seq += start2;
 					printFsa_ptr(header2, qseq2, compressor, out);
 					qseq2->seq -= start2;


=====================================
runinput.h
=====================================
@@ -23,9 +23,9 @@
 /* pointers determining how to deliver the input */
 extern void (*printFsa_ptr)(Qseqs*, Qseqs*, CompDNA*, FILE*);
 extern void (*printFsa_pair_ptr)(Qseqs*, Qseqs*, Qseqs*, Qseqs*, CompDNA*, FILE*);
-long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out);
-long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out);
-long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int fiveClip, int kmersize, char *trans, FILE *out);
+long unsigned run_input(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out);
+long unsigned run_input_PE(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out);
+long unsigned run_input_INT(char **inputfiles, int fileCount, int minPhred, int fiveClip, int minlen, char *trans, FILE *out);
 void bootFsa(Qseqs *header, Qseqs *qseq, CompDNA *compressor, FILE *out);
 void printFsa(Qseqs *header, Qseqs *qseq, CompDNA *compressor, FILE *out);
 void printFsa_pair(Qseqs *header, Qseqs *qseq, Qseqs *header_r, Qseqs *qseq_r, CompDNA *compressor, FILE *out);


=====================================
runkma.c
=====================================
@@ -134,7 +134,7 @@ char * nameLoad(Qseqs *name, FILE *infile) {
 	return (char *) name->seq;
 }
 
-int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
+int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
 	
 	int i, j, tmp_template, tmp_tmp_template, file_len, bestTemplate, tot;
 	int template, bestHits, t_len, start, end, aln_len, status, rand, sparse;
@@ -343,6 +343,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 		alnThread->points = seedPoint_init(1024, rewards);
 		alnThread->NWmatrices = NWmatrices;
 		alnThread->kmersize = kmersize;
+		alnThread->minlen = minlen;
 		alnThread->template_lengths = template_lengths;
 		alnThread->templates_index = templates_index;
 		alnThread->mq = mq;
@@ -416,6 +417,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 	alnThread->points = points;
 	alnThread->NWmatrices = NWmatrices;
 	alnThread->kmersize = kmersize;
+	alnThread->minlen = minlen;
 	alnThread->mq = mq;
 	alnThread->sam = (sam == 1) ? 1 : 0;
 	alnThread->scoreT = scoreT;
@@ -483,6 +485,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 		thread->thread_num = thread_num;
 		thread->template = -2;
 		thread->mq = mq;
+		thread->minlen = minlen;
 		thread->scoreT = scoreT;
 		thread->evalue = evalue;
 		thread->bcd = bcd;
@@ -1139,6 +1142,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 	thread->thread_num = thread_num;
 	thread->template = 0;
 	thread->mq = mq;
+	thread->minlen = minlen;
 	thread->scoreT = scoreT;
 	thread->evalue = evalue;
 	thread->bcd = bcd;
@@ -1216,7 +1220,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 						updateMatrix(matrix_out, thread->template_name, templates_index[template]->seq, matrix, t_len);
 					}
 					if(extendedFeatures) {
-						getExtendedFeatures(thread->template_name, matrix, templates_index[template]->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+						getExtendedFeatures(thread->template_name, matrix, templates_index[template]->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], aligned_assem->fragmentCountAln, aligned_assem->readCountAln, extendedFeatures_out);
 					}
 					if(vcf) {
 						updateVcf(thread->template_name, templates_index[template]->seq, evalue, t_len, matrix, vcf, vcf_out);
@@ -1237,7 +1241,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 						fprintf(res_out, "%-12s\t%8ld\t%8u\t%8d\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%4.1e\n",
 							thread->template_name, read_score, (unsigned) expected, t_len, 0.0, cover, 0.0, q_cover, (double) depth, (double) q_value, p_value);
 						if(extendedFeatures) {
-							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], 0, 0, extendedFeatures_out);
 						}
 					}
 				} else {
@@ -1289,7 +1293,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
 	return status;
 }
 
-int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
+int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
 	
 	/* runKMA_MEM is a memory saving version of runKMA,
 	   at the cost it chooses best templates based on kmers
@@ -2156,6 +2160,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 		thread->num = i;
 		thread->thread_num = thread_num;
 		thread->mq = mq;
+		thread->minlen = minlen;
 		thread->scoreT = scoreT;
 		thread->evalue = evalue;
 		thread->bcd = bcd;
@@ -2215,6 +2220,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 	thread->num = 0;
 	thread->thread_num = thread_num;
 	thread->mq = mq;
+	thread->minlen = minlen;
 	thread->scoreT = scoreT;
 	thread->evalue = evalue;
 	thread->bcd = bcd;
@@ -2319,7 +2325,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 						updateMatrix(matrix_out, thread->template_name, thread->template_index->seq, matrix, t_len);
 					}
 					if(extendedFeatures) {
-						getExtendedFeatures(thread->template_name, matrix, thread->template_index->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+						getExtendedFeatures(thread->template_name, matrix, thread->template_index->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], aligned_assem->fragmentCountAln, aligned_assem->readCountAln, extendedFeatures_out);
 					}
 					if(vcf) {
 						updateVcf(thread->template_name, thread->template_index->seq, evalue, t_len, matrix, vcf, vcf_out);
@@ -2353,7 +2359,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
 						fprintf(res_out, "%-12s\t%8ld\t%8u\t%8d\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%4.1e\n",
 							thread->template_name, read_score, (unsigned) expected, t_len, 0.0, cover, 0.0, q_cover, (double) depth, (double) q_value, p_value);
 						if(extendedFeatures) {
-							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], 0, 0, extendedFeatures_out);
 						}
 					}
 				} else {


=====================================
runkma.h
=====================================
@@ -26,5 +26,5 @@
 unsigned char * ustrdup(unsigned char *src, size_t n);
 int load_DBs_KMA(char *templatefilename, long unsigned **alignment_scores, long unsigned **uniq_alignment_scores, int **template_lengths, unsigned shm);
 char * nameLoad(Qseqs *name, FILE *infile);
-int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);
-int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);
+int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);
+int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);


=====================================
savekmers.c
=====================================
@@ -3663,13 +3663,17 @@ int save_kmers_HMM(const HashMapKMA *templates, const Penalties *rewards, int *b
 			}
 			
 			/* set minLen */
-			minLen = template_lengths[1] / 2;
-			for(i = 1; i < templates->DB_size; ++i) {
-				if(minLen > (template_lengths[i] / 2)) {
-					minLen = template_lengths[i] / 2;
+			/*
+			minLen = *template_lengths >> 1;
+			i = templates->DB_size;
+			while(--i) {
+				if(minLen > ++*template_lengths >> 1)) {
+					minLen = *template_lengths >> 1;
 				}
 			}
-			minLen = MIN(minLen, templates->kmersize);
+			minLen = exhaustive < minLen ? minLen : exhaustive;
+			*/
+			minLen = exhaustive;
 		}
 		return 1;
 	} else if((DB_size = templates->DB_size) < USHRT_MAX) {
@@ -4771,7 +4775,7 @@ void ankerAndClean_MEM(int *regionTemplates, int *Score, int *Score_r, char *inc
 int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int *bestTemplates, int *bestTemplates_r, int *Score, int *Score_r, CompDNA *qseq, CompDNA *qseq_r, const Qseqs *header, int *extendScore, const int exhaustive, volatile int *excludeOut, FILE *out) {
 	
 	/* return 0 for match, 1 otherwise */
-	static int *Sizes = 0;
+	static int *Sizes = 0, minlen = 0;
 	static double coverT = 0.5, mrs = 0.5;
 	static KmerAnker **tVF_scores, **tVR_scores;
 	static SeqmentTree **tSeqments;
@@ -4798,6 +4802,7 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 			/* set coverT */
 			coverT = *((double *)(bestTemplates_r));
 			mrs = *((double *)(Score));
+			minlen = exhaustive;
 			i = *bestTemplates;
 			Sizes = smalloc(i * sizeof(int));
 			tVF_scores = smalloc(2 * i * sizeof(KmerAnker *));
@@ -5322,7 +5327,8 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 				cStart = tmp_score->start;
 				cover = queSeqmentTree(chainSeqments->root, cStart, best_score->end);
 				/* verify chain */
-				if(cover < coverT * (best_score->end - cStart) && mrs * (best_score->end - cStart) <= best_score->score) {
+				len = best_score->end - cStart;
+				if(minlen <= len && cover < coverT * len && mrs * len <= best_score->score) {
 					/* get chain */
 					rc = 1;
 				} else {
@@ -5340,7 +5346,8 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 						cover = queSeqmentTree(chainSeqments->root, cStart, best_score->end);
 						
 						/* verify chain */
-						if(cover < coverT * (best_score->end - cStart) && mrs * (best_score->end - cStart) <= best_score->score) {
+						len = best_score->end - cStart;
+						if(minlen <= len && cover < coverT * len && mrs * len <= best_score->score) {
 							/* get chain */
 							rc = 1;
 						} else {
@@ -5361,7 +5368,8 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 				cStart_r = tmp_score->start;
 				cover = queSeqmentTree(chainSeqments->root, cStart_r, best_score_r->end);
 				/* verify chain */
-				if(cover < coverT * (best_score_r->end - cStart_r) && mrs * (best_score_r->end - cStart_r) <= best_score_r->score) {
+				len = best_score_r->end - cStart_r;
+				if(minlen <= len && cover < coverT * len && mrs * len <= best_score_r->score) {
 					/* get chain */
 					rc |= 2;
 				} else {
@@ -5379,7 +5387,8 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 						cover = queSeqmentTree(chainSeqments->root, cStart_r, best_score_r->end);
 						
 						/* verify chain */
-						if(cover < coverT * (best_score_r->end - cStart_r) && mrs * (best_score_r->end - cStart_r) <= best_score_r->score) {
+						len = best_score_r->end - cStart_r;
+						if(minlen <= best_score_r->end - cStart_r && cover < coverT * len && mrs * len <= best_score_r->score) {
 							/* get chain */
 							rc |= 2;
 						} else {
@@ -5438,7 +5447,7 @@ int save_kmers_chain(const HashMapKMA *templates, const Penalties *rewards, int
 int save_kmers_sparse_chain(const HashMapKMA *templates, const Penalties *rewards, int *bestTemplates, int *bestTemplates_r, int *Score, int *Score_r, CompDNA *qseq, CompDNA *qseq_r, const Qseqs *header, int *extendScore, const int exhaustive, volatile int *excludeOut, FILE *out) {
 	
 	/* return 0 for match, 1 otherwise */
-	static int *Sizes = 0;
+	static int *Sizes = 0, minlen = 0;
 	static double coverT = 0.5, mrs = 0.5;
 	static KmerAnker **tVF_scores;
 	static SeqmentTree **tSeqments;
@@ -5464,6 +5473,7 @@ int save_kmers_sparse_chain(const HashMapKMA *templates, const Penalties *reward
 			/* set coverT */
 			coverT = *((double *)(bestTemplates_r));
 			mrs = *((double *)(Score));
+			minlen = exhaustive;
 			i = *bestTemplates;
 			Sizes = smalloc(i * sizeof(int));
 			tVF_scores = smalloc(i * sizeof(KmerAnker *));
@@ -5865,7 +5875,7 @@ int save_kmers_sparse_chain(const HashMapKMA *templates, const Penalties *reward
 					len = best_score->end - start;
 					
 					/* verify chain */
-					if(coverT * len <= cover || best_score->score < mrs * len) {
+					if(len < minlen || coverT * len <= cover || best_score->score < mrs * len) {
 						/* silence anker */
 						best_score->score = 0;
 					}


=====================================
spltdb.c
=====================================
@@ -394,7 +394,7 @@ unsigned get_ankers_spltDB(int *infoSize, int *out_Tem, CompDNA *qseq, Qseqs *he
 	return num;
 }
 
-int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename, int argc, char **argv, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
+int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename, int argc, char **argv, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num) {
 	
 	/* https://www.youtube.com/watch?v=LtXEMwSG5-8 */
 	
@@ -1396,6 +1396,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 		thread->num = i;
 		thread->thread_num = thread_num;
 		thread->mq = mq;
+		thread->minlen = minlen;
 		thread->scoreT = scoreT;
 		thread->evalue = evalue;
 		thread->bcd = bcd;
@@ -1455,6 +1456,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 	thread->num = 0;
 	thread->thread_num = thread_num;
 	thread->mq = mq;
+	thread->minlen = minlen;
 	thread->scoreT = scoreT;
 	thread->evalue = evalue;
 	thread->bcd = bcd;
@@ -1524,7 +1526,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 	}
 	templatefilename[file_len] = 0;
 	if(extendedFeatures == 2) {
-		getExtendedFeatures(templatefilename, 0, 0, 0, 0, 0, 0, extendedFeatures_out);
+		getExtendedFeatures(templatefilename, 0, 0, 0, 0, 0, 0, 0, 0, extendedFeatures_out);
 	}
 	if(assembly_KMA_Ptr == &skip_assemble_KMA) {
 		alignLoadPtr = &alignLoad_skip;
@@ -1573,7 +1575,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 			seq_seeker = 0;
 			index_seeker = 0;
 			if(extendedFeatures == 2) {
-				getExtendedFeatures(templatefilename, 0, 0, 0, 0, 0, 0, extendedFeatures_out);
+				getExtendedFeatures(templatefilename, 0, 0, 0, 0, 0, 0, 0, 0, extendedFeatures_out);
 			}
 		} else if(w_scores[template] > 0) {
 			
@@ -1642,7 +1644,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 						updateMatrix(matrix_out, thread->template_name, thread->template_index->seq, matrix, t_len);
 					}
 					if(extendedFeatures) {
-						getExtendedFeatures(thread->template_name, matrix, thread->template_index->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+						getExtendedFeatures(thread->template_name, matrix, thread->template_index->seq, t_len, aligned_assem, fragmentCounts[template], readCounts[template], aligned_assem->fragmentCountAln, aligned_assem->readCountAln, extendedFeatures_out);
 					}
 					if(vcf) {
 						updateVcf(thread->template_name, thread->template_index->seq, evalue, t_len, matrix, vcf, vcf_out);
@@ -1669,7 +1671,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
 						fprintf(res_out, "%-12s\t%8ld\t%8u\t%8d\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%8.2f\t%4.1e\n",
 							thread->template_name, read_score, (unsigned) expected, t_len, 0.0, 0.0, 0.0, 0.0, (double) depth, (double) q_value, p_value);
 						if(extendedFeatures) {
-							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], extendedFeatures_out);
+							getExtendedFeatures(thread->template_name, 0, 0, 0, aligned_assem, fragmentCounts[template], readCounts[template], 0, 0, extendedFeatures_out);
 						}
 					}
 				} else {


=====================================
spltdb.h
=====================================
@@ -35,4 +35,4 @@ struct spltDBbuff {
 int print_ankers_spltDB(int *out_Tem, CompDNA *qseq, int rc_flag, const Qseqs *header, const int flag, FILE *out);
 int print_ankers_Sparse_spltDB(int *out_Tem, CompDNA *qseq, int rc_flag, const Qseqs *header, const int flag, FILE *out);
 unsigned get_ankers_spltDB(int *infoSize, int *out_Tem, CompDNA *qseq, Qseqs *header, int *flag, FILE *inputfile);
-int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename, int argc, char **argv, int ConClave, int kmersize, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);
+int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename, int argc, char **argv, int ConClave, int kmersize, int minlen, Penalties *rewards, int extendedFeatures, double ID_t, int mq, double scoreT, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int sam, int nc, int nf, unsigned shm, int thread_num);


=====================================
stdstat.c
=====================================
@@ -142,6 +142,18 @@ double p_chisqr(long double q) {
 	return 1 - 1.772453850 * erf(sqrt(0.5 * q)) / tgamma(0.5);
 }
 
+double power(double x, unsigned n) {
+	
+	double y;
+	
+	if(n) {
+		y = power(x, n >> 1);
+		return (n & 1) ? y * y * x : y * y;
+	}
+	
+	return 1.0;
+}
+
 double binP(int n, int k, double p) {
 	
 	int i, j, nk;
@@ -153,9 +165,11 @@ double binP(int n, int k, double p) {
 	
 	q = 1 - p;
 	if(k == 0) {
-		return pow(q, n);
+		P = power(q, n);
+		return P != 0.0 ? P : 1.0e-308;
 	} else if(n == k) {
-		return pow(p, n);
+		P = power(p, n);
+		return P != 0.0 ? P : 1.0e-308;
 	} else if(p == 0 || q == 0) {
 		return 0.0;
 	}
@@ -175,9 +189,9 @@ double binP(int n, int k, double p) {
 	}
 	
 	if(nk < k) {
-		P *= pow(p, k - nk);
+		P *= power(p, k - nk);
 	} else if(k < nk) {
-		P *= pow(q, nk - k);
+		P *= power(q, nk - k);
 	}
 	
 	return P != 0.0 ? P : 1.0e-308;


=====================================
stdstat.h
=====================================
@@ -25,5 +25,6 @@ int cmp_or(int t, int q);
 int cmp_and(int t, int q);
 double fastp(long double q);
 double p_chisqr(long double q);
+double power(double x, unsigned n);
 double binP(int n, int k, double p);
 unsigned minimum(unsigned *src, unsigned n);


=====================================
vcf.c
=====================================
@@ -91,6 +91,7 @@ void initialiseVcf(FileBuff *fileP, char *templateFilename) {
 	
 	fileP->next = (unsigned char *) update;
 	fileP->bytes = avail;
+	
 }
 
 void updateVcf(char *template_name, long unsigned *template_seq, double evalue, int t_len, AssemInfo *matrix, int filter, FileBuff *fileP) {
@@ -165,6 +166,7 @@ void updateVcf(char *template_name, long unsigned *template_seq, double evalue,
 		if(bestScore) {
 			/* determine base at current position */
 			bestNuc = baseCall(bestNuc, nuc, bestScore, depthUpdate, evalue, &assembly[pos]);
+			
 			/* discard unimportant changes */
 			if(nuc != toupper(bestNuc)) {
 				/* INFO */
@@ -173,6 +175,7 @@ void updateVcf(char *template_name, long unsigned *template_seq, double evalue,
 				AF = (double) AD / DP;
 				RAF = (double) bestScore / DP;
 				DEL = assembly[pos].counts[5];
+				
 				/* FORMAT */
 				Q = pow(depthUpdate - (bestScore << 1), 2) / depthUpdate;
 				P = p_chisqr(Q);
@@ -232,6 +235,7 @@ void updateVcf(char *template_name, long unsigned *template_seq, double evalue,
 					*update++ = bestNuc;
 					--avail;
 				}
+				
 				check = sprintf(update, "\t%d\t%s\tDP=%d;AD=%d;AF=%.2f;RAF=%.2f;DEL=%d;", QUAL, *FILTER_ptr, DP, AD, AF, RAF, DEL);
 				update += check; avail -= check;
 				check = sprintf(update, "AD6=%d,%d,%d,%d,%d,%d\t", assembly[pos].counts[0], assembly[pos].counts[1], assembly[pos].counts[2], assembly[pos].counts[3], assembly[pos].counts[4], assembly[pos].counts[5]);
@@ -253,6 +257,7 @@ void updateVcf(char *template_name, long unsigned *template_seq, double evalue,
 			update += check; avail -= check;
 			check = sprintf(update, "Q:P:FT\t%.2f:%4.1e:%s\n", 0.0, 1.0, FILTER);
 			update += check; avail -= check;
+		
 		}
 	} while((pos = assembly[pos].next) != 0);
 	


=====================================
version.h
=====================================
@@ -17,4 +17,4 @@
  * limitations under the License.
 */
 
-#define KMA_VERSION "1.2.16"
+#define KMA_VERSION "1.2.19"



View it on GitLab: https://salsa.debian.org/med-team/kma/compare/d6256546f2607e48f3b102a45fbb211d8a491620...8e34dd7cd271ba9e5268983e3d6fd67672666176

-- 
View it on GitLab: https://salsa.debian.org/med-team/kma/compare/d6256546f2607e48f3b102a45fbb211d8a491620...8e34dd7cd271ba9e5268983e3d6fd67672666176
You're receiving this email because of your account on salsa.debian.org.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20191227/182eed4e/attachment-0001.html>


More information about the debian-med-commit mailing list