Foro - Perl en Español

por **Alfumao** » 2012-10-10 07:49 @367

Hola a todos,

Estoy buscando un módulo que me permita pasar valores de variables desde Perl a scripts de R.

Encontré algunos pero no daban muy buena impresión en ciertas páginas dedicadas Perl, decían que no eran muy fiables en general...

Al final me he quedado con este, pero no entiendo muy bien en el manual cómo se haría la interpolación de variables de Perl a R.

¿Alguien ha utilizado este módulo alguna vez y me puede poner un ejemplo?

Muchísimas gracias por adelantado.

#!usr/bin/perl
use warnings;
use strict;
use Getopt::Long;
use Statistics::R;
 
#usage: perl /home/alfumao/scripts/PerleR_ToC.pl -g [ref_gemone] -p [/path_to_Table] -t [Tablename]
 
 
my ( $path, $genome, $table );
GetOptions(
    'path=s'         => \$path,
    'genome=s'        => \$genome,
    'table=s'        => \$table,
    );
 
chdir $path or die "ERROR: Unable to enter $path: $!\n";
print "current dir =$path\n";
# Create a communication bridge with R and start R
my $R = Statistics::R->new();
$R->set('genome', "$genome");
$R->set('file', "$table");
 
my $out_genome = $R->get('genome');
print "genome = $out_genome\n";
my $out_table = $R->get('file');
print "file = $out_table\n";
 
my$ToC=$R->run(
q`library(Rsamtools)`,
q`library(GenomicFeatures)`,
q`bamlist=list()`,
q`src_files=list.files(pattern="*.bam")`,
q`for(filename in src_files)`,
q`{`,
q`tmp=readBamGappedAlignments(filename)`,
q`bamlist[[length(bamlist)+1]]=GRanges(seqnames=rname(tmp),`,
q`ranges=IRanges(start=start(tmp),end=end(tmp)),`,
q`strand=rep("*",length(tmp)))`,
q`}`,
q`names(bamlist)=src_files`,
q`library(GenomicFeatures)`,
q`txdb=makeTranscriptDbFromUCSC(genome,tablename='knownGene')`,
q`tx_by_gene=transcriptsBy(txdb,'gene')`,
q`ex_by_gene=exonsBy(txdb,'gene')`,
q`toc=data.frame(rep(NA,length(tx_by_gene)))`,
q`for(i in 1:length(bamlist))`,
q`{`,
q`toc[,i]=countOverlaps(tx_by_gene,bamlist[[i]])`,
q`}`,
q`rownames(toc)=names(tx_by_gene)`,
q`colnames(toc)=names(bamlist)`,
q`dim(toc)`,
q`head(toc)`,
q`write.table(toc, file ,quote=FALSE,sep="\t")`,
q`print("ok")`
);
 
 
$R->stop();
Coloreado en 0.003 segundos,  usando GeSHi 1.0.8.4

pero me sale el siguiente error:

Error stopping R: 256

El código de R está testado y me funciona en R.

¿Alguna idea?

por **explorer** » 2012-10-15 06:11 @299

Habría que saber qué significa el error 256 para R.

Sobre el código,

No son necesarias las comillas dobles en las líneas 21 y 22
Al ser un código R tan largo, es mucho más cómodo escribirlo en modo here-doc (documento incrustado). En la documentación de run(), en la página de manual del módulo, tienes un ejemplo

Edito: el mensaje de error lo genera el método stop() del módulo Statistics::R; el valor de 256 es lo que devuelve el intento de parada de R.

por **Alfumao** » 2012-10-17 03:36 @192

Gracias, explorer, he hecho los cambios y revisado la posibilidad del here-doc, como me has recomendado, y funciona perfectamente.

¡Gracias a Perl que me permite hacer tantas cosas interesantes para mi trabajo!

Porque aunque sea una estupidez, aún me resisto a ponerme a aprender Python...

por **Alfumao** » 2012-10-18 03:32 @189

Era todo demasiado idílico...

Al cambiar el código del here-doc han empezado a saltar errores que no entiendo.

Yo creía que el here-doc era una forma de introducir código de R y que Perl no iba a interferir con dicho código, solo enviarlo a R y hacer lo que fuera pertinente, pero me encuentro que Perl tiene problemas con el código pidiéndome que declare variables que hay dentro del here-doc (o eso creo entender)

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#!usr/bin/perl
use warnings;
use strict;
use Getopt::Long;
use Statistics::R;
 
#usage: perl /storage/Runs/CIC/ANALYSIS/JL_scripts/PerleR_ToC.pl -g [ref_gemone] -p [/path_to_Table] -t [DatabaseTablename] -c [path/ToC]
 
 
my ( $path, $genome, $table, $toc, $libType );
GetOptions(
    'path=s'         => \$path,
    'genome=s'       => \$genome,
    'tablename=s'    => \$table,
    'counttable=s'   => \$toc,
    #'sconditions=s'   => \@conds,
    'libtype=s'      => \$libType,  #solo para un tipo de librerias
    );
 
chdir $path or die "ERROR: Unable to enter $path: $!\n";
print "current dir =$path\n";
 
# Create a communication bridge with R and start R
my $R = Statistics::R->new();
 
$R->set('genome', $genome);
my $out_genome = $R->get('genome');
print "genome = $out_genome\n";
 
$R->set('tablename', "$table");
my $out_table = $R->get('tablename');
print "database table = $out_table\n";
 
$R->set('file', "$toc");
my $out_toc = $R->get('file');
print "ToC = $out_toc\n";
 
#$R->set('conds', "@conds");
#my @out_conds = $R->get('conds');
#print "Conditions = @out_conds\n";
 
#$R->set('libType', "$libType");
#my $out_libType = $R->get('libType');
#print "Libraries type = $out_libType\n";
 
# Here-doc with multiple R commands:
my $cmds = <<EOF;
 
library(DESeq)
 
#personalizar el nombre de la tabla
countsTable <- read.delim( "/solexa/storage/Runs/CIC/ANALYSIS/JA-01-mRNA/Sequences/FASTQ_original/mm9_JA01_ToC.tsv", header=TRUE, row.names = 1, stringsAsFactors=TRUE )
head (countsTable )
 
 
ToCDesign = data.frame(
row.names = colnames( countsTable ),
condition = c( "untreated", "treated" ),
libType = c( "single-end", "single-end" ) )
ToCDesign
 
singleSamples = ToCDesign$libType == "single-end"
countTable = countsTable [ , singleSamples ]
condition = ToCDesign$condition[ singleSamples ]
head(countTable)
condition
condition = factor( c( "untreated", "treated" ) )
 
cds = newCountDataSet( countTable, condition )
cds = estimateSizeFactors( cds )
sizeFactors( cds )
head( counts( cds, normalized=TRUE ) )
cds = estimateDispersions( method ="blind", cds )
 
str( fitInfo(cds) )
plotDispEsts( cds )
head( fData(cds) )
 
#2)DATA QUALITY ASSESMENT
 
#Heatmap of the count table
cdsFullBlind = estimateDispersions( cds, method = "blind" )
vsdFull = varianceStabilizingTransformation( cdsFullBlind )
library("RColorBrewer")
library("gplots")
 
select = order(rowMeans(counts(cds)), decreasing=TRUE)[1:30]
hmcol = colorRampPalette(brewer.pal(9, "GnBu"))(100)
 
heatmap.2(exprs(vsdFull)[select,], col = hmcol, trace="none", margin=c(10, 6))
heatmap.2(counts(cds)[select,], col = hmcol, trace="none", margin=c(10,6))
 
dists = dist ( t( exprs(vsdFull) ) )
mat = as.matrix( dists )
 
rownames(mat) = colnames(mat) = with(pData(cdsFullBlind), paste(condition, libType, sep=" : "))
heatmap.2(mat, trace="none", col = rev(hmcol), margin=c(13, 13))
dev.off()
 
print('ok')
EOF
my $out2 = $R->run($cmds);
 
 
#$R->stop();
Coloreado en 0.002 segundos,  usando GeSHi 1.0.8.4

Los errores que me da son los siguientes (al poner línea 47 asumo que es el here-doc total):

Global symbol "@condition" requires explicit package name at /home/scripts/PerleR_DE_Qual.pl line 47.

Bareword "singleSamples" not allowed while "strict subs" in use at /home/scripts/PerleR_DE_Qual.pl line 47.

Execution of /home/scripts/PerleR_DE_Qual.pl aborted due to compilation errors.

¿Cómo se puede arreglar esto?

por **explorer** » 2012-10-18 06:36 @317

El problema está en la línea 64:

condition = ToCDesign$condition[ singleSamples ]

Perl "cree" que $condition[ singleSamples ] es el acceso a un elemento de un array @condition, y como tienes activada la opción 'strict', Perl te avisa de que es una variable que no has declarado antes (ya se ve entonces lo útil que es poner 'strict').

Lo que pasa es que, por defecto, los here-doc se comportan como si fueran un texto doblemente entrecomillado, por lo que Perl buscará por cosas que se parezcan a variables Perl, para interpolarlas.

En cambio, si se trata de un texto literal, y no queremos que Perl haga ninguna interpolación, solo tenemos que indicarlo a Perl, de la siguiente manera:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

my $cmds = <<'EOF';
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

Tienes más información en perlop, en la sección <<EOF:

Comillas simples
Las comillas simples indican que el texto se va a tratar literalmente, sin interpolación de su contenido.

por **Alfumao** » 2012-10-18 08:06 @379

Muchas gracias una vez más explorer

Foro - Perl en Español

Módulo Statistics-R

Módulo Statistics-R

Publicidad

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

Re: Módulo Statistics-R

¿Quién está conectado?