websoft/bin/concise2html

   1 #!/usr/bin/perl
   2
   3 use strict;
   4 use warnings;
   5 use Jpred;
   6 use Getopt::Long;
   7
   8 my ($in, $out);
   9 my $format = 'seq';
  10 GetOptions(
  11    "in=s"    => \$in,
  12    "out=s"   => \$out,
  13    "format=s" => \$format
  14 );
  15
  16 if (!$in) { $in = "-"; }
  17 if (!$out) { $out = "-"; }
  18
  19 # Read in the concise CSV file. Data is stored in hash's with the key
  20 # being the name of the type of data and the value being the sequence
  21 # data, the order of the alignments is also stored for display in the
  22 # order with which they are aligned.
  23
  24 my (@align, @seq, %predict, %prob);
  25
  26 my @predictions = qw(Lupas_14 Lupas_21 Lupas_28 jnetpred JNETHMM JNETALIGN JNETPSSM JNETJURY JNETSOL25 JNETSOL5 JNETSOL0 JNETCONF);
  27
  28 my @prob = qw(JNETPROPH JNETPROPB JNETPROPC);
  29
  30 # Arrays that holds the order in which things are printed
  31 my @print_pred = qw(jnetpred JNETALIGN JNETHMM JNETPSSM);
  32 my @coils = qw(Lupas_14 Lupas_21 Lupas_28);
  33 my @solvent = qw(JNETSOL25 JNETSOL5 JNETSOL0);
  34
  35 # Hash that relates data names to printed names
  36 my %convert = (
  37         JNETALIGN => "jalign",
  38         JNETHMM => "jhmm",
  39         jnetpred => "Jnet",
  40         JNETPSSM => "jpssm",
  41
  42         Lupas_14 => "Lupas 14",
  43         Lupas_21 => "Lupas 21",
  44         Lupas_28 => "Lupas 28",
  45
  46         JNETSOL25 => "Jnet_25",
  47         JNETSOL5 => "Jnet_5",
  48         JNETSOL0 => "Jnet_0",
  49         JNETCONF => "Jnet Rel"
  50 );
  51
  52 # Read in the data
  53 open(IN, "<$in") or die($!);
  54 while (<IN>) {
  55         if (/:/) {
  56                 my ($id, $seq) = split(":", $_);
  57
  58                 if ($id =~ /;/) {       # Find sequence alignments
  59                         my ($align, $seqid) = split(";", $id);
  60                         push @align, $seqid;
  61                         $seq =~ s/,//g;
  62                         chomp($seq);
  63                         push @seq, $seq;
  64                         }
  65
  66                 foreach (@predictions) {        # Find the predictions
  67                         if ($id eq $_) {
  68                                 if ($predict{$id}) { die("Same ID ($id) found twice in concise file!"); }
  69                                 $seq =~ s/,//g;
  70                                 chomp($seq);
  71                                 $predict{$id} = $seq;
  72                                 }
  73                         }
  74
  75                 foreach (@prob) {
  76                         if ($id eq $_) {
  77                                 if ($prob{$id}) { die("Same ID ($id) found twice in concise file!"); }
  78                                 chomp($seq);
  79                                 @_ = split(",", $seq);
  80                                 foreach (@_) {
  81                                         if ($_ >= 0.85) { $_ = 9; }
  82                                         else { $_ = int(($_ * 10) + 0.5); }
  83                                         }
  84                                 my $temp = join("", @_);
  85                                 $prob{$id} = $temp;
  86                                 }
  87                         }
  88                 }
  89         }
  90 close(IN);
  91
  92 ############################################
  93 #
  94 # Sets up the references for prediction
  95 # methods and deffn methods
  96 #
  97
  98 my @predname = qw(jnetpred JNETALIGN JNETPSSM JNETHMM);
  99 foreach (@predname) {
 100         if (!$predict{$_}) { warn("\"$_\" data not present\n"); next; }
 101         $predict{$_} =~ s/[TCYWXZ]/-/g;
 102         }
 103
 104 ######################################################################
 105 # Print out the sequences to HTML
 106 ######################################################################
 107
 108 open(HTML, ">$out") or die($!);
 109 print HTML "<!DOCTYPE html PUBLIC \"-//W3C//DTD HTML 3.2//EN\"><html><head><meta http-equiv=\"Content-Type\" content=\"text/html; charset=iso-8859-1\"><title>Jpred secondary structure prediction results</title></head><body bgcolor=\"#ffffff\"><pre><code>";
 110
 111 #
 112 # Prints the sequence passed to it as a reference to an array, second
 113 # array contains the "conservation", which varies between 0 and 1 for
 114 # each position.
 115 #
 116
 117 sub printseq {
 118         my $seqref = shift @_;
 119         my $scoreref = shift @_;
 120
 121         unless ($scoreref) {
 122                 foreach (0..$#{$seqref}) {
 123                         print HTML ${$seqref}[$_];
 124                 }
 125                 return;
 126         }
 127
 128 #       foreach (0..$#{$seqref} - 1) {
 129         foreach (0..$#{$seqref}) {
 130                 my $colour = "#";       # Work out the colour
 131                 if (${$scoreref}[$_] == 1) {
 132                         $colour .= "0000ff";
 133                 }
 134                 else {
 135                         $colour .= (sprintf "%-2.2lx", ${$scoreref}[$_] * 255).(sprintf "%-2.2lx", 0).(sprintf "%-2.2lx", 0);
 136                 }
 137
 138                 print HTML "<font color=\"$colour\">${$seqref}[$_]</font>";
 139         }
 140 }
 141
 142 # Prints the gap at the begining of the line between the sequence ID
 143 # and the ":". First argument is the length of the ID and the second is
 144 # the maximimum length.
 145
 146 sub printgap {
 147         my ($strlen, $max_len) = @_;
 148         foreach ($strlen .. $max_len) {
 149                 print HTML " ";
 150                 }
 151         print HTML ": ";
 152         }
 153
 154 # Works out a consensus score according to Zvelebil et al. (1987)
 155 # J. Mol. Biol. 195, pp957, and does the smoothing
 156
 157 sub cons {
 158         my %conservation = (
 159                 I => [1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
 160                 L => [1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
 161                 V => [1, 0, 0, 0, 0, 1, 0, 1, 0, 0],
 162                 C => [1, 0, 0, 0, 0, 1, 0, 0, 0, 0],
 163                 A => [1, 0, 0, 0, 0, 1, 1, 0, 0, 0],
 164                 G => [1, 0, 0, 0, 0, 1, 1, 0, 0, 0],
 165                 M => [1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
 166                 F => [1, 0, 0, 0, 0, 0, 0, 0, 1, 0],
 167                 Y => [1, 0, 0, 1, 0, 0, 0, 0, 1, 0],
 168                 W => [1, 0, 0, 1, 0, 0, 0, 0, 1, 0],
 169                 H => [1, 1, 0, 1, 1, 0, 0, 0, 1, 0],
 170                 K => [1, 1, 0, 1, 1, 0, 0, 0, 0, 0],
 171                 R => [0, 1, 0, 1, 1, 0, 0, 0, 0, 0],
 172                 E => [0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
 173                 Q => [0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
 174                 D => [0, 0, 1, 1, 0, 1, 0, 0, 0, 0],
 175                 N => [0, 0, 0, 1, 0, 1, 0, 0, 0, 0],
 176                 S => [0, 0, 0, 1, 0, 1, 1, 0, 0, 0],
 177                 T => [1, 0, 0, 1, 0, 1, 0, 0, 0, 0],
 178                 P => [0, 0, 0, 0, 0, 1, 0, 0, 0, 1],
 179                 Z => [0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
 180                 '-' => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 181                 '_' => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 182                 '.' => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 183                 X => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
 184                 B => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], # Not mentioned, so use X
 185                 U => [1, 0, 0, 0, 0, 1, 0, 0, 0, 0], # Not mentioned, so use C
 186                 );
 187
 188         my $aref = shift @_;
 189         my (@cons, @align);
 190
 191         # Get each sequence
 192         foreach (@{$aref}) {
 193                 push @align, $_;
 194         }
 195
 196         my $seq_len = length($align[0]);
 197
 198         # Foreach position
 199         foreach my $pos (0..$seq_len - 1) {
 200                 my @temp;
 201
 202                 # Foreach sequence, get the residue for this position
 203                 foreach (@align) { push @temp, substr($_, $pos, 1); }
 204
 205                 my ($p, $i, $j);
 206                 for ($i = 0; $i < $#temp; $i++) {
 207                         my $one = uc($temp[$i]);
 208                         for ($j = $i + 1; $j < $#temp; $j++) {
 209                                 my $two = uc($temp[$j]);
 210                                 foreach (0..9) {
 211                                         # This debugging is in here
 212                                         # incase a residue occurs that we
 213                                         # haven't got in the hash above
 214                                         unless (defined $conservation{$one}[$_]) {
 215 #                                               warn "first $one, $_\n";
 216 #                                               die "'$conservation{$one}[$_]'\n";
 217                                                 next;
 218                                         }
 219                                         unless (defined $conservation{$two}[$_]) {
 220 #                                               warn "second $two, $_\n";
 221 #                                               die "'$conservation{$two}[$_]'\n";
 222                                                 next;
 223                                         }
 224                                         if ($conservation{$one}[$_] != $conservation{$two}[$_]) { $p++; }
 225                                 }
 226                         }
 227                 }
 228
 229                 # This loop looks to be v. slow
 230                 if ($p) {
 231                         if (($p = 0.9 - 0.1 * $p) < 0) { $p = 0; }
 232                 }
 233                 else {
 234                         for ($i = 1; $i < $#temp; $i++) {
 235                                 if ($temp[$i] ne $temp[0]) {
 236                                         $p++;
 237                                         last;
 238                                 }
 239                         }
 240                         if ($p) { $p = 0.9; }
 241                         else { $p = 1; }
 242                 }
 243
 244                 push @cons, $p;
 245         }
 246
 247         # Average the result over 3 residues
 248 #       my (@final, $i);
 249 #       push @final, $cons[0];
 250 #       for ($i = 1; $i < ($#cons); $i++) {
 251 #               push @final, ($cons[$i - 1] + $cons[$i] + $cons[$i + 1]) / 3;
 252 #               }
 253 #       push @final, $cons[$#cons + 1];
 254
 255         return @cons;
 256 }
 257
 258 ############################################################
 259 # Uncomment this to get the conservation scoring worked out
 260 ############################################################
 261 #my @score = &cons(\@seq);
 262
 263 #
 264 # Print the sequences
 265 #
 266
 267 # Find the longest identifier and set the space at the start of the line
 268
 269 my $max_len = 0;
 270 foreach (@align) {
 271         if (length($_) > $max_len) { $max_len = length($_); }
 272         }
 273 if ($max_len <= 6) { $max_len = 11; }
 274
 275 foreach (0..$#align) {
 276         # Begining of line
 277         if ($_ == 0) {
 278                 print HTML "$align[$_]";
 279                 }
 280         else {
 281            if ($format eq 'seq') {
 282               ## differentiate between Uniprot and Uniparc (UniRef90_UNI\w+) IDs as they require
 283               ## different SRS query strings
 284               if ($align[$_] =~ /UniRef90_(UPI\w+)/) {
 285                  print HTML "<a href=$SRSSERVER?db=allebi&query=$1>$align[$_]</a>";
 286               } else {
 287                  print HTML "<a href=$SRSSERVER?db=allebi&query=$align[$_]>$align[$_]</a>";
 288               }
 289            } else {
 290               print HTML "$align[$_]";
 291            }
 292                 }
 293         &printgap(length($align[$_]) ,$max_len);
 294
 295         # Print the sequence
 296         @_ = split("", $seq[$_] );
 297 #       &printseq([@_], [@score]);
 298         &printseq([@_]);
 299
 300         # End of the line
 301         if ($_ == 0) {
 302                 print HTML " : $align[$_] \n";
 303                 }
 304         else {
 305            if ($format eq 'seq') {
 306               ## differentiate between Uniprot and Uniparc IDs as they require
 307               ## different SRS query strings
 308               if ($align[$_] =~ /UniRef90_(UPI\w+)/) {
 309                  print HTML " : <a href=$SRSSERVER?db=allebi&query=$1>$align[$_]</a>\n";
 310               } else {
 311                  print HTML " : <a href=$SRSSERVER?db=allebi&query=$align[$_]>$align[$_]</a>\n";
 312               }
 313            } else {
 314               print HTML " : $align[$_]\n";
 315            }
 316                 }
 317         }
 318 print HTML "\n";
 319
 320 ######################################################################
 321
 322 #
 323 # Create a guide sequence of the original query
 324 # (of the form 1----------11--- etc.)
 325 #
 326
 327 my ($guide, $i);
 328 my $seq_len = length($seq[0]);
 329
 330 for ($i = 1; $i < $seq_len; $i += 10) {
 331         my $dash;
 332         my $no_dash = 9 - length($i);
 333         foreach (0..$no_dash) { $dash .= "-"; }
 334         $guide .= $i.$dash;
 335         }
 336
 337 while (length($guide) > $seq_len) {     # Trim spare guide away
 338         chop($guide);
 339         }
 340
 341 printgap(0, $max_len);
 342 print HTML $guide." :\n";
 343
 344 print HTML "OrigSeq";
 345 printgap(7, $max_len);
 346 print HTML $seq[0];
 347 print HTML " : OrigSeq\n\n";
 348
 349 # Print out secondary structure predictions, goes through in order of
 350 # @print_pred so comes out right
 351
 352 foreach (@print_pred) {
 353         if (!$predict{$_}) { next; }
 354         if ($_ eq "JNETALIGN") {        #  Add a line before JNETALIGN which will put jnetpred on its own
 355                 print HTML "\n";
 356                 }
 357
 358         print HTML $convert{$_};
 359         printgap(length($convert{$_}), $max_len);
 360         @_ = split("", $predict{$_});
 361         foreach (@_) {
 362                 if ($_ eq "E") {
 363                         print HTML "<font color=ffa800>$_</font>";
 364                         }
 365                 elsif ($_ eq "H") {
 366                         print HTML "<font color=e90055>$_</font>";
 367                         }
 368                 elsif ($_ eq "-") {
 369                         print HTML $_;
 370                         }
 371                 }
 372         print HTML " : $convert{$_}\n";
 373         }
 374 print HTML "\n";
 375
 376 foreach (@coils) {
 377         if (!$predict{$_}) { next; }
 378         print HTML $convert{$_};
 379         printgap(length($convert{$_}), $max_len);
 380         @_ = split("", $predict{$_});
 381         foreach (@_) {
 382                 if ($_ eq "C") { print HTML "<font color=00aa00>$_</font>"; }
 383                 else { print HTML $_; }
 384                 }
 385         print HTML " : $convert{$_}\n";
 386         }
 387 print HTML "\n";
 388
 389 foreach (@solvent) {
 390         if (!$predict{$_}) { next; }
 391         print HTML $convert{$_};
 392         printgap(length($convert{$_}), $max_len);
 393         @_ = split("", $predict{$_});
 394         foreach (@_) {
 395                 if ($_ eq "B") {
 396                         print HTML "<font color=aa0000>$_</font>";
 397                         }
 398                 else { print HTML $_; }
 399                 }
 400         print HTML " : $convert{$_}\n";
 401         }
 402
 403 printf HTML $convert{"JNETCONF"};
 404 printgap(8, $max_len);
 405
 406 @_ = split("", $predict{"JNETCONF"});
 407 foreach (@_) {
 408         if ($_ >= 7) { print HTML "<font color=00aa00>$_</font>"; }
 409         else { print HTML "$_"; }
 410         }
 411 print HTML " : ".$convert{"JNETCONF"}."\n";
 412
 413 print HTML "</code></pre>";
 414 print HTML "<hr>";
 415
 416 #
 417 # Key to data
 418 #
 419
 420 print HTML "<pre><code>";
 421 print HTML "<b>\nNotes
 422 Key:</b>
 423 Colour code for alignment:
 424 Blue            - Complete identity at a position
 425 Shades of red   - The more red a position is, the higher the level of
 426                   conservation of chemical properties of the amino acids
 427 Jnet            - Final secondary structure prediction for query
 428 jalign          - Jnet alignment prediction
 429 jhmm            - Jnet hmm profile prediction
 430 jpssm           - Jnet PSIBLAST pssm profile prediction
 431
 432 Lupas           - Lupas Coil prediction (window size of 14, 21 and 28)
 433
 434 Note on coiled coil predictions  - = less than 50% probability
 435                                  c = between 50% and 90% probability
 436                                  C = greater than 90% probability
 437
 438 Jnet_25         - Jnet prediction of burial, less than 25% solvent accesibility
 439 Jnet_5          - Jnet prediction of burial, less than 5% exposure
 440 Jnet_0          - Jnet prediction of burial, 0% exposure
 441 Jnet Rel        - Jnet reliability of prediction accuracy, ranges from 0 to 9, bigger is better.
 442 ";
 443
 444 print HTML "</code></pre>$JPREDFOOT</body></html>";
 445 close(HTML);