Commit 0cf35de8 authored by Mathieu Giraud's avatar Mathieu Giraud
Browse files

Merge branch 'cdr3' - experimental '-3' option

parents 7c8f49e6 9f824625
......@@ -26,6 +26,8 @@
#include "tools.h"
#include "affectanalyser.h"
#include <sstream>
#include <string>
#include <math.h>
Segmenter::~Segmenter() {}
......@@ -518,6 +520,9 @@ FineSegmenter::FineSegmenter(Sequence seq, Germline *germline, Cost segment_c)
Dend=0;
segment_cost=segment_c;
CDR3start = -1;
CDR3end = -1;
// TODO: factoriser tout cela, peut-etre en lancant deux segmenteurs, un +, un -, puis un qui chapote
// Strand +
......@@ -669,7 +674,6 @@ FineSegmenter::FineSegmenter(Sequence seq, Germline *germline, Cost segment_c)
info = string_of_int(Vend + FIRST_POS) + " " + string_of_int(Jstart + FIRST_POS) ;
finishSegmentation();
}
......@@ -781,6 +785,64 @@ void FineSegmenter::FineSegmentD(Germline *germline){
}
}
void FineSegmenter::findCDR3(){
string str = getSequence().sequence;
list<string> codon_start;
codon_start.push_back("TGT");
codon_start.push_back("TGC");
list<string> codon_end;
codon_end.push_back("TTT");
codon_end.push_back("TTC");
codon_end.push_back("TGG");
list<int> p_start;
list<int> p_end;
size_t loc;
std::list<string>::const_iterator it;
for (it = codon_start.begin(); it != codon_start.end(); ++it) {//filter 1 : start codon must be in V
loc = 0;
while ( loc != string::npos & loc < Vend){
loc = str.find(*it, loc+3);
if (loc != string::npos & loc < Vend) {
p_start.push_front(loc);
}
}
}
for (it = codon_end.begin(); it != codon_end.end(); ++it) {//filter 2 : end codon must be in J
loc = Jstart;
while ( loc != string::npos){
loc = str.find(*it, loc+3);
if (loc != string::npos) {
p_end.push_back(loc);
}
}
}
CDR3start = -1;
CDR3end = -1;
std::list<int>::const_iterator it1;
for (it1 = p_start.begin(); it1 != p_start.end(); ++it1) {
std::list<int>::const_iterator it2;
for (it2 = p_end.begin(); it2 != p_end.end(); ++it2) {
if ( (*it2-*it1)%3 == 0){ //filter 3 : start/stop codon must be seprated by a multiple of 3
if ( fabs((*it2-*it1)-36 ) < fabs((CDR3end-CDR3start)-36) ){ //filter 4 : cdr3 length must be close to 12 AA
CDR3start = *it1;
CDR3end = *it2;
}
}
}
}
}
JsonList FineSegmenter::toJsonList(Germline *germline){
JsonList result;
......@@ -801,6 +863,15 @@ JsonList FineSegmenter::toJsonList(Germline *germline){
seg.add("3", germline->rep_3.label(best_J));
seg.add("3start", Jstart);
if (CDR3start >= 0)
{
JsonList *json_cdr;
json_cdr=new JsonList();
json_cdr->add("start", CDR3start);
json_cdr->add("stop", CDR3end);
seg.add("cdr3", *json_cdr);
}
result.add("seg", seg);
}
......
......@@ -46,6 +46,7 @@ protected:
string sequence;
int Vend, Jstart;
int Dstart, Dend;
int CDR3start, CDR3end;
bool reversed, segmented, dSegmented;
string removeChevauchement();
......@@ -208,6 +209,7 @@ class FineSegmenter : public Segmenter
* @param germline: germline used
*/
void FineSegmentD(Germline *germline);
void findCDR3();
JsonList toJsonList(Germline *germline);
......
......@@ -190,6 +190,7 @@ void usage(char *progname)
<< "Fine segmentation options (second pass, see warning in doc/algo.org)" << endl
<< " -f <string> use custom Cost for fine segmenter : format \"match, subst, indels, homo, del_end\" (default "<<VDJ<<" )"<< endl
<< " -3 CDR3 detection (experimental)" << endl
<< endl
<< "Additional clustering" << endl
......@@ -262,7 +263,7 @@ int main (int argc, char **argv)
int minPts = DEFAULT_MINPTS ;
Cost cluster_cost = DEFAULT_CLUSTER_COST ;
Cost segment_cost = DEFAULT_SEGMENT_COST ;
bool detect_CDR3 = false;
int save_comp = 0;
int load_comp = 0;
......@@ -304,7 +305,7 @@ int main (int argc, char **argv)
//$$ options: getopt
while ((c = getopt(argc, argv, "AhaiIg:G:V:D:J:k:r:vw:e:C:f:l:c:m:M:N:s:b:Sn:o:L%:y:z:uU")) != EOF)
while ((c = getopt(argc, argv, "AhaiIg:G:V:D:J:k:r:vw:e:C:f:l:c:m:M:N:s:b:Sn:o:L%:y:z:uU3")) != EOF)
switch (c)
{
......@@ -494,7 +495,10 @@ int main (int argc, char **argv)
break;
// Fine segmentation
case '3':
detect_CDR3 = true;
break;
case 'f':
segment_cost=strToCost(optarg, VDJ);
break;
......@@ -1120,7 +1124,11 @@ int main (int argc, char **argv)
if (segmented_germline->rep_4.size())
seg.FineSegmentD(segmented_germline);
if (detect_CDR3)
seg.findCDR3();
// Output representative, possibly segmented...
// to stdout, CLONES_FILENAME, and CLONE_FILENAME-*
cout << seg << endl ;
......@@ -1362,6 +1370,10 @@ int main (int argc, char **argv)
{
if (germline->rep_4.size())
s.FineSegmentD(germline);
if (detect_CDR3)
s.findCDR3();
cout << s << endl;
segmented = true ;
break ;
......
......@@ -549,26 +549,30 @@ Segment.prototype = {
findPotentialField : function () {
result = [""];
var clone = this.m.clone(1);
for (var i in this.m) {
if (this.isDNA(this.m[i])){
result.push(i)
if (this.isDNA(this.m[i]) ||this.isDNA(m[i]) ){
if (result.indexOf(i) == -1) result.push(i);
}
}
for (var i in clone) {
if (this.isDNA(clone[i])){
result.push(i)
for (var j=0; (j<10 & j<this.m.clones.length) ; j++){
var clone = this.m.clone(j);
console.log(clone)
for (var i in clone) {
if (this.isDNA(clone[i]) || this.isPos(clone[i]) ){
if (result.indexOf(i) == -1) result.push(i);
}
}
}
for (var i in clone) {
if (this.isDNA(clone.seg[i])){
result.push(i)
if (typeof clone.seg != 'undefined'){
for (var i in clone.seg) {
if (this.isDNA(clone.seg[i]) || this.isPos(clone.seg[i]) ){
if (result.indexOf(i) == -1) result.push(i);
}
}
}
}
return result;
},
......@@ -585,6 +589,16 @@ Segment.prototype = {
}
},
isPos : function (obj) {
if (obj == null) {
return false;
}else if (typeof obj.start != 'undefined' && typeof obj.stop != 'undefined') {
return true;
}else{
return false;
}
},
isAA : function (string) {
if (string == null) {
return false;
......@@ -652,21 +666,36 @@ Sequence.prototype = {
},
computeAAseq : function () {
var start = 0;
var stop = this.seq.length;
//for nikos we can use his cdr3 to find the correct start
var start = -1;
var stop = -1;
var clone = this.m.clone(this.id);
if (typeof clone["<option>_sequence.JUNCTION.raw nt seq</option>"] != "undefined"){
start = clone.sequence.indexOf(clone["<option>_sequence.JUNCTION.raw nt seq</option>"]) % 3;
if (typeof clone.seg != "undefined" && typeof clone.seg["cdr3"] != "undefined"){
start = clone.seg["cdr3"].start
stop = clone.seg["cdr3"].stop
}else if (typeof clone["_sequence.JUNCTION.raw nt seq"] != "undefined"){
// .clntab. TODO: move this to fuse.py
var junc = clone["_sequence.JUNCTION.raw nt seq"][0]
start = clone.sequence.indexOf(junc)
stop = start + junc.length
}
var i=start;
for (var i=0; i<this.seq.length; i++) this.seqAA[i] =this.seq[i]; // "&nbsp";
var i = 0
while (i<this.seq.length){
if (i < start || i > stop)
{
i++
continue
}
var code = "";
var pos;
while (code.length<3 & i<this.seq.length){
while (code.length<3 & i<=stop){
if (this.seq[i] != "-") {
code += this.seq[i];
this.seqAA[i] = "&nbsp";
......@@ -715,7 +744,12 @@ Sequence.prototype = {
var jColor = "";
if (this.m.colorMethod == "J") jColor = "style='color : " + clone.colorJ + "'";
var window_start = this.pos[clone.sequence.indexOf(clone.id)]
var window_start = this.pos[clone.sequence.indexOf(clone.id)];
if (typeof clone.seg != "undefined" && typeof clone.seg["cdr3"] != "undefined"){
window_start = clone.seg["cdr3"].start;
}else if (typeof clone["<option>_sequence.JUNCTION.raw nt seq</option>"] != "undefined"){
window_start = clone.sequence.indexOf(clone["<option>_sequence.JUNCTION.raw nt seq</option>"]);
}
var highlights = [];
for (var i in segment.highlight){
......@@ -773,7 +807,7 @@ Sequence.prototype = {
if (typeof clone[field] != 'undefined'){
p = clone[field]; //check clone meta-data
}else if (typeof clone[field] != 'undefined'){
}else if (typeof clone.seg != 'undefined' &&typeof clone.seg[field] != 'undefined'){
p = clone.seg[field]; //check clone seg data
}else if (typeof this.m[field] != 'undefined'){
p = this.m[field]; //check model
......@@ -809,11 +843,11 @@ tableAA = {
'TCG' : 'S',
'TAT' : 'Y',
'TAC' : 'Y',
'TAA' : '#',
'TAG' : '#',
'TAA' : '*',
'TAG' : '*',
'TGT' : 'C',
'TGC' : 'C',
'TGA' : '#',
'TGA' : '*',
'TGG' : 'W',
'CTT' : 'L',
'CTC' : 'L',
......
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