Mercurial > repos > greg > snpeff_v2_from_pablo
diff snpEff_2_1a/scripts/smoothScatter.pl @ 0:f8eaa3f8194b default tip
Uploaded snpEff_v2_1a_core.tgz from Pablo Cingolani
author | greg |
---|---|
date | Fri, 20 Apr 2012 14:47:09 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/snpEff_2_1a/scripts/smoothScatter.pl Fri Apr 20 14:47:09 2012 -0400 @@ -0,0 +1,82 @@ +#!/usr/bin/perl +#------------------------------------------------------------------------------- +# +# Plot a smooth scatter plot +# Data is feed as two column of numbers +# +# Note: Any line that does not match a numeric regular expression, is filtered out). +# +# Pablo Cingolani +#------------------------------------------------------------------------------- + +#------------------------------------------------------------------------------- +# Main +#------------------------------------------------------------------------------- + +# Parse command line option (file base name) +$base = 'smoothScatter'; +if( $ARGV[0] ne '' ) { $base = $ARGV[0]; } + +$pngFile = "$base.png"; +$txtFile = "$base.txt"; + +# Read STDIN and create an R table +open TXT, "> $txtFile" or die "Cannot open output file '$txtFile'\n"; +print TXT "x\ty\n"; +for( $ln = 0 ; $l = <STDIN> ; ) { + chomp $l; + ($x, $y) = split /\t/, $l; + + # Does the string contain exactly one number? (can be float) + if(( $x =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ ) && ( $y =~ /^[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?$/ )) { print TXT "$x\t$y\n"; } +} +close TXT; + +#--- +# Create an R program, save histogram plot as PNG image +#--- + +open R, "| R --vanilla --slave " or die "Cannot open R program\n"; +print R <<EOF; + +smoothLowess <- function( x, y, title, q=1.0 ) { + # Show only this part of the data + xmin <- quantile( x, 1-q ) + xmax <- quantile( x, q ) + + ymin <- quantile( y, 1-q ) + ymax <- quantile( y, q ) + + keep <- (x >= xmin) & (x <= xmax) & (y >= ymin) & (y <= ymax); + qx <- x[ keep ] + qy <- y[ keep ] + + smoothScatter(qx, qy, main=title, ylab='Y (column 2)', xlab='X (column 1)'); + lines( lowess(qx,qy), col='orange' ); +} + +png('$pngFile', width = 1024, height = 1024); +par( mfrow=c(2,1) ); + +data <- read.csv("$txtFile", sep='\\t', header = TRUE); +x <- data\$x +y <- data\$y + +smoothLowess(x, y, "Smooth scatter plot and Lowess", 1.0); +smoothLowess(x, y, "Smooth scatter plot and Lowess: Quantile [2% - 98%]", 0.98); + +dev.off(); +quit( save='no' ) +EOF + +close R; + +#--- +# Show figure +#--- + +$os = `uname`; +$show = "eog"; +if( $os =~ "Darwin" ) { $show = "open"; } +`$show $pngFile`; +