#!/usr/bin/env perl
use strict;
#use warnings;

=head1 NAME

    MycPermCheck

=head1 VERSION

    1.1

=cut

=head1 LICENSE

    This work is licensed under a Creative Commons Attribution-NoDerivatives 4.0 International License.
    http://creativecommons.org/licenses/by-nd/4.0/

=cut

=head1 DESCRIPTION
    
    MycPermCheck predicts the permeability probability of small molecules against the Mycobacterium tuberculosis cell wall.

=head1 SYNOPSIS
    
    ./mycpermcheck1.1 -i <input file> [-s <sort mode> -c <Calculation Mode>]
    
    sort modes can be n (by name), o (off, order is unchanged) or p (by probability, default)

    calculation mode can be a (all, default), f (first molecule of identically named), m (average of identically named molecule)

=cut

### Check parameters
my %PAR = @ARGV;
if( !exists $PAR{-i} ) {
    print "MycPermCheck 1.1\nUsage: ./mycpermcheck1.1 -i <input file> [-s <sort mode> -c <calculation mode>]\n\nsort modes can be n (by name), o (off, order is unchanged) or p (by probability, default)\n\ncalculation mode can be a (all, default), f (first) or m (mean of isomeric forms)\n\nFor more information use \"perldoc mycpermcheck1.1\"\n\n" and exit;
}

my $isoform = 'a';
$isoform = $PAR{-c} if(exists $PAR{-c} );
if( $isoform ne 'a' and $isoform ne 'f' and $isoform ne 'm' ){
    print "MycPermCheck 1.1\nUsage: ./mycpermcheck1.1 -i <input file> [-s <sort mode> -c <calculation mode>]\n\nsort modes can be n (by name), o (off, order is unchanged) or p (by probability, default)\n\ncalculation mode can be a (all, default), f (first) or m (mean of isomeric forms)\n\nFor more information use \"perldoc mycpermcheck1.1\"\n\n" and exit;
}


my $fname2 = $PAR{-i};
if( $fname2 eq '') {
    print "Please specify file!\n\n" and exit;
}

my $sort = '';
$sort = $PAR{-s} if( exists $PAR{-s} );
if( $sort ne 'p' and $sort ne 'n' and $sort ne 'o')  {
    $sort = 'p';
}


### Define descriptor coefficients and centers
my %COEF_qp = ('PISA'      => [7.831499e-05, 250.7896295],
	       'FOSA'      => [-5.137187e-05, 236.1673680],
	       'QPlogPo/w' => [6.060713e-03, 3.2621865],
	       'accptHB'   => [-3.181890e-03, 5.6655750],
	       'glob'      => [-7.429219e-02, 0.8191964]);

my %COEF_pad = ('C2SP2'    => [-2.951990e-03, 8.1930000],
		'XLogP'    => [-5.353284e-03, 1.8380930],
		'TPSA'     => [5.160086e-05, 92.8078396],
		'HybRatio' => [5.197880e-02, 0.2705913],
		'LOBMAX'   => [-4.130163e-03, 2.1456477]);

=head1 INTERNAL FUNCTIONS

=cut

=over

=item logreg_qp()
    
    Method logreg_qp() applies the MycPermCheck logistic regression function for a given PC1 coordinate input for the QikProp based model.

=cut

sub logreg_qp {
    my $pc1 = $_[0];
    my $prob = 1/(1+exp(-45.187*$pc1));
    return $prob;
}

=item logreg_pad()
    
    Method logreg_pad() applies the MycPermCheck logistic regression function for a given PC1 coordinate input for the PaDEL based model.

=cut

sub logreg_pad {
    my $pc1 = $_[0];
    my $prob = 1/(1+exp(52.943*$pc1));
    return $prob;
}

=item loadandformat()

    Method loadandformat() reads and checks the input file for file format and correctness.

=cut

sub loadandformat {
    my $fname2 = $_[0];

    open( CHECK, "$fname2") or die "File not found!\n\n";
    my $headline = <CHECK>;
    my @HEADER = split(',', $headline);
    
    my $fileformat = "qikprop";
    if ( $headline =~ /Name/ ) {
	$fileformat = "padel";
    }
    elsif( $headline !~ /molecule,/ ) {
	print "$fname2 has wrong file format! Please use QikProp or PaDEL CSV files only.\n" and exit;
    }
    close CHECK or die "Can't close file! Something went wrong.\n\n";
    
    ### Check and save indices of wanted descriptors
    my @DESC = ();
    @DESC = qw(FOSA QPlogPo/w PISA accptHB glob) if( $fileformat eq 'qikprop');
    @DESC = qw(LOBMAX TPSA C2SP2 HybRatio XLogP) if( $fileformat eq 'padel');

    my @DESC_ind = ();
    foreach my $desc ( @DESC ) {
	die "descriptor $desc is missing!\n" if( $headline !~ /$desc/ );
	push(@DESC_ind, grep { $HEADER[$_] =~ /$desc/ } 0..$#HEADER);
    }

    return ($fileformat,\@DESC_ind,\@HEADER,\@DESC);
}

=item noiso()

    Perform MycPermCheck without averaging of isomeric forms.

=cut

sub noiso {
    my $fname2 = $_[0];
    my $fileformat = $_[1];
    my @DESC_ind = @{$_[2]};
    my @HEADER = @{$_[3]};
    
    my @OUTPUT = ();
    my %NAMES = ();
    my %DATA_centered = ();
        
    open( INPUT, "$fname2") or die "File not found!\n\n";
    while( <INPUT> ) {
	next if ($_ =~ /molecule/);
	next if ($_ =~ /Name/);
	my @LINE = split(',', $_);

        if ( $isoform eq 'f' ) {
	    next if ( exists $NAMES{$LINE[0]} );
	}

	next if ( $LINE[1] eq '' );
	
	my %DATA = ($HEADER[0] => $LINE[0],
		    $HEADER[$DESC_ind[0]] => $LINE[$DESC_ind[0]],
		    $HEADER[$DESC_ind[1]] => $LINE[$DESC_ind[1]],
		    $HEADER[$DESC_ind[2]] => $LINE[$DESC_ind[2]],
		    $HEADER[$DESC_ind[3]] => $LINE[$DESC_ind[3]],
		    $HEADER[$DESC_ind[4]] => $LINE[$DESC_ind[4]]);
	
	%DATA_centered = ($HEADER[0] => $LINE[0]);
	foreach( keys %DATA ) {
	    next if ($_ eq 'molecule');
	    next if ($_ eq 'Name');
	    $DATA_centered{$_} = $DATA{$_}-$COEF_qp{$_}[1] if( $fileformat eq 'qikprop');
	    $DATA_centered{$_} = $DATA{$_}-$COEF_pad{$_}[1] if( $fileformat eq 'padel');
	}
	
	######################################################################    
	### MycPermCheck:
	### Calculate PC1 coordinate and Calculate logistic regression
	my $probability = '';
	if( $fileformat eq 'qikprop' ) {
	    my $PC1 = 0.1536226 * 0.6510271 * 9.99875 * ($COEF_qp{PISA}[0]*$DATA_centered{PISA} + $COEF_qp{FOSA}[0]*$DATA_centered{FOSA} + $COEF_qp{'QPlogPo/w'}[0]*$DATA_centered{'QPlogPo/w'} + $COEF_qp{accptHB}[0]*$DATA_centered{accptHB} + $COEF_qp{glob}[0]*$DATA_centered{glob});
	    
	    $probability = &logreg_qp($PC1);
	}
	elsif( $fileformat eq 'padel' ) {
	    my $PC1 = 0.1649452 * 0.6063332 * 9.99875 * ($COEF_pad{C2SP2}[0]*$DATA_centered{C2SP2} + $COEF_pad{XLogP}[0]*$DATA_centered{XLogP} + $COEF_pad{TPSA}[0]*$DATA_centered{TPSA} + $COEF_pad{HybRatio}[0]*$DATA_centered{HybRatio} + $COEF_pad{LOBMAX}[0]*$DATA_centered{LOBMAX});
	    
	    $probability = &logreg_pad($PC1);
	}
	######################################################################
	
	
	### Save result
	push(@OUTPUT, [$LINE[0],$probability,"$LINE[$DESC_ind[0]],$LINE[$DESC_ind[1]],$LINE[$DESC_ind[2]],$LINE[$DESC_ind[3]],$LINE[$DESC_ind[4]]"]);
	
	### Memorize molecule names
	$NAMES{$LINE[0]} = 1;
    }
    close INPUT or die "Can't close file! Something went wrong.\n\n";
    return @OUTPUT;
}

=item iso()

    Perform MycPermCheck with averaging of isomeric forms.

=cut

sub iso {
    my $fname2 = $_[0];
    my $fileformat = $_[1];
    my @DESC_ind = @{$_[2]};
    my @HEADER = @{$_[3]};
    
    my @OUTPUT = ();
    my @NAMES = ();
    my %DATA = ();
    my %DATA_avg = ();
    my %DATA_centered = ();

    open( INPUT, "$fname2") or die "File not found!\n\n";
    while( <INPUT> ) {
	next if ($_ =~ /molecule/);
	next if ($_ =~ /Name/);
	my @LINE = split(',', $_);
	next if ( $LINE[1] eq '' );
	
	my $name = $LINE[0];
	push(@NAMES, $name);
	push(@{$DATA{$name}}, {$HEADER[$DESC_ind[0]] => $LINE[$DESC_ind[0]],
			       $HEADER[$DESC_ind[1]] => $LINE[$DESC_ind[1]],
			       $HEADER[$DESC_ind[2]] => $LINE[$DESC_ind[2]],
			       $HEADER[$DESC_ind[3]] => $LINE[$DESC_ind[3]],
			       $HEADER[$DESC_ind[4]] => $LINE[$DESC_ind[4]]});
 
    }
    close INPUT or die "Can't close file! Something went wrong.\n\n";	
    
    ### Averaging
    foreach my $ligand ( keys %DATA ) {
	my $d1 = 0;
	my $d2 = 0;
	my $d3 = 0;
	my $d4 = 0;
	my $d5 = 0;
	my $nr = (@{$DATA{$ligand}});

	foreach my $iso ( @{$DATA{$ligand}} ) {
	    $d1 += $iso->{$HEADER[$DESC_ind[0]]};
	    $d2 += $iso->{$HEADER[$DESC_ind[1]]};
	    $d3 += $iso->{$HEADER[$DESC_ind[2]]};
	    $d4 += $iso->{$HEADER[$DESC_ind[3]]};
	    $d5 += $iso->{$HEADER[$DESC_ind[4]]};
	}
	$DATA_avg{$ligand} = {$HEADER[$DESC_ind[0]] => $d1/$nr,
			      $HEADER[$DESC_ind[1]] => $d2/$nr,
			      $HEADER[$DESC_ind[2]] => $d3/$nr,
			      $HEADER[$DESC_ind[3]] => $d4/$nr,
			      $HEADER[$DESC_ind[4]] => $d5/$nr};
    }

    my %BEENTHERE = ();

    foreach my $name ( @NAMES ) {
	next if ( exists $BEENTHERE{$name} );
	foreach ( keys %{$DATA_avg{$name}} ) {
	    $DATA_centered{$name}{$_} = $DATA_avg{$name}{$_}-$COEF_qp{$_}[1] if( $fileformat eq 'qikprop');
	    $DATA_centered{$name}{$_} = $DATA_avg{$name}{$_}-$COEF_pad{$_}[1] if( $fileformat eq 'padel');
	}

	my $probability = '';
	if( $fileformat eq 'qikprop' ) {
	    my $PC1 = 0.1536226 * 0.6510271 * 9.99875 * ($COEF_qp{PISA}[0]*$DATA_centered{$name}{PISA} + $COEF_qp{FOSA}[0]*$DATA_centered{$name}{FOSA} + $COEF_qp{'QPlogPo/w'}[0]*$DATA_centered{$name}{'QPlogPo/w'} + $COEF_qp{accptHB}[0]*$DATA_centered{$name}{accptHB} + $COEF_qp{glob}[0]*$DATA_centered{$name}{glob});
	    
	    $probability = &logreg_qp($PC1);
	}
	elsif( $fileformat eq 'padel' ) {
	    my $PC1 = 0.1649452 * 0.6063332 * 9.99875 * ($COEF_pad{C2SP2}[0]*$DATA_centered{$name}{C2SP2} + $COEF_pad{XLogP}[0]*$DATA_centered{$name}{XLogP} + $COEF_pad{TPSA}[0]*$DATA_centered{$name}{TPSA} + $COEF_pad{HybRatio}[0]*$DATA_centered{$name}{HybRatio} + $COEF_pad{LOBMAX}[0]*$DATA_centered{$name}{LOBMAX});
	    
	    $probability = &logreg_pad($PC1);
	}
	######################################################################
	
	### Save result
	push(@OUTPUT, [$name,$probability,"$DATA_avg{$name}{$HEADER[$DESC_ind[0]]},$DATA_avg{$name}{$HEADER[$DESC_ind[1]]},$DATA_avg{$name}{$HEADER[$DESC_ind[2]]},$DATA_avg{$name}{$HEADER[$DESC_ind[3]]},$DATA_avg{$name}{$HEADER[$DESC_ind[4]]}"]);

	$BEENTHERE{$name} = 1;
    }

    return @OUTPUT;
}

=item sorting()

    Provides the selectable sorting methods

=cut

sub sorting {
    if( $sort eq 'p' ) {
	$b->[1] <=> $a->[1] or lc($a->[0]) cmp lc($b->[0])
    }
    elsif( $sort eq 'n' ) {
	lc($a->[0]) cmp lc($b->[0]) or $a->[1] <=> $b->[1]
    }
}


##########################################################################################
### Start main program

my ($fileformat,$DESC_ind,$HEADER,$DESC) = &loadandformat($fname2);
if( $isoform eq 'a' or $isoform eq 'f') {
    my @OUTPUT = &noiso($fname2, $fileformat, $DESC_ind, $HEADER);
    print "Name,Probability," . join(',',@{$DESC}) . "\n";
    foreach( sort sorting @OUTPUT ) {
	print "$_->[0]," . sprintf("%.3f", $_->[1]) . ",$_->[2]" . "\n";
    }
}
elsif( $isoform eq 'm' ) {
    my @OUTPUT = &iso($fname2, $fileformat, $DESC_ind, $HEADER) if( $isoform eq 'm');
    print "Name,Probability," . join(',',@{$DESC}) . "\n";
    foreach( sort sorting @OUTPUT ) {
	print "$_->[0]," . sprintf("%.3f", $_->[1]) . ",$_->[2]" . "\n";
    }
}

##########################################################################################

=back

=head1 TODO

    - "Sort mode off" does not work with "Calculate Mean of Isomeric Forms": FIXED in 1.1

=head1 CITE

    Merget et al. (2013) MycPermCheck: The Mycobacterium tuberculosis permeability prediction tool for small molecules, Bioinformatics, 29(1): 62-68.

=cut

=head1 AUTHOR

    Benjamin Merget, 2014
    
=cut