Mercurial > repos > greg > snpeff_v2_from_pablo

comparison snpEff_2_1a/scripts/qqplot.pl @ 0:f8eaa3f8194b default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded snpEff_v2_1a_core.tgz from Pablo Cingolani
author greg
date Fri, 20 Apr 2012 14:47:09 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f8eaa3f8194b
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