Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.
http://perlenespanol.com/foro/

Módulo Statistics-R

http://perlenespanol.com/foro/modulo-statistics-r-t7587.html

Página 1 de 1

Módulo Statistics-R

NotaPublicado: 2012-10-10 07:49 @367
por Alfumao
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.

Re: Módulo Statistics-R

NotaPublicado: 2012-10-10 09:55 @455
por explorer
Por lo que veo en el manual, con set() se le puede pasar a R un escalar o una referencia a un array, y con get(), obtenerlos.

Re: Módulo Statistics-R

NotaPublicado: 2012-10-11 02:10 @132
por Alfumao
Gracias, explorer, estaba mirando otra versión del modulo (u otro módulo que se llama muy parecido) y era bastante más parco en explicaciones y más complicado de descifrar :wink:

Re: Módulo Statistics-R

NotaPublicado: 2012-10-15 03:43 @197
por Alfumao
Hola de nuevo, ahora vienen los problemas usando el módulo.
He intentado ejecutar este código:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!usr/bin/perl
  2. use warnings;
  3. use strict;
  4. use Getopt::Long;
  5. use Statistics::R;
  6.  
  7. #usage: perl /home/alfumao/scripts/PerleR_ToC.pl -g [ref_gemone] -p [/path_to_Table] -t [Tablename]
  8.  
  9.  
  10. my ( $path, $genome, $table );
  11. GetOptions(
  12.     'path=s'         => \$path,
  13.     'genome=s'        => \$genome,
  14.     'table=s'        => \$table,
  15.     );
  16.  
  17. chdir $path or die "ERROR: Unable to enter $path: $!\n";
  18. print "current dir =$path\n";
  19. # Create a communication bridge with R and start R
  20. my $R = Statistics::R->new();
  21. $R->set('genome', "$genome");
  22. $R->set('file', "$table");
  23.  
  24. my $out_genome = $R->get('genome');
  25. print "genome = $out_genome\n";
  26. my $out_table = $R->get('file');
  27. print "file = $out_table\n";
  28.  
  29. my$ToC=$R->run(
  30. q`library(Rsamtools)`,
  31. q`library(GenomicFeatures)`,
  32. q`bamlist=list()`,
  33. q`src_files=list.files(pattern="*.bam")`,
  34. q`for(filename in src_files)`,
  35. q`{`,
  36. q`tmp=readBamGappedAlignments(filename)`,
  37. q`bamlist[[length(bamlist)+1]]=GRanges(seqnames=rname(tmp),`,
  38. q`ranges=IRanges(start=start(tmp),end=end(tmp)),`,
  39. q`strand=rep("*",length(tmp)))`,
  40. q`}`,
  41. q`names(bamlist)=src_files`,
  42. q`library(GenomicFeatures)`,
  43. q`txdb=makeTranscriptDbFromUCSC(genome,tablename='knownGene')`,
  44. q`tx_by_gene=transcriptsBy(txdb,'gene')`,
  45. q`ex_by_gene=exonsBy(txdb,'gene')`,
  46. q`toc=data.frame(rep(NA,length(tx_by_gene)))`,
  47. q`for(i in 1:length(bamlist))`,
  48. q`{`,
  49. q`toc[,i]=countOverlaps(tx_by_gene,bamlist[[i]])`,
  50. q`}`,
  51. q`rownames(toc)=names(tx_by_gene)`,
  52. q`colnames(toc)=names(bamlist)`,
  53. q`dim(toc)`,
  54. q`head(toc)`,
  55. q`write.table(toc, file ,quote=FALSE,sep="\t")`,
  56. q`print("ok")`
  57. );
  58.  
  59.  
  60. $R->stop();
Coloreado en 0.004 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?

Re: Módulo Statistics-R

NotaPublicado: 2012-10-15 06:11 @299
por explorer
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.

Re: Módulo Statistics-R

NotaPublicado: 2012-10-17 03:36 @192
por Alfumao
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...

Re: Módulo Statistics-R

NotaPublicado: 2012-10-18 03:32 @189
por Alfumao
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
  1. #!usr/bin/perl
  2. use warnings;
  3. use strict;
  4. use Getopt::Long;
  5. use Statistics::R;
  6.  
  7. #usage: perl /storage/Runs/CIC/ANALYSIS/JL_scripts/PerleR_ToC.pl -g [ref_gemone] -p [/path_to_Table] -t [DatabaseTablename] -c [path/ToC]
  8.  
  9.  
  10. my ( $path, $genome, $table, $toc, $libType );
  11. GetOptions(
  12.     'path=s'         => \$path,
  13.     'genome=s'       => \$genome,
  14.     'tablename=s'    => \$table,
  15.     'counttable=s'   => \$toc,
  16.     #'sconditions=s'   => \@conds,
  17.     'libtype=s'      => \$libType,  #solo para un tipo de librerias
  18.     );
  19.  
  20. chdir $path or die "ERROR: Unable to enter $path: $!\n";
  21. print "current dir =$path\n";
  22.  
  23. # Create a communication bridge with R and start R
  24. my $R = Statistics::R->new();
  25.  
  26. $R->set('genome', $genome);
  27. my $out_genome = $R->get('genome');
  28. print "genome = $out_genome\n";
  29.  
  30. $R->set('tablename', "$table");
  31. my $out_table = $R->get('tablename');
  32. print "database table = $out_table\n";
  33.  
  34. $R->set('file', "$toc");
  35. my $out_toc = $R->get('file');
  36. print "ToC = $out_toc\n";
  37.  
  38. #$R->set('conds', "@conds");
  39. #my @out_conds = $R->get('conds');
  40. #print "Conditions = @out_conds\n";
  41.  
  42. #$R->set('libType', "$libType");
  43. #my $out_libType = $R->get('libType');
  44. #print "Libraries type = $out_libType\n";
  45.  
  46. # Here-doc with multiple R commands:
  47. my $cmds = <<EOF;
  48.  
  49. library(DESeq)
  50.  
  51. #personalizar el nombre de la tabla
  52. 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 )
  53. head (countsTable )
  54.  
  55.  
  56. ToCDesign = data.frame(
  57. row.names = colnames( countsTable ),
  58. condition = c( "untreated", "treated" ),
  59. libType = c( "single-end", "single-end" ) )
  60. ToCDesign
  61.  
  62. singleSamples = ToCDesign$libType == "single-end"
  63. countTable = countsTable [ , singleSamples ]
  64. condition = ToCDesign$condition[ singleSamples ]
  65. head(countTable)
  66. condition
  67. condition = factor( c( "untreated", "treated" ) )
  68.  
  69. cds = newCountDataSet( countTable, condition )
  70. cds = estimateSizeFactors( cds )
  71. sizeFactors( cds )
  72. head( counts( cds, normalized=TRUE ) )
  73. cds = estimateDispersions( method ="blind", cds )
  74.  
  75. str( fitInfo(cds) )
  76. plotDispEsts( cds )
  77. head( fData(cds) )
  78.  
  79. #2)DATA QUALITY ASSESMENT
  80.  
  81. #Heatmap of the count table
  82. cdsFullBlind = estimateDispersions( cds, method = "blind" )
  83. vsdFull = varianceStabilizingTransformation( cdsFullBlind )
  84. library("RColorBrewer")
  85. library("gplots")
  86.  
  87. select = order(rowMeans(counts(cds)), decreasing=TRUE)[1:30]
  88. hmcol = colorRampPalette(brewer.pal(9, "GnBu"))(100)
  89.  
  90. heatmap.2(exprs(vsdFull)[select,], col = hmcol, trace="none", margin=c(10, 6))
  91. heatmap.2(counts(cds)[select,], col = hmcol, trace="none", margin=c(10,6))
  92.  
  93. dists = dist ( t( exprs(vsdFull) ) )
  94. mat = as.matrix( dists )
  95.  
  96. rownames(mat) = colnames(mat) = with(pData(cdsFullBlind), paste(condition, libType, sep=" : "))
  97. heatmap.2(mat, trace="none", col = rev(hmcol), margin=c(13, 13))
  98. dev.off()
  99.  
  100. print('ok')
  101. EOF
  102. my $out2 = $R->run($cmds);
  103.  
  104.  
  105. #$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?

Re: Módulo Statistics-R

NotaPublicado: 2012-10-18 06:36 @317
por explorer
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
  1. 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.

Re: Módulo Statistics-R

NotaPublicado: 2012-10-18 08:06 @379
por Alfumao
Muchas gracias una vez más explorer ;)
Todos los horarios son UTC - 6 horas [ DST ]
Página 1 de 1
Powered by phpBB® Forum Software © phpBB Group
https://www.phpbb.com/