-
Notifications
You must be signed in to change notification settings - Fork 2
/
Rlite.pm
98 lines (62 loc) · 1.8 KB
/
Rlite.pm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
package Rlite;
## light interface to some useful R functions.
use strict;
use warnings;
use Carp;
use Resource;
####
sub G_test_matrix {
my ($data_rows_aref, $correction) = @_; # ( [ [1, 5], [11, 14], ...]) contingency table, one row_aref at a time.
my @data_rows = @$data_rows_aref;
## correction can be "williams" or "yates"
if ($correction && ! ($correction eq 'williams' || $correction eq 'yates')) {
confess "invalid correction specification";
}
my $Rsub_file = "Rlite/Rsubs/g.test.r";
my $Rsub_dir = &Resource::find_resource($Rsub_file) or confess "Error, cannot lcoate file: $Rsub_file";
my $ncol = scalar @{$data_rows[0]};
## build the Rscript
my $rscript = "source (\"$Rsub_dir/$Rsub_file\")\n"
. "x = c(";
my @vals;
foreach my $row (@data_rows) {
push (@vals, @$row);
}
$rscript .= join (",", @vals) . ")\n";
$rscript .= "m = matrix(x, ncol=$ncol, byrow=TRUE)\n";
$rscript .= "result=g.test(m";
if ($correction) {
$rscript .= " ,correct=\"$correction\"";
}
$rscript .= ")\n";
$rscript .= "print(result\$statistic)\n";
$rscript .= "print(result\$p.value)\n";
my $result = &_Run_Rscript($rscript);
my @txt_rows = split (/\n/, $result);
pop @txt_rows;
my $pvalue = pop @txt_rows;
pop @txt_rows;
pop @txt_rows;
my $Gval = pop @txt_rows;
$Gval =~ s/^\s+//;
return ($Gval, $pvalue);
}
################################################
## private:
####
sub _Run_Rscript {
my($rscript) = @_;
my $tmpfile = "R.tmp.$$.script";
open (my $fh, ">$tmpfile") or confess "Error, cannot write tmpfile $rscript";
print $fh $rscript;
close $fh;
## run R, catpure the output
my $command = "R --vanilla < $tmpfile";
my $result = `$command`;
if ($?) {
confess "Error running R command:\n$command\nret($?) ";
}
unlink($tmpfile);
return ($result);
}
1;