Commit 079f87c1 authored by sallet's avatar sallet
Browse files

start a method to find uncomatible gene

parent eea1b74e
......@@ -331,6 +331,106 @@ bool OneAltEst :: CompatibleWith(Prediction *pred)
return false;
}
Gene* OneAltEst :: GetUncompatibleGene(Prediction *pred)
{
int idxf, idxe=0;
Gene *g;
bool UncompatibleGeneFound;
//locate overlapping gene
if (PAR.getI("AltEst.strandSpecific"))
g = pred->FindGene(start,end, strand);
else
g = pred->FindGene(start,end);
// no overlapping gene
if (g == NULL)
{
UncompatibleGeneFound = false;
// return false;
return NULL;
}
else
{
fprintf(stdout, "J'ai trouve un gene ! \n ");
g->Print();
// stugy g to be sure it is incompatible
// check first exon start is in transcribed matured region
int nbFeature = g->nbFea();
for (idxf = 0; idxf < nbFeature; idxf++)
{
if ((g->vFea[idxf]->start-1 <= vi_ExonStart[idxe]) &&
(g->vFea[idxf]->end-1 >= vi_ExonStart[idxe]))
{
if (! g->vFea[idxf]->IsTranscribedAndUnspliced() )
{
UncompatibleGeneFound = true;
return g;
}
else break;
}
}
idxf++;
bool firstOk = false;
// test des frontieres suivantes
while ((idxe < exonsNumber-1) && (idxf < nbFeature))
{
if (!firstOk && (g->vFea[idxf]->start-1 == vi_ExonEnd[idxe]+1))
firstOk = true;
if (!firstOk && ( ! g->vFea[idxf]->IsTranscribedAndUnspliced()) ) // IG or intron: broken
{
UncompatibleGeneFound = true;
return g;
//return false;
}
if (firstOk && (g->vFea[idxf]->end == vi_ExonStart[idxe+1]))
{
if ( !g->vFea[idxf]->IsTranscribedAndUnspliced() ) // IF or intron
{
idxe++;
idxf++;
firstOk = false;
continue;
}
else
{
UncompatibleGeneFound = true;
return g;
//return false;
}
}
idxf++;
}
// test de la fin
for (; idxf < nbFeature; idxf++)
{
if ((g->vFea[idxf]->start-1 <= vi_ExonEnd[idxe]) &&
(g->vFea[idxf]->end-1 >= vi_ExonEnd[idxe]))
{
//return (g->vFea[idxf]->IsTranscribedAndUnspliced());
if (!g->vFea[idxf]->IsTranscribedAndUnspliced())
{
return g;
}
else
{
return NULL;
}
}
}
}
//return false;
fprintf(stdout, "je suis a la fin\n");
return NULL;
}
// --------------------------------------
// Penalize according to the EST
// --------------------------------------
......
......@@ -64,6 +64,7 @@ class OneAltEst
int minEstLen, int maxEstLen);
bool IsInconsistentWith(OneAltEst*);
bool CompatibleWith(Prediction *pred);
Gene* GetUncompatibleGene(Prediction *pred);
void Print ();
inline void UpdateBoundaries() { start = vi_ExonStart[0]; end = vi_ExonEnd[vi_ExonEnd.size()-1]; };
inline char* GetId() { return id; };
......
......@@ -306,8 +306,126 @@ std::vector<Prediction*> AllAltPredict (DNASeq* TheSeq, int fromPos, int toPos,
if (AltEstDB->voae_AltEst[altidx].IsToRemove()) continue;
if (AltEstDB->voae_AltEst[altidx].CompatibleWith(pred)) continue;
// positions around the AltEst hit
localFrom = Max(fromPos, AltEstDB->voae_AltEst[altidx].GetStart()-RepredictMargin);
localTo = Min(toPos, AltEstDB->voae_AltEst[altidx].GetEnd()+RepredictMargin);
// DynaProg end at the lastpos + 1 to account for final signals.
int FirstNuc = (Forward ? localFrom : localTo+1);
int LastNuc = (Forward ? localTo+1 : localFrom);
Dag->Init(FirstNuc-Dir, LastNuc+Dir);
Dag->WeightThePrior();
for (int nuc = FirstNuc; nuc != LastNuc+Dir; nuc += Dir)
{
// recuperation des infos
MSensor->GetInfoAt(TheSeq, nuc, &Data);
AltEstDB->Penalize(altidx,nuc,&Data);
if (Forward)
Dag->ShortestPathAlgoForward(nuc,Data);
else
Dag->ShortestPathAlgoBackward(nuc,Data);
if (nuc && (nuc % GCLatency == 0)) Dag->MarkAndSweep(nuc,GCVerbose,GCLatency);
}
Dag->WeightThePrior();
Dag->BuildPrediction(localFrom, localTo, Forward);
Dag->pred->TrimAndUpdate(TheSeq);
AltPred = Dag->pred;
if (AltPred)
{
if ( (AltPred->vGene[0]->cdsStart == -1) || (AltPred->vGene[0]->cdsEnd == -1))
{
delete Dag->pred;
Dag->Clean();
continue;
}
// Delete the gene of the alt prediction which doesn't overlap the EST
AltPred->DeleteOutOfRange(AltEstDB->voae_AltEst[altidx].GetStart(),AltEstDB->voae_AltEst[altidx].GetEnd(), AltEstDB->voae_AltEst[altidx].GetStrand());
// If genes overlapping the EST was found and if the prediction is original
if ( (AltPred->nbGene > 0) && (AltPred->IsOriginal(pred,vPred,ExonBorderMatchThreshold)) )
{
fprintf(stderr,"Optimal path length = %.4f\n",- AltPred->optimalPath);
baseGene = pred->FindGene(AltPred->vGene[0]->trStart,AltPred->vGene[0]->trEnd, AltPred->vGene[0]->GetStrand());
if (baseGene)
{
baseGene->hasvariant++;
AltPred->vGene[0]->isvariant = true;
AltPred->vGene[0]->hasvariant = baseGene->hasvariant;
AltPred->vGene[0]->geneNumber = baseGene->geneNumber;
baseGene->tuStart = ( baseGene->tuStart ) ? Min(baseGene->tuStart,AltPred->vGene[0]->trStart)
: Min(baseGene->trStart,AltPred->vGene[0]->trStart);
baseGene->tuEnd = ( baseGene->tuEnd ) ? Max(baseGene->tuEnd,AltPred->vGene[0]->trEnd)
: Max(baseGene->trEnd,AltPred->vGene[0]->trEnd);
}
else
{
fprintf(stderr,"New gene predicted by alternative spliced gene prediction.\n");
AltPred->vGene[0]->geneNumber = pred->nbGene + newGene++;
}
vPred.push_back(AltPred);
}
else
{
delete Dag->pred;
}
}
Dag->Clean();
}
return vPred;
}
// -------------------------------------------------------------------------
// Compute alternative Predictions based on EST
// -------------------------------------------------------------------------
std::vector<Prediction*> AllAltPredict2 (DNASeq* TheSeq, int fromPos, int toPos, MasterSensor* MSensor,
AltEst *AltEstDB, Prediction *pred)
{
int ExonBorderMatchThreshold = PAR.getI("AltEst.ExonBorderMatchThreshold");
int RepredictMargin = PAR.getI("AltEst.RepredictMargin");
int GCVerbose = PAR.getI("EuGene.VerboseGC");
int GCLatency = PAR.getI("EuGene.GCLatency");
DATA Data;
int Forward = 1;//PAR.getI("Sense");
int Dir = (Forward ? 1 : -1);
int localFrom,localTo;
std::vector <Prediction*> vPred;
Prediction* AltPred;
Gene* baseGene;
int newGene = 0; // if a splice variant has no base gene, it is a "new" gene. counter needed for gene number
// Active the appropriated tracks according to the eugene mode Prokaryote or Eukaryote
InitActiveTracks(0);
// Create a dag which would be reinitialze for each altEst
DAG* Dag = new DAG(TheSeq, PAR);
Dag->LoadDistLength();
for (int altidx = 0; altidx < AltEstDB->totalAltEstNumber; altidx++)
{
if (AltEstDB->voae_AltEst[altidx].IsToRemove()) continue;
fprintf (stdout, "\nAltEst %i - %i %c CompatibleWith : %i", AltEstDB->voae_AltEst[altidx].GetStart(),
AltEstDB->voae_AltEst[altidx].GetEnd(), AltEstDB->voae_AltEst[altidx].GetStrand(), AltEstDB->voae_AltEst[altidx].CompatibleWith(pred));
fprintf (stdout, "GetUncompatibleGene :\n");
Gene* ggg= AltEstDB->voae_AltEst[altidx].GetUncompatibleGene(pred);
if (ggg)
{
ggg->Print();
}
else {
fprintf (stdout, "NULL\n");
}
if (AltEstDB->voae_AltEst[altidx].CompatibleWith(pred)) continue;
// positions around the AltEst hit
localFrom = Max(fromPos, AltEstDB->voae_AltEst[altidx].GetStart()-RepredictMargin);
localTo = Min(toPos, AltEstDB->voae_AltEst[altidx].GetEnd()+RepredictMargin);
// DynaProg end at the lastpos + 1 to account for final signals.
int FirstNuc = (Forward ? localFrom : localTo+1);
int LastNuc = (Forward ? localTo+1 : localFrom);
......@@ -653,7 +771,7 @@ int main (int argc, char * argv [])
if (true)
{
vPred = AllAltPredict (TheSeq, fromPos, toPos, MS, AltEstDB, pred);
vPred = AllAltPredict2 (TheSeq, fromPos, toPos, MS, AltEstDB, pred);
}
else
{
......
Supports Markdown
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