Skip to content
Commits on Source (4)
Show whitespace changes
Inline Side-by-side
alnfrags.c
View file @ 8e34dd7c
......@@ -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;
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;
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;
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;
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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -42,6 +42,8 @@ struct assem {
unsigned len;
unsigned aln_len;
unsigned size;
unsigned fragmentCountAln;
unsigned readCountAln;
};
struct assembly {
......@@ -62,6 +64,7 @@ struct assemble_thread {
int file_count;
int spin;
int mq;
int minlen;
int bcd;
int sam;
int thread_num;
......
debian/changelog
View file @ 8e34dd7c
kma (1.2.19-1) unstable; urgency=medium
* Team upload.
* New upstream version
-- Steffen Moeller <moeller@debian.org> Fri, 27 Dec 2019 20:12:38 +0100
kma (1.2.16-1) unstable; urgency=medium
* Team upload.
......
ef.c
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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
View file @ 8e34dd7c
......@@ -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;
......