Interaction des langages R et Perl pour les statistiques

Le but de cet article est de vous présenter comment faire interagir les langages R et Perl.
8 commentaires Donner une note à l'article (5)

Article lu   fois.

L'auteur

Profil ProSite personnel

Liens sociaux

Viadeo Twitter Facebook Share on Google+   

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.

Image non disponible
Extrait du tableau de données

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 :

  1. 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 ;
  2. 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 :

Sans pont partagé
Sélectionnez
use Statistics::R;
# Création du pont entre Perl et R et démarrage de R
my $R = Statistics::R->new();
…
Pont partagé
Sélectionnez
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).

r_bin
Sélectionnez
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 :

 
Sélectionnez
# 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) :

 
Sélectionnez
$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().

 
Sélectionnez
my $fichier_commandes = 'C:/temp/Data.r';
$R->run_from_file($fichier_commandes);
print $R->get('c');
Data.r
Sélectionnez
a <- 2
b <- 5
c <- a * b

On 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 :
 
Sélectionnez
#!/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).

data frame non fonctionnel dans Perl
Sélectionnez
#!/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;
Résultat
Sélectionnez
$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 :
Création d'une image
Sélectionnez
#!/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
Chargement d'un fichier de données
Sélectionnez
#!/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');
Résultat
Sélectionnez
      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
Sélectionnez
# Chargement d'une bibliothèque
$R->run(q`library(FactoMineR)`);
  • Sauvegarde des objets :
Sauvegarde des objets
Sélectionnez
# 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
Sélectionnez
# Suppression des objets
$R->run(
    q`rm(v1,v2,data)`,
    q`ls<-ls()`
);
  • Chargement du fichier .Rdata
Chargement du fichier .RData
Sélectionnez
# 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 :

Procédure ACP
Sélectionnez
#!/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).
Extrait du fichier de résultat
Sélectionnez
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 !

Script sur les données du decathlon
Sélectionnez
#!/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 Donner une note à l'article (5)

Vous avez aimé ce tutoriel ? Alors partagez-le en cliquant sur les boutons suivants : Viadeo Twitter Facebook Share on Google+   

  

Les sources présentées sur cette page sont libres de droits et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une œuvre intellectuelle protégée par les droits d'auteur. Copyright © 2011-2014 stoyak. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.