cov.pl

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

###############################################################################
###############################################################################
### This script calculates the variance, Pearson product-moment correlation ###
### coefficient and r^2.                                                    ###
###############################################################################
###############################################################################

use strict;
use Getopt::Long;

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

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

&GetOptions(
 "help!"     => \$help,    # print this help
 "in=s"      => \$infile,  # calculate covariance and correlation coefficient within this file
 "col=s"     => \$col,     # column number of values to calculate cov and r2, seperated by colon (starts from 0)
 "del:s"     => \$del,     # delimeter of columns[default is \t][optional]
 "gnuplot=s" => \$gnuplot, # path to gnuplot if scatter plot is desired
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

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

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

$del ||= "\t";

##############################################################################
### 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 -in string3d_human_exp_keggTest_bitScore_alignColsSum_calib_kegg_histo.txt -col 0:5\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 -in=../../databases/locuslink/locuslink_031007.txt ";
}

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

sub createPNG{

  (my $in, my $col1, my $col2) = @_;
  my $outfile = $in;
  $col1++;
  $col2++;
  $outfile =~ s/\.txt/\.png/;

  open (GP, "|$gnuplot") or die "no gnuplot";
  # force buffer to flush after each write
  use FileHandle;
  GP->autoflush(1);
  my $script = <<EOF;

  set terminal postscript eps color
  set output 'temp.eps'
  plot "$in" using $col1:$col2
EOF
  print GP $script;
  close GP;

  system("convert -density 144x144 temp.eps $outfile");
  unlink("temp.eps");
}

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


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

if (defined($col)) {
  (my $col1, my $col2, my $tmp) = split(/\:/,$col);
  &createPNG($infile, $col1, $col2) if(defined $gnuplot);

  my $avg1 = 0;
  my $avg2 = 0;
  my $i = 0;
  my @lines = ();
  while(<>) {
    if(!/^\#/){
      $i++;
      my $lineBuffer = $_;
      chomp($lineBuffer);
      my @lineArray = split($del, $lineBuffer);
      $avg1 += $lineArray[$col1];
      $avg2 += $lineArray[$col2];
      push(@lines, $lineBuffer);
    }
  }

  if($i==0){
    print "$infile\tnan\tnan\tnan\tnan\n"
  }
  $avg1 = $avg1/$i;
  $avg2 = $avg2/$i;

  my $cov12 = 0;
  my $var1 = 0;
  my $var2 = 0;
  foreach my $line (@lines) {
    if($line !~ /^\#/){
      my @lineArray = split($del, $line);
      $cov12  += ($lineArray[$col1]-$avg1)*($lineArray[$col2]-$avg2);
      $var1 += ($lineArray[$col1]-$avg1)**2;
      $var2 += ($lineArray[$col2]-$avg2)**2;
    }
  }

  $cov12 = $cov12/($i-1);
  $var1 = $var1/($i-1);
  $var2 = $var2/($i-1);

  my $std1 = sqrt($var1);
  my $std2 = sqrt($var2);

  # Pearson product-moment correlation coefficient
  my $r = $cov12/($std1*$std2);
  # correlation with itself is obviously 1
  my $r11 = 1;
  my $r22 = 1;

  my $r2 = $r**2;
  my $r2_11 = 1;
  my $r2_22 = 1;
  
  printf "%s\tcov:(%.3f,%.3f;%.3f,%.3f)\tcor(%.3f,%.3f;%.3f,%.3f)\tr^2(%.3f,%.3f;%.3f,%.3f)\n",$infile,$var1,$cov12,$cov12,$var2,$r11,$r,$r,$r22,$r2_11,$r2,$r2,$r2_22;
}
else {
  die "\nTry \"$0 -h\" for a complete list of options\n\n";
}

## end of main
## end of script