Commit b06a678c authored by flothoni's avatar flothoni Committed by Mathieu Giraud

clone.js; Add a function that allow to found the best matching sequence

take an array of sequences
parent a88b5521
......@@ -281,6 +281,69 @@ Clone.prototype = {
},
/**
* Return the best matching sequence from a list of sequence
* Can be used to find the best matching primer from a list of primers
* If no perfect match, can search in germline sequence (longer) and use alignment option
* @param {Array} sequences Array of sequence to search in clone
* @param {Array} extend Array of germline sequence to extend. Available option are 5 and 3
* @return {Array} [best sequence found, use extension]
*/
getBestMatchingSequence: function(sequences, extend){
var best_seq = []
var best_rst = []
var best_score = 0
// Look for perfect match in clone sequence
for (var seq_pos = 0; seq_pos < sequences.length; seq_pos++) {
var sequence = sequences[seq_pos]
if (this.sequence.indexOf(sequence) != -1) {
best_seq.push(sequence)
}
}
if (best_seq.length > 0){
// TODO; if multiple perfect match, return longer sequence
return [best_seq[0], false]
}
// No perfect match, look in germline sequence with alignement tool
var genes = (extend != undefined) ? extend : [5, 3]
// Look if germline sequence is available
for (var g = 0; g < genes.length; g++) {
var gene_way = genes[g]
var gene = this.getGene(gene_way);
if (gene !== undefined) {
var germName = this.germline.substring(0, 3);
if (this.m.germline[germName] != undefined) {
germseq_found = this.m.findGermlineFromGene(gene);
if (germseq_found != undefined){
// Get germline sequence; look for best sequence
germseq = germseq_found.toUpperCase().replace(/\./g, '');
for (seq_pos = 0; seq_pos < sequences.length; seq_pos++) {
var sequence = sequences[seq_pos]
var rst = bsa_align(true, germseq, sequence, [1, -2], [-2, -1])
if (rst[0] > best_score){
best_seq = [sequence]
best_rst = [rst]
best_score = rst[0]
} else if (rst[0] == best_score){
best_seq.push(sequence)
best_rst.push(rst)
}
}
}
}
}
}
if (!best_seq.length) {
return [undefined, false]
}
return [best_seq[0], true]
},
/**
* Get the amino-acid sequence of the provided field (in the seg part)
*/
......
......@@ -735,3 +735,69 @@ QUnit.test("clonedb", function(assert) {
assert.equal(m.clones[1].numberInCloneDB(), undefined, "");
assert.equal(m.clones[1].numberSampleSetInCloneDB(), undefined, "");
});
QUnit.test("primer detection align", function(assert) {
var m = new Model();
m.parseJsonData(json_data, 100);
m.populatePrimerSet()
// Clone with seg start and stop position
var clone8 = {
"germline": "IGH",
"id": "id7",
"name": "clone with exact position",
"reads": [0,0,0,0],
"seg": {
"3": {
"delLeft": 0,
"name": "IGHJ1*01",
"start": 297
},
"5": {
"delRight": 0,
"name": "IGHV1-2*01",
"stop": 296
}
},
"sequence": "init",
}
var clone = new Clone(clone8, m, 0, c_attributes);
m.initClones();
primer5_ecngs_igh = m.primersSetData.ecngs.IGH.primer5
primer3_ecngs_igh = m.primersSetData.ecngs.IGH.primer3
primer3_biomed_igh = m.primersSetData.biomed2.IGH.primer5
primer3_biomed_igh = m.primersSetData.biomed2.IGH.primer3
primer5_ecngs_trd = m.primersSetData.ecngs.TRD.primer5
primer3_ecngs_trd = m.primersSetData.ecngs.TRD.primer3
//"IGHV1-2*01 // IGHJ1*01",
// var igh_V1_2_J1_full
clone.sequence = "CAGGTGCAGCTGGTGCAGTCTGGGGCTGAGGTGAAGAAGCCTGGGGCCTCAGTGAAGGTCTCCTGCAAGGCTTCTGGATACACCTTCACCGGCTACTATATGCACTGGGTGCGACAGGCCCCTGGACAAGGGCTTGAGTGGATGGGACGGATCAACCCTAACAGTGGTGGCACAAACTATGCACAGAAGTTTCAGGGCAGGGTCACCAGTACCAGGGACACGTCCATCAGCACAGCCTACATGGAGCTGAGCAGGCTGAGATCTGACGACACGGTCGTGTATTACTGTGCGAGAGAGCTGAATACTTCCAGCACTGGGGCCAGGGCACCCTGGTCACCGTCTCCTCAG"
assert.deepEqual( clone.getBestMatchingSequence(primer5_ecngs_igh, false), ["CTGGGTGCGACAGGCCCCT", false], "seq full; primer5, no extend" )
assert.deepEqual( clone.getBestMatchingSequence(primer5_ecngs_igh, [5]), ["CTGGGTGCGACAGGCCCCT", false], "seq full; primer5, extend 5" )
// not full match, but if extend, have bas_align call
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, []), [undefined, false],"seq full; primer3, no extend")
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, [5]), ["GGTCACCGTCTCCTCAGGTAAG", true],"seq full; primer3, extend 5" )
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, [3]), ["GGTCACCGTCTCCTCAGGTAAG", true],"seq full; primer3, extend 3" )
// Sequence with 5' and 3' deletion
clone.seg[5].stop = 176
clone.seg[3].start = 177
clone.sequence = "CCTGGACAAGGGCTTGAGTGGATGGGACGGATCAACCCTAACAGTGGTGGCACAAACTATGCACAGAAGTTTCAGGGCAGGGTCACCAGTACCAGGGACACGTCCATCAGCACAGCCTACATGGAGCTGAGCAGGCTGAGATCTGACGACACGGTCGTGTATTACTGTGCGAGAGAGCTGAATACTTCCAGCACTGGGGC"
assert.deepEqual( clone.getBestMatchingSequence(primer5_ecngs_igh, false), [undefined, false], "seq short; primer5, no extend" )
assert.deepEqual( clone.getBestMatchingSequence(primer5_ecngs_igh, [5]), ["CTGGGTGCGACAGGCCCCT", true], "seq short; primer5, extend 5" )
// not full match, but if extend, have bas_align call
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, []), [undefined, false], "seq short; primer3, no extend")
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, [5]), ["GGTCACCGTCTCCTCAGGTAAG", true], "seq short; primer3, extend 5" )
assert.deepEqual( clone.getBestMatchingSequence(primer3_ecngs_igh, [3]), ["GGTCACCGTCTCCTCAGGTAAG", true], "seq short; primer3, extend 3" )
});
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