[med-svn] [Git][med-team/kma][upstream] New upstream version 1.3.10
Nilesh Patra
gitlab at salsa.debian.org
Sat Feb 6 07:24:58 GMT 2021
Nilesh Patra pushed to branch upstream at Debian Med / kma
Commits:
339b0789 by Nilesh Patra at 2021-02-06T12:48:20+05:30
New upstream version 1.3.10
- - - - -
17 changed files:
- KMAspecification.pdf
- align.c
- alnfrags.c
- alnfrags.h
- assembly.c
- assembly.h
- kma.c
- mt1.c
- mt1.h
- nw.c
- nw.h
- runkma.c
- runkma.h
- sam.c
- spltdb.c
- spltdb.h
- version.h
Changes:
=====================================
KMAspecification.pdf
=====================================
Binary files a/KMAspecification.pdf and b/KMAspecification.pdf differ
=====================================
align.c
=====================================
@@ -200,6 +200,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
} else {
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.pos = 0;
aligned->s[0] = 0;
@@ -214,6 +215,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
if(aligned->mapQ < mq || score < kmersize) {
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.pos = 0;
aligned->s[0] = 0;
@@ -225,6 +227,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
/* initialize */
Stat.len = 0;
Stat.score = 0;
+ Stat.match = 0;
Stat.gaps = 0;
value = points->tStart[start] - 1;
Stat.pos = value;
@@ -278,6 +281,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
Stat.pos -= (NWstat.len - NWstat.gaps);
Stat.score = NWstat.score;
Stat.len = NWstat.len;
+ Stat.match = NWstat.match;
Stat.gaps = NWstat.gaps;
} else {
aligned->start = q_s;
@@ -294,6 +298,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
memset(aligned->s + Stat.len, '|', end);
memcpy(aligned->q + Stat.len, qseq + q_s, end);
Stat.len += end;
+ Stat.match += end;
end = points->qEnd[start];
for(i = points->qStart[start]; i < end; ++i) {
nuc = qseq[i];
@@ -334,6 +339,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
/* gap is too big to give a positive score */
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
aligned->s[0] = 0;
aligned->len = 0;
@@ -354,6 +360,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
memcpy(aligned->q + Stat.len, Frag_align->q, NWstat.len);
Stat.score += NWstat.score;
Stat.len += NWstat.len;
+ Stat.match += NWstat.match;
Stat.gaps += NWstat.gaps;
}
} else {
@@ -410,6 +417,7 @@ AlnScore KMA(const HashMapCCI *template_index, const unsigned char *qseq, int q_
memcpy(aligned->q + Stat.len, Frag_align->q, NWstat.len);
Stat.score += NWstat.score;
Stat.len += NWstat.len;
+ Stat.match += NWstat.match;
Stat.gaps += NWstat.gaps;
} else {
Frag_align->end = 0;
@@ -560,6 +568,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
} else {
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.pos = 0;
points->len = 0;
@@ -574,6 +583,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
if(mapQ < mq || score < kmersize) {
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.pos = 0;
points->len = 0;
@@ -583,6 +593,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
/* initialize */
Stat.len = 0;
Stat.score = 0;
+ Stat.match = 0;
Stat.gaps = 0;
value = points->tStart[start] - 1;
Stat.pos = value;
@@ -615,6 +626,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
Stat.pos -= (NWstat.len - NWstat.gaps);
Stat.score = NWstat.score;
Stat.len = NWstat.len;
+ Stat.match = NWstat.match;
Stat.gaps = NWstat.gaps;
}
}
@@ -626,6 +638,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
q_s = points->qStart[start];
end = points->qEnd[start] - q_s;
Stat.len += end;
+ Stat.match += end;
end = points->qEnd[start];
for(i = points->qStart[start]; i < end; ++i) {
nuc = qseq[i];
@@ -665,6 +678,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
/* gap is too big to give a positive score */
Stat.score = 0;
Stat.len = 1;
+ Stat.match = 0;
Stat.gaps = 0;
points->len = 0;
return Stat;
@@ -678,6 +692,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
}
Stat.score += NWstat.score;
Stat.len += NWstat.len;
+ Stat.match += NWstat.match;
Stat.gaps += NWstat.gaps;
}
} else {
@@ -711,6 +726,7 @@ AlnScore KMA_score(const HashMapCCI *template_index, const unsigned char *qseq,
}
Stat.score += NWstat.score;
Stat.len += NWstat.len;
+ Stat.match += NWstat.match;
Stat.gaps += NWstat.gaps;
}
points->len = 0;
=====================================
alnfrags.c
=====================================
@@ -35,12 +35,12 @@
#include "threader.h"
#include "updatescores.h"
-int (*alnFragsPE)(HashMapCCI**, int*, int*, int, double, int, CompDNA*, CompDNA*, CompDNA*, CompDNA*, unsigned char*, unsigned char*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*) = alnFragsUnionPE;
+int (*alnFragsPE)(HashMapCCI**, int*, int*, int, double, double, int, CompDNA*, CompDNA*, CompDNA*, CompDNA*, unsigned char*, unsigned char*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*) = alnFragsUnionPE;
-int alnFragsSE(HashMapCCI **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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsSE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
int t_i, template, read_score, bestHits, aln_len, start, end, W1, arc, rc;
- int q_start, q_end, *qBoundPtr;
+ int q_start, q_end, t_len, *qBoundPtr;
double score, bestScore;
AlnScore alnStat;
@@ -118,14 +118,15 @@ int alnFragsSE(HashMapCCI **templates_index, int *matched_templates, int *templa
aln_len = alnStat.len;
start = alnStat.pos;
end = start + aln_len - alnStat.gaps;
- if(template_lengths[abs(template)] < end) {
+ t_len = template_lengths[abs(template)];
+ if(template_lengths[abs(template)] < end && ((t_len < q_len) ? mrc * t_len : mrc * q_len) <= alnStat.match) {
end -= template_lengths[abs(template)];
}
/* penalty for non complete mapping */
read_score = alnStat.score;
/* full gene award */
- if((start == 0) && (end == template_lengths[abs(template)])) {
+ if((start == 0) && (end == t_len)) {
read_score += abs(W1);
}
//read_score += (((start != 0) + (end != template_lengths[abs(template)])) * W1);
@@ -174,10 +175,10 @@ int alnFragsSE(HashMapCCI **templates_index, int *matched_templates, int *templa
return 1;
}
-int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
int t_i, template, read_score, W1;
- int compScore, bestHits, bestHits_r, aln_len, start, end, arc, rc;
+ int compScore, bestHits, bestHits_r, aln_len, start, end, t_len, arc, rc;
double score;
AlnScore alnStat;
@@ -256,10 +257,11 @@ int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *t
/* get read score */
aln_len = alnStat.len;
read_score = alnStat.score;
- if(minlen <= aln_len && 0 < read_score) {
+ t_len = template_lengths[abs(template)];
+ if(minlen <= aln_len && 0 < read_score && ((t_len < qseq_comp->seqlen) ? mrc * t_len : mrc * qseq_comp->seqlen) <= alnStat.match) {
start = alnStat.pos;
end = alnStat.pos + alnStat.len - alnStat.gaps;
- if(start == 0 && end == template_lengths[abs(template)]) {
+ if(start == 0 && end == t_len) {
read_score += abs(W1);
}
score = 1.0 * read_score / aln_len;
@@ -302,11 +304,11 @@ int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *t
/* get read score */
aln_len = alnStat.len;
read_score = alnStat.score;
- if(minlen <= aln_len && 0 < read_score) {
+ if(minlen <= aln_len && 0 < read_score && ((t_len < qseq_r_comp->seqlen) ? mrc * t_len : mrc * qseq_r_comp->seqlen) <= alnStat.match) {
start = alnStat.pos;
end = alnStat.pos + alnStat.len - alnStat.gaps;
- if(start == 0 && end == template_lengths[abs(template)]) {
+ if(start == 0 && end == t_len) {
read_score += abs(W1);
}
score = 1.0 * read_score / aln_len;
@@ -339,7 +341,7 @@ int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *t
}
read_score += bestTemplates[t_i];
- if(best_start_pos[t_i] == 0 && best_end_pos[t_i] == template_lengths[abs(template)]) {
+ if(best_start_pos[t_i] == 0 && best_end_pos[t_i] == t_len) {
read_score += abs(W1);
}
if(compScore < read_score) {
@@ -501,10 +503,10 @@ int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *t
return 3;
}
-int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
int t_i, template, read_score, W1, PE;
- int compScore, bestHits, bestHits_r, aln_len, start, end, arc, rc;
+ int compScore, bestHits, bestHits_r, aln_len, start, end, t_len, arc, rc;
double score;
AlnScore alnStat;
@@ -584,11 +586,12 @@ int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int
/* get read score */
aln_len = alnStat.len;
read_score = alnStat.score;
- if(minlen <= aln_len && 0 < read_score) {
+ t_len = template_lengths[abs(template)];
+ if(minlen <= aln_len && 0 < read_score && ((t_len < qseq_comp->seqlen) ? mrc * t_len : mrc * qseq_comp->seqlen) <= alnStat.match) {
start = alnStat.pos;
end = alnStat.pos + alnStat.len - alnStat.gaps;
- if(start == 0 && end == template_lengths[abs(template)]) {
+ if(start == 0 && end == t_len) {
read_score += abs(W1);
}
score = 1.0 * read_score / aln_len;
@@ -629,11 +632,11 @@ int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int
/* get read score */
aln_len = alnStat.len;
read_score = alnStat.score;
- if(minlen <= aln_len && 0 < read_score) {
+ if(minlen <= aln_len && 0 < read_score && ((t_len < qseq_r_comp->seqlen) ? mrc * t_len : mrc * qseq_r_comp->seqlen) <= alnStat.match) {
start = alnStat.pos;
end = alnStat.pos + alnStat.len - alnStat.gaps;
- if(start == 0 && end == template_lengths[abs(template)]) {
+ if(start == 0 && end == t_len) {
read_score += abs(W1);
}
score = 1.0 * read_score / aln_len;
@@ -825,9 +828,10 @@ int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int
return 3;
}
-int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB) {
+int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_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, start, end, arc, rc;
+ int t_len;
double score, bestScore;
AlnScore alnStat, alnStat_r;
@@ -902,7 +906,8 @@ int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *t
alnStat = KMA_score(templates_index[template], qseq, qseq_comp->seqlen, 0, qseq_comp->seqlen, qseq_comp, mq, scoreT, points, NWmatrices);
}
- if(0 < alnStat.score && minlen <= alnStat.len) {
+ t_len = template_lengths[abs(template)];
+ if(0 < alnStat.score && minlen <= alnStat.len && ((t_len < qseq_comp->seqlen) ? mrc * t_len : mrc * qseq_comp->seqlen) <= alnStat.match) {
if(arc) {
if(rc < 0) {
/* rc */
@@ -936,7 +941,7 @@ int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *t
}
read_score = alnStat.score + alnStat_r.score;
- if(start == 0 && end == template_lengths[abs(template)]) {
+ if(start == 0 && end == t_len) {
read_score += abs(W1);
}
score = 1.0 * read_score / aln_len;
@@ -948,7 +953,7 @@ int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *t
}
/* save best match(es) */
- if(read_score > kmersize && score >= scoreT) {
+ if(read_score > kmersize && score >= scoreT && ((t_len < qseq_r_comp->seqlen) ? mrc * t_len : mrc * qseq_r_comp->seqlen) <= alnStat.match) {
if(score > bestScore) { // save as best match
bestScore = score;
*best_read_score = read_score;
@@ -1016,7 +1021,7 @@ void * alnFrags_threaded(void * arg) {
long *seq_indexes;
long unsigned *alignment_scores, *uniq_alignment_scores;
unsigned char *qseq_fr, *qseq_rr;
- double scoreT;
+ double scoreT, mrc;
FILE *inputfile, *frag_out_raw;
FileBuff *frag_out_all;
CompDNA *qseq_comp, *qseq_r_comp, *qseq_fr_comp, *qseq_rr_comp;
@@ -1051,6 +1056,7 @@ void * alnFrags_threaded(void * arg) {
mq = thread->mq;
sam = thread->sam;
scoreT = thread->scoreT;
+ mrc = thread->mrc;
template_lengths = thread->template_lengths;
templates_index = thread->templates_index;
@@ -1103,9 +1109,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, minlen, qseq_comp, qseq_r_comp, qseq_fr_comp, qseq_rr_comp, qseq->seq, qseq_r->seq, qseq_fr, qseq_rr, 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, seq_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
+ unmapped = alnFragsPE(templates_index, matched_templates, template_lengths, mq, scoreT, mrc, minlen, qseq_comp, qseq_r_comp, qseq_fr_comp, qseq_rr_comp, qseq->seq, qseq_r->seq, qseq_fr, qseq_rr, 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, seq_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
} else { // SE
- 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, seq_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
+ unmapped = alnFragsSE(templates_index, matched_templates, template_lengths, mq, scoreT, mrc, 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, seq_indexes, frag_out_raw, points, NWmatrices, excludeOut, excludeDB);
}
} else {
unmapped = 0;
=====================================
alnfrags.h
=====================================
@@ -47,6 +47,7 @@ struct aln_thread {
int mq;
int sam;
double scoreT;
+ double mrc;
CompDNA *qseq_comp;
CompDNA *qseq_r_comp;
Qseqs *qseq;
@@ -60,9 +61,9 @@ struct aln_thread {
};
#define ALNTHREAD 1;
#endif
-extern int (*alnFragsPE)(HashMapCCI**, int*, int*, int, double, int, CompDNA*, CompDNA*, CompDNA*, CompDNA*, unsigned char*, unsigned char*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*);
-int alnFragsSE(HashMapCCI **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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
-int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+extern int (*alnFragsPE)(HashMapCCI**, int*, int*, int, double, double, int, CompDNA*, CompDNA*, CompDNA*, CompDNA*, unsigned char*, unsigned char*, unsigned char*, unsigned char*, Qseqs*, Qseqs*, int, int*, int*, long unsigned*, long unsigned*, int*, int*, int*, int*, int*, int*, int, long*, FILE*, AlnPoints*, NWmat*, volatile int*, volatile int*);
+int alnFragsSE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsUnionPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsPenaltyPE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
+int alnFragsForcePE(HashMapCCI **templates_index, int *matched_templates, int *template_lengths, int mq, double scoreT, double mrc, int minlen, CompDNA *qseq_comp, CompDNA *qseq_r_comp, CompDNA *qseq_fr_comp, CompDNA *qseq_rr_comp, unsigned char *qseq, unsigned char *qseq_r, unsigned char *qseq_fr, unsigned char *qseq_rr, 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, long *seq_indexes, FILE *frag_out_raw, AlnPoints *points, NWmat *NWmatrices, volatile int *excludeOut, volatile int *excludeDB);
void * alnFrags_threaded(void * arg);
=====================================
assembly.c
=====================================
@@ -283,7 +283,7 @@ void * assemble_KMA_threaded(void *arg) {
long unsigned depth, depthVar;
short unsigned *counts;
const char bases[6] = "ACGTN-";
- double score, scoreT, evalue;
+ double score, scoreT, mrc, evalue;
unsigned char bestNuc;
FILE **files, *file, *xml_out;
AlnScore alnStat;
@@ -314,6 +314,7 @@ void * assemble_KMA_threaded(void *arg) {
mq = thread->mq;
minlen = thread->minlen;
scoreT = thread->scoreT;
+ mrc = thread->mrc;
evalue = thread->evalue;
bcd = thread->bcd;
sam = thread->sam;
@@ -433,9 +434,9 @@ void * assemble_KMA_threaded(void *arg) {
start = alnStat.pos;
end = start + aln_len - alnStat.gaps;
- /* Get normed score */
+ /* Get normed score check read coverage */
read_score = alnStat.score;
- if(minlen <= aln_len) {
+ if(minlen <= aln_len && ((t_len < qseq->len) ? mrc * t_len : mrc * qseq->len) <= alnStat.match) {
score = 1.0 * read_score / aln_len;
} else {
read_score = 0;
@@ -781,7 +782,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
long unsigned depth, depthVar;
short unsigned *counts;
const char bases[6] = "ACGTN-";
- double score, scoreT, evalue;
+ double score, scoreT, mrc, evalue;
unsigned char bestNuc;
FILE **files, *file, *xml_out;
AlnScore alnStat;
@@ -812,6 +813,7 @@ void * assemble_KMA_dense_threaded(void *arg) {
mq = thread->mq;
minlen = thread->minlen;
scoreT = thread->scoreT;
+ mrc = thread->mrc;
evalue = thread->evalue;
bcd = thread->bcd;
sam = thread->sam;
@@ -986,9 +988,9 @@ void * assemble_KMA_dense_threaded(void *arg) {
start = alnStat.pos;
end = start + aln_len - alnStat.gaps;
- /* Get normed score */
+ /* Get normed score check read coverage */
read_score = alnStat.score;
- if(minlen <= aln_len) {
+ if(minlen <= aln_len && ((t_len < qseq->len) ? mrc * t_len : mrc * qseq->len) <= alnStat.match) {
score = 1.0 * read_score / aln_len;
} else {
read_score = 0;
@@ -1622,7 +1624,7 @@ void * assemble_KMA(void *arg) {
int thread_num, mq, bcd, start, end, q_start, q_end;
int stats[5], buffer[8], *qBoundPtr;
short unsigned *counts;
- double score, scoreT, evalue;
+ double score, scoreT, mrc, evalue;
long double var, nucHighVar;
char *s, *s_next;
unsigned char *t, *q, *t_next, *q_next;
@@ -1654,6 +1656,7 @@ void * assemble_KMA(void *arg) {
delta = qseq->size;
mq = thread->mq;
minlen = thread->minlen;
+ mrc = thread->mrc;
scoreT = thread->scoreT;
evalue = thread->evalue;
bcd = thread->bcd;
@@ -1857,9 +1860,9 @@ void * assemble_KMA(void *arg) {
start = alnStat.pos;
end = start + aln_len - alnStat.gaps;
- /* Get normed score */
+ /* Get normed score check read coverage */
read_score = alnStat.score;
- if(minlen <= aln_len) {
+ if(minlen <= aln_len && ((t_len < qseq->len) ? mrc * t_len : mrc * qseq->len) <= alnStat.match) {
score = 1.0 * read_score / aln_len;
} else {
read_score = 0;
=====================================
assembly.h
=====================================
@@ -79,6 +79,7 @@ struct assemble_thread {
int kmersize;
int thread_num;
double scoreT;
+ double mrc;
double evalue;
char *template_name;
FILE **files;
=====================================
kma.c
=====================================
@@ -158,7 +158,8 @@ static void helpMessage(int exeStatus) {
fprintf(helpOut, "#\t-and\t\tBoth mrs and p_value thresholds\n#\t\t\thas to reached to in order to\n#\t\t\treport a template hit.\t\tor\n");
fprintf(helpOut, "#\t-mq\t\tMinimum mapping quality\t\t0\n");
fprintf(helpOut, "#\t-mrs\t\tMinimum alignment score,\n#\t\t\tnormalized to alignment length\t0.50\n");
- fprintf(helpOut, "#\t-mct\t\tMax overlap between templates\t0.50\n");
+ fprintf(helpOut, "#\t-mrc\t\tMinimum read coverage\t\t0.10\n");
+ fprintf(helpOut, "#\t-mct\t\tMax overlap between templates\t0.10\n");
fprintf(helpOut, "#\t-reward\t\tScore for match\t\t\t1\n");
fprintf(helpOut, "#\t-penalty\tPenalty for mismatch\t\t-2\n");
fprintf(helpOut, "#\t-gapopen\tPenalty for gap opening\t\t-3\n");
@@ -206,7 +207,7 @@ int kma_main(int argc, char *argv[]) {
static unsigned xml, nc, nf, shm, exhaustive, verbose;
static char *outputfilename, *templatefilename, **templatefilenames;
static char **inputfiles, **inputfiles_PE, **inputfiles_INT, ss;
- static double ID_t, scoreT, coverT, evalue;
+ static double ID_t, scoreT, coverT, mrc, evalue;
static Penalties *rewards;
int i, j, args, exe_len, fileCount, size, escape, tmp, step1, step2;
unsigned totFrags;
@@ -258,7 +259,8 @@ int kma_main(int argc, char *argv[]) {
mq = 0;
bcd = 1;
scoreT = 0.5;
- coverT = 0.5;
+ coverT = 0.1;
+ mrc = 0.1;
ID_t = 1.0;
one2one = 0;
ss = 'q';
@@ -686,6 +688,15 @@ int kma_main(int argc, char *argv[]) {
exit(1);
}
}
+ } else if(strcmp(argv[args], "-mrc") == 0) {
+ ++args;
+ if(args < argc) {
+ mrc = strtod(argv[args], &exeBasic);
+ if(*exeBasic != 0) {
+ fprintf(stderr, "Invalid argument at \"-mrc\".\n");
+ exit(1);
+ }
+ }
} else if(strcmp(argv[args], "-mct") == 0) {
++args;
if(args < argc) {
@@ -1274,7 +1285,7 @@ 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, minlen, rewards, ID_t, mq, scoreT, evalue, bcd, Mt1, ref_fsa, print_matrix, vcf, xml, sam, nc, nf, thread_num);
+ runKMA_Mt1(myTemplatefilename, outputfilename, strjoin(argv, argc), kmersize, minlen, rewards, ID_t, mq, scoreT, mrc, evalue, bcd, Mt1, ref_fsa, print_matrix, vcf, xml, sam, nc, nf, thread_num);
free(myTemplatefilename);
fprintf(stderr, "# Closing files\n");
} else if(step2) {
@@ -1293,11 +1304,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, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
+ status |= runKMA_spltDB(templatefilenames, targetNum, outputfilename, argc, argv, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, mrc, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
} else if(mem_mode) {
- status |= runKMA_MEM(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
+ status |= runKMA_MEM(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, mrc, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
} else {
- status |= runKMA(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
+ status |= runKMA(myTemplatefilename, outputfilename, exeBasic, ConClave, kmersize, minlen, rewards, extendedFeatures, ID_t, mq, scoreT, mrc, evalue, bcd, ref_fsa, print_matrix, print_all, vcf, xml, sam, nc, nf, shm, thread_num, verbose);
}
free(myTemplatefilename);
fprintf(stderr, "# Closing files\n");
=====================================
mt1.c
=====================================
@@ -80,7 +80,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, 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 xml, 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 mrc, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int xml, 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;
@@ -290,6 +290,7 @@ void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
@@ -354,6 +355,7 @@ void runKMA_Mt1(char *templatefilename, char *outputfilename, char *exePrev, int
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
=====================================
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, 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 xml, 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 mrc, double evalue, int bcd, int Mt1, int ref_fsa, int print_matrix, int vcf, int xml, int sam, int nc, int nf, int thread_num);
=====================================
nw.c
=====================================
@@ -49,11 +49,13 @@ AlnScore NW(const long unsigned *template, const unsigned char *queryOrg, int k,
if(t_len == 0 || q_len == 0) {
if(t_len == q_len) {
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = 0;
aligned->s[0] = 0;
} else if(t_len == 0) {
Stat.len = q_len;
+ Stat.match = 0;
Stat.gaps = q_len;
Stat.score = W1 + (q_len - 1) * U;
memset(aligned->s, '_', q_len);
@@ -62,6 +64,7 @@ AlnScore NW(const long unsigned *template, const unsigned char *queryOrg, int k,
memcpy(aligned->q, query, q_len);
} else {
Stat.len = t_len;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = W1 + (t_len - 1) * U;
memset(aligned->s, '_', t_len);
@@ -248,6 +251,7 @@ AlnScore NW(const long unsigned *template, const unsigned char *queryOrg, int k,
E_ptr = E + (m * (q_len + 1));
nuc_pos = m + t_s;
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
while(E_ptr[n] != 0) {
if(nuc_pos == template_length) {
@@ -257,6 +261,7 @@ AlnScore NW(const long unsigned *template, const unsigned char *queryOrg, int k,
aligned->t[Stat.len] = getNuc(template, nuc_pos);
aligned->q[Stat.len] = query[n];
aligned->s[Stat.len] = (aligned->t[Stat.len] == aligned->q[Stat.len]) ? '|' : '_';
+ ++Stat.match;
++nuc_pos;
E_ptr += (q_len + 1);
++n;
@@ -324,11 +329,13 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
if(t_len == 0 || q_len == 0) {
if(t_len == q_len) {
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = 0;
aligned->s[0] = 0;
} else if(t_len == 0) {
Stat.len = q_len;
+ Stat.match = 0;
Stat.gaps = q_len;
Stat.score = W1 + (q_len - 1) * U;
memset(aligned->s, '_', q_len);
@@ -337,6 +344,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
memcpy(aligned->q, query, q_len);
} else {
Stat.len = t_len;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = W1 + (t_len - 1) * U;
memset(aligned->s, '_', t_len);
@@ -567,6 +575,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
E_ptr = E + (m * (bq_len + 1));
nuc_pos = m + t_s;
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
while(E_ptr[n] != 0) {
if(nuc_pos == template_length) {
@@ -576,6 +585,7 @@ AlnScore NW_band(const long unsigned *template, const unsigned char *queryOrg, i
aligned->t[Stat.len] = getNuc(template, nuc_pos);
aligned->q[Stat.len] = query[q_pos];
aligned->s[Stat.len] = (aligned->t[Stat.len] == aligned->q[Stat.len]) ? '|' : '_';
+ ++Stat.match;
++nuc_pos;
E_ptr += (bq_len + 1);
++q_pos;
@@ -643,14 +653,17 @@ AlnScore NW_score(const long unsigned *template, const unsigned char *queryOrg,
if(t_len == 0 || q_len == 0) {
if(t_len == q_len) {
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = 0;
} else if(t_len == 0) {
Stat.len = q_len;
+ Stat.match = 0;
Stat.gaps = q_len;
Stat.score = W1 + (q_len - 1) * U;
} else {
Stat.len = t_len;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = W1 + (t_len - 1) * U;
}
@@ -827,12 +840,14 @@ AlnScore NW_score(const long unsigned *template, const unsigned char *queryOrg,
n = pos[1];
nuc_pos = m + t_s;
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
while(E_ptr[n] != 0) {
if(nuc_pos == template_length) {
nuc_pos = 0;
}
if((E_ptr[n] & 7) == 1) {
+ ++Stat.match;
++nuc_pos;
E_ptr += (q_len + 1);
++n;
@@ -883,14 +898,17 @@ AlnScore NW_band_score(const long unsigned *template, const unsigned char *query
if(t_len == 0 || q_len == 0) {
if(t_len == q_len) {
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = 0;
} else if(t_len == 0) {
Stat.len = q_len;
+ Stat.match = 0;
Stat.gaps = q_len;
Stat.score = W1 + (q_len - 1) * U;
} else {
Stat.len = t_len;
+ Stat.match = 0;
Stat.gaps = 0;
Stat.score = W1 + (t_len - 1) * U;
}
@@ -1110,12 +1128,14 @@ AlnScore NW_band_score(const long unsigned *template, const unsigned char *query
E_ptr = E + (m * (bq_len + 1));
nuc_pos = m + t_s;
Stat.len = 0;
+ Stat.match = 0;
Stat.gaps = 0;
while(E_ptr[n] != 0) {
if(nuc_pos == template_length) {
nuc_pos = 0;
}
if((E_ptr[n] & 7) == 1) {
+ ++Stat.match;
++nuc_pos;
E_ptr += (bq_len + 1);
++q_pos;
=====================================
nw.h
=====================================
@@ -37,6 +37,7 @@ struct alnScore {
int score;
int len;
int pos;
+ int match;
int gaps;
};
=====================================
runkma.c
=====================================
@@ -135,7 +135,7 @@ char * nameLoad(Qseqs *name, FILE *infile) {
return (char *) name->seq;
}
-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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
int i, j, tmp_template, tmp_tmp_template, file_len, bestTemplate, tot;
int template, bestHits, t_len, start, end, aln_len, status, rand, sparse;
@@ -281,6 +281,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
fprintf(stderr, "# Running KMA.\n");
t0 = clock();
+
/* allocate stuff */
i = 1;
alnThreads = 0;
@@ -423,6 +424,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
}
kmaPipe(0, 0, inputfile, &i);
status |= i;
+
i = 0;
sfwrite(&i, sizeof(int), 1, frag_out_raw);
fflush(frag_out_raw);
@@ -488,6 +490,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
@@ -1150,6 +1153,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
@@ -1293,7 +1297,7 @@ int runKMA(char *templatefilename, char *outputfilename, char *exePrev, int ConC
}
/* Close files */
- fclose(seq_in);
+ close(seq_in_no);
fclose(res_out);
if(alignment_out) {
fclose(alignment_out);
@@ -1322,7 +1326,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, 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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
/* runKMA_MEM is a memory saving version of runKMA,
at the cost it chooses best templates based on kmers
@@ -2227,6 +2231,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
@@ -2291,6 +2296,7 @@ int runKMA_MEM(char *templatefilename, char *outputfilename, char *exePrev, int
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
=====================================
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, 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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
-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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
=====================================
sam.c
=====================================
@@ -98,6 +98,7 @@ char * makeCigar(Qseqs *Cigar, const Aln *aligned) {
void saminit(Qseqs *template_name, FILE *name_file, int *template_lengths, int DB_size) {
+ fprintf(stdout, "@HD\tVN:1.6\tGO:reference\n");
while(--DB_size) {
fprintf(stdout, "@SQ\tSN:%s\tLN:%d\n", nameLoad(template_name, name_file), *++template_lengths);
}
=====================================
spltdb.c
=====================================
@@ -395,7 +395,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, 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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose) {
/* https://www.youtube.com/watch?v=LtXEMwSG5-8 */
@@ -1429,6 +1429,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
@@ -1493,6 +1494,7 @@ int runKMA_spltDB(char **templatefilenames, int targetNum, char *outputfilename,
thread->mq = mq;
thread->minlen = minlen;
thread->scoreT = scoreT;
+ thread->mrc = mrc;
thread->evalue = evalue;
thread->bcd = bcd;
thread->sam = sam;
=====================================
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, 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 xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
+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 mrc, double evalue, int bcd, int ref_fsa, int print_matrix, int print_all, int vcf, int xml, int sam, int nc, int nf, unsigned shm, int thread_num, int verbose);
=====================================
version.h
=====================================
@@ -17,4 +17,4 @@
* limitations under the License.
*/
-#define KMA_VERSION "1.3.9"
+#define KMA_VERSION "1.3.10"
View it on GitLab: https://salsa.debian.org/med-team/kma/-/commit/339b078952ca08927f0b733f18a91a69676c6ff7
--
View it on GitLab: https://salsa.debian.org/med-team/kma/-/commit/339b078952ca08927f0b733f18a91a69676c6ff7
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/20210206/640eeda2/attachment-0001.html>
More information about the debian-med-commit
mailing list