#!/usr/local/bin/perl

#############################################################################
# Copyright (C) 2008, Matrix Science Limited.                               #
#                                                                           #
# This script is free software. You can redistribute and/or                 #
# modify it under the terms of the GNU General Public License               #
# as published by the Free Software Foundation; either version 2            #
# of the License or, (at your option), any later version.                   #
#                                                                           #
# These modules are distributed in the hope that they will be useful,       #
# but WITHOUT ANY WARRANTY; without even the implied warranty of            #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the              #
# GNU General Public License for more details.                              #
#############################################################################

  use strict;

# requires Mascot Parser 2.2 or later
  use lib "../bin";
  use msparser 2.002;

# validate command line arguments  
  if (!defined($ARGV[4])) { 
    &usage();
    exit 1;
  }
  my $filename = $ARGV[0];
  my $thresh_type = lc($ARGV[1]);
  unless ($thresh_type eq "identity" || $thresh_type eq "homology") {
    die "thresh_type must be either identity or homology";
  }
  my $start_val = $ARGV[2] + 0;
  my $end_val = $ARGV[3] + 0;
  my $num_vals = int($ARGV[4] + 0);
  if ($start_val < 1E-18) {
    die "start cannot be less than 1E-18";
  }
  if ($end_val > $start_val) {
    unless ($num_vals > 1) {
      die "num_vals must be at least 2";
    }
  } elsif ($end_val == $start_val) {
    $num_vals = 1;
  } else {
    die "start cannot be greater than end";
  }

# create Mascot Parser results summary object   
  my $objResFile = new msparser::ms_mascotresfile($filename);
  unless ($objResFile->isValid) {
    die $objResFile->getErrorNumber(1);
  }
  unless ($objResFile->anyPeptideSummaryMatches($msparser::ms_mascotresfile::SEC_PEPTIDES)) {
    die "this script only works for MS/MS";
  }
  my $msresFlags = $msparser::ms_mascotresults::MSRES_MAXHITS_OVERRIDES_MINPROB
              | $msparser::ms_mascotresults::MSRES_MUDPIT_PROTEIN_SCORE;
  my $objSummary = new msparser::ms_peptidesummary(
    $objResFile, 
    $msresFlags,
    0.05,
    1,
    "",
    0,
    5
  );
  unless ($objResFile->isValid) {
    die $objResFile->getErrorNumber(1);
  }

# create a Mascot Parser searchparams object
  my $objParams = new msparser::ms_searchparams($objResFile);

# tabulate values at geometrically spaced intervals
  my(@expect,           # expect value
    @score_diff,        # delta between threshold score and p<0.05 value
    @Ftpeps,            # count of matches found in target database only
    @Rtpeps,            # count of matches found in decoy database only
    @Btpeps,            # count of matches found in both databases (concatenated only)
    );
  my $factor = 1;
  if ($num_vals > 1) {
    $factor = ($end_val / $start_val) ** (1 / ($num_vals - 1));
  }
  for (my $i = $num_vals; $i >= 1; $i--) {
    my $ex = $end_val / $factor ** ($num_vals - $i);
    push @expect, $ex;
    push @score_diff, 10 * log($ex / 0.05) / log(10);
  }

# loop through all queries (spectra)
  my $matchCount = 0;   # count of matched queries (only for target database if Mascot automatic decoy)
  for (my $queryNum = 1; $queryNum <= $objResFile->getNumQueries(); $queryNum++) {
    my $objQuery = new msparser::ms_inputquery($objResFile, $queryNum);
    my($threshold, $score);
    if ($thresh_type eq "identity") {
      $threshold = $objSummary->getPeptideIdentityThreshold($queryNum, 20);
    } elsif ($thresh_type eq "homology") {
      $threshold = $objSummary->getHomologyThreshold($queryNum, 20);
      unless ($threshold) {
      # if no homology threshold, then use identity threshold
        $threshold = $objSummary->getPeptideIdentityThreshold($queryNum, 20);
      }
    }
    my $objPeptide = $objSummary->getPeptide($queryNum, 1);   # only consider top-ranked matches
    if (defined($objPeptide) && $objPeptide->getAnyMatch()) {
      $matchCount++;
      $score = $objPeptide->getIonsScore();
      if (($score + $score_diff[0]) >= $threshold) {
        if ($objParams->getDECOY() > 0) {
        # Mascot automatic decoy
        # have a match at the highest expect value in target database
        # increment table rows for each expect value
          for (my $i = 0; $i <= $#expect; $i++) {
            if (($score + $score_diff[$i]) >= $threshold) {
              $Ftpeps[$i]++;
            }
          }
        } else {
        # manual concatenated decoy
        # have a match at the highest expect value, get list of proteins 
          my $tempString = $objResFile->getSectionValueStr($msparser::ms_mascotresfile::SEC_PEPTIDES, "q$queryNum\_p1");
          my($left, $right) = split(/;/, $tempString, 2);
          my @temp_prots = $right =~ /"(.+?)"/g;
        # tally up whether match is to target, decoy, or both
          my $numF = 0;
          my $numR = 0;
          foreach my $accession (@temp_prots) {
            if ($accession =~ /_REVERSE/ 
              || $accession =~ /###REV###/
              || $accession =~ /_RANDOM/
              || $accession =~ /###RND###/) {
              $numR++;
            } else {
              $numF++;
            }
          }
        # increment table rows for each expect value
          for (my $i = 0; $i <= $#expect; $i++) {
            if (($score + $score_diff[$i]) >= $threshold) {
              if ($numR && $numF) {
              # match found in both target and decoy
                $Btpeps[$i]++;
              } else {
                if ($numR) {
                  $Rtpeps[$i]++;
                } else {
                  $Ftpeps[$i]++;
                }
              }
            }
          }
        }
      }
    }
  }

# get decoy matches from Mascot automatic decoy
  if ($objParams->getDECOY() > 0) {
    $msresFlags = $msparser::ms_mascotresults::MSRES_MAXHITS_OVERRIDES_MINPROB
                | $msparser::ms_mascotresults::MSRES_MUDPIT_PROTEIN_SCORE
                | $msparser::ms_mascotresults::MSRES_DECOY;
    $objSummary = new msparser::ms_peptidesummary(
      $objResFile, 
      $msresFlags,
      0.05,
      1,
      "",
      0,
      5
    );
    unless ($objResFile->isValid) {
      die $objResFile->getErrorNumber(1);
    }
    for (my $queryNum = 1; $queryNum <= $objResFile->getNumQueries(); $queryNum++) {
      my $objQuery = new msparser::ms_inputquery($objResFile, $queryNum);
      my($threshold, $score);
      if ($thresh_type eq "identity") {
        $threshold = $objSummary->getPeptideIdentityThreshold($queryNum, 20);
      } elsif ($thresh_type eq "homology") {
        $threshold = $objSummary->getHomologyThreshold($queryNum, 20);
        unless ($threshold) {
        # if no homology threshold, then use identity threshold
          $threshold = $objSummary->getPeptideIdentityThreshold($queryNum, 20);
        }
      }
      my $objPeptide = $objSummary->getPeptide($queryNum, 1);   # only consider top-ranked matches
      if (defined($objPeptide) && $objPeptide->getAnyMatch()) {
        $score = $objPeptide->getIonsScore();
        for (my $i = 0; $i <= $#expect; $i++) {
          if (($score + $score_diff[$i]) >= $threshold) {
            $Rtpeps[$i]++;
          } else {
            last;
          }
        }
      }
    }
  }
    
# print table in ascending expect order
  print "\n\"FDR from target-decoy search for matches exceeding $thresh_type threshold\"\n";
  print "\"($matchCount matches with non-zero score from " . $objResFile->getNumQueries() . " spectra)\"\n\n";
  print "\"Calc. Expect.\"\t\"Actual FDR\"\t\"Target\"\t\"Decoy\"\t\"Both\"\n";
  for (my $i = $#expect; $i >= 0; $i--) {
    my $FDR;
    if ($objParams->getDECOY() > 0) {
      $FDR = $Rtpeps[$i] / $Ftpeps[$i];
    } else {
      $FDR = 2 * $Rtpeps[$i] / ($Ftpeps[$i] + $Rtpeps[$i]);
    }
    print sprintf("%.4g", $expect[$i]) 
      . "\t" 
      . sprintf("%.4f", $FDR)
      . "\t$Ftpeps[$i]\t$Rtpeps[$i]\t$Btpeps[$i]\n";
  }
      
  exit 0;

sub usage() {
  print "\nTabulate target-decoy FDR data for plotting ROC curves, etc.\n\n";
  print "If using a concatenated database, decoy accession strings\n";
  print "must include _REVERSE or ###REV### or _RANDOM or ###RND###\n\n";
  print "The program must be run from the mascot cgi directory\n\n";
  print "Usage:   fdr_table.pl file thresh_type start end num_vals\n\n";
  print "Example: fdr_table.pl ../data/20060503/F123456.dat identity 0.01 10 20\n\n";
  print "         file is the path to a Mascot results file\n";
  print "         thresh_type is either identity or homology\n";
  print "         start is the lowest calculated expectation value\n";
  print "         end is the highest calculated expectation value\n";
  print "         num_vals is the number of rows to output, incl. start and end\n";
}
