#!/usr/bin/perl

use strict;
use lib "/net/isi-backup/caleb/new_node_stuff/my_perl_stuff";
use cwstats_fgu202;

#what are the distributions we need?
#those genes within KEGG pathways|those genes in this pathway|those KEGG pathways genes in these CNVs|those in this pathway in these CNVs

if(@ARGV!=4){ die "do_stats.perl [MR gene file] [outfile] [sig cut off for FDR] [sig cut off for hypergeo]\n"; }

open(DOUT,">$ARGV[1]")||die;

#so, the parolap files contain the NCBI genes overlapping the CNVs so we can use those as long as we restrict to
#those KEGG pathway genes
my %kpathgenes;
open(KGENES,"/net/isi-backup/caleb/new_node_stuff/kegg/map_kegg_genes2ncbi-genesandpos.out.paths")||die;
while(<KGENES>){
   if(/path/){
      my ($gene,$paths)=split(/\t/);
      $gene =~ s/hsa://;
      $kpathgenes{$gene}=1;
      }
    }

#now get the genes in the neuropaths
my @files = </net/isi-backup/caleb/new_node_stuff/kegg/neurodegenerative_pathways/*genes>;
my %npathgenes;
foreach my $file (@files){
   (my $path = $file) =~ s/\/net\/isi-backup\/caleb\/new_node_stuff\/kegg\/neurodegenerative_pathways\/|.genes$//g;
   #die "$path";
   open(GLIST,$file)||die;
   while(<GLIST>){
      my ($gene,$ps)=split(/\t/);
      $gene =~ s/hsa://;
      #die "path $path gene $gene";
      if(!exists $kpathgenes{$gene}){ die "problem path $path gene $gene"; }
      $npathgenes{$path}->{$gene}=1;
      $npathgenes{all}->{$gene}=1;
      }
   close GLIST;
   }


#now get the genes from the MR regions and see if any of these are neuro genes
my %ncounts;
my $kpathgcnt=0;
my %npathgcnt;
my %npathgcntgenes;
my $MRgcnt=0;
open(MRIN,$ARGV[0])||die;
while(<MRIN>){
   if(/^REG/){}
      chomp;
      $MRgcnt++;
      my @line = split(/\t/);
      my $gene = $line[0];
      #is this a KEGG gene? Is this a neuro gene?
      if(exists $kpathgenes{$gene}){
         $kpathgcnt++;
         foreach my $npath (keys %npathgenes){
            #die $npath;
            if(exists $npathgenes{$npath}->{$gene}){
               $npathgcnt{$npath}++;
               $npathgcntgenes{$npath}->{$gene}=1;
               }
            } #go through each npath
         } #closes if kpathgene
   }



#now need to do the stats
#(num of genes with MP,num of genes with this MP,num of MP genes in reg, num of genes in reg with this MP)
#So this is asking, Of those genes with phenotypes in our regions, how over represented is a particular phenotype, given all genes
#with phenotypes amongst which there are so many with this particular phenotype
#the question for me is what should the background be? I'm going to take all genes with phenotypes, but perhaps I need to ask the
#question of whether genes giving phenotypes are over represented as a whole in these regions.


print DOUT " *************  ALL TERMS **************\n\n";
my @stats;
foreach my $path (sort keys %npathgenes){
   if(!defined $npathgenes{$path}){ $npathgenes{$path}=0; }
   my $allkpathgcnt = keys %kpathgenes;
   my $allnpathgcnt = keys %{$npathgenes{$path}};
   #die "$path do_hyperGlogcalc_this_or_greater($allkpathgcnt,$allnpathgcnt,$kpathgcnt,$npathgcnt{$path}";
   my $pvalG = sprintf("%2.4f",do_hyperGlogcalc_this_or_greater($allkpathgcnt,$allnpathgcnt,$kpathgcnt,$npathgcnt{$path}));
   my $pvalF = sprintf("%2.4f",do_hyperGlogcalc_this_or_fewer($allkpathgcnt,$allnpathgcnt,$kpathgcnt,$npathgcnt{$path}));
   my $statcomment;
   if($pvalG < $pvalF){
      push @stats,[$pvalG,"MANY",$path,$allkpathgcnt,$allnpathgcnt,$kpathgcnt,$npathgcnt{$path}];
      if($pvalG < $ARGV[3]){ $statcomment="SIG MANY"; }
      }
   else{
      push @stats,[$pvalF,"FEW",$path,$allkpathgcnt,$allnpathgcnt,$kpathgcnt,$npathgcnt{$path}];
      if($pvalF < $ARGV[3]){ $statcomment="SIG FEW"; }
      }
   if(($pvalG < $ARGV[3])||($pvalF < $ARGV[3])){
      print DOUT "$path Gpv $pvalG Fpv $pvalF ToTG $allkpathgcnt ToTP $allnpathgcnt RegG $kpathgcnt RegP $npathgcnt{$path} $statcomment\n";
      print DOUT "GENES hit for $path - ".join(" : ", sort {$a <=> $b} keys %{$npathgcntgenes{$path}})."\n";
      print DOUT "*******************************************************\n";
      }
   }

my ($pvalcut,$numsigterms,$fdr) = do_FDR(\@stats,$ARGV[2]);
print DOUT "\n\n *************     FDR    **************\n\n";
print DOUT "FDR pvalcut $pvalcut numsigterms $numsigterms fdr $fdr\n";
print DOUT "Prob\tRepresent\tPhID\tDes\tToTG\tToTPh\tRegG\tRegPh\n";
foreach my $ph (@stats){
   if($ph->[0] <= $pvalcut){ print DOUT join("\t",@{$ph})."\n"; }
   }
print DOUT "\n\n *************     FDR - END TOP LEVEL   **************\n\n";
#
#my ($pvalcutall,$numsigtermsall,$fdrall) = do_FDR(\@statsall,$ARGV[2]);
#print DOUT " *************  FDR - ALL **************\n\n";
#print DOUT "FDR pvalcut $pvalcutall numsigterms $numsigtermsall fdr $fdrall\n";
#print DOUT "Prob\tRepresent\tPhID\tDes\tToTG\tToTPh\tRegG\tRegPh\n";
#foreach my $ph (@statsall){
#   if($ph->[0] <= $pvalcut){ print DOUT join("\t",@{$ph})."\n"; }
#   }
#print DOUT "\n\n *************     FDR - END ALL   **************\n\n";
#
#
#
