Commit a0d13502 authored by Mathieu Giraud's avatar Mathieu Giraud

vidjil.cpp, -c segment: uses multigermline

The CMD_SEGMENT scenario now tries the germlines in the multigermline (until the read is segmented).
parent d4034f77
...@@ -653,7 +653,6 @@ int main (int argc, char **argv) ...@@ -653,7 +653,6 @@ int main (int argc, char **argv)
MultiGermline *multigermline = new MultiGermline(); MultiGermline *multigermline = new MultiGermline();
if (command == CMD_GERMLINES || command == CMD_CLONES || command == CMD_WINDOWS)
{ {
cout << "Build Kmer indexes" << endl ; cout << "Build Kmer indexes" << endl ;
...@@ -677,7 +676,9 @@ int main (int argc, char **argv) ...@@ -677,7 +676,9 @@ int main (int argc, char **argv)
germline = new Germline(germline_system, 'X', germline = new Germline(germline_system, 'X',
rep_V, rep_D, rep_J, rep_V, rep_D, rep_J,
delta_min, delta_max); delta_min, delta_max);
germline->new_index(seed);
if (command == CMD_WINDOWS || command == CMD_CLONES)
germline->new_index(seed);
multigermline->insert(germline); multigermline->insert(germline);
} }
...@@ -1315,30 +1316,39 @@ int main (int argc, char **argv) ...@@ -1315,30 +1316,39 @@ int main (int argc, char **argv)
//reads = OnlineFasta(f_reads, 1, " "); //reads = OnlineFasta(f_reads, 1, " ");
Fasta rep_V(f_rep_V, 2, "|", cout);
Fasta rep_D(f_rep_D, 2, "|", cout);
Fasta rep_J(f_rep_J, 2, "|", cout);
Germline *germline;
germline = new Germline(germline_system, 'X',
rep_V, rep_D, rep_J,
delta_min, delta_max);
while (reads->hasNext()) while (reads->hasNext())
{ {
reads->next(); reads->next();
FineSegmenter s(reads->getSequence(), germline, segment_cost);
if (s.isSegmented()) { Sequence seq = reads->getSequence() ;
if (segment_D) bool segmented = false ;
s.FineSegmentD(germline);
cout << s << endl; for (list<Germline*>::const_iterator it = multigermline->germlines.begin(); it != multigermline->germlines.end(); ++it)
} else { {
cout << "Unable to segment" << endl; Germline *germline = *it ;
cout << reads->getSequence();
cout << endl << endl; FineSegmenter s(seq, germline, segment_cost);
}
} if (s.isSegmented())
{
if (segment_D && germline->rep_4.size())
s.FineSegmentD(germline);
cout << s << endl;
segmented = true ;
break ;
}
else
{
if (verbose)
cout << "# " << germline->code << ": unable to segment" << endl;
}
}
if (!segmented)
{
seq.label += " unsegmented";
cout << seq << endl;
}
}
} else { } else {
cerr << "Ooops... unknown command. I don't know what to do apart from exiting!" << endl; cerr << "Ooops... unknown command. I don't know what to do apart from exiting!" << endl;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment