|
0
|
1 #!/usr/bin/perl
|
|
|
2 #-------------------------------------------------------------------------------
|
|
|
3 #
|
|
|
4 # Plot a QQ plot (using R)
|
|
|
5 # Data is feed as a 1 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 = 'QQ-plot';
|
|
|
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 vector
|
|
|
24 open TXT, "> $txtFile" or die "Cannot open output file '$txtFile'\n";
|
|
|
25 print TXT "x\n";
|
|
|
26 for( $ln = 0 ; $l = <STDIN> ; ) {
|
|
|
27 chomp $l;
|
|
|
28
|
|
|
29 # Does the string contain exactly one number? (can be float)
|
|
|
30 if( $l =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ ) { print TXT "$l\n"; }
|
|
|
31 }
|
|
|
32 close TXT;
|
|
|
33
|
|
|
34 #---
|
|
|
35 # Create an R program, save QQ-plot as PNG image
|
|
|
36 #---
|
|
|
37
|
|
|
38 open R, "| R --vanilla --slave " or die "Cannot open R program\n";
|
|
|
39 print R <<EOF;
|
|
|
40
|
|
|
41 qqplot <- function( x, title ) {
|
|
|
42 keep <- (x > 0) & (x <= 1) & ( ! is.na(x) );
|
|
|
43 x <- x[keep]
|
|
|
44 s <- sort(x);
|
|
|
45 ly <- -log10(s);
|
|
|
46
|
|
|
47 n <- length(s);
|
|
|
48 lx <- -log10( (1:n) / (n+1) )
|
|
|
49
|
|
|
50 # Show auto range
|
|
|
51 #par( mfrow=c(2,1) );
|
|
|
52 #plot( lx, ly, main=title, xlab="-Log[ rank / (N+1) ]", ylab="-Log[ p ]" );
|
|
|
53 #abline( 0 , 1 , col='red');
|
|
|
54
|
|
|
55 # Show full range in both plots
|
|
|
56 range <- c(0 , max(lx, ly) );
|
|
|
57 plot( lx, ly, xlim=range, ylim=range, main=title, xlab="-Log[ rank / (N+1) ]", ylab="-Log[ p ]" );
|
|
|
58 abline( 0 , 1 , col='red');
|
|
|
59 }
|
|
|
60
|
|
|
61 png('$pngFile', width = 1024, height = 1024);
|
|
|
62
|
|
|
63 data <- read.csv("$txtFile", sep='\t', header = TRUE);
|
|
|
64 qqplot( data\$x, "$base" );
|
|
|
65
|
|
|
66 dev.off();
|
|
|
67 quit( save='no' )
|
|
|
68 EOF
|
|
|
69
|
|
|
70 close R;
|
|
|
71
|
|
|
72 #---
|
|
|
73 # Show figure
|
|
|
74 #---
|
|
|
75
|
|
|
76 $os = `uname`;
|
|
|
77 $show = "eog";
|
|
|
78 if( $os =~ "Darwin" ) { $show = "open"; }
|
|
|
79 `$show $pngFile`;
|
|
|
80
|
