silico.biotoul.fr
 

InfoBio TD bioperl

From silico.biotoul.fr

Jump to: navigation, search

L'objet du TD est de se familiariser avec le langage Perl et le module BioPerl.

Perl (pour Pratictal Extraction and Reporting Language) est un langage de programmation/script. Il met en oeuvre le paradigme de la programmation itérative.

Contents

Types de données

  • nombre entier
  • nombre réels
  • chaîne de caractères. Comme dans un script shell, dans une chaîne délimitée par des simples quote ' , les variables ne sont pas remplacées par leur valeur contrairement aux chaînes délimitées par des double quote ".
  • descripteur de fichier
  • référence (pointeur)

Structures de données

  • les scalaires (SCALAR) utilisant la notation $. Par exemple $filename = 'sequences.fasta' contiendra typiquement le nom (ou chemin) d'un fichier.
  • les listes ou tableaux (ARRAY) utilisant la notation @. Par exemple @filenames = ['comE_sequences.fasta', 'comD_sequences.fasta']. Pour accéder à un élément, on utilise un indice : $filenames[3] pour le quatrième élément.
  • les tableaux associatifs (HASH) avec la notation %. Par exemple %gene_annotations = { 'comD'=>'sensor', 'comE'=>'response regulator' } ou bien %gene_annotations = [ 'comD', 'sensor', 'comE', 'response regulator']. Pour accéder à un élément : $gene_annotation{'comD'}.

Pour plus de détails, cf. http://perldoc.perl.org/perldata.html

Quelques fonctions essentielles

L'ensemble des fonctions disponibles : http://perldoc.perl.org/perlfunc.html#Perl-Functions-by-Category

Opérations sur les listes

Les listes peuvent se manipuler comme des piles (LIFO), files (FIFO) ou autres grâce aux opérations :

  • $elem = pop @list; récupère (et supprime de la liste) le dernier élément
  • push @list, $elem; opération duale de pop : ajoute l'élément en fin de liste
  • $elem = shift @list; récupère (et supprime) le premier élément de la liste
  • unshift @list, $elem; dévinez !
  • $nb_elem = scalar @list; interprète la liste dans un contexte scalaire ce qui correspond à récupérer le nombre d'éléments.
  • sort @list trie la liste par ordre lexicographique. Si l'on veut trier par ordre numérique croissant (11>2 alors que "11"<"2") il faut spécifier la routine de tri : sort {$a <=> $b} @list.
  • itérations sur chaque éléments: foreach my $elm (@list) { print $elm; }

Opérations sur les tableaux associatifs

  • @list = keys %hash;
  • les clés triées : @list = sort keys %hash;
  • les clés triées en fonction des valeurs en ordre numérique décroissant : @list = sort { $hash{$b} <=> $hash{$a} } keys %hash;

cf. http://perldoc.perl.org/functions/sort.html pour plus d'exemples de tri.

  • les valeurs : @list = values %hash;
  • itérations sur chaque élément : foreach my $key (keys %hash) { print $hash{$key}: }
  • test de l'existence d'une clé : defined $hash{$key}

Quelques particularités

Certaines variables sont prédéfinies :

  • $_ et @_ correspondent à la variable par défaut ou variable courante.

Exemple : 

foreach (@list) {
   print $_;
}
  • <> correspond au fichier passé en paramètre au script.

Exemple : lacement du script suivant avec la commande: sh> perl le_script.pl un_fichier_texte.txt 

#!/usr/bin/perl
 
while (<>) {
   chomp;
   print $. . $_;
}

Commentaires : la boucle while itère sur chaque ligne du fichier. La fonction chomp supprime le caractère de fin de ligne (ici, le fait de l'appeler sans paramètre l'applique à la variable par défaut $_ qui est une ligne du fichier texte). La variable spéciale $. correspond au numéro de ligne du fichier. Le '.' qui suit correspond à la concaténation de chaînes de caractères.

Manipulation de fichiers

  • Lecture d'un fichier 
open(F, '<filename');
while (<F>) {
  ...
}
close(F);
  • Ecriture dans un fichier 
open(F,'>output_file');
print F 'some content';
close(F);

Branchements conditionnels 

if ($condition) {
 ...
}
# commentaire : on peut mettre n'importe quoi dans la condition, elle sera évaluée à vrai du moment que la valeur est différent de 0 ou null
 
# avec un else
if (open(F,"$filename")) {
 ...
}
else {
 ...
}
 
# avec un else if
if ($condition) {
 ...
}
elsif ($autre_condition} {
 ...
}
# le if en fin de statement
print $some_debuging if $verbose;
 
# l'inverse : le ''à moins que'' (unless)
print "Success" unless $failure;

Boucles 

# WHILE
my $i=0;
while ($i<$max_iteration && $do_not_stop) {
  ...
   $i++;
}
 
# FOR
for (my $i=0 ; $i < scalar @list ; $i++) {
  $list[$i] *= 2;
}

Fonctions 

# définition
sub nom_fonction {
 # corps fonction
}
 
# on peut ensuite l'appeler
nom_fonction;
 
# passage des paramètres
# le passage se fait toujours sous forme d'une liste
####################################################
sub fonction_2 {
  my @params = @_;
}
 
# appel
fonction_2($param1, $param2);
 
 
sub fonction_3 {
  my ($param_1, $param_2) = @_; # on s'intéresse seulement aux 2 premiers
}
 
# appel
fonction_3(@params);
 
# si l'on souhaite passer plusieurs listes ou un hash, on passe par des références (pointeurs)
sub fonction_4 {
  my $h_ref = shift; # récupération d'une référence sur un hash (1er paramètre)
  my $array_ref = shift; # récupération du 2ème paramètre qui est une référence sur un tableau
 
  # pour déréférencer :
  my %h = %$h_ref;
  my @a = @$array_ref;
 
  # sinon on y accède avec un cast
  my $troisieme_element = @{$array_ref}[2];
 
}
 
# appel avec passage par référence
fonction_4(\%hash, \@list);

Il existe encore pas mal de possibilités mais on se contentera de celles-ci pour aujourd'hui.

Expression régulières

Pour terminer cette courte introduction, nous allons voir que perl est particulièrement adapté à l'exploitation des expressions régulières. Les possibilités sont assez nombreuses et nous ne verrons qu'une infime partie. Pour plus de détails, vous pouvez vous référer au tutoriel http://perldoc.perl.org/perlretut.html 

# test de l'occurrence d'une expression régulière dans une chaîne
while (<>) {
   next if $_ =~ /^#/; # skip comments (lines beginning with '#')
   ...
}
 
# récupération du contenu des occurrences d'une expression régulière
my ($heure, $minute, $seconde) = $time =~ /(\d+):(\d+):(\d+)/; # match "...dernier accès 9:40:20 pour le fichier ..."

Ainsi, une expression régulière est délimitée par des '/'. Les parenthèses permettent de séparer les différents match pour les affecter à des variables. On peut matcher des chiffres, lettres, espaces avec :

  • \d : digit
  • \s : white space
  • \S : tout sauf white space
  • \w : word character
  • [a-z] : minuscule
  • [a-zA-Z] : lettre
  • ^ : début de chaîne
  • $ : fin de chaîne
  • ... cf. le tutoriel plus haut

Pour le nombre d'occurrences :

  •  ? : 0 ou 1 fois
  • * : 0 ou plus
  • + : 1 ou plus
  • {4} : 4 fois
  • {3,} : 3 ou plus
  • {2,5} : au moins 2 et au maximum 5

BioPerl

BioPerl est un module perl dédié à la biologie. Il fournit notamment un ensemble de routines pour lire/écrire des fichiers de séquences dans différents formats (EMBL, GenBank, FASTA, ...), de récupérer une ou des séquences dans une banque de données, de lancer une analyse (par exemple une recherche par similarité de séquences) locale ou distante, d'encapsuler un appel à un programme de la suite EMBOSS.

Pour l'utiliser, il faut sélectionner le(s) module(s) qui nous intéressent. Par exemple, pour lire un fichier au format EMBL : 

#!/usr/bin/perl
 
use strict;
use Bio::SeqIO;
use Data::Dumper;
 
my $filename = shift;
 
my $seqin = Bio::SeqIO->new(
	-file   => "<$filename",
	#~ -format => $informat,
);
 
my $seq = $seqin->next_seq;
print Dumper($seq);

Exercice 1 : Récupérez la fiche de CAA65843 et faire un script qui prend en paramètre un fichier dans ce format et qui extrait pour chaque domaine son nom et sa position (début et fin). Les faire afficher. Voici la sortie attendue: 

[user@host TP-bioperl]$ ./extract_domains.pl CBEL.gp 
>CAA65843.CBM_1.23-55
PSFGNCGSDAAGVSCCQSTQYCQPWNANYYQCL
>CAA65843.APPLE_Factor_XI_like.57-129
LPAKCAQQFPNVDFNGDDIQTIYGIQPGECCTRCSETAGCKAYTFVNSNPGQPACYLKSGTGTRTPSVGAVSG
>CAA65843.CBM_1.162-191
YGSCGSSNGATCCPSGYYCQPWNDSFYQCI
>CAA65843.APPLE_Factor_XI_like.193-265
PPAKCSKQLTDKDYYGNDIKTVYVSLPSLCCDACASTAGCKAYTYINNNPGQPVCYLKSAAGTATTKIGAVSG

Pour cela vous aurez besoin des fonctions fournies par la librairie bioperl : 

# Récupérer l'identifiant de la séquence
my $accession = $seq->accession();
 
# Récupérer la liste de Features contenues dans la fiche
my @features = $seq->get_SeqFeatures();
 
# Type de feature (source, Protein, sig_peptide, Region, ...)
my $feature_type = $feature->primary_tag();
 
# Dans une feature, récupérer la liste de valeurs pour une sous-feature (region_name, db_xref, note, ...)
my @tag_values = $feature->get_tag_values('region_name');
 
# Récupérer le début et la fin (positions) de la feature sur la séquence
my ($start, $end) = ( $feature->start(), $feature->end() );
 
# Extraire une région de la séquence
my $subseq = $seq->subseq( $start, $end );


Exercice 2 : Utiliser le script en ligne à l'adresse http://prosite.expasy.org/mydomains/ sur le site de PROSITE pour récupérer un dessin de la protéine et de ses domaines. Il s'agira de modifier le script précédent pour qu'il construise l'URL permettant de tracer le schéma et qu'il récupère l'image à cette adresse.

Pour connaître la longueur de la séquence : 

my $length = $seq->length();

Pour récupérer un fichier à une URL (par exemple http://silico.biotoul.fr) : 

bash> wget "http://silico.biotoul.fr" -O nom_du_fichier_sauvegardé

Pour lancer une commande externe depuis un script perl : 

my $url = "http://prosite.expasy.org/cgi-bin/prosite/PSImage.cgi?hscale=1.0&type=1&len=268&hit=23,55,CBM_1,CBM_1+57,129,APPLE_Fact,APPLE_Fact+162,191,CBM_1,CBM_1+193,265,APPLE_Fact,APPLE_Fact";
my $output = "CAA65843.domains.png";
# Exécution de wget
`wget "$url -O $output`;
Retrieved from "http://silico.biotoul.fr/p/InfoBio_TD_bioperl"