countStats.pl

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

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

use strict;
use warnings;
use Getopt::Long;
no autovivification;
use List::Util qw[min max sum];
use Scalar::Util qw(looks_like_number);  

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

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

&GetOptions(
    "help!"              => \$help,             # print this help
    "missIntInDTUfile=s" => \$missIntInDTUfile, # missingInteractionsInDTU_filt.tsv.gz
) 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";
}
##############################################################################
sub median {
    my @values = sort{$a<=>$b} @_;

    my $med = 0;
    # if odd number of elements take middle element. Array size is given in size-1
    if(@values%2==1){
	$med = $values[int(@values/2)];
    }
    # if even number of elements take average of middle two elements
    else{
	my $n= @values;	
	$med = ($values[int(@values/2)-1] + $values[int(@values/2)]) / 2;
    }
    return $med;
}
##############################################################################
sub mean {
    my $avg = sum(@_)/@_;
    return $avg;
}

##############################################################################
sub std {
    my $avg = &mean(@_);

    my $std = 0;
    foreach my $v (@_){
	$std += ($v-$avg)**2;
    }
    @_ = (1,1) if(@_ == 1);
    $std = sqrt($std/(@_-1));
    
    return $std;
}
##############################################################################
# median absolute deviation
sub mad {
    my $med = &median(@_);
    
    my @stdDiffs = ();
    foreach my $v (@_){
	push(@stdDiffs, abs($v-$med));
    }
        
    my $mad = &median(@stdDiffs);
    return $mad;
}
##############################################################################
### END OF SUBROUTINES########################################################
############################################################################## 


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

my %dutEnst = ();
open(F, "zcat $missIntInDTUfile|") or die "\nERROR: Failed to open $missIntInDTUfile: $!\n\n";
while(<F>) {
    chomp($_);
    my @a = split(/\t/);
    $dutEnst{$a[3]} = 1;
}
close(F);

my %data = ();
my %dut = ();
my %colNames = ();
print STDERR "Reading in data ... ";
my $header = <>;
my @colNames = split(/\t/, $header);
while(<>) {
    my @a = split(/\t/);
    for(my $i = 15; $i < @a-1; $i++) {
	if(!exists $data{$i}) {
	    my @v = ($a[$i]);
	    $data{$i} = \@v;
	    $dut{$i}  = \@v if(exists $dutEnst{$a[1]});
	} else {
	    push(@{$data{$i}}, $a[$i]);
	    push(@{$dut{$i}},  $a[$i]) if(exists $dutEnst{$a[1]});
	}
    }
}
print STDERR "done\n";

print STDERR "Calculating stats ... ";
print "#Type\tSampleId\tMin\tMax\tMean\tMedian\tStd\tMad\tMeanMedianDiff\tRelMeanMedianDiff\n";
foreach my $col (sort {$a<=>$b} keys %data) {
    my @values = @{$data{$col}};
    my $mean = &mean(@values);
    my $med = &median(@values);
    my $std = &std(@values);
    my $mad = &mad(@values);
    my $min = min(@values);
    my $max = max(@values);
    printf("ALL\t%s\t%i\t%i\t%.3f\t%i\t%.3f\t%i\t%.3f\t%.3f\n", $colNames[$col], $min, $max, $mean, $med, $std, $mad, abs($mean-$med), (min($mean, $med) == 0 ? "99999999" : (abs($mean-$med)/min($mean, $med))));
}
foreach my $col (sort {$a<=>$b} keys %dut) {
    my @values = @{$dut{$col}};
    my $mean = &mean(@values);
    my $med = &median(@values);
    my $std = &std(@values);
    my $mad = &mad(@values);
    my $min = min(@values);
    my $max = max(@values);
    printf("DUT\t%s\t%i\t%i\t%.3f\t%i\t%.3f\t%i\t%.3f\t%.3f\n", $colNames[$col], $min, $max, $mean, $med, $std, $mad, abs($mean-$med), (min($mean, $med) == 0 ? "99999999" : (abs($mean-$med)/min($mean, $med))));
}
print STDERR "done\n";
print STDERR "ALL DONE\n";