[med-svn] [Git][med-team/kma][upstream] New upstream version 1.2.19
Steffen Möller
gitlab at salsa.debian.org
Fri Dec 27 19:14:32 GMT 2019
Steffen Möller pushed to branch upstream at Debian Med / kma
Commits:
a28edffe by Steffen Moeller at 2019-12-27T19:12:34Z
New upstream version 1.2.19
- - - - -
22 changed files:
- alnfrags.c
- alnfrags.h
- assembly.c
- assembly.h
- 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;
=====================================
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/commit/a28edffe0a2b4a8ddcdea4e7912fc595876a1d1a
--
View it on GitLab: https://salsa.debian.org/med-team/kma/commit/a28edffe0a2b4a8ddcdea4e7912fc595876a1d1a
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/1dfb6d31/attachment-0001.html>
More information about the debian-med-commit
mailing list