Commit aba59213 authored by Mathieu Giraud's avatar Mathieu Giraud
Browse files

segment.cpp: SEG_METHOD_ONE, KmerSegmenter

This is evil. Vidjil was designed to handle recombinations,
not to be a read mapper. See #1724.

Moreover, we should extract a window in a more reproducible way.
See #2643.
parent f0b48068
......@@ -20,7 +20,8 @@ enum SEGMENTATION_METHODS {
SEG_METHOD_53, // Regular or incomplete germlines, 5'-3'
SEG_METHOD_543, // Regular or incomplete germlines, 5'-3', with an additional middle gene (such a D gene)
SEG_METHOD_MAX12, // Pseudo-germline, most two frequent kmer affectations (-2)
SEG_METHOD_MAX1U // Pseudo-germline, most frequent kmer affection and unknwon affectation (-4)
SEG_METHOD_MAX1U, // Pseudo-germline, most frequent kmer affection and unknwon affectation (-4)
SEG_METHOD_ONE // Map a read onto a genomic region, without recombination. Evil.
} ;
......
......@@ -444,6 +444,41 @@ KmerSegmenter::KmerSegmenter(Sequence seq, Germline *germline, double threshold,
reversed = (nb_strand[0] > nb_strand[1]) ;
if (germline->seg_method == SEG_METHOD_ONE) {
KmerAffectAnalyser kaa(*(germline->index), sequence);
KmerAffect kmer = KmerAffect(germline->affect_4, 1, germline->seed.size());
int c = kaa.count(kmer);
// E-value
double pvalue = kaa.getProbabilityAtLeastOrAbove(kmer, c);
evalue = pvalue * multiplier ;
if (evalue >= threshold)
{
because = UNSEG_TOO_FEW_ZERO ;
return ;
}
segmented = true ;
because = reversed ? SEG_MINUS : SEG_PLUS ;
int pos = sequence.size() / 2 ;
info = "=" + string_of_int(c) + " @" + string_of_int(pos) ;
// getJunction() will be centered on pos
box_V->end = pos;
box_J->start = pos;
finishSegmentation();
return ;
}
if ((germline->seg_method == SEG_METHOD_MAX12)
|| (germline->seg_method == SEG_METHOD_MAX1U))
{ // Pseudo-germline, MAX12 and MAX1U
......
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