#!/usr/bin/perl -w
#
# Note: All decent programs, including Pymol or Jmol can color PDB files
#   according to B factor. It would be natural to generate a PDB file with
#   variability indicated in the B factor field. However, for some residues
#   the variability information may not be available. What color should be
#   used for such cases? The color of highly variable residues? We chose
#   a different color at the cost of more complicated code.
#
# The script will determine similarities or entropies on fly from a supplied
#   multiple sequence alignment. It will not use .gprf file, because it 
#   contains only ONE sequence and there is less flexibility when generating
#   alignment with the PDB sequence. With the supplied alignment the user has the
#   freedom to choose which of the sequences best resembles the PDB sequence.
#   (This problem was encountered with flavivirus E proteins.)
#

use FindBin;
use lib "$FindBin::RealBin";
use PDB;
use Alignments;
use color;

$params = parse_cmd_line_args();
$values = create_mapping($params);
$colors = determine_colors($params,$values);
output_jmol($params,$colors);
output_pymol($params,$colors);
output_molmol($params,$colors);
output_colorscale($params) unless exists($$params{'color_motifs'});

exit;


#-----------------------------------------

sub usage
{
    my (@msg) = @_;

    if (scalar @msg)
    {
        print STDERR join('',@msg);
        exit;
    }
    print STDERR 
        "About: The script colors a pdb file based on multiple sequence alignment.\n",
        "   The program will calculate the entropies or similarities\n",
        "   on fly and will need some of the parameters required also by motifmaker.\n",
        "Usage: variability-plot [OPTIONS]\n",
        "Options:\n",
        "   -a, --alignment <file.aln>                      Multiple sequence alignment to determine entropies/similarities from. With the -M option,\n",
        "                                                       the alignment should be the one used to generate the .prf file.\n",
        "   -c, --chain <letter>                            What chain of the PDB file should be aligned.\n",
        "   -C, --colors <#c8c8c8,#ff0000,..>               Use these colors to color the motifs: background,motif1,motif2,...\n",
        "       --color-scale <#ff0000,#ffffff,#0000ff>     Use the specified color scale for variability plots.\n",
        "       --color-scale-vals <min,avg,max>            All values smaller than min will colored red, etc.\n",
        "   -e, --entropies                                 Color according to the total entropies.\n",
        "   -f, --floating-colors                           Let minimum value correspond to red, maximum to blue and average to white.\n",
        "   -i, --id <string>                               What sequence from the alignment should be the PDB aligned with.\n",
        "   -l, --labels                                    Generate labels.\n",
        "   -L, --list <int,int,..>                         Color only selected motifs, indexed from 0. (See also the -M option.)\n",
        "   -m, --unscaled-max                              Use the unscaled-maximum approach.\n",
        "   -M, --motifs <file.prf>                         Color only residues which are part of a motif as listed in the .prf file.\n",
        "   -p, --prefix <string>                           Prefix of the output files.\n",
        "   -P, --pdb-file <file.pdb>                       The PDB file to be colored.\n",
        "   -s, --similarities                              Color according to similarities. (Default.)\n",
        "   -S, --scaled-sum                                Use the scaled-sum approach.\n",
        "       --substitute <list>                         If the alignment file contains letters which are not permitted, such as 'X',\n",
        "                                                       for example, replace the letter. The format is a comma separated list of\n",
        "                                                       expressions ORI:NEW. (For example 'X:-,.:-,\\ :-,*:-')\n",
        "   -u, --unscaled-sum                              Use the unscaled-sum approach. (This is the default.)\n",
        "   -V, --vector-lib <file.pm>                      Use the specified vectorlib.pm (see also create-vectorlib)\n",
        "   -v, --verbose                                   Print some usefull information.\n",
        "\n";
    exit;

}


sub parse_cmd_line_args
{
    my $cmd_args = join(' ',@ARGV);

    my $vectorlib = 'vectorlib-swissprot-56.0.pm';
    my $total_entropy_method = 'unscaled_sum';
    my $params;

    $$params{'default_color'} = 'C8C8C8';
    @default_colors = ('FAFF5C','5757FF','8CFF80','FF6E4A');
    my @motif_colors = @default_colors;

    $$params{'color_min'} = "#ff0000";
    $$params{'color_avg'} = "#ffffff";
    $$params{'color_max'} = "#0000ff";
    $$params{'color_min_value'} = 0.0;
    $$params{'color_avg_value'} = 0.5;
    $$params{'color_max_value'} = 1.0;

    $$params{'substitute'} = $Alignments::common_substitutions;
    $$params{'pdb_opts'}->{'uc'} = 1;
    $$params{'pdb_opts'}->{'chain'} = ' ';
    while (my $arg=shift(@ARGV))
    {
        if ( $arg eq '-a' || $arg eq '--alignment' ) { $$params{'aln_file'}=shift(@ARGV); next; }
        if ( $arg eq '-c' || $arg eq '--chain' ) { $$params{'pdb_opts'}->{'chain'}=shift(@ARGV); next; }
        if ( $arg eq '-C' || $arg eq '--colors' ) { $$params{'colors'}=shift(@ARGV); next; }
        if ( $arg eq '--color-scale-vals' ) 
        { 
            my $colors = shift(@ARGV);
            my ($min,$avg,$max) = split(/,/,$colors);
            $$params{'color_min_value'} = $min unless $min eq '';
            $$params{'color_avg_value'} = $avg unless $avg eq '';
            $$params{'color_max_value'} = $max unless $max eq '';
            next; 
        }
        if ( $arg eq '--color-scale' ) 
        { 
            my $colors = shift(@ARGV);
            my ($min,$avg,$max) = split(/,/,$colors);
            $$params{'color_min'} = $min unless $min eq '';
            $$params{'color_avg'} = $avg unless $avg eq '';
            $$params{'color_max'} = $max unless $max eq '';
            next; 
        }
        if ( $arg eq '-L' || $arg eq '--list' ) { $$params{'motif_list'}=shift(@ARGV); next; }
        if ( $arg eq '-e' || $arg eq '--entropies' ) { $$params{'plot_entropies'}=1; next; }
        if ( $arg eq '-f' || $arg eq '--floating-colors' ) { $$params{'floating_colors'}=1; next; }
        if ( $arg eq '-i' || $arg eq '--id' ) { $$params{'refseq_id'}=shift(@ARGV); next; }
        if ( $arg eq '-l' || $arg eq '--labels' ) { $$params{'labels'}=1; next; }
        if ( $arg eq '-m' || $arg eq '--unscaled-max' ) { $total_entropy_method = 'unscaled_max'; next; }
        if ( $arg eq '-M' || $arg eq '--motifs' ) 
        { 
            $$params{'prf_file'}     = shift(@ARGV); 
            $$params{'color_motifs'} = 1;
            next; 
        }
        if ( $arg eq '-p' || $arg eq '--prefix' ) { $$params{'prefix'}=shift(@ARGV); next; }
        if ( $arg eq '-P' || $arg eq '--pdb-file' ) { $$params{'pdb_file'}=shift(@ARGV); next; }
        if ( $arg eq '-s' || $arg eq '--similarities' ) { $$params{'plot_similarities'}=1; next; }
        if ( $arg eq '-S' || $arg eq '--scaled-sum' ) { $total_entropy_method = 'scaled_sum'; next; }
        if ( $arg eq '-u' || $arg eq '--unscaled-sum' ) { $total_entropy_method = 'unscaled_sum'; next; }
        if ( $arg eq '-V' || $arg eq '--vector-lib' ) { $vectorlib = shift(@ARGV); next; }
        if (                 $arg eq '--substitute' ) { $$params{'substitute'} = shift(@ARGV); next; }
        if ( $arg eq '-h' || $arg eq '--help' || $arg eq '-?' ) { usage(); }
        usage("Unknown argument: \"$arg\".\n");
    }

    if ( !exists($$params{'pdb_file'}) ) { usage("Missing the -P option.\n"); }
    if ( !exists($$params{'aln_file'}) ) { usage("Missing the -a option.\n"); }
    if ( !exists($$params{'refseq_id'}) ) 
    { 
        my $aln = Alignments::parse_alignment($$params{'aln_file'});
        my @seqs = ();
        for my $seq (@$aln) { push @seqs,$$seq{id}; }
        usage("No reference sequence was specified.\nThe alignment contains the following sequences:\n\t", join(",\n\t",@seqs),"\n");
    }
    if ( exists($$params{'colors'}) )
    {
        $$params{'colors'} =~ s/\#//g;

        @motif_colors = split(/,/,$$params{'colors'});
        $$params{'default_color'} = shift(@motif_colors);
    }
    if ( exists($$params{'motif_list'}) )
    {
        # Respect the order given by the user: the motifs do not have to
        #   be sorted, therefore rearrange colors accordingly.
        #
        my $colors = {};
        my @list = split(/,/,$$params{'motif_list'});
        for (my $i=0; $i<scalar @list; $i++)
        {
            my $motif = $list[$i];
            $$params{'motif_list'}{$motif} = 1;
            $$colors{$motif} = $motif_colors[$i];
        }
        $$params{'motif_colors'} = $colors;
    }

    require $vectorlib;
    require utils;
    if ( exists($$params{'aln_file'}) ) { Alignments::select_total_entropy_method($total_entropy_method); }

    $$params{'prefix'} = Utils::init_result_dir($$params{'prefix'});
    $$params{'pdb_outfile'} = $$params{'prefix'} . '-colored.pdb';
    $$params{'colorscale_file'} = $$params{'prefix'} . '-colorscale.png';

    open(LOG,">>$$params{'prefix'}.log") || usage("$$params{'prefix'}.log: $!");
    print LOG "variability-plot args: $cmd_args\n";
    close(LOG);

    return $params;
}



# The hash $values will be returned
#       $i                          .. The index of the residue in the PDB sequence (relative to chain)
#       $$values{$i}{'value'}        .. entropy or similarity value
#       $$values{$i}{'aa_pdb'}       .. amino acid one-letter code in the PDB sequence
#       $$values{$i}{'aa_ref'}       .. amino acid one-letter code in the reference sequence from the alignment
#       $$values{$i}{'refseq_pos'}   .. Position relative to reference sequence from the alignment
#       $$values{$i}{'absaln_pos'}   .. Absolute position of the residue in the alignment
#       $$values{$i}{'motif'}        .. Is this residue a part of one of the requested motifs?
#       $$values{$i}{'pdb_residue'}  .. PDB residue, as returned by PDB::read_pdb
#
sub create_mapping
{
    my ($params) = @_;

    my $pdb = PDB::read_pdb($$params{'pdb_file'});
    $pdb    = PDB::strip_pdb($pdb, $$params{'pdb_opts'});
    my $pdb_seq = PDB::get_sequence($pdb,$$params{'pdb_opts'});
    if ( !$pdb_seq || !length($pdb_seq) ) 
    { 
        usage("Could not determine the PDB sequence from the file $$params{'pdb_file'}, chain id '$$params{pdb_opts}->{chain}'.\n") 
    }

    # MolMol does not recognise the chain IDs listed in PDB, it works with
    #   sequential numbering instead.
    $$params{'pdb_opts'}{'chain_num'} = '';
    if ( $$params{'pdb_opts'}{'chain'} ne ' ' ) 
    {
        my $chain_numbering = PDB::get_chain_numbering($pdb);
        my $chain = $$params{'pdb_opts'}{'chain'};
        if ( !exists($$chain_numbering{$chain}) )
        {
            usage("The PDB structure does not contain the chain '$chain'.");
        }
        $$params{'pdb_opts'}{'chain_num'} = $$chain_numbering{$chain} + 1;
    }

    # Calculate the entropies or similarities on fly.
    my $aln = Alignments::parse_alignment($$params{'aln_file'},{'substitute'=>$$params{'substitute'}});
    my ($consensus,$entropies,$similarities) = Alignments::create_consensus_sequence($aln);
    $consensus = ''; # To keep compiler happy, we will not need the consensus.

    # Align the PDB sequence with the reference sequence from the alignment. To do this, get the reference sequence.
    my $ref_seq = Alignments::get_sequence($aln,$$params{'refseq_id'});
    if ( !$ref_seq || !length($ref_seq) ) 
    { 
        my @seqs = ();
        for my $seq (@$aln) { push @seqs,$$seq{id}; }
        usage(
                "No such sequence '$$params{refseq_id}' in the alignment $$params{'aln_file'}.\n",
                "The alignment contains the following sequences:\n\t", join(",\n\t",@seqs), "\n"
             ) ;
    }

    # Mapping of the indexes from relative to absolute position in the alignment from the user. (Will be needed later.)
    #   This is required when gaps are present.
    my $refseq_rel_to_abs = Alignments::relative_to_absolute($ref_seq);

    # Now run the alignment with the PDB sequence.
    my @seqs = ();
    push @seqs, { id=>'pdb',    seq=>$pdb_seq };
    push @seqs, { id=>'refseq', seq=>$ref_seq };
    my ($refseq_aln) = Alignments::run_clustal(\@seqs,{delete=>0,prefix=>$$params{'prefix'}.'-pdb2ref'});

    # Parse the alignment and determine mapping of indexes between the two sequences to treat gaps properly.
    # !!! NOTE:
    #   If making any changes here, make sure that you test the new code THOROUGHLY. These relative-to-absolute
    #   and relative-to-relative mappings are likely to give you a real headache.
    #
    my $pdb_aligned_seq = Alignments::get_sequence($refseq_aln,'pdb');
    my $ref_aligned_seq = Alignments::get_sequence($refseq_aln,'refseq');
    if ( !$pdb_aligned_seq || !$ref_aligned_seq ) { usage("Could not determine results of the PDB to reference sequence alignment.\n"); }
    my $pdb_to_ref_rel_rel = Alignments::relative_to_relative($pdb_aligned_seq,$ref_aligned_seq);

    my $motif_residues = {};
    if ( $$params{'color_motifs'} && exists($$params{'prf_file'}) )
    {
        my $prf = Utils::read_prf_file_motifs($$params{'prf_file'});
        for my $motif (keys %$prf)
        {
            # Omit motifs which were not requested. Without the -L option, all motifs
            #   will be printed.
            if ( exists($$params{'motif_list'}) && !exists($$params{'motif_list'}{$motif}) ) { next; }
            
            my @residues = @{$$prf{$motif}{'residues'}};
            my $from = $residues[0]->{'aln_column'} - 1;
            my $to   = $residues[-1]->{'aln_column'} - 1;
            for ($idx=$from; $idx<=$to; $idx++)
            {
                $$motif_residues{$idx} = $motif;
            }
        }
    }

    my $values = {};

    # Finally, generate a hash of values.
    # Go through all indexes of the rel-to-rel mapping.
    for (my $iref_pdb=0; $iref_pdb<scalar @$pdb_to_ref_rel_rel; $iref_pdb++)
    {
        my $iref = $$pdb_to_ref_rel_rel[$iref_pdb];     # iref is an index in the reference sequence with no gaps present
        my $ialn = $$refseq_rel_to_abs[$iref];          # ialn is an absolute index in the *-pdb2ref alignment

        if ( $iref==-1 || $ialn==-1 ) { next }          # Negative indices indicate a gap in the alignment
    
        my $pdb_letter = substr($pdb_seq,$iref_pdb,1);

        # A sanity check: The $ref_letter taken from the user-supplied alignment should match
        #   the $pdb_letter:
        #
        #   my $ref_letter = substr($ref_seq,$ialn,1);      
        #   print STDERR "$pdb_letter:$iref_pdb .. $ref_letter:$ialn .. $$entropies[$ialn], $$similarities[$ialn]\n";
    
        $$values{$iref_pdb}{'value'}  = exists($$params{'plot_entropies'}) ? $$entropies[$ialn] : $$similarities[$ialn];
        $$values{$iref_pdb}{'motif'}  = exists($$motif_residues{$ialn}) ? $$motif_residues{$ialn} : -1;
        $$values{$iref_pdb}{'aa_pdb'} = $pdb_letter;
        $$values{$iref_pdb}{'aa_ref'} = substr($ref_seq,$ialn,1);
        $$values{$iref_pdb}{'refseq_pos'}  = $iref;
        $$values{$iref_pdb}{'absaln_pos'}  = $ialn;
        $$values{$iref_pdb}{'pdb_residue'} = $$pdb[$iref_pdb];
    }

    return $values;
}




sub determine_colors
{
    my ($params,$values) = @_;

    my @colors = ();

    if ( exists($$params{'color_motifs'}) )
    {
        my $ncolors = scalar @default_colors;
        
        for my $idx (sort {$a<=>$b} keys %$values)
        {
            if ( $$values{$idx}{'motif'} == -1 ) { next; }

            my $motif = $$values{$idx}{'motif'};
            my $color = exists($$params{'motif_colors'}{$motif}) ? $$params{'motif_colors'}{$motif} : $default_colors[ $motif % $ncolors ];
            $$values{$idx}{'color'} = $color;
            push @colors, $$values{$idx};
        }
        return \@colors;
    }

    my $min = $$params{'color_min_value'};
    my $avg = $$params{'color_avg_value'};
    my $max = $$params{'color_max_value'};

    if ( exists($$params{'floating_colors'}) && $$params{'floating_colors'}==1 )
    {
        my $n = 0;
        for my $idx (keys %$values)
        {
            my $value = $$values{$idx}{'value'};
            if ( !$n ) 
            { 
                $min = $value;
                $max = $value;
                $avg = 0;
            }
            if ( $min>$value ) { $min=$value }
            if ( $max<$value ) { $max=$value }
            $avg += $value;
            $n++;
        }
        $avg /= $n;
    }

    $$params{'min'} = $min;
    $$params{'max'} = $max;
    $$params{'avg'} = $avg;

    #   if ( $verbose )
    #   {
    #       printf STDERR "min,avg,max = %.2f,%.2f,%.2f\n", $min,$avg,$max;
    #   }

    for my $idx (sort {$a<=>$b} keys %$values)
    {
        my $color = Color::get_color_string($$params{'color_min'},$$params{'color_avg'},$$params{'color_max'}, $$values{$idx}{'value'}, $min,$avg,$max);
        $color =~ s/^\#//;
        $$values{$idx}{'color'} = $color;
        push @colors, $$values{$idx};
    }
    
    return \@colors;
}



sub output_colorscale
{
    my ($params) = @_;

    require Image::Magick;

    my $img_width  = 500;
    my $img_height = 40; 
    my $font_size  = 12;
    
    my $min_color = $$params{'color_min'};
    my $avg_color = $$params{'color_avg'};
    my $max_color = $$params{'color_max'};
    my $min_value = $$params{'color_min_value'};
    my $avg_value = $$params{'color_avg_value'};
    my $max_value = $$params{'color_max_value'};

    my $img=Image::Magick->new(size=>"${img_width}x$img_height");
    $img->ReadImage('xc:white');
    
    my @text_params = $img->QueryFontMetrics(text=>'1.0', pointsize=>$font_size);
    my $text_width  = $text_params[4];  # tex_ppem, y_ppem, ascender, descender, width, height, max_advance
    my $text_height = $text_params[5];

    my $margins = $text_width*0.5;
    my $padding = $text_height*0.3;
    my $colorbar_width  = $img_width - 2*$margins;
    my $colorbar_height = $img_height - $padding - $text_height;

    my $nsteps = 100;
    my $step   = $colorbar_width*1.0 / $nsteps;

    my $woffset = $margins;
    for (my $i=0; $i<$nsteps; $i++)
    {
        my $value  = $i /$nsteps;
        my $color  = Color::get_color_string($min_color,$avg_color,$max_color,$value,$min_value,$avg_value,$max_value);
        
        my $w_from = sprintf "%d", $woffset + $i*$step;
        my $w_to   = sprintf "%d", $w_from + $step;
        my $h_from = sprintf "%d", 0;
        my $h_to   = sprintf "%d", $h_from + $colorbar_height;
        
        $img->Draw(primitive=>'Rectangle',points=>"$w_from,$h_from $w_to,$h_to",fill=>$color,bordercolor=>$color,stroke=>$color, antialias=>'false');
    }
    for (my $i=0; $i<=1.0; $i+=0.1)
    {
        my $text = sprintf "%.1f", $i;
        
        my $w = sprintf "%d", $woffset + $i*$colorbar_width;
        my $w_text = sprintf "%d", $w - $text_width*0.5;
        my $h_from = sprintf "%d", $colorbar_height;
        my $h_to   = sprintf "%d", $h_from + $padding;
        my $h_text = sprintf "%d", $h_from + $text_height;

        if ( $w<$woffset ) { $w = sprintf "%d", $woffset }
        if ( $w>=$woffset+$colorbar_width-1 ) { $w=sprintf "%d", $woffset+$colorbar_width-1 }
        
        $img->Draw(primitive=>'Line',points=>"$w,$h_from $w,$h_to", antialias=>'false');
        $img->Annotate(text=>$text, x=>$w_text, y=>$h_text, antialias=>'true', pointsize=>$font_size);
    }
    
    $img->Write($$params{'colorscale_file'});
}


sub output_jmol
{
    my ($params,$colors) = @_;
    open(JMOL,">$$params{'prefix'}.jmol") || usage("$$params{'prefix'}.jmol: $!");

    print JMOL << "EOT";
var defaultConf = "language=\\\"en_US\\\"; background [230,230,230];"
+ "cpk on; cartoon off; wireframe off; label off; hide hetero;"
+ "select all; color [x$$params{'default_color'}];"
EOT

    my $chain = $$params{'pdb_opts'}->{'chain'} ne ' ' ? ":$$params{'pdb_opts'}->{'chain'}" : '';
    for my $col (@$colors)
    {
        my $resi  = $$col{'pdb_residue'}{'res_number'};
        my $color = $$col{'color'};
        print JMOL "+ \"select $resi$chain;color [x$color]; label off;\"\n";
    }
    close(JMOL);
}



sub output_molmol
{
    my ($params,$colors) = @_;
    open(MOLMOL,">$$params{'prefix'}.mol") || usage("$$params{'prefix'}.mol: $!");

    my (@color) = Color::parse_color($$params{'default_color'});
    my $default_color = sprintf "%.2f %.2f %.2f", @color;

    print MOLMOL << "EOT";
StyleAtom invisible
StyleBond invisible
AutoScale
SelectAtom ''
ColorAtom $default_color
SelectAtom ''
XMacStand ribbon.mac
PaintRibbon atom
EOT

    my $chain = $$params{'pdb_opts'}->{'chain_num'};
    for my $col (@$colors)
    {
        my $resi = $$col{'pdb_residue'}{'res_number'};
        (@color) = Color::parse_color($$col{'color'});
        print MOLMOL  "SelectAtom '#$chain:$resi'\n";
        printf MOLMOL "ColorAtom %.2f %.2f %.2f\n", @color;
    }
    close(MOLMOL);
}





sub output_pymol
{
    my ($params,$colors) = @_;
    open(PYMOL,">$$params{'prefix'}.pymol") || usage("$$params{'prefix'}.pymol: $!");

    my $chain = $$params{'pdb_opts'}->{'chain'} ne ' ' ? "& c. $$params{'pdb_opts'}->{'chain'}" : '';
    my @items = ();
    my $pymol_colors = {};  # Used only when coloring motifs. Otherwise, gradient is generated by pymol on fly.
    for my $col (@$colors)
    {
        my $resi  = $$col{'pdb_residue'}{'res_number'};
        my $value = $$col{'value'};
        if ( $$params{'color_motifs'} )
        {
            my $ncolors = scalar keys %$pymol_colors;
            if ( !exists($$pymol_colors{$$col{'color'}}) )
            {
                my ($r,$g,$b) = Color::parse_color($$col{'color'});

                $$pymol_colors{$$col{'color'}}->{'name'}  = "col_$ncolors";
                $$pymol_colors{$$col{'color'}}->{'value'} = "($r,$g,$b)";
            }
            
            $value = "'$$pymol_colors{$$col{'color'}}->{'name'}'";
        }
        push @items, "[selection + \" i. $resi $chain\", $value]";
    }
    my $to_be_colored = join(",\n\t", @items);
    my $labels = exists($$params{'labels'}) && $$params{'labels'} ? "cmd.label(sele+\" & n. ca\",'\"%s%s\" % (resn,resi)')" : '';

    my ($r,$g,$b) = Color::parse_color($$params{'default_color'});
    my $col_unknown = "$r,$b,$b";

    print PYMOL '# Modify colors and selection and then run this script from'                   . "\n";
    print PYMOL '# the pymol console as'                                                        . "\n";
    print PYMOL '#   run /path/to/this-script.pymol'                                            . "\n";
    print PYMOL '#   '                                                                          . "\n";
    print PYMOL '    '                                                                          . "\n";
    print PYMOL 'from pymol import cmd    '                                                     . "\n";
    print PYMOL                                                                                   "\n";
    print PYMOL 'selection   = ""             # When working with multiple files at once,   '   . "\n";
    print PYMOL '                             #   set this e.g. to "MyPDB" to prevent     '     . "\n";
    print PYMOL '                             #   coloring of other PDBs.    '                  . "\n";
    print PYMOL                                                                                   "\n";
    print PYMOL "col_unknown = ($col_unknown) # Light yellow                                       \n";

    if ( $$params{'color_motifs'} )
    {
        print PYMOL 'cmd.set_color("col_unknown",col_unknown)      ' . "\n";
        for my $py_col (keys %$pymol_colors)
        {
            print PYMOL "cmd.set_color(\"$$pymol_colors{$py_col}{'name'}\", $$pymol_colors{$py_col}{'value'})\n";
        }

        print PYMOL "\n";
        print PYMOL "def color_selection(items):\n";
        print PYMOL "    for i in range(len(items)):\n";
        print PYMOL "        sele  = items[i][0]\n";
        print PYMOL "        color = items[i][1]\n";
        print PYMOL "        $labels\n";
        print PYMOL "        cmd.color(color,sele)\n";
        print PYMOL "\n";
        print PYMOL 'cmd.color("col_unknown", selection)           ' . "\n";
        print PYMOL                                                    "\n";
        print PYMOL 'if len(selection)!=0: selection += " & "  '     . "\n";
        print PYMOL                                                    "\n";
        print PYMOL 'color_selection([         '                     . "\n";
        print PYMOL "    $to_be_colored                                 \n";
        print PYMOL '    ])            '                             . "\n";
    }
    else
    {
        print PYMOL << "EOT";
col_min     = (1.0,0.1,0.1)  # Reddish
col_avg     = (0.9,0.9,0.9)  # Very light gray
col_max     = (0.1,0.1,1.0)  # Blueish
nbins       = 15             # Number of color shades
dvalue      = 1.0/nbins      # The length of color intervals


# ------------------------------------------------------------
# No changes should be required after this point
# 

def value_to_rgb(value,min=$$params{'min'},avg=$$params{'avg'},max=$$params{'max'}):
    if value<0: value=0
    if value<avg:
        r = (value/avg)*(col_avg[0]-col_min[0]) + col_min[0]
        g = (value/avg)*(col_avg[1]-col_min[1]) + col_min[1]
        b = (value/avg)*(col_avg[2]-col_min[2]) + col_min[2]
        return (r,g,b)

    if value>1: value=1
    r = ((value-avg)/(max-avg))*(col_max[0]-col_avg[0]) + col_avg[0]
    g = ((value-avg)/(max-avg))*(col_max[1]-col_avg[1]) + col_avg[1]
    b = ((value-avg)/(max-avg))*(col_max[2]-col_avg[2]) + col_avg[2]
    return (r,g,b)


def init_colors(min=$$params{'min'},avg=$$params{'avg'},max=$$params{'max'}):
    cmd.set_color("col_unknown",col_unknown)
    colors = []
    for ibin in range(nbins+1):
        value = dvalue*ibin
        color = "col_" + str(ibin)
        colors.append(color)
        cmd.set_color(color,value_to_rgb(value,min,avg,max))
    return colors


def color_selection(items):
    for i in range(len(items)):
        sele  = items[i][0]
        simil = items[i][1]

        $labels

        if simil==-1:
            cmd.color("col_unknown",sele)
        else:
            idx = int(simil/dvalue)
            if idx<0: idx = 0
            if idx>nbins: idx = nbins
            cmd.color(colors[idx],sele)



colors = init_colors()
cmd.color("col_unknown", selection)

if len(selection)!=0: selection += " & "

color_selection([
    $to_be_colored
    ])
EOT
    }

    close(PYMOL);
}


