#!/usr/bin/perl

#plan. Read in the data from Supp3.GTExResults.txt and get list of genes;
#Then go to the HPO file in /net/isi-backup/steve/GENCODYS/NijmegenID/SystematicWithoutSheikh/dataOut, identify the patients with the gene
#Then go to /net/isi-backup/steve/GENCODYS/NijmegenID/SystematicWithoutSheikh/01.getThePatients.dn.out to get the HPO phenotypes for those patients contributing genes

use strict;

my %GTExdata;
my $id = 0;
open(GTExDATA,"GTEx_0.7_Enrichments.out")||die;
my $GTExhead = <GTExDATA>;
my @pvalues = ();
while(<GTExDATA>){
   $id++;
   chomp;
   my @line = split(/\t/);
   my $HPOID=$line[1];
   if ($HPOID eq "") {$HPOID = "blank-$id";}
   my $GTExID = $id;
   $GTExdata{$HPOID}->{$GTExID}->{HPOName} = undef;
   $GTExdata{$HPOID}->{$GTExID}->{Patients} = undef;
   $GTExdata{$HPOID}->{$GTExID}->{Pvalue} = $line[0];
   push (@pvalues, $line[0]);
   $GTExdata{$HPOID}->{$GTExID}->{PatientsContributing} = undef;
   #if(($HPOID eq "HP0004279")&&($GTExID eq "GTEx:0003700 - transcription factor activity")){ die "$GTExdata{$HPOID}->{$GTExID}->{Patients} = $line[2]"; }
   (my $genes = $line[2]) =~ s/"//g;
   my @genes = split(/,/,$genes);
   #@{$GTExdata{$HPOID}->{Genes}} = @genes;
   foreach my $gene (@genes){
      $GTExdata{$HPOID}->{$GTExID}->{Genes}->{$gene}=1;
      }
   $GTExdata{$HPOID}->{$GTExID}->{NoGenes} = scalar(@genes);
   #die "$GTExdata{$HPOID}->{Genes}->[1]";
   }

my $FDRthresh = Benjamini_Hochberg_05(@pvalues);

#my %pats2phenos;
#open(PATPHENOS,"/net/isi-backup/steve/GENCODYS/NijmegenID/SystematicWithoutSheikh/01.getThePatients.dn.out")||die;
#while(<PATPHENOS>){
#   chomp;
#   my @line = split(/\t/);
#   my @subphenos = split(/,/,$line[1]);
#   my $pats = split(/,/,$line[2]);
#   foreach my $pat (@pats){
#      $pats2phenos{$pat}->{
#not going to use this as I'm confused about the child terms Steve lists. I'm going to go straight from the data, I presume from Anneke.

#/net/isi-backup/steve/GENCODYS/NijmegenID/qry250Karray_clinical_data.txt
#Read in the patient phenotypes
my %pats2phenos;
open(PATPHENOS,"/net/isi-backup/steve/GENCODYS/NijmegenID/qry250Karray_clinical_data.txt")||die;
my $head = <PATPHENOS>;
while(<PATPHENOS>){
   chomp;
   my @line = split;
   $pats2phenos{$line[0]}->{$line[1]}=1;
   }

#Now, for each HPO in the enrichments, open up the CNV file, and associate the genes with the patients
#DID STEVE USE COMPLETE OVERLAPS??
#I think he used partials minus Shaikhs. I'm going to read the gene files and remove patients that don't have any genes in the final gene file
#my @GENEFILES = </net/isi-backup/steve/GENCODYS/NijmegenID/SystematicPartials/dataOut/DenovoCNVs_all_hg18.ensembl*partial_noPartialShaikh>;
#my @GENEFILES = </net/isi-backup/steve/GENCODYS/NijmegenID/Systematic/dataOut/
#dataOut? dataOut2?

#my @CNVGFILES = </net/isi-backup/steve/GENCODYS/NijmegenID/Systematic/dataOut/DenovoCNVs_all_hg18.ensembl.*complete_noCompleteShaikh>;
my @CNVGFILES = </net/isi-scratch/steve/SystematicOnlySmallPartials/dataOut/DenovoCNVs_all_hg18.ensembl.*partial_noCompleteShaikh>;
#first, let's go through the complete_noCompleteShaikh files to get the genes that were selected after final filtering
my %CNVgenes;
foreach my $file (@CNVGFILES){
   #if ($file =~ /complete/){ next; }
   my ($rest1,$hpo,$rest2) = split(/.ensembl.|.partial/,$file);
   #print "$file :: $hpo\n";
   #die;
   open(GENES,$file)||die;
   while(<GENES>){
      my @line = split(/\t/);
      $CNVgenes{$hpo}->{$line[0]}=1;
      #die $line[0];
      }
   close(GENES);
   }
 

my @CNVFILES = </net/isi-scratch/steve/SystematicOnlySmallPartials/dataOut/DenovoCNVs_all_hg18.ensembl.*>;
my %CNVGENE2PAT;
my %PAT2CNVGENE;
my %HPO_CNVPATS;
foreach my $file (@CNVFILES){
   if ($file =~ /partial|linksToHPOGenes/){ next; }
   my ($rest1,$hpo,$rest2) = split(/.ensembl.|.complete/,$file);
   #print "$file :: $hpo\n";
   open(CNVS,$file)||die;
   my $pat="";
   while(<CNVS>){
      chomp;
      if(/^ENS/){
         my @line = split(/_/);
         if($pat eq ""){ die "$file\nno pat :: $_"; }
         if(!exists $CNVgenes{$hpo}->{$line[0]}){ next; }
         $CNVGENE2PAT{$line[0]}->{$pat}=1;
         $PAT2CNVGENE{$pat}->{$line[0]}=1;
         #die "$file\nPAT2CNVGENE{$pat}->{$line[0]}";
         }
      else{
         my @line = split(/\t/);
         $pat = $line[4];
         $HPO_CNVPATS{$hpo}->{$pat}=1;
         #die "HPO_CNVPATS{$hpo}->{$pat}";
         }
      }
   close(CNVS);
   }

#OK, all data in so now form sets of patients with the HPO that contribute to the enrichment and those that don't
#@{$GTExdata{$HPOID}->{$GTExID}->{Genes}}
#$pats2phenos{$PAT}->{HPO}
#$HPO_CNVPATS{$hpo}->{$pat}
#$CNVGENE2PAT{ensg}->{$pat}

#I also want to note how many genes in the enrichment each patient has
foreach my $HPO (keys %GTExdata){
   #go through each GTEx enrichment for this HPO term
   foreach my $GTEx (keys %{$GTExdata{$HPO}}){
      #go through each of the patients considered
      my @pats;
      my $CONTRIBPATS = 0;
      if (scalar(keys(%{$HPO_CNVPATS{$HPO}})) < 3) {next;}
      if (!($GTExdata{$HPO}->{$GTEx}->{Pvalue} < $FDRthresh)) {next;}
      foreach my $pat (keys %{$HPO_CNVPATS{$HPO}}){
         my $pathasnogenes=0;
         #print "$pat\n";
         #check this pat has CNV genes
         my $numgenesforpat = keys %{$PAT2CNVGENE{$pat}};
         #print "$numgenesforpat\n";
         if(!$numgenesforpat){
            #print "PAT HAS NO GENES :: HPO $HPO :: pat $pat\n";
            $pathasnogenes=1;
            #next;
            }
         #how many genes in the enrichment does this patient's CNV(s) hit?
         my $ENRICH_GCOUNT=0;
         my @ENRICH_GENES;
         foreach my $patgene (keys %{$PAT2CNVGENE{$pat}}){
            #die "GTExdata{$HPO}->{$GTEx}->{Genes}->{$patgene}";
            if(exists $GTExdata{$HPO}->{$GTEx}->{Genes}->{$patgene}){
               $ENRICH_GCOUNT++;
               push @ENRICH_GENES,$patgene;
               }
            }
         #add this patient
         if($pathasnogenes){ $ENRICH_GCOUNT = -1; }
         push @pats,"$pat($ENRICH_GCOUNT;@ENRICH_GENES)";
         if($ENRICH_GCOUNT>0){ $CONTRIBPATS++; }
         }
#      if($CONTRIBPATS != $GTExdata{$HPO}->{$GTEx}->{PatientsContributing}){ die "GTExdata{$HPO}->{$GTEx}->{PatientsContributing} $GTExdata{$HPO}->{$GTEx}->{PatientsContributing} CONTRIBPATS $CONTRIBPATS"; }
 #     my $numpats = @pats;
 #     if($numpats != $GTExdata{$HPO}->{$GTEx}->{Patients}){ die "GTExdata{$HPO}->{$GTEx}->{Patients} $GTExdata{$HPO}->{$GTEx}->{Patients} numpats $numpats\n@pats "; }
      print "$HPO\t$GTEx\t";
      print join(",",@pats)."\n";
      } #closes GTEx
   } #closes HPO


# 1995 paper
sub Benjamini_Hochberg_05 {
	my @pvals = sort{$a <=> $b} @_;
	my $m = scalar(@pvals);
	my $thresh = 0;
	for (my $i = 0; $pvals[$i] <= $i/$m*0.05; $i++) {
		$thresh = $pvals[$i];
	}
	return ($thresh);
}
