I. Introduction▲
Le langage « R », distribué gratuitement, est très utilisé dans le domaine des statistiques et de l'analyse des données, notamment grâce à une bibliothèque large et complète. Perl est un langage de script qui permet de manipuler aisément processus et fichiers textes. Les mathématiciens et bio-informaticiens l'apprécient également pour sa gestion des expressions régulières. Il peut donc être très intéressant de conjuguer ces deux langages, afin de profiter de leurs avantages respectifs.
Statistics::R (Controls the R interpreter through Perl) est l'un des modules qui permet de faire interagir les langages R et Perl. La mise à jour de cet article utilise la version 0.32 du module.
À la première publication de ce tutoriel, la version 0.09 du module était détailllée. Le nouvel auteur qui l'a pris en charge l'a complètement réécrit.
Il est ainsi possible dans un script Perl de faire appel à R et à ses bibliothèques, de récupérer les objets ainsi créés et de les exploiter dans le pipeline d'analyse Perl.
Les étapes principales d'utilisation du module sont les suivantes :
- déclaration de l'objet Perl $R ;
- ouverture du pont entre R et Perl ;
- chargement des bibliothèques et/ou création des objets R et/ou des calculs statistiques ;
- fermeture du pont.
On prendra ici l'exemple de l'analyse de données transcriptomiques en utilisant la bibliothèque FactoMineRFactoMineR. L'un des scripts ci-dessous pourra donc être utile aux bio-informaticiens débutants !
II. Contexte▲
L'exemple est basé sur l'étude de données transcriptomiques. Le fichier de soumission contient les données d'expression de 20562 gènes pour deux conditions A et B, chacune représentée par quatre échantillons (de A1 à A4 et de B1 à B4).
Ce tableau de données contient donc huit échantillons (huit individus ou huit colonnes) et 20652 gènes (variables quantitatives). Une variable qualitative sera ajoutée pour décrire chaque échantillon.
L'objectif de cette étude est de résumer et décrire ce jeu de données et d'identifier des relations avec les conditions A et B.
Classiquement, le fichier de données représente les gènes en lignes et les individus en colonnes. En effet, en transcriptomique, le nombre de gènes est très largement supérieur au nombre d'échantillons. Ce format permet ainsi d'être supporté par la plupart des tableurs et des éditeurs de texte.
Pour l'utilisation de la fonction ACP (analyse en composantes principales) disponible dans la bibliothèque FactoMineR, ce tableau sera donc importé et transposé.
III. Explication par l'exemple▲
III-A. Méthodes obsolètes▲
Dans le but de garder une trace des anciennes options dans le cas où certains lecteurs utiliseraient l'ancienne version du module, voici une liste des méthodes devenues obsolètes :
- log_dir : espace de travail. Il représentait le répertoire de travail où le pont était créé entre les deux langages. R et Perl et devaient avoir les droits de lecture et d'écriture.
- startR : ouverture du pont entre R et Perl ;
- start_sharedR : ouverture du pont ou utilisation d'une communication existante ;
- stopR : fermeture du pont entre R et Perl ;
- Rbin : retour du chemin de l'exécutable R ;
- send ($CMD) : envoi des commandes à exécuter par R ;
- is_started : TRUE si l'interpréteur R est démarré ;
- clean_up : nettoyage de l'environnement en supprimant tous les objets ;
- error : retour du dernier message d'erreur.
III-B. Retour à notre exemple▲
Pour créer un lien entre votre programme Perl et R, vous devez faire appel au constructeur new du module. Ce dernier comprend deux options :
- r_bin : afin de préciser le chemin complet du répertoire bin de R. Il n'est utile de spécifier ce chemin que si le module ne parvient pas à le trouver automatiquement ;
- shared : permet de créer un pont partagé par plusieurs instances de R que vous pouvez créer dans votre programme.
Voici un exemple de code :
use Statistics::R;
# Création du pont entre Perl et R et démarrage de R
my $R = Statistics::R->new();
…use Statistics::R;
# Création du pont entre Perl et R et démarrage de deux instances de R
# pouvant communiquer entre elles.
my $R1 = Statistics::R->new( shared => 1);
my $R2 = Statistics::R->new( shared => 1);Supposons que l'installation de « R » soit faite dans un répertoire non habituel (sous Windows ou Linux/Mac).
use Statistics::R;
# Création du pont entre Perl et R et démarrage de R
my $R = Statistics::R->new( r_bin => '/home/logiciel/R/bin' );Remarque : il est bien sûr nécessaire que R soit préalablement installé !
Pour l'envoi des commandes R, on utilise la méthode run. Cette méthode prend en argument un scalaire comprenant la ou les commandes R ou bien un tableau pouvant contenir une liste de commandes.
Les commandes R sont susceptibles de contenir des apostrophes, des guillemets, des dollars… : des caractères ayant une signification en Perl. Il faut ainsi faire très attention à l'interpolation. De ce fait, nous utiliserons les opérateurs q et qq (apostrophe et guillemet) permettant d'avoir une meilleure lisibilité du code.
Exemples :
# Pas d'interpolation, j'affecte la valeur 2 à x.
$R->run(q`x<-2`);
# Interpolation, j'affecte la valeur 4 à y. Cette valeur est issue d'une variable Perl.
my $valeur = 4;
$R->run(qq`y<-$valeur`);Voici différentes manières de faire appel à plusieurs commandes R dans une seule méthode run (vous pourriez également utiliser plusieurs méthodes run) :
$R->run(
q`x<-2`,
qq`y<-$valeur`,
);
# OU
my $commandeR = <<"CMD";
x<-2
y<-$valeur
CMD
$R->run($commandeR);
# OU
$R->run(<<"CMD";
x<-2
y<-$valeur
CMD
);N.B. L'écriture <<"CMD" implique une interpolation du fait des guillemets. Si vous ne souhaitez pas interpoler votre commande, il faudra écrire de la manière suivante : <<'CMD'.
Vous avez la possibilité d'appeler plusieurs commandes R présentes dans un fichier séparé via la méthode run_from_file().
my $fichier_commandes = 'C:/temp/Data.r';
$R->run_from_file($fichier_commandes);
print $R->get('c');a <- 2
b <- 5
c <- a * bOn obtient : 10.
Les méthodes run et run_from_file retournent ce qui s'affiche dans la console R si vous aviez fait un print dans le code R. C'est toujours utile si vous souhaitez exploiter cette sortie.
- Lecture des données R dans le script Perl :
#!/usr/bin/perl
use warnings;
use strict;
use Statistics::R;
# Communication Perl et R
my $R = Statistics::R->new();
# Affiche [1] 3
my $resultat = $R->run(
q`x<-1`,
q`y<-2`,
q`sum<-x+y`,
q`print(sum)`,
);
my $x = $R->get('x');
my $y = $R->get('y');
my $somme = $R->get('sum');
# Affiche
# $resultat = [1] 3
# $x = 1
# $y = 2
# $somme = 3
print "\$resultat = $resultat\n";
print "\$x = $x\n";
print "\$y = $y\n";
print "\$somme = $somme\n";Le module met la méthode get à disposition pour lire les variables scalaires et les vecteurs R. De la même façon, il existe une méthode set pour définir une variable scalaire ou un vecteur R.
N.B. Vous ne pouvez extraire que des variables scalaires et des vecteurs. Si vous cherchez à récupérer une matrice ou un tableau (data.frame), vous n'obtiendrez dans Perl qu'une liste « aplatie » et non un tableau à plusieurs dimensions (version 0.32 du module).
#!/usr/bin/perl
use warnings;
use strict;
use Statistics::R;
use Data::Dumper;
# Communication Perl et R
my $R = Statistics::R->new();
# Affiche [1] 3
my $resultat = $R->run(
q`n = c(1,2,3)`,
q`m = c("a","b","c")`,
q`df = data.frame(n, m)`,
);
my $x = $R->get('n');
my $y = $R->get('m');
my $df = $R->get('df');
print Dumper $x;
print Dumper $y;
print Dumper $df;$VAR1 = [
'1',
'2',
'3'
];
$VAR1 = [
'a',
'b',
'c'
];
$VAR1 = [
'n',
'm',
'1',
'1',
'a',
'2',
'2',
'b',
'3',
'3',
'c'
];On peut constater que les variables $x, $y et $df sont des références de tableaux. Autant $x et $y sont correctes pour notre exemple, autant, $df n'est pas exploitable. Il ne contient qu'une liste « aplatie », un tableau à une seule dimension. Cela signifie qu'il n'est pas possible d'extraire un data.frame dans une variable Perl via la méthode get. Espérons que le module évolue afin de nous fournir à l'avenir un tableau à plusieurs dimensions !
- Création d'une image :
#!/usr/bin/perl
use warnings;
use strict;
use Statistics::R;
# Communication Perl et R
my $R = Statistics::R->new();
# Création d'une image
my $image = 'C:/temp/test.png';
my $commandeR = <<"CMD";
png(file ="$image", bg = "transparent",)
v1<-c(1,2,3)
v2<-c(4,5,6)
plot(v1,v2,type="l",main="exemple de plot")
dev.off()
CMD
$R->run($commandeR);- chargement d'un fichier de données (« Data.txt ») qui contient les informations suivantes :
| Data | Col1 | Col2 | Col3 |
| Row1 | A | 1 | z |
| Row2 | B | 2 | y |
| Row3 | C | 3 | x |
| Row4 | D | 4 | w |
#!/usr/bin/perl
use warnings;
use strict;
use Statistics::R;
use Data::Dumper;
# Communication Perl et R
my $R = Statistics::R->new();
# Chargement d'un fichier
my $fichier = 'C:/temp/Data.txt';
my $commandeR = <<"CMD";
data<-read.table("$fichier", header=T, row.names=1, sep="\t")
CMD
$R->run($commandeR);
print $R->run('data'),"\n";
print Dumper $R->get('data'); Col1 Col2 Col3
Row1 A 1 z
Row2 B 2 y
Row3 C 3 x
Row4 D 4 w
$VAR1 = [
'Col1',
'Col2',
'Col3',
'Row1',
'A',
'1',
'z',
'Row2',
'B',
'2',
'y',
'Row3',
'C',
'3',
'x',
'Row4',
'D',
'4',
'w'
];Attention : Liste « aplatie » !
Chargement d'une bibliothèque :
# Chargement d'une bibliothèque
$R->run(q`library(FactoMineR)`);- Sauvegarde des objets :
# Sauvegarde des objets
my $r_save = 'C:/Documents and Settings/stoyak/Article_Dvp/Save.RData';
$R->run(qq`save(v1,v2,data, file="$r_save")`);- Suppression des objets :
# Suppression des objets
$R->run(
q`rm(v1,v2,data)`,
q`ls<-ls()`
);- Chargement du fichier .Rdata
# Chargement du fichier .RData
print $R->run(
qq`load("$r_save")`,
q`ls<-ls()`
);- Fermeture du pont
La fonction stop permet d'arrêter une instance R, mais si vous ne le faites pas, Perl s'en charge à la fin du programme, son utilisation n'est donc pas obligatoire ($R->stop();).
IV. Exemples▲
IV-A. Exemple de données transcriptomiques▲
Pour les bio-informaticiens, cet exemple consiste en l'étude de données transcriptomiques, en utilisant la méthode ACP de la bibliothèque FactoMineR.
La procédure ACP ci-dessous consiste en ces quelques étapes :
- déclaration des variables Perl et de l'objet Statistics::R ;
- ouverture du pont entre R et Perl ;
- chargement de la bibliothèque FactoMineR qu'il faudra installer au préalable dans une console R (install.packages("FactoMineR")) ;
- import et transposition du fichier de données ;
- ajout de la variable qualitative « condition » ;
- appel de la fonction PCA et création du fichier avec les graphes sur les variables et les individus ;
- calcul du pourcentage d'inertie pour les quatre premiers axes ;
- description des axes ;
- lecture des objets créés par R dans Perl ;
- sauvegarde des objets dans un fichier « .Rdata » ;
- fermeture du pont entre R et Perl.
Voici le script complet :
#!/usr/bin/perl
use strict;
use Carp;
use warnings;
use File::Basename;
use Statistics::R;
#===========================================
# Auteur : Stoyak - Developpez.com
# But : Statistics-R demo
# Date : 17/07/2014
#===========================================
#===========================================
# Définition des paramètres
#===========================================
# Fichier de données
my $file = 'C:/temp/Data.txt';
my %parameters = (
'file' => $file,
);
#===========================================
# Appel de la procédure ACP
#===========================================
ACP( \%parameters );
#===========================================
# Procédure ACP
#===========================================
sub ACP {
my ($ref_params) = @_;
my ( $file_name, $dir_name, $extension ) = fileparse( $ref_params->{file}, qr/\.[^.]*/ );
# Fichiers de résultats ('.txt' et '.pdf')
my $pdf_acp = $dir_name . '/PCA.pdf';
my $result_file = $dir_name . '/Result.txt';
open my $fh, '>', $result_file || die "Impossible d'écrire dans le fichier $result_file\n";
# Déclaration de l'objet
my $R = Statistics::R->new() or die "Problem with R : $!\n";
# Chargement de la bibliothèque FactoMineR (http://factominer.free.fr/)
# q` pour le chargement de la bibliothèque
$R->run(q`library(FactoMineR)`);
# Transposition du fichier original
# qq` pour l'interpolation des variables Perl
$R->run(<<"CMD"
data<-read.table("$ref_params->{file}", header=T, row.names=1, sep="\t")
data<-as.data.frame(t(data))
condition<-as.factor(c(rep("A",4),rep("B",4)))
data<-cbind.data.frame(condition,data)
colnames(data)[1]<-"condition"
CMD
);
# Fichier avec les graphes sur les variables et les individus
$R->run(<<"CMD"
res.pca<-PCA(data,ncp=6,graph=FALSE,quali.sup=1)
pdf("$pdf_acp")
plot(res.pca, choix="ind", axes=1:2, new.plot=F)
plot(res.pca, choix="var", axes=1:2, new.plot=F)
plot(res.pca, choix="ind", axes=2:3, new.plot=F)
plot(res.pca, choix="var", axes=2:3, new.plot=F)
dev.off()
CMD
);
# Pourcentage d'inertie pour les quatre premiers axes
# qq` pour la récupération de l'objet round par Perl
$R->run(q`round<-round(res.pca$eig[1:4,],2)`);
my $inertie = $R->get('round');
# Description des axes
my $dimdesc = $R->run(<<'CMD'
dimdesc<-dimdesc(res.pca)
print(dimdesc)
CMD
);
# Affichage des objets "round" et "dimdesc" dans le fichier de résultat
print {$fh} "Inertia:\n$inertie\n\nDimdesc:\n$dimdesc\n";
# Sauvegarde des objets "data'" et "res.pca" (fichier .RData)
my $r_save = $dir_name . '/Data.RData';
$R->run(qq`save(data,res.pca, file="$r_save")`);
# Fermeture du pont
$R->stop();
close $fh;
return;
}On obtient donc trois fichiers :
- « Data.RData » ;
- « PCA.pdf » (exemple du graphe des individus selon les deux premiers axes) ;
- « Result.txt » (extrait du fichier de résultat).
Inertia:
eigenvalue percentage of variance cumulative percentage of variance
comp 1 6508.60 31.65 31.65
comp 2 5249.58 25.53 57.18
comp 3 2841.99 13.82 71.01
comp 4 2036.85 9.91 80.91
Dimdesc:
$Dim.1
$Dim.1$quanti
correlation P-value
S213 0.9975904 3.491288e-08
S7134 0.9957727 1.882513e-07
S9211 0.9954746 2.309098e-07
S15923 0.9937002 6.221021e-07
S4168 0.9931895 7.857080e-07
.../...
$Dim.1$quali
P-value
condition 9.77511e-05
$Dim.1$category
Estimate P-value
B 77.91396 9.77511e-05
A -77.91396 9.77511e-05
$Dim.2
$Dim.2$quanti
correlation P-value
S10368 0.9911161 1.741237e-06
S6443 0.9910109 1.803658e-06
S5688 0.9907445 1.968439e-06
S442 0.9900412 2.450838e-06
S1041 0.9892404 3.089000e-06
.../...
$Dim.3
$Dim.3$quanti
correlation P-value
S6447 0.9858095 7.068008e-06
S12141 0.9729702 4.837528e-05
S16799 0.9705853 6.223013e-05
S193 0.9672014 8.605207e-05
S17531 0.9661212 9.475946e-05
S11034 0.9645758 1.082002e-04
.../...IV-B. Exemple de la bibliothèque FactoMineR▲
Pour les non-bio-informaticiens, ou pour ceux qui souhaiteraient tester ces méthodes sans avoir de fichier à disposition, un exemple est disponible dans la bibliothèque FactoMineR. Le jeu de données consiste en un tableau qui contient les performances réalisées par des athlètes lors de deux compétitions. Le fichier est disponible. Le script présenté ci-dessous reprend les codesles codes qui vous sont détaillés par les auteurs de la bibliothèque.
Les données :
Le tableau de données contient 41 lignes et 13 colonnes.
Les colonnes 1 à 12 sont des variables continues: les dix premières colonnes correspondent aux performances des athlètes pour les dix épreuves du décathlon et les colonnes 11 et 12 correspondent respectivement au rang et au nombre de points obtenus. La dernière colonne est une variable qualitative correspondant au nom de la compétition (Jeux Olympiques de 2004 ou Décastar 2004).
Et voici donc le script !
#!/usr/bin/perl
use strict;
use Carp;
use warnings;
use File::Basename;
use Statistics::R;
use Data::Dumper;
#===========================================
# Auteur : Stoyak - Developpez.com
# But : Statistics-R demo
# Date : 17/07/2014
#===========================================
#===========================================
# Définition des paramètres
#===========================================
# Répertoire de sortie
my $dir = 'C:/temp';
my %parameters = ( 'dir' => $dir, );
#===========================================
# Appel de la procédure ACP
#===========================================
ACP( \%parameters );
#===========================================
# Procédure ACP
#===========================================
sub ACP {
my ($ref_params) = @_;
# Fichiers de résultats ('.txt' et '.pdf')
my $pdf_acp = $ref_params->{dir} . '/DecathPCA.pdf';
my $result_file = $ref_params->{dir} . '/DecathResult.txt';
open my $fh, '>', $result_file || die "Impossible d'écrire dans le fichier $result_file\n";
# Déclaration de l'objet
my $R = Statistics::R->new() or die "Problem with R : $!\n";
# Chargement de la bibliothèque FactoMineR (http://factominer.free.fr/) et des données
$R->run( q`library(FactoMineR)`, q`data(decathlon)` );
# ACP sur les individus et les variables
$R->run(
<<"CMD"
pdf("$pdf_acp")
res.pca<-PCA(decathlon, scale.unit=TRUE, ncp=5, quanti.sup=c(11: 12), quali.sup=13, graph=FALSE)
plot(res.pca, choix="ind", axes=1:2, new.plot=F, habillage=13)
plot(res.pca, choix="var", axes=1:2, new.plot=F, habillage=13)
dev.off()
CMD
);
# Description des axes
my $dimdesc = $R->run(
<<"CMD"
dimdesc<-dimdesc(res.pca, axes=c(1,2))
print(dimdesc)
CMD
);
# Print de l'objet 'dimdesc' dans le fichier de résultat
print {$fh} "Dimdesc:\n$dimdesc\n";
close $fh;
return;
}Après exécution, vous devez avoir dans votre répertoire de sortie les deux fichiers « DecathPCA.pdf » et « DecathResult.txt » dont les contenus sont décrits par les auteurs.
V. Remerciements▲
Je tiens à remercier djibril de m'avoir suggéré d'écrire cet article et de m'avoir introduit dans la communauté des rédacteurs de developpez.com. Je le remercie pour m'avoir gentiment (!) grondée pour effectuer la mise à jour.
Je tiens également à remercier _Max_, ClaudeLELOUP et dourouc05 pour leurs corrections.
Je tiens à citer de nouveau la bibliothèque FactoMineR grâce à laquelle j'ai exposé des exemples plus complexes et complets.
N'hésitez pas à laissez des commentaires, des remarques ou propositions pour l'amélioration de ce tutoriel : 8 commentaires 