TP5 : Utilisations des programmes d'alignement en ligne de commande

L'objectif du TP est d'apprendre à installer puis utiliser les logiciels d'alignement en ligne de commande sur son ordinateur. Un des avantages des lignes de commande est la possibilités d'utililiser toutes les options pour un logiciel, ce qui n'est pas forcément possible par une interface web à cause du trop grand nombre d'options.
Ces logiciels peuvent être utiliser séparémment pour chaque analyse, mais il est également possible d'utiliser des scripts (shell) afin de pouvoir créer des chaines de traitement. Ces chaines de traitement permettent d'exécuter de manière automatique plusieurs fois le meme logiciels sur des données différentes et/ou une suite de logiciels pour effectuer un traitement complexe. Dans ce TP nous verrons comment créer de telles scripts.

Nous utiliserons les logiciels suivant : Une partie du TP a été réalisé grâce à des exercices fournis par le site web d'Olivier Elemento et grâce au livre Bioinformatique - Principes d'utilisation des outils (Jean-Loup Risler, Denis Tagu, 2010, Éditions Quae)

I. Installation des logiciels utilisés

Sous Linux (Ubuntu), de nombreux paquets logiciels existent déjà pour la bioinformatique (un paquet logiciel se réfère à un logiciel contenu dans un fichier d'archive qui, manipulé par un système de gestion de paquets logiciels, installe ou supprime le logiciel en question dans un ordinateur). Il vous suffit donc de les installer via un logiciel de gestionnaire de paquet (comme Synaptic, apt-get, la logithèque Ubuntu, ...). Vous trouverez ci-dessous les paquets à installer pour chaque logiciel.
Particularité de Muscle
Muscle est régulièrement mis à jour, ainsi Seaview peut ne pas disposer de la dernière version. Donc si vous souhaitez avoir la version la plus récente de Muscle, il faut l'installer à part.
Pour cela, il faut aller sur le site de Muscle, ici, puis dans 'Download', et télécharger la version qui correspond à votre environnement matériel. Le fichier téléchargé contient l'éxecutable (binaire) du logiciel Muscle.

Le format .tar.gz est un format de données archivées (tar) et compressées (gz=gzip). Pour extraire les données du fichier téléchargé, il faut utiliser la commande tar -zxvf fichier .
Cette commande permet d'extraire les données du fichier.
Si vous le souhaitez, installez la dernière version de Muscle (Muscle 3.8)
Flèche vers le haut

II. Exécution des logiciels séparément

II-1.Comment lancer / exécuter un logiciel

Pour lancer une commande ou éxécuter un programme, il suffit de taper son nom dans un shel (un terminal). Dans le cas d'un programme, si celui-ci ne se trouve pas dans le répertoire courant, il faudra que la variable d'environnement $PATH précise le répertoire dans lequel il se trouve. Sous Linux par défaut, les exécutables (binaires) sont dans les répertoires /bin ou /usr/bin ou /usr/local/bin et ces répertoires sont listés dans la variable d'environnement $PATH.
Afin de confirmer cela, vous pouvez voir la valeur de la variable d'environnement $PATH, en tapant la commande echo $PATH dans votre terminal.
Pour tous les programmes installés en utilisant un logiciel de gestionnaire de paquet, les exécutables sont automatiquement mis dans le répertoire /usr/bin.
Sinon, il y a deux possibilités :
  • vous pouvez mettre ou copier le fichier éxécutable souhaité dans le répertoire /usr/bin
  • vous pouvez ajouter un répertoire au $PATH, il suffit d'utiliser la commande export (par exemple, pour rajouter le répertoire /home/arigon/Bureau/ProgMuscle au $PATH : export PATH=$PATH:/home/arigon/Bureau/ProgMuscle)
Attention, le 2eme choix implique que votre éxécutable doit toujours rester au meme endroit pour qu'il puisse être éxécuter.
Pour tous les logiciels installés, vérifier que vous pouvez bien les lancer depuis votre terminal quelques soit le répertoire où vous êtes, et le cas échéant faire les modifications nécessaires (copie de l'éxécutable dans /usr/bin ou modification de la variable $PATH)

La plupart des logiciels installés propose deux modes lorsqu'on les exécute sur un terminal :
  • un mode interactif où les paramètres sont demandés à l'utilisateur au fur et à mesure (fichier en entrée, en sortie, options)
  • un mode en ligne de commande où tous les paramètres sont spécificiés aprés le nom du programme au moment de l'exécution
En tapant dans le terminal le nom d'un des programmes, vous lancez le mode interactif, s'il existe, sinon il en résulte une erreur.
Flèche vers le haut

II-2. Utilisation de Readseq

Il existe une multitude de formats, et passer de l'un a l'autre peut se faire grace au programme Readseq. Il s'utilise en ligne de commande de la manière suivante :
readseq fichierEntree -all -format=fmt -output=fichierSortie

fmt est une chaine de caractère précisant le format de sortie, à choisir. Beaucoup de formats sont possibles dont : GenBank|gb, EMBL|em, Pearson|Fasta|fa, Clustal, MSF, XML, ...

Pour avoir plus de précisions sur les options, vous pouvez taper readseq -help ou man readseq dans le terminal ou regarder les sites web suivant : aide readseq

Par exemple, pour obtenir un fichier au format FASTA à partir d'un fichier au format PHYLIP :
readseq sequences.phy -all -format=Pearson/Fasta -output=sequences.fasta

Aller chercher des séquences au format Genbank sur le site du NCBI et convertisser ses séquences au format FASTA à l'aide de Readseq (essayer avec une puis plusieurs séquences).
Readseq est également disponible en mode interactif en tapant la commande readseq dans le terminal.
Faites une conversion de format en utilisant le mode interactif.
Vous pouvez également l'utiliser au travers d'une interface web (ici).
Flèche vers le haut

II-3. Utilisation de BLAST et formatdb

BLAST travaille sur des bases de données pré-indéxées, c'est-à-dire traitées de manière à respecter une structure de données précise. L'indexation d'une base de données se fait au moyen de la commande formatdb (remplacée par makeblastdb pour les dernière version de Blast, fin 2012). Sur le site du NCBI, il est néanmoins possible de télécharger un certain nombre de bases de données pré-indexées.

Construction d'une banque BLAST

La première étape lors de l'utilisation de BLAST consiste à construire une base de séquences utilisable par BLAST. Le programme formatdb, fourni dans le package blast2 (ou makeblastdb en fonction de la version du package blast utilisé), permet de créer une telle base de données, à partir de séquences stockées au format FASTA.
Quelques options de formatdb:
-i nom du fichier en entrée (au format FASTA)
-p type des séquence (par défaut T pour protéique, sinon F pour nucléique)
-n nom de la base de données résultat (par défaut nom du fichier en entrée)
-o option pour parser les SeqID et créer des indexes (par défaut F, T si les séquences sont au format NCBI)
Pour avoir plus de précisions sur les options, vous pouvez taper formatdb -help ou man formatdb dans le terminal ou regarder le site web suivant : aide formatdb. (Information sur la commande makeblastdb ici)

Par exemple, supposons que l'on souhaite formater le fichier ecoli.nt contenant un grand ensemble de séquences d'E. coli. La commande nécéssaire est :
formatdb -i ecoli.nt -p F -n maBdBlastEcoli
A partir des séquences de génomes complets d'Escherichia coli (400 séquences) que vous trouverez ici, construire une banque Blast. Ces séquences sont au format NCBI, donc vous pouvez utiliser l'option -o T : essayer avec et sans cette option.
Remarque :
Si l'option -o T a été utilisée lors de la construction de la banque BLAST, il est possible de récupérer une séquence stockée dans une base de données BLAST (créée avec formatdb ou makeblastdb). Pour cela, il faut utiliser la commnade fastacmd (remplacée par blastdbcmd pour les dernière version de Blast, fin 2012). Le package blast2 fournit le programme fastacmd (ou blastdbcmd en fonction de la version du package blast utilisé):
fastacmd -d base_de_donnée -s id_sequence
Pour avoir plus de précisions sur les options, vous pouvez taper man fastacmd dans le terminal ou regarder le site web suivant : aide fastacmd. (Information sur la commande blastdbcmd ici)
Essayer de retrouver une des séquences de la banque en utilisant la commande fastacmd, par exemple la séquence dont le numéro d'accession est AE000432.
Exécution d'un BLAST

Il est ensuite possible de faire des requètes sur la base ainsi créée. Pour cela, il est possible d'tiliser le programme blastall fourni dans le package blast2.
Quelques options :
-p type de BLAST (blastn, blastp, blastx, tblastn, tblastx)
-d noms des banques BLAST à utiliser (sans les extensions de fichier)
-i fichier en entrée au format fasta
-o fichier de sortie (par défaut tdout)
-e sueil de la E-value (par défaut=10.0)
Pour avoir plus de précisions sur les options, vous pouvez taper blastall ou man formatdb dans le terminal ou regarder les sites web suivants : aide blastall 1 ou aide blastall 2.

Par exemple, supposons que notre séquence inconnue soit contenue dans le fichier test.txt, la commande permettant de la comparer aux séquences de la banque maBdBlastEcoli est :
blastall -p blastn -d maBdBlastEcoli -i test.txt -o test.out

Supposons maintenent que nous recherchons les homologues d'une séquence nucléotidique contenue dans le fichier K03455seqref.fas (format FASTA), avec un seuil égal à 0.15 (E-value < 0.15), dans la base GenBank (supposons celle-ci installée en local et nommée nr). Pour cela, on souhaite afficher aux environs de 5000 résultats (-v), sans aucun alignements (-b 0). La commande nécéssaire est :
blastall -p blastn -d nr -i K03455seqref.fas -e 0.15 -v 5000 -b 0

Notons que blastall ne permet pas de choisir l'espèce (l'interface Web du NCBI le permet, au travers d'une requète Entrez).
Aller chercher une séquence de génome complet sur le site du NCBI en utilisant le mot de clé "ecoli" afin de prendre une séquence proche de celles utilisées pour contruire la banque BLAST, puis lancer un BLAST afin de comparer cette séquence avec les séquences de la banque précédemment créée.
Flèche vers le haut

II-4. Utilisation de Strecher / Matcher

Pour visualiser les alignements obtenus avec Stretcher et Matcher, vous pouvez utiliser Seaview.

Strecher

Stretcher s'utilise en ligne de commande avec 3 paramètres obligatoires : les noms des 2 fichiers des séquences à comparer et le nom du fichier de sortie.
Si le nom du fichier de sortie n'est pas précisé, Stretcher le demandera à l'utilsiateur en en proposant un par défaut.

La commande permettant de comparer la séquence du fichier s1 et la séquence du fichier s2 en mettant le résultat dans le fichier nommé fichierSortie est :
stretcher [-asequence] s1 [-bsequence] s2 [-output] fichierSortie
([-paramètre] = pas nécessaire d'indiquer le nom du paramètre)

Pour avoir plus de précisions sur les autres options, vous pouvez taper stretcher -help ou man stretcher dans le terminal ou regarder le site web suivant : aide stretcher.
En utilisant la commande stretcher, comparer ces 2 séquences nucléiques : xl-actin-DNA.fa et xl-actin-mRNA.fa (cf. TP1)
Stretcher peut également être utilisé en mode intéractif.
Faites à nouveau la comparaison en utilisant le mode interactif.
Matcher

Matcher s'utilise de la même manière que Stretcher.
Pour avoir plus de précisions sur les paramètres et options, vous pouvez taper matcher -help ou man matcher dans le terminal ou regarder le site web suivant : aide matcher.
En utilisant la commande matcher, puis en mode intéractif, comparer ces 2 séquences protéiques, les facteurs de transcription krox 24 et sp1 : krox24.seq et sp1.seq (cf. TP1)
Flèche vers le haut

II-5. Utilisation de ClustalW

Pour visualiser les alignements obtenus avec ClustalW, vous pouvez utiliser Seaview.

ClustalW propose un mode interactif et un mode en ligne de commande.

Mode interactif

Dans le mode interactif, un menu avec sous-menu apparait et permet d'indiquer le(s) fichier(s) de séquences à aligner, de modifier les paramètres d'alignement et de lancer la constrution de l'alignement souhaité.
Aligner les séquences protéiques partageant le motif bZip en utilisant les paramètres par défaut : fichier au format fasta (cf. TP4)
Essayer maintenant de ré-aligner les memes séquences en modifiant les paramètres : méthode d'alignement par paire (précise/approximée), paramètres de la méthode d'alignement par paire (pénalité de gap, matrice,...), paramètres de la méthode d'alignement multiple (pénalité de gap, matrice,...), format de sortie, ...
ClustalW permet d'ajouter une séquence à un alignement déjà connu ou d'aligner 2 alignements déjà connus en utilisant l'alignement de profile.
Ainsi, essayer d'aligner cette séquence avec l'alignement précédemment obtenu (vous pouvez modifier les paramètres par défaut si vous le souhaitez).
Vous pouvez également récupérer d'autres séquences partageant le meme motif (via le site PROSITE et NCBI). Une fois que vous avez ces séquences, il y a 2 méthodes :
  • soit les aligner entre elles puis aligner cette alignement avec le précédent
  • soit demander à ClustalW d'aligner les séquences (récupérées non alignées) avec l'alignement précédent
Comparer les alignements obtenus avec les 2 méthodes.

Mode ligne de commande

En mode ligne de commande, un seul paramètre est obligatoire : le nom du fichier de séquence à aligner La commande permettant d'aligner les séquences d'un fichier avec les paramètres par défaut est :
clustalw -INFILE=monFichier.fasta

Tous les paramètres/options modifiables en mode interactif sont modifiables en ligne de commande.
Quelques exemples d'options :
  • Réglage général rapide : option -QUICKTREE
  • Réglage encore plus rapide (plus de 1000 séquences): options -QUICKTREE -KTUPLE=2 -TOPDIAGS=3 -WINDOW=3
  • Options de fichier de sortie : -OUTFILE=nomFichierSortie -OUTPUT=formatSortie (GCG, GDE, PHYLIP, PIR or NEXUS)
  • Option -MAXDIV : à chaque itération la structure de guidagedétermine les groupes de séqunces à aligner; Cependant, si les groupes partagent moins de 30% d'identité alors l'alignement est reporté (mot-clé DELAYED utilisé lors de l'exécution). Les alignements reportés sont traités à la fin de l'agorithme. Pour changer cette valeur on utilise l'otpion -MAXDIV=%identité
  • Pour ajouter une séquence à un alignement déjà connu ou aligner 2 alignements déjà connus (alignement de profils) : options -PROFILE1=file.ext -PROFILE2=file.ext
Pour avoir plus de précisions sur les paramètres et options, vous pouvez taper clustalw -help ou man clustalw dans le terminal ou regarder le site web suivant : aide clustalw.
Refaire les alignements construits en mode intéractif en utilisant cette fois-ci le mode en ligne de commande.
Remarque :
CulstalW propose de convertir des sequences en utilisant l'option -CONVERT
Essayer de convertir une séquence d'un format vers un autre : utiliser une(des) séquence(s) vu dans un exercice précédent ou aller chercher des séquences sur le site du NCBI
Flèche vers le haut

II-6. Utilisation de Muscle

Pour visualiser les alignements obtenus avec Muscle, vous pouvez utiliser Seaview.

Muscle n'a pas de mode interactif, il faut l'utiliser en mode ligne de commande uniquement.
Le seul paramètre obligatoire est le nom du fichier de séquence à aligner. Cependant si le nom de fichier de sortie n'est pas spécifié, Muscle affiche l'alignement résultat dans le terminal. Il est donc préférable de spécifier le fichier de sortie. La commande permettant d'aligner les séquences d'un fichier avec les paramètres par défaut est :
muscle -in monFichierEntree.fasta -out monFichierSortie.fasta

L'algorithme de Muscle est constitué de 3 étapes : 1 étape d'alignement progressif, 1 étape d'amélioration de l'alignement progressif et 1 étape de raffinement itératif. Il est possible de n'exécuter que les 2 premières étapes (version muscle-prog) ou que la 1ere étape (version muscle-fast) :
  • version rapide (muscle-prog) = étape 1 et 2 : dans le cas où il y a beaucoup de séquences ou des séquences longues à aligner
    option -maxiters 2
  • version encore plus rapide (muscle-fast) = étape 1 : dans le cas où il y a beaucoup de séquences à aligner et que ces séquences sont proches
    proteique : options -maxiters 1 -diags -sv -distance1 kbit20_3
    nucléique : options -maxiters 1 -diags
Pour ajouter une séquence à un alignement déjà connu ou aligner 2 alignements déjà connus (alignement de profils) : options -profile -in1 existing_aln.fasta -in2 new_seq.fasta ou -profile -in1 file1.fasta -in2 file2.fasta

D'autres options sont paramétrables, tel que les pénalités de gap, les matrices de substitution utilisées ,...
Pour avoir plus de précisions sur les paramètres et options, vous pouvez taper muscle ou muscle -help ou man muscle dans le terminal ou regarder le site web suivant : aide muscle.
Refaire les alignements construits avec ClustalW en utilisant les différentes versions de Muscle et comparer les résultats.
Flèche vers le haut

III. Lancement des logiciels via un script shell

Un script shell permet d'automatiser une série d'opérations. Il se présente sous la forme d'un fichier contenant une ou plusieurs commandes qui seront exécutées de manière séquentielle.

Le fichier d'un script shell a comme extension .sh .

La première chose à faire dans un script shell est d'indiquer quel shell est utilisé. En effet, la syntaxe du langage change un peu selon si on utilise "sh", "bash", "ksh", etc. Par défaut, le shell utilisé par Ubuntu est le bash (Bourne Again Shell). Il faut donc commencer le fichier par :
#!/bin/bash.
On écrit ensuite les différentes instructions à exécuter.

Pour exécuter un script à partir du dossier où est placé le script, on peut : Pour exécuter un script à partir de n'importe quel dossier, il faut rendre exécutable le fichier puis modifier le PATH (voir section II-1)

Voila quelques sites pour vous aider à créer votre script shell (du moins détaillé au plus détaillé) :
Tutoriel script shell
Introduction aux scripts shell
Guide avancé d'écriture des scripts Bash
Il est possible d'appeler un script shell avec des paramètres/arguments (exemple : monScript param1 param2). Plusieurs variables existent pour récupérer dans le script les informations concernant les paramètres/arguments : Remarque :
Il faut savoir qu'il existe en shell des commandes trés pratiques qui peuvent être utilisées dans des scripts ou directement dans le terminal : Vous pouvez également utiliser le langage perl pour parser des fichiers tel que des sorties Blast, ....

De nombreux sites web fournissent des aides pour écrire des shell, construire des expessions régulières, utiliser les différentes commande UNIX dont sed, awk, ... : exemples de site web d'aide 1, site web d'aide 2

Premier script : traitement Blast
Créer un premier script qui prend en paramètre 2 noms de fichiers. Ce script devra convertir l'ensemble des séquences du premier fichier en un fichier de séquences au format fasta, puis créer un banque blast à partir de ces séquences, puis lancer un blast de la séquence requete contenue dans le 2eme fichier (en paramètre) sur la banque créée précédemment. A partir du résultat de Blast, récupérer les noms des séquences correspondant aux 10 meilleurs hits Blast (pour cela, vous pouvez utiliser la commande awk dans votre script) puis aligner ces séquences avec la séquence requete (contenue dans le 2eme fichier en paramètre).
Exécuter votre script avec en paramètre les 2 fichiers suivants : fichier des sequences de génomes complets d'Escherichia coli (400 séquences) pour la banque Blast au format GenBank et fichier d'une sequence de Salmonella à comparer.
Améliorations du premier script :
Second script : alignements à la chaine
Créer un second script qui prend en paramètre un nom de répertoire et qui construit les alignements multiples des séquences contenus dans chaque fichier de ce répertoire. Il y aura donc en résultat autant de fichiers d'alignement que de fihciers de séquences présents dans le répertoire.
Créer un répertoire dans lequel vous mettrez ces fichiers de séquences protéiques au format fasta (fichier de séquences partageant le motif bZip, et des fichiers de séquences de la banque de données HOBACGEN du PBIL fichier2, fichier3, fichier4) et exécuter votre script sur ce répertoire afin de calculer les alignements multiples pour chacun des fichiers de ce répertoire.

Exemple de scripts :