#!/usr/bin/perl

use strict;
use warnings;
use Data::Dumper;
use Getopt::Long;

my @header  = ();
my %data_p  = ();
my %data_s  = ();
my %data_n  = ();
my %data_nn = ();
my $sd_th   = 2; # cutoff to remove flattly distributed genes

while(<>){
  chomp;
  
  if ($.==1){
    my @F=split"\t",$_;
    @header=@F;
    next;
  }
 
  my @F=split"\t",$_;
  my $last_column=$#F;
  foreach my $p (2 .. $last_column){
    $data_p{$F[0]}->{$header[$p]}=$F[$p];
    $data_s{$header[$p]}+=$F[$p];
  }

}

foreach my $gene (keys %data_p) {
  foreach my $sample (keys %{$data_p{$gene}}){
 
    if($data_s{$sample} && $data_p{$gene}->{$sample}){
      $data_n{$gene}->{$sample} = $data_p{$gene}->{$sample};
    }
    else{
      $data_n{$gene}->{$sample}=0;
    }
  }
}

my $removed = 0;
my $totals = 0;
foreach my $gene (keys %data_n) {
  my @values = ();

  foreach my $sample (keys %{$data_n{$gene}}){
    push @values, $data_n{$gene}->{$sample};
  }

  #remove genes with a flat distribution over all samples
  $totals++;
  if( &sd(@values)<= $sd_th ) {
    $removed++;
    next;
  }  
  
  foreach my $sample (keys %{$data_n{$gene}}){
    $data_nn{$gene}->{$sample} = ( $data_n{$gene}->{$sample} - &mean(@values) )/( &sd(@values) );
  }
}

print STDERR "$removed genes from $totals removed due to SD < $sd_th\n";
my $h_c = 1;
foreach my $gene (sort keys %data_nn){
  if($h_c){
    my @h_val = ();
    foreach my $sample (sort keys %{$data_nn{$gene}}){
      push @h_val, $sample;
    }
    print join("\t",'ID', 'Desc', @h_val)."\n";
    $h_c=0;
  }
  
  print "$gene\t$gene\t";
  my @out_val = ();
  foreach my $sample (sort keys %{$data_nn{$gene}}){
    push @out_val, $data_nn{$gene}->{$sample};
  }
  print join("\t",@out_val)."\n";
}


# sub routines
sub mean{
  my @val=@_;
  if(@val == 1){
    return 'NA';
  }
  my $total = 0;
  foreach (@val) {
    $total += $_;
  }
  return $total/@val;
}

sub sd {
  my(@val) = @_;
  if(@val == 1){
    return 'NA';
  }
  my $average = &mean(@val);
  my $sqtotal = 0;
  foreach(@val) {
    $sqtotal += ($average-$_) ** 2;
  }
  my $std = ($sqtotal / (@val-1)) ** 0.5;
  return $std;
}
