Je réalise un projet de bioinformatique dans un cadre scolaire, mais ne possède pas de connaissances en biologie. Heureusement, je m'y connais mieux en informatique... Pour résumer, mon projet consiste à construire un arbre phylogénétique de différents virus grâce aux séquences ADN du NCBI. En fait, mon intérêt se porte davantage sur l'informatique que sur la bioinformatique donc je voudrais réaliser un programme qui crée un tel arbre sans me soucier des séquences qu'il doit traiter ni d'une interprétation épidiomologique des résultats. Je dois cependant obtenir un arbre fiable.
Ma première étape que je m'étais fixée était la recherche de régions codantes. J'ai cru comprendre qu'elle était facultative, mais je dois avouer ne pas être en mesure de justifier clairement son intérêt... J'obtiens donc, en concaténant mes régions codantes, une nouvelle séquence ADN sur laquelle je vais travailler.
Ensuite, j'utilise l'algorithme de Needleman-Wunsch pour calculer les distances de Levensthein entre chacune de mes séquences et créer une matrice des distances. Je crois avoir compris le fonctionnement de l'algorithme mais vous propose tout de même mon code au cas où :
int needlemanWunsch(Nucleotide sequence1[], Nucleotide sequence2[], int tailleSequence1, int tailleSequence2)
{
int* lignes[2] = { NULL, NULL };
int* colonnes[2] = { NULL,NULL };
lignes[0] = malloc(sizeof(int) * (tailleSequence1+1));
lignes[1] = malloc(sizeof(int) * (tailleSequence1+1));
colonnes[0] = malloc(sizeof(int) * (tailleSequence2+1));
colonnes[1] = malloc(sizeof(int) * (tailleSequence2+1));
int const COUTS_SUBSTITUTIONS[4][4] = { {0,TRANSVERSION,TRANSITION,TRANSVERSION},
{TRANSVERSION,0,TRANSVERSION,TRANSITION},
{TRANSITION,TRANSVERSION,0,TRANSVERSION},
{TRANSVERSION,TRANSITION,TRANSVERSION,0} };
if (lignes[0] == NULL || lignes[1] == NULL || colonnes[0] ==NULL || colonnes[1] ==NULL)
{
printf("echec (needlemanWunsch) 0");
exit(EXIT_FAILURE);
}
for (int i = 0; i < tailleSequence1; i++)
{
lignes[0][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i < tailleSequence1; i++)
{
lignes[1][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i < tailleSequence2; i++)
{
colonnes[0][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i < tailleSequence2; i++)
{
colonnes[1][i] = i * COUT_INSERTION_DELETION;
}
int k = 1;
while (k <= tailleSequence1 && k <= tailleSequence2)
{
lignes[k % 2][k] = min(lignes[(k + 1) % 2][k - 1] + COUTS_SUBSTITUTIONS[sequence1[k-1]][sequence2[k-1]],
min(lignes[(k + 1) % 2][k - 1] + COUT_INSERTION_DELETION,
colonnes[(k + 1) % 2][k - 1] + COUT_INSERTION_DELETION));
colonnes[k % 2][k] = lignes[k % 2][k];
for (int i = k; i < tailleSequence1; i++)
{
lignes[k % 2][i] = min(lignes[(k + 1)%2][i - 1] + COUTS_SUBSTITUTIONS[sequence1[i]][sequence2[k-1]],
min(lignes[(k + 1) % 2][i-1] + COUT_INSERTION_DELETION,
colonnes[(k + 1) % 2][i - 1] + COUT_INSERTION_DELETION));
}
for (int j = k; j < tailleSequence2; j++)
{
colonnes[k % 2][j] = min(colonnes[(k + 1) % 2][j - 1] + COUTS_SUBSTITUTIONS[sequence1[k-1]][sequence2[j]],
min(colonnes[(k + 1) % 2][j - 1] + COUT_INSERTION_DELETION,
lignes[(k + 1) % 2][j - 1] + COUT_INSERTION_DELETION));
}
k++;
}
int cout = lignes[(k + 1) % 2][tailleSequence1 - 1];
if (tailleSequence1 > tailleSequence2)
{
cout = colonnes[(k + 1) % 2][tailleSequence2 - 1];
}
free(lignes[0]);
free(lignes[1]);
free(colonnes[0]);
free(colonnes[1]);
return cout;
}
Hélas, je n'obtiens pas de résultat satisfaisant et je ne sais pas si le problème est algorithmique ou provient d'une subtilité biologique (qui pourrait ne pas être si subtile que ça pour les initiés). Voici les résultats que j'ai obtenus :
- avec UPGMA (sans recherche des régions codantes) :
- avec Neighbor Joining (sans recherche des régions codantes) :
- avec UPGMA (avec recherche des régions codantes) :
- avec Neighbor Joining (avec recherche des régions codantes) :
Voilà, si quelqu'un a le courage de m'expliquer les erreurs sans doute nombreuses qui conduisent à de tels arbres...
N'hésitez pas à demander si vous pensez qu'il y a besoin de plus d'informations ou de lignes de code pour comprendre le problème (même si je suis sûr que vous n'hésiterez pas).
Merci d'avance de prendre le temps d'aider un boulet comme moi (mais un gentil boulet) !
Juste en lisant en diagonale le code que tu as publié je remarque la longueur excessive des fonctions (visuellement). Dis-toi que chaque commentaire qui explique le bloc de code qui le suis devrait certainement donner naissance à une fonction à part entière.
J'ai vu un peu trop de malloc qui indique que tu n'as pas assez isolé de sdd sans doute.
Edit: ah oui, as-tu déjà utilisé un debuger ?
- Edité par White Crow 30 décembre 2021 à 12:57:55
Ça va être difficile de t'aider sans connaître les algorithmes que tu utilises sans parler des notions de biologie...
Je vais essayer quand même. Si j'ai bien compris, le problème est le suivant : tu calcules des distances ; ces distances ne sont pas correctes. C'est bien ça ?
J'ai une petite expérience du calcul numérique, dans ce cas voici ce que je ferais :
Construire un exemple simple sur un bout de papier, et calculer les distances à la main en écrivant toutes les étapes. Ce sera la référence.
Mettre dans le programme des points d'arrêt correspondant aux différentes étapes pour consulter les valeurs des variables correspondant au calcul de référence (on peut utiliser des 'printf' − c'est ce que je fais, j'ai pas honte !)
Suivre le déroulement du programme en comparant les valeurs affichées par le programme et les valeurs du calcul de référence. L'endroit où ça diverge indique l'endroit où il y a une erreur. En général, une fois qu'on a isolé l'instruction qui donne un résultat faux, il devient relativement facile de trouver la cause de l'erreur.
Corriger l'erreur et recommencer (au cas où il y en aurait d'autres).
Je me souviens d'une ancienne discussion où quelqu'un avait programmé la méthode du pivot sauf qu'elle ne donnait pas de bons résultats. S'en était suivie une discussion de plus d'une semaine. Malgré mes conseils (je lui avais conseillé la même chose) il n'a pas cherché à piste l'erreur ligne après ligne. Un soir, j'en ai eu marre, j'ai fait son boulot, ça m'a pris une heure. Avec plein de 'printf' partout j'ai trouvé l'erreur. Elle était subtile et je ne l'aurais jamais remarquée en examinant le listing. Une soirée à construire un cas simple à la main puis examiner le programme pas à pas, c'est quand même bien plus rapide qu'une semaine pour rien.
Juste en lisant en diagonale le code que tu as publié je remarque la longueur excessive des fonctions (visuellement). Dis-toi que chaque commentaire qui explique le bloc de code qui le suis devrait certainement donner naissance à une fonction à part entière.
J'ai vu un peu trop de malloc qui indique que tu n'as pas assez isolé de sdd sans doute.
Edit: ah oui, as-tu déjà utilisé un debuger ?
- Edité par White Crow il y a environ 1 heure
Je ne sais pas ce que signifie isoler un sdd et n'ai pas utilisé de debuger... Sinon, j'essayerai d'appliquer le conseil pour rendre mon code plus lisible, en espérant trouver du temps car j'ai pas mal de choses à faire à côté...
robun a écrit:
Ça va être difficile de t'aider sans connaître les algorithmes que tu utilises sans parler des notions de biologie...
Je vais essayer quand même. Si j'ai bien compris, le problème est le suivant : tu calcules des distances ; ces distances ne sont pas correctes. C'est bien ça ?
J'ai une petite expérience du calcul numérique, dans ce cas voici ce que je ferais :
Construire un exemple simple sur un bout de papier, et calculer les distances à la main en écrivant toutes les étapes. Ce sera la référence.
Mettre dans le programme des points d'arrêt correspondant aux différentes étapes pour consulter les valeurs des variables correspondant au calcul de référence (on peut utiliser des 'printf' − c'est ce que je fais, j'ai pas honte !)
Suivre le déroulement du programme en comparant les valeurs affichées par le programme et les valeurs du calcul de référence. L'endroit où ça diverge indique l'endroit où il y a une erreur. En général, une fois qu'on a isolé l'instruction qui donne un résultat faux, il devient relativement facile de trouver la cause de l'erreur.
Corriger l'erreur et recommencer (au cas où il y en aurait d'autres).
Je me souviens d'une ancienne discussion où quelqu'un avait programmé la méthode du pivot sauf qu'elle ne donnait pas de bons résultats. S'en était suivie une discussion de plus d'une semaine. Malgré mes conseils (je lui avais conseillé la même chose) il n'a pas cherché à piste l'erreur ligne après ligne. Un soir, j'en ai eu marre, j'ai fait son boulot, ça m'a pris une heure. Avec plein de 'printf' partout j'ai trouvé l'erreur. Elle était subtile et je ne l'aurais jamais remarquée en examinant le listing. Une soirée à construire un cas simple à la main puis examiner le programme pas à pas, c'est quand même bien plus rapide qu'une semaine pour rien.
Je m'y prends parfois de cette manière également et en effet, cette méthode est souvent utile. Cependant, le problème est plus de savoir si j'ai bien appliqué la marche à suivre. Je vérifierai quand même mes résultats sur un exemple pour assurer le coup, mais même si j'obtiens des résultats corrects, je ne saurais dire si je les utilise correctement... Pire encore, la bioinformatique ne fournit que des prédictions et en théorie, je devrais bel et bien avoir des erreurs. Il serait donc possible que mon programme soit correct mais pas optimisé, conduisant à un nombre trop grand d'erreurs, et donc à des arbres phylogénétiques à côté de la plaque. Autrement dit, il faudrait peut-être ajouter une étape intermédiaire pour obtenir quelque chose de plus satisfaisant (sans forcément obtenir un résultat parfait).
Je ne sais pas si je suis vraiment clair, mais le problème pourrait peut-être être résumé ainsi : mes résultats sont-ils mauvais à cause d'un choix d'algorithmes désuets ou à cause d'une mauvaise application/implémentation des algorithmes ?
SDD = Structure De Données ... subtil, ce White Crow ...
Phylogénétique: Branche de la génétique traitant des modifications génétiques au sein des espèces animales ou végétales.
Tu as beaucoup de malloc sans vérification. Je vais sortir de la boule à mites ma petite fonction pour malloc (qui vérifie les erreurs): void *getArea(char *name, int size) { void *area = malloc(size); if(area) return area; perror(name); printf("Required: %d\n", size); exit(EXIT_FAILURE); }
- Edité par PierrotLeFou 30 décembre 2021 à 15:00:44
Le Tout est souvent plus grand que la somme de ses parties.
Cependant, le problème est plus de savoir si j'ai bien appliqué la marche à suivre. Je vérifierai quand même mes résultats sur un exemple pour assurer le coup, mais même si j'obtiens des résultats corrects, je ne saurais dire si je les utilise correctement... Pire encore, la bioinformatique ne fournit que des prédictions et en théorie, je devrais bel et bien avoir des erreurs. Il serait donc possible que mon programme soit correct mais pas optimisé, conduisant à un nombre trop grand d'erreurs, et donc à des arbres phylogénétiques à côté de la plaque. Autrement dit, il faudrait peut-être ajouter une étape intermédiaire pour obtenir quelque chose de plus satisfaisant (sans forcément obtenir un résultat parfait).
2. Le programme est-il optimisé par rapport au besoin spécifique ?
3. Quelles étapes faut-il ajouter ?
On commence par lequel ? Ici c'est un forum de programmation, donc à mon avis c'est le 1. Il s'agit donc de vérifier des calculs de distance. Ne parlons pas du reste, ça va paralyser la vérification des calculs.
Question 1 : est-ce que les distances que tu obtiens sont correctes ?
A - Je ne sais pas. Dans ce cas il faut absolument que tu construises un exemple simple pour lequel tu sais faire les calculs.
B - Je m'en fiche, peu importe, ce n'est pas le problème. Dans ce cas on va passer au problème suivant.
C - Oui. Dans ce cas on va passer au problème suivant.
D - Non. Dans ce cas il faut suivre l'exécution du programme pas à pas comme je le préconisais.
Ta dernière réponse où tu mélanges les problèmes m'a bien embrouillé, j'espère que toi, par contre, tu as les idées claires !
SDD = Structure De Données ... subtil, ce White Crow ...
Phylogénétique: Branche de la génétique traitant des modifications génétiques au sein des espèces animales ou végétales.
Tu as beaucoup de malloc sans vérification. Je vais sortir de la boule à mites ma petite fonction pour malloc (qui vérifie les erreurs): void *getArea(char *name, int size) { void *area = malloc(size); if(area) return area; perror(name); printf("Required: %d\n", size); exit(EXIT_FAILURE); }
- Edité par PierrotLeFou il y a environ 23 heures
D'accord, je vois ce qu'est une structure de données. Je vais étudier cette fonction getArea de plus près.
robun a écrit:
Quentinus Prime a écrit:
Cependant, le problème est plus de savoir si j'ai bien appliqué la marche à suivre. Je vérifierai quand même mes résultats sur un exemple pour assurer le coup, mais même si j'obtiens des résultats corrects, je ne saurais dire si je les utilise correctement... Pire encore, la bioinformatique ne fournit que des prédictions et en théorie, je devrais bel et bien avoir des erreurs. Il serait donc possible que mon programme soit correct mais pas optimisé, conduisant à un nombre trop grand d'erreurs, et donc à des arbres phylogénétiques à côté de la plaque. Autrement dit, il faudrait peut-être ajouter une étape intermédiaire pour obtenir quelque chose de plus satisfaisant (sans forcément obtenir un résultat parfait).
2. Le programme est-il optimisé par rapport au besoin spécifique ?
3. Quelles étapes faut-il ajouter ?
On commence par lequel ? Ici c'est un forum de programmation, donc à mon avis c'est le 1. Il s'agit donc de vérifier des calculs de distance. Ne parlons pas du reste, ça va paralyser la vérification des calculs.
Question 1 : est-ce que les distances que tu obtiens sont correctes ?
A - Je ne sais pas. Dans ce cas il faut absolument que tu construises un exemple simple pour lequel tu sais faire les calculs.
B - Je m'en fiche, peu importe, ce n'est pas le problème. Dans ce cas on va passer au problème suivant.
C - Oui. Dans ce cas on va passer au problème suivant.
D - Non. Dans ce cas il faut suivre l'exécution du programme pas à pas comme je le préconisais.
Ta dernière réponse où tu mélanges les problèmes m'a bien embrouillé, j'espère que toi, par contre, tu as les idées claires !
- Edité par robun il y a environ 23 heures
Oui, il semblerait que tu exprimes mieux que moi le problème... Pour la question 1, il s'agissait de la réponse A, mais j'ai corrigé mon algorithme et me voilà soulagé car il donne de meilleurs résultats. Voici le nouveau :
int needlemanWunsch(Nucleotide sequence1[], Nucleotide sequence2[], int tailleSequence1, int tailleSequence2)
{
int* lignes[2] = { NULL, NULL };
int* colonnes[2] = { NULL,NULL };
lignes[0] = malloc(sizeof(int) * (tailleSequence1+1));
lignes[1] = malloc(sizeof(int) * (tailleSequence1+1));
colonnes[0] = malloc(sizeof(int) * (tailleSequence2+1));
colonnes[1] = malloc(sizeof(int) * (tailleSequence2+1));
int const COUTS_SUBSTITUTIONS[4][4] = { {0,TRANSVERSION,TRANSITION,TRANSVERSION},
{TRANSVERSION,0,TRANSVERSION,TRANSITION},
{TRANSITION,TRANSVERSION,0,TRANSVERSION},
{TRANSVERSION,TRANSITION,TRANSVERSION,0} };
if (lignes[0] == NULL || lignes[1] == NULL || colonnes[0] ==NULL || colonnes[1] ==NULL)
{
exit(EXIT_FAILURE);
}
for (int i = 0; i <= tailleSequence1; i++)
{
lignes[0][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i <= tailleSequence1; i++)
{
lignes[1][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i <= tailleSequence2; i++)
{
colonnes[0][i] = i * COUT_INSERTION_DELETION;
}
for (int i = 0; i <= tailleSequence2; i++)
{
colonnes[1][i] = i * COUT_INSERTION_DELETION;
}
int k = 0;
int tailleMini = min(tailleSequence1, tailleSequence2);
while (k < tailleMini)
{
lignes[k % 2][k + 1] = min(lignes[(k + 1) % 2][k] + COUTS_SUBSTITUTIONS[sequence1[k]][sequence2[k]], /* \ */
min(lignes[(k + 1) % 2][k + 1] + COUT_INSERTION_DELETION, /* | */
colonnes[(k + 1) % 2][k + 1] + COUT_INSERTION_DELETION)); /* - */
colonnes[k % 2][k + 1] = lignes[k % 2][k + 1];
for (int i = k + 1; i < tailleSequence1; i++)
{
lignes[k % 2][i + 1] = min(lignes[(k + 1) % 2][i] + COUTS_SUBSTITUTIONS[sequence1[i]][sequence2[k]], /* \ */
min(lignes[(k + 1) % 2][i + 1] + COUT_INSERTION_DELETION, /* | */
lignes[k % 2][i] + COUT_INSERTION_DELETION)); /* - */
}
for (int i = k + 1; i < tailleSequence2; i++)
{
colonnes[k % 2][i + 1] = min(colonnes[(k + 1) % 2][i] + COUTS_SUBSTITUTIONS[sequence1[k]][sequence2[i]], /* \ */
min(colonnes[k % 2][i] + COUT_INSERTION_DELETION, /* | */
colonnes[(k + 1) % 2][i + 1] + COUT_INSERTION_DELETION)); /* - */
}
k++;
}
int cout = lignes[(k + 1) % 2][tailleSequence1];
if (tailleSequence1 > tailleSequence2)
{
cout = colonnes[(k + 1) % 2][tailleSequence2];
}
free(lignes[0]);
free(lignes[1]);
free(colonnes[0]);
free(colonnes[1]);
return cout;
}
Maintenant, je devrais mieux m'en sortir. Je vais donc mettre le sujet en résolu, et si je ne trouve toujours pas la réponse aux questions 2 et 3, je ferai peut-être un nouveau sujet mais plus spécifique pour moins embrouiller les gens.
Merci, je vais y jeter un oeil, ça pourra toujours m'aider, même si j'ai fini par m'en sortir. Et pas d'inquiétude pour le Python, je parle Python presque couramment !
Problèmes de bioinformatique
× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
Le Tout est souvent plus grand que la somme de ses parties.
Le Tout est souvent plus grand que la somme de ses parties.
Le Tout est souvent plus grand que la somme de ses parties.