createBoxplotRscripts.pl

#!/usr/bin/perl -w
# Author: Abdullah Kahraman
# Date: 06.01.2015

###############################################################################
###############################################################################
### bla bla bla bla bla bla bla bla bla bla bla bla bla bla bla bla         ###
###############################################################################
###############################################################################

use strict;
use warnings;
use Getopt::Long;
use File::Basename;

my (
    # variable for parameters which are read in from commandline
    $help,
   );

##############################################################################
### read all needed parameters from commandline ##############################

&GetOptions(
    "help!"      => \$help,     # print this help
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

##############################################################################

# help
if ($help) {printHelp(); exit}

##############################################################################
### SETTINGS #################################################################
##############################################################################

##############################################################################
### SUBROUTINES ##############################################################
############################################################################## 

###############################################################################
sub printHelp {
###############################################################################
    # prints a help about the using and parameters of this scripts 
    # (execute if user types commandline parameter -h)
    # param:  no paramaters
    # return: no return value

    my (
	$usage,
	$sourceCode,
	@rows,
	$row,
	$option,
	$scriptInfo,
	$example,
       );

    $usage = "$0\n";


    print "\nUsage: " .  $usage . "\n";

    print "Valid options are:\n\n";
    open(MYSELF, "$0") or
      die "Cannot read source code file $0: $!\n";
    $sourceCode .= join "", <MYSELF>;
    close MYSELF;
    $sourceCode =~ s/^.+?\&GetOptions\(\n//s;
    $sourceCode =~ s/\n\).+$//s;
    @rows = split /\n/, $sourceCode;
    foreach $row (@rows){
        $option = $row;
	$option =~ s/\s+\"//g;
	$option =~ s/\"\s.+\#/\t\#/g;
	$option =~ s/=./\t<value> [required]/;
	$option =~ s/:./\t<value> [optional]/;
	$option =~ s/!/\t<non value> [optional]/;

	$row =~ s/^.*//;
	print "\t";
	printf("%-1s%-30s%-30s\n", "-",$option,$row);

    } # end of foreach $row (@rows)
    print "\n";
    print "Options may be abreviated, e.g. -h for --help\n\n";

    $example  = "$0";
}

##############################################################################
### END OF SUBROUTINES########################################################
############################################################################## 


############
### MAIN ###
############

for (my $i=15; $i<=38; $i++) {
    my $input = "string3d_human/scores_".$i."_hit.txt.gz";
    my $output = "scores_".$i."_hit_plots";

    my $script = <<SCRIPT;
library(vioplot)
header = c("isHit", "KEGGpathway", "stringId1", "stringId2", "qUniprotAc1", "qUniprotAc2", "probisUniprotId1", "probisUniprotId2", "probisPdbId1", "probisPdbId2", "uniprotId1", "uniprotId2", "qPdbId-chainId2", "bs1-bs2", "probisEvalue1", "probisRmsd1", "probisZvalue1", "probisTotal1", "probisHit1", "probisFalse1", "probisRatio1", "probisEvalue2", "probisRmsd2", "probisZvalue2", "probisTotal2", "probisHit2", "probisFalse2", "probisRatio2", "hhsearchProb1", "hhsearchIdent1", "hhsearchInsideAlign1", "hhsearchOutsideAlign1", "hhsearchRatio1", "hhsearchProb2", "hhsearchIdent2", "hhsearchInsideAlign2", "hhsearchOutsideAlign2", "hhsearchRatio2");
i=$i
rowsTop = read.table(gzfile("$input"), nrows = 5)
classes = sapply(rowsTop, class)
data = read.table(gzfile("$input"),colClasses = classes, comment.char = "")
SCRIPT

    if($i==30 or $i==38) { 
	$script .= 'data[,2] = as.numeric(sub("%","",data[,2]))/100'."\n";
    }
    # for scores that indicate better performance with increasing value, we need to sort backwards
    if($i==15 or $i==16 or $i==20 or $i==22 or $i==23 or $i==27 or $i==32 or $i==37) {
	$script .= 'data = data[order(data[,2]),]'."\n";
    } else {
	$script .= 'data = data[order(-data[,2]),]'."\n";
    }

    $script .= <<SCRIPT;
data.dens <- density(data[,2])
png("$output.png", width=1400, height=1400)
par(mfrow=c(2,2))
vioplot(data[,2], names=c(header[i]))
data.box = boxplot(data[,2])
data.hist = hist(data[,2], breaks=100, main=paste("Histogram of",header[i],sep=" "))
plot(data.dens)
dev.off()
cat(i, "\\t", header[i], "\\t", data.hist\$breaks, "\\t", data.hist\$mids, "\\t", data.hist\$counts, "\\t", data.hist\$density, "\\n", data.box\$stats, "\\n", file="$output.txt")
SCRIPT

    open(O, ">$output.R") or die "\nERROR: Failed to create $output.R: $!\n\n";
    print O $script;
    close(O);
}