Condel.pm

=head1 LICENSE
                                                                                                                     
 Copyright (c) 1999-2011 The European Bioinformatics Institute and                                                   
 Genome Research Limited.  All rights reserved.                                                                      
                                                                                                                     
 This software is distributed under a modified Apache license.                                                       
 For license details, please see

   http://www.ensembl.org/info/about/code_licence.html                                                               
                                                                                                                     
=head1 CONTACT                                                                                                       

 Graham Ritchie <grsr@ebi.ac.uk>
    
=cut

=head1 NAME

 Condel

=head1 SYNOPSIS

 mv Condel.pm ~/.vep/Plugins
 perl variant_effect_predictor.pl -i variations.vcf --plugin Condel,/path/to/config/Condel/config,b

=head1 DESCRIPTION

 This is a plugin for the Ensembl Variant Effect Predictor (VEP) that calculates
 the Consensus Deleteriousness (Condel) score (1) for a missense mutation 
 based on the pre-calculated SIFT (2) and PolyPhen-2 (3) scores from the Ensembl 
 API (4). It adds one new entry class to the VEP's Extra column, Condel which is
 the calculated Condel score. Condel is developed by the Biomedical Genomics Group 
 of the Universitat Pompeu Fabra, at the Barcelona Biomedical Research Park 
 (bg.upf.edu/Condel). The code in this plugin is based on a script provided by this 
 group and slightly reformatted to fit into the Ensembl API.

 The plugin takes 3 command line arguments, the first is the path to a Condel 
 configuration directory which contains cutoffs and the distribution files etc., 
 the second is either "s", "p", or "b" to output the Condel score, prediction or 
 both (the default is both), and the third argument is either 1 or 2 to use the 
 original version of Condel (1), or the newer version (2) - 2 is the default and 
 is recommended to avoid false positive predictions from Condel in some 
 circumstances.

 An example Condel configuration file and a set of distribution files can be found 
 in the config/Condel directory in this repository. You should edit the 
 config/Condel/config/condel_SP.conf file and set the 'condel.dir' parameter to
 the full path to the location of the config/Condel directory on your system.

 References:

 (1) Gonzalez-Perez A, Lopez-Bigas N. 
     Improving the assessment of the outcome of non-synonymous SNVs with a Consensus deleteriousness score (Condel)
     Am J Hum Genet 88(4):440-449 (2011)
     doi:10.1016/j.ajhg.2011.03.004

 (2) Kumar P, Henikoff S, Ng PC.
     Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm
     Nature Protocols 4(8):1073-1081 (2009)
     doi:10.1038/nprot.2009.86
 
 (3) Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, Kondrashov AS, Sunyaev SR. 
     A method and server for predicting damaging missense mutations
     Nature Methods 7(4):248-249 (2010)
     doi:10.1038/nmeth0410-248
 
 (4) Flicek P, et al.
     Ensembl 2012
     Nucleic Acids Research (2011)
     doi: 10.1093/nar/gkr991

=cut

package Condel;

use strict;
use warnings;

use Bio::EnsEMBL::Variation::Utils::BaseVepPlugin;

use base qw(Bio::EnsEMBL::Variation::Utils::BaseVepPlugin);

sub version {
    return '2.4';
}

sub feature_types {
    return ['Transcript'];
}

sub get_header_info {
    return {
        Condel => "Consensus deleteriousness score for an amino acid substitution based on SIFT and PolyPhen-2",
    };
}

sub new {
    my $class = shift;
    
    my $self = $class->SUPER::new(@_);

    # parse the config file and distribution files 

    my $config_dir = $self->params->[0];
    $self->{output} = $self->params->[1] || 'b';
    $self->{version} = $self->params->[2] || 2;
    
    my $config_file = "$config_dir/condel_SP.conf";

    open(CONF, "<$config_file") || die "Could not open $config_file";

    my @conf = <CONF>;
    
    my $safe_conf = 0;
    
    for (my $i = 0; $i < @conf; $i++){
        if ($conf[$i] =~ /condel\.dir=\'(\S+)\'/){
            $self->{config}->{'condel.dir'} = $1;
            $safe_conf++;
        }
        elsif ($conf[$i] =~ /(cutoff\.HumVar\.\w+)=\'(\S+)\'/){
            $self->{config}->{$1} = $2;
            $safe_conf++;
        }
        elsif ($conf[$i] =~ /(max\.HumVar\.\w+)=\'(\S+)\'/){
            $self->{config}->{$1} = $2;
            $safe_conf++;
        }
    }

    if ($safe_conf < 3){
        die "Malformed config file!!!\n\n";
    }

    open(SIFT, $self->{config}->{'condel.dir'}."/methdist/sift.data");
    my @sift = <SIFT>;
    close SIFT;

    for (my $i = 0; $i < @sift; $i++){
        if ($sift[$i] =~ /(\S+)\s+(\S+)\s+(\S+)/){
            $self->{sift}->{tp}->{$1} = $2;
            $self->{sift}->{tn}->{$1} = $3;
        }
    }

    open(POLYPHEN, $self->{config}->{'condel.dir'}."/methdist/polyphen.data");
    my @polyphen = <POLYPHEN>;
    close POLYPHEN;

    for (my $i = 0; $i < @polyphen; $i++){
        if ($polyphen[$i] =~ /(\S+)\s+(\S+)\s+(\S+)/){
            $self->{polyphen}->{tp}->{$1} = $2;
            $self->{polyphen}->{tn}->{$1} = $3;
        }
    }

    return $self;
}

sub run {
    my ($self, $tva) = @_;
    
    my $pph_score   = $tva->polyphen_score;
    my $pph_pred    = $tva->polyphen_prediction;
    my $sift_score  = $tva->sift_score;

    my $results = {};

    if (defined $pph_score && defined $sift_score && $pph_pred ne 'unknown') {

        my ($condel_pred, $condel_score) = $self->compute_condel($sift_score, $pph_score);

        $condel_score = sprintf "%.3f", $condel_score;
            
        my $output = "$condel_pred($condel_score)";
        
        $output = $condel_pred if ($self->{output} =~ /^p/);
        $output = $condel_score if ($self->{output} =~ /^s/);

        $results = {
            Condel => $output,
        };
    }

    return $results;
}

sub compute_condel {

    my ($self, $sift_score, $polyphen_score) = @_;

    my $USE_V2 = $self->{version} == 2;

    my $class;
    
    my $base = 0;
    my $int_score = 0;

    $sift_score     = sprintf("%.3f", $sift_score);
    $polyphen_score = sprintf("%.3f", $polyphen_score);
    
    if ($sift_score <= $self->{config}->{'cutoff.HumVar.sift'}){
        $int_score += sprintf("%.3f", (1 - $sift_score/$self->{config}->{'max.HumVar.sift'})*(1-$self->{sift}->{tn}->{$sift_score}));
        $base += $USE_V2 ? 1 : 1-$self->{sift}->{tn}->{$sift_score};
    }
    else {
        $int_score += sprintf("%.3f", (1 - $sift_score/$self->{config}->{'max.HumVar.sift'})*(1-$self->{sift}->{tp}->{$sift_score}));
        $base += $USE_V2 ? 1 : 1-$self->{sift}->{tp}->{$sift_score};
    }
    
    if ($polyphen_score >= $self->{config}->{'cutoff.HumVar.polyphen'}){
        $int_score += sprintf("%.3f", $polyphen_score/$self->{config}->{'max.HumVar.polyphen'}*(1-$self->{polyphen}->{tn}->{$polyphen_score}));
        $base += $USE_V2 ? 1 : 1-$self->{polyphen}->{tn}->{$polyphen_score};
    }
    else {
        $int_score += sprintf("%.3f", $polyphen_score/$self->{config}->{'max.HumVar.polyphen'}*(1-$self->{polyphen}->{tp}->{$polyphen_score}));
        $base += $USE_V2 ? 1 : 1-$self->{polyphen}->{tp}->{$polyphen_score};
    }

    if ($base == 0){
        $int_score = -1;
        $class = 'not_computable_was';
    }
    else {
        $int_score = sprintf("%.3f", $int_score/$base);
    }

    if ($int_score >= 0.469){
        $class = 'deleterious';
    }
    elsif ($int_score >= 0 && $int_score < 0.469) {
        $class = 'neutral';
    }

    # if the user wants an array, return the class and score, otherwise just return the class
        
    return ($class, $int_score);
}

1;