websoft/bin/prettify.pl

   1 #!/usr/bin/perl
   2
   3 =head1 NAME
   4
   5 prettify.pl
   6
   7 =head1 SYNOPSYS
   8
   9 prettify.pl [--blc I<path>] [--alscript I<path>] --postscript I<path>
  10
  11 =head1 DESCRIPTION
  12
  13 This program uses a BLC file to create an Alscript command file.
  14
  15 =head1 TODO
  16
  17 Add more options, make prettier Alscript command files.
  18
  19 =head1 BUGS
  20
  21 Alscript can't count properly and so we have to add one to some of the indicies to make it work.
  22
  23 =head1 AUTHORS
  24
  25 Started by Patrick Audley <paudley@compbio.dundee.ac.uk>
  26 Extended by Jonathan Barber <jon@compbio.dundee.ac.uk>
  27
  28 =cut
  29
  30 use strict;
  31 use warnings;
  32 use Getopt::Long;
  33 use BLC;
  34
  35 my $prefix = "seq.protein";
  36 my $bin = "/software/submit_all";
  37
  38 my $blc = '-';
  39 my $alscript = '-';
  40 my $psout;
  41
  42 GetOptions(
  43         "blc:s" => \$blc,
  44         "alscript:s" => \$alscript,
  45         "postscript=s" => \$psout,
  46 ) or die "$!\n";
  47
  48 die "--postscript option required\n" unless $psout;
  49
  50 my $data = new BLC;
  51 $data->read_file($blc);
  52
  53 # Find which entries for secondary structure are in the sequence data
  54 my @jpred = map { /^(?:jnet|jhmm|jalign|jpssm|lupas_14|lupas_21|lupas_28)$/ or 0 } $data->get_seq_ids;
  55
  56 # Replace space characters with '-' for the FSM
  57 my @sequences = map { my $c = $_; $c =~  s/ /-/g; $c } $data->get_sequences;
  58
  59 # Find the positions of the secondary structure in the BLC file
  60 my ($top, $bottom, $alignEnd);
  61 foreach (0 .. $#jpred) { $top = $_ and last if $jpred[$_] }
  62 foreach (reverse(0 .. $#jpred)) { $bottom = $_ and last if $jpred[$_] }
  63 my $end = length $sequences[0];
  64
  65 # ... and then add 1 to them as Alscript can't count properly
  66 $alignEnd = $top;
  67 $top++;
  68 $bottom++;
  69
  70 # alscript bits
  71 open ALSOUT, ">$alscript" or die "$alscript: Can't open the Alscript command file\n";
  72
  73 # Commands for Alscript
  74 # Add one as Alscript can't count, and make sure it's over a minimum otherwise
  75 # it cores...
  76 my $foo = @sequences > 50 ? @sequences + 1 : 50;
  77 my $MAX_NSEQ = "MAX_NSEQ ".$foo;
  78 # Buffer size for reading the BLC file needs to be at least the size of alignment +2
  79 # This also needs to be larger than the longest line in the Alscript command file
  80 # If any lines are longer than MAX_INPUT_LEN then Alscript is likely to fail with 'success'
  81 my $MAX_INPUT_LEN = "MAX_INPUT_LEN ".($foo > 50 ? $foo + 2 : 50);
  82 my $BLOCK_FILE = '#';
  83 unless ($blc eq '-') { $BLOCK_FILE = "BLOCK_FILE $blc"; }
  84 my $OUTPUT_FILE = "OUTPUT_FILE $psout";
  85
  86 # Print the header
  87 print ALSOUT <<END;
  88 SILENT_MODE
  89 $BLOCK_FILE
  90 $MAX_NSEQ
  91 $MAX_INPUT_LEN
  92 $OUTPUT_FILE
  93 PORTRAIT
  94 POINTSIZE 8
  95 IDENT_WIDTH 12
  96 X_OFFSET 2
  97 Y_OFFSET 2
  98 DEFINE_FONT 0 Helvetica         DEFAULT
  99 DEFINE_FONT 1 Helvetica         REL             0.75
 100 DEFINE_FONT 7 Helvetica         REL             0.6
 101 DEFINE_FONT 3 Helvetica-Bold    DEFAULT
 102 DEFINE_FONT 4 Times-Bold        DEFAULT
 103 DEFINE_FONT 5 Helvetica-BoldOblique     DEFAULT
 104 #
 105 DEFINE_COLOUR 3  1 0.62 0.67    # Turquiose
 106 DEFINE_COLOUR 4  1 1 0          # Yellow
 107 DEFINE_COLOUR 5  1 0 0          # Red
 108 DEFINE_COLOUR 7  1 0 1          # Purple
 109 DEFINE_COLOUR 8  0 0 1          # Blue
 110 DEFINE_COLOUR 9  0 1 0          # Green
 111 DEFINE_COLOUR 10 0.41 0.64 1.00 # Pale blue
 112 DEFINE_COLOUR 11 0.41 0.82 0.67 # Pale green
 113 DEFINE_COLOUR 50 0.69 0.18 0.37 # Pink (helix)
 114 DEFINE_COLOUR 51 1.00 0.89 0.00 # Gold (strand)
 115 NUMBER_INT 10
 116 SETUP
 117 #
 118 # Highlight specific residues.
 119 # Avoid highlighting Lupas 'C' predictions by
 120 # limiting the highlighting to the alignments
 121 Scol_CHARS      C 1 1 $end $alignEnd   4
 122 Ccol_CHARS      H ALL    5
 123 Ccol_CHARS      P ALL    8
 124 SURROUND_CHARS LIV   ALL
 125 #
 126 # Replace known structure types with whitespace
 127 SUB_CHARS 1 $top $end $bottom H SPACE
 128 SUB_CHARS 1 $top $end $bottom E SPACE
 129 SUB_CHARS 1 $top $end $bottom - SPACE
 130 END
 131
 132 # Find the range of each of the secondary structure elements
 133 foreach (0..$#jpred) {
 134         if ($jpred[$_]) { # Then it's a secondary structure definition
 135                 my $fsm = fsm($sequences[$_], qw(- E H c C));
 136                 my $line = $_ + 1;
 137
 138                 foreach (@{$fsm->{'E'}}) {
 139                         my ($start, $end) = @{$_};
 140                         my $colour = 51;
 141                         print ALSOUT "STRAND $start $line $end\n";
 142                         print ALSOUT "COLOUR_TEXT_REGION $start $line $end $line $colour\n";
 143                 }
 144
 145                 foreach (@{$fsm->{'H'}}) {
 146                         my ($start, $end) = @{$_};
 147                         my $colour = 50;
 148                         print ALSOUT "HELIX $start $line $end\n";
 149                         print ALSOUT "COLOUR_TEXT_REGION $start $line $end $line $colour\n";
 150                 }
 151         }
 152 }
 153
 154 # Munge out the jpred predictions using a basic FSM
 155 # Pass the sequence and the states that are present.
 156 # Returns HoA's, keys are the states, arrays are start and end
 157 # points of the states.
 158 # Idea taken from Patrick Audley
 159 sub fsm {
 160         my ($sequence, @states) = @_;
 161         my $state;      # The state we're in
 162         my $start = 0;  # The start position of the current state
 163         my %res;        # Hash to store results in
 164         my $pos = 0;    # The position in the FSM
 165
 166         $state = substr $sequence, 0, 1;
 167         foreach my $char (split //, $sequence) {
 168                 $pos++;
 169                 next if $char eq $state; # Same state
 170
 171                 # Otherwise we've changed states
 172                 my $flag = 0;
 173                 foreach (@states) {
 174                         if ($_ eq $char) {
 175                                 push @{$res{$state}}, [ $start, $pos - 1 ];
 176                                 $state = $char;
 177                                 $start = $pos;
 178                                 $flag = 1;
 179                                 last;
 180                         }
 181                 }
 182                 unless ($flag) {
 183                         warn "'$char' not a defined state at position $pos.\n";
 184                 }
 185         }
 186
 187         return \%res;
 188 }