0
|
1 #!/usr/bin/perl
|
|
2 #-------------------------------------------------------------------------------
|
|
3 #
|
|
4 # Plot a smooth scatter plot
|
|
5 # Data is feed as two column of numbers
|
|
6 #
|
|
7 # Note: Any line that does not match a numeric regular expression, is filtered out).
|
|
8 #
|
|
9 # Pablo Cingolani
|
|
10 #-------------------------------------------------------------------------------
|
|
11
|
|
12 #-------------------------------------------------------------------------------
|
|
13 # Main
|
|
14 #-------------------------------------------------------------------------------
|
|
15
|
|
16 # Parse command line option (file base name)
|
|
17 $base = 'smoothScatter';
|
|
18 if( $ARGV[0] ne '' ) { $base = $ARGV[0]; }
|
|
19
|
|
20 $pngFile = "$base.png";
|
|
21 $txtFile = "$base.txt";
|
|
22
|
|
23 # Read STDIN and create an R table
|
|
24 open TXT, "> $txtFile" or die "Cannot open output file '$txtFile'\n";
|
|
25 print TXT "x\ty\n";
|
|
26 for( $ln = 0 ; $l = <STDIN> ; ) {
|
|
27 chomp $l;
|
|
28 ($x, $y) = split /\t/, $l;
|
|
29
|
|
30 # Does the string contain exactly one number? (can be float)
|
|
31 if(( $x =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ ) && ( $y =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ )) { print TXT "$x\t$y\n"; }
|
|
32 }
|
|
33 close TXT;
|
|
34
|
|
35 #---
|
|
36 # Create an R program, save histogram plot as PNG image
|
|
37 #---
|
|
38
|
|
39 open R, "| R --vanilla --slave " or die "Cannot open R program\n";
|
|
40 print R <<EOF;
|
|
41
|
|
42 smoothLowess <- function( x, y, title, q=1.0 ) {
|
|
43 # Show only this part of the data
|
|
44 xmin <- quantile( x, 1-q )
|
|
45 xmax <- quantile( x, q )
|
|
46
|
|
47 ymin <- quantile( y, 1-q )
|
|
48 ymax <- quantile( y, q )
|
|
49
|
|
50 keep <- (x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax);
|
|
51 qx <- x[ keep ]
|
|
52 qy <- y[ keep ]
|
|
53
|
|
54 smoothScatter(qx, qy, main=title, ylab='Y (column 2)', xlab='X (column 1)');
|
|
55 lines( lowess(qx,qy), col='orange' );
|
|
56 }
|
|
57
|
|
58 png('$pngFile', width = 1024, height = 1024);
|
|
59 par( mfrow=c(2,1) );
|
|
60
|
|
61 data <- read.csv("$txtFile", sep='\\t', header = TRUE);
|
|
62 x <- data\$x
|
|
63 y <- data\$y
|
|
64
|
|
65 smoothLowess(x, y, "Smooth scatter plot and Lowess", 1.0);
|
|
66 smoothLowess(x, y, "Smooth scatter plot and Lowess: Quantile [2% - 98%]", 0.98);
|
|
67
|
|
68 dev.off();
|
|
69 quit( save='no' )
|
|
70 EOF
|
|
71
|
|
72 close R;
|
|
73
|
|
74 #---
|
|
75 # Show figure
|
|
76 #---
|
|
77
|
|
78 $os = `uname`;
|
|
79 $show = "eog";
|
|
80 if( $os =~ "Darwin" ) { $show = "open"; }
|
|
81 `$show $pngFile`;
|
|
82
|