#!/usr/bin/perl
#
# Name: screen.pl
# Purpose: to select shRNAs for a gene set based on the oligo design
#          output produced by Picky. The shRNAs must satisfy the basic
#          shRNA design parameters outlined below. They must also
#          satisfy the user specified uphill free energy ladder $range
#          from 5'- to 3'-end on the siRNA targeting site to facilitate
#          dsRNA unzip and RISC incorporation. The output from this
#          program are cherry-picked shRNAs that can be found to
#          uniquely target some genes in a gene set. Note that not all
#          genes may get an shRNA to target it using this design
#          approach. You may either relax the Picky oligo design and
#          the $range parameters, or to use the pair of Perl programs
#          unscreened_shRNA.pl and merge_screened_shRNA.pl to design
#          shRNAs specifically for a particular target gene.
# Created: 12/24/2006
# Copyright: (c) 2006-2009 Hui-Hsien Chou <hhchou@iastate.edu>
# License: See Picky License.pdf or visit <www.complex.iastate.edu>

use strict;
use warnings;
use Getopt::Long;

my $usage = <<'USAGE';
Usage: screen.pl {optional arguments} PickyOligos.picky > output.html
  -text: generate text output instead; not fancy HTML table
  -loopsize: the shRNA loop length beyond the two 3'-overhang bases
  -range: the range of free energy ladder from 5'- to 3'-end on the target
  -hygro: a DNA fragment unique to your host genome for shRNA loop design
  -seed: use the same seed will generate the same set of loops for shRNAs
USAGE

# default values of the command line arguments; the DNA fragment is
# from a hygromycin oligo designed with C. elegans gene set as
# nontargets and found to be unique to C. elegans; make sure you use
# Picky to find a DNA fragment that's unique to your host genome
my ($text, $range, $loopsize, $hygro, $seed) =
    (0, 12, 7, 'CCGCAAGGAATCGGTCAATACACTACATGGCGTGATTTCATATGCGCGATTGCT', time^$$);

# read in command line parameters
die $usage unless GetOptions("text" => \$text, "range=i" => \$range,
    "loopsize=i" => \$loopsize, "hygro=s" => \$hygro, "seed=i" => \$seed );

# sanity check of arguments
die $usage if @ARGV<1;
die "one Picky output file at a time" if @ARGV>1;
die "valid range is between 5 to 18" if ($range<5 || $range>18);
die "valid loopsize is between 5 to 10" if ($loopsize<5 || $loopsize>10);
die "DNA fragment should be at least 30 bases" if length($hygro)<30;
$hygro = uc($hygro); # change $hygro to uppercase
die "invalid DNA fragment content" if ($hygro=~tr/ACGT/ACGT/)!=length($hygro);

# if the user supplies a $seed value, use it to initialize random
# number so the same set of loop constructs will be repeated if all
# other arguments and the input are the same; this will repeat the
# shRNA designed. The program generates a random value by default.
srand($seed);

# this subroutine formats outlink URL to a sequence database based on
# the input sequence name pattern; you very likely will need to modify
# this subroutine to match the specific database related to your gene
# sets; examples given here are for the Wormbase links.
sub outlink {
  if ($_[0] =~ /^(.+)\|(.+)$/) {
    return "<a href=\"http://www.wormbase.org/db/gene/gene?name=$1;class=Gene\"",
           " target=\"wormbase\">$1</a>,\n",
	   "<a href=\"http://www.wormbase.org/db/seq/sequence?name=$2;class=Transcript\"",
	   " target=\"wormbase\">$2</a>";
  } else {
    return $_[0];
  }
}

# Gibbs free energy nearest-neighbor model parameters from Xia
# et. al. Thermodynamic parameters for an expanded nearest-neighbor
# model for formation of RNA duplexes with Watson-Crick base
# pairs. Biochemistry. 1998 Oct 20; 37(42):14719-35.
my %energies = (
"AA" => -0.93,	#    "AA/TT"
"AC" => -2.24,	#    "AC/TG"
"AG" => -2.08,	#    "AG/TC"
"AT" => -1.10,	#    "AT/TA"
"CA" => -2.11,	#    "CA/GT"
"CC" => -3.26,	#    "CC/GG"
"CG" => -2.36,	#    "CG/GC"
"CT" => -2.08,	#    "CT/GA"
"GA" => -2.35,	#    "GA/CT"
"GC" => -3.42,	#    "GC/CG"
"GG" => -3.26,	#    "GG/CC"
"GT" => -2.24,	#    "GT/CA"
"TA" => -1.33,	#    "TA/AT"
"TC" => -2.35,	#    "TC/AG"
"TG" => -2.11,	#    "TG/AC"
"TT" => -0.93,	#    "TT/AA"
);

# to complement DNA fragment
sub complement {
  my $s = $_[0];
  $s =~ tr/ACGT/TGCA/;
  return $s;
}

# to reverse and complement DNA fragment
sub reverse_complement {
  my $s = reverse($_[0]);
  $s =~ tr/ACGT/TGCA/;
  return $s;
}

# subroutine to find a loop forming segment from the $hygro template;
# the input to this subroutine is the two base overhang at 3'-end on an
# siRNA target site; this routine tries to avoid its reverse complement
# to ensure forming the shRNA loop beginning with it; if no suitable
# loop can be found this routine returns a poly-A loop segment.
sub find_loop {
  my $s = &reverse_complement($_[0]);

  my ($i, $j) = (0, 0);
  for ( ; $i<length($hygro); $i++) {
    $j = index($hygro, $s, $i);
    $j = length($hygro) if $j<0;
    last if $j-$i>=$loopsize;
    $i = $j;
  }
  if ($j-$i>=$loopsize) {
    return $_[0].substr($hygro, $i+int(rand($j-$loopsize-$i)), $loopsize);
  } else {
    return $_[0].('A'x$loopsize);
  }
}

# the following subroutines check some preferred shRNA features
my ($temp_diff, $at_count, $tail_eng, $bind_eng, $preferred_features);
sub check_temp_diff {
  return $temp_diff > 15.0;
}
sub check_AT_count {
  return $at_count > 4;
}
sub check_tail {
  return $tail_eng > -8.0;
}
sub check_bind {
  return $bind_eng < -40.0;
}
sub enough_preferred_features {
  return $preferred_features > 3;
}
my @check_list = ([\&check_temp_diff, \$temp_diff], [\&check_AT_count, \$at_count],
                  [\&check_tail, \$tail_eng], [\&check_bind, \$bind_eng]);

# print the HTML header part which defines the CSS style and Javascript code
print <<HEADER unless $text;
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>shRNA choices</title>
<script type="text/javascript" src="jquery-1.3.2.min.js"></script>
<script type="text/javascript" src="jquery.tablesorter.min.js"></script>
<script type="text/javascript" src="jquery.tablesorter.pager.js"></script>
<script type="text/javascript">
// call the tablesorter plugin
\$(document).ready(function() {
  \$("#shRNAs")
    .tablesorter({
      sortList: [[0,0]],
      headers: { 1: {sorter: false}, 2: {sorter: false}, 3: {sorter: false},
                 4: {sorter: false}, 12: {sorter: false} }
      })
   .tablesorterPager({container: \$("#pager"), positionFixed: false});
});
</script>
<style type="text/css">
<!--
th.header {
  background-image: url(unsorted.gif);
  cursor: pointer;
  background-repeat: no-repeat;
  background-position: bottom center;
  padding-bottom: 10px;
}
th.headerSortUp {
  background-image: url(decreasing.gif);
  background-color: lightyellow;
}
th.headerSortDown {
  background-image: url(increasing.gif);
  background-color: lightyellow;
}
a:link, a:visited, a:hover {
  text-decoration: none;
}
.header {
  font-weight: bold;
  font-size: 10pt;
  text-align: left;
}
.entry {
  font-size: 9pt;
}
.sequence {
  font-family: "courier";
  font-size: 8pt;
}
.sense {
  color: red;
}
.antisense {
  color: green;
}
.loop {
  color: black;
}
.terminator {
  color: brown;
}
.linker {
  color: blue;
}
.bottom {
  color: magenta;
}
img.first, img.prev, img.next, img.last {
  vertical-align: middle;
}
-->
</style>
</head>
<body>
<h1>shRNA choices</h1>
<table border=0>
<tr valign="top"><td>
<table border=0 cellpadding=0>
<tr bgcolor="lightgrey"><td>Grey shRNAs are recommended</td></tr>
<tr bgcolor="lightblue"><td>Light blue indicates preferred values</td></tr>
<tr bgcolor="lightyellow"><td><img src="unsorted.gif"/>Click on these symbols to sort</td></tr>
<tr bgcolor="lightyellow"><td><img src="unsorted.gif"/>Shift-click to sort multiple fields</td></tr>
<tr bgcolor="lightgreen"><td>Regeneration seed: $seed</td</tr>
</table>
</td><td width="50px"></td><td>
<span class="sense">siRNA sense strand</span><br>
<span class="loop">siRNA loop area</span><br>
<span class="antisense">siRNA anti-sense strand</span><br>
<span class="linker">Sticky vector linkers (use yours)</span><br>
<span class="terminator">RNA polymerase terminator (use yours)</span><br>
<span class="bottom">Strand not transcribed</span><br>
</td></tr>
</table>
<br>
<table id="shRNAs" class="tablesorter" border="2" cellpadding="5">
<thead>
<tr class="header">
<th>ID</th>
<th>Target Gene(s)<br>(Click for details)</th>
<th>Target<br>Loc(s)</th>
<th>Nontarget Gene<br>(Click for details)</th>
<th>Nontarget<br>Loc</th>
<th>Target<br>Temp</th>
<th>Nontarget<br>Temp</th>
<th>Temp<br>Diff</th>
<th>3&apos; AT<br>Count</th>
<th>3&apos; Open<br>Energy</th>
<th>Binding<br>Energy</th>
<th>Score<br>(0-4)</th>
<th>shRNA&nbsp;Schematic&nbsp;Drawings&nbsp;(Oligos&nbsp;hidden&nbsp;under&nbsp;URL&nbsp;links)</th>
</tr>
</thead>
<tbody>
HEADER

# loop through Picky designed oligos to find valid shRNA candidates
my ($count, $header_skipped) = ( 1, 0 );
while (<>) {
  # skip the first header line in a .picky file
  unless ($header_skipped) {
    $header_skipped = 1;
    next;
  }
  # split each Picky output line into fields
  my ($oligo, @fs) = split /\t/;
  # scan through each oligo sequence to find shRNA target sites that
  # match the following basic shRNA design parameters:
  #   1) G base at 5'-end of the target site
  #   2) A base at the 3'-end of the target site
  #   3) T base at base position 10
  #   4) None G base at position 13
  # The regular expression pattern below also captures the two overhang
  # bases on both ends of the target site; they are needed to produce
  # the correct shRNA construct.
  while ($oligo =~ /(?=(..G.{8}T..[^G].{5}A..))/g) {
    # check if the shRNA candidate satisfies the monotonic energy ladder
    # requirement that it is increasing toward 3'-end within the $range
    my ($i, $range_sum) = (2, 0.0);
    for (my $j=0; $j<$range; $i++, $j++) {
      $range_sum += $energies{substr($1, $i, 2)};
    }
    for ( ; $i<20; $i++) {
      my $eng = $energies{substr($1, $i, 2)};
      last if $eng<$range_sum/$range;
      $range_sum -= $energies{substr($1, $i-$range, 2)};
      $range_sum += $eng;
    }
    next if $i<20; # if not, skip this shRNA candidate
    # save shRNA candidate match location on the oligo
    my $loc = length($`);
    # The following four additional values are then computed for
    # each target site:
    #   1) The Gibbs free energy at 37 C when the antisense siRNA binds to
    #      the target site. Note that we didn't include the initiation
    #      and A/U end compensation energy terms since they will be the
    #      same for all candidates matching the shRNA design pattern
    #      thus they do not change the shRNA preference order.
    #   2) The Gibbs free energy at 37 C for the last six bases on the
    #      3'-end of the target site to unwind
    #   3) The AT base count at the last six bases of the 3'-end of the
    #      target site
    #   4) The Picky computed melting temperature difference between
    #      the oligo target(s) and closest nontarget. This value tells
    #      how specific a shRNA candidate chosen from the oligo will be.
    # These values are output and can be used to rate the shRNAs.

    $temp_diff = $fs[1]-$fs[5];     # melting temperature difference
    my $tail = substr($1, -8, -2);  # last six bases on the target site;
                                    # the two base overhang at 3'-end
                                    # does not participate in binding
    $at_count = ($tail =~ tr/AT/AT/);  # count AT bases
    # calculate tail unwinding energy
    $tail_eng = 0.0;
    for ($i=0; $i<length($tail)-1; $i++) {
      $tail_eng += $energies{substr($tail, $i, 2)};
    }
    # calculate siRNA/target binding energy; note the two 5' overhang
    # bases are participating in the binding
    $bind_eng = 0.0;
    for ($i=0; $i<20; $i++) {
      $bind_eng += $energies{substr($1, $i, 2)};
    }
    # construct the shRNA parts
    my $sens = substr($1, 2, 19);
    my $loop = &find_loop(substr($1, -2));
    my $anti = &reverse_complement(substr($1, 0, 21));
    my $oligo1 = $sens.$loop.$anti;
    my $oligo2 = &reverse_complement($sens.$loop.$anti);
    # generate output
    if ($text) {
      # text output in FASTA format; target(s) listed on the header line
      print ">shRNA_${count}_S";
      for ($i=6; $i<@fs; $i+=4) {
	print " $fs[$i]:", $fs[$i+2]+$loc;
      }
      print "\n$oligo1\n";
      print ">shRNA_${count}_A\n$oligo2\n";
    } else {
      # HTML output in table format
      print "<tr class=\"entry\"><td>\n";
      # shRNA ID
      print "$count</td><td>\n";
      # target gene list
      for ($i=6; $i<@fs; $i+=4) {
	print outlink($fs[$i]), "<br>\n";
      }
      # target locations
      print "</td><td>\n";
      for ($i=6; $i<@fs; $i+=4) {
	print $fs[$i+2]+$loc, "<br>\n";
      }
      # closest nontarget gene (there is only one)
      print "</td><td>\n", outlink($fs[2]), "\n";
      # approximate nontarget location
      print "</td><td>\n", $fs[4]+$loc, "\n";
      # target melting temperature
      print "</td><td>\n", $fs[1], "\n";
      # nontarget melting temperature
      print "</td><td>\n", $fs[5], "\n</td>\n";
      # check and print the following shRNA features
      #   melting temperature difference
      #   3'-end AT count
      #   3'-end tail unwind energy
      #   target binding energy
      $preferred_features = 0;
      for (@check_list) {
	if (&{$_->[0]}) {
	  print '<td bgcolor="lightblue">'; $preferred_features++;
	} else {
	  print '<td>';
	}
	print ${$_->[1]}, "</td>\n";
      }
      # preferred feature score
      print '<td>', $preferred_features, "</td>\n";
      if (&enough_preferred_features) {
	print "<td class=\"sequence\" bgcolor=\"lightgrey\">\n";
      } else {
	print "<td class=\"sequence\">\n";
      }
      print "<a href=\"http://$oligo1\"><span class=\"linker\">5&apos;-</span>\n";
      print '<span class="sense">', $sens, "</span>\n";
      print '<span class="loop">', $loop, "</span>\n";
      print '<span class="antisense">', $anti, "</span>\n";
      print '<span class="linker">-3&apos;-</span><span class="terminator">X</span></a><br>', "\n";
      print "<a href=\"http://$oligo2\"><span class=\"linker\">3&apos;-</span>\n";
      print '<span class="bottom">', &complement($sens), "</span>\n";
      print '<span class="bottom">', &complement($loop), "</span>\n";
      print '<span class="bottom">', &complement($anti), "</span>\n";
      print '<span class="linker">-5&apos;-</span><span class="terminator">Y</span></a>', "\n";
      print "</td></tr>\n";
    }
    $count++; # increases the shRNA counter
  }
}

# reduce counter to actual total
$count--;

# get the current date and time
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
my $now = sprintf "%4d-%02d-%02d %02d:%02d:%02d",$year+1900,$mon+1,$mday,$hour,$min,$sec;

# print the HTML footer part
print <<FOOTER unless $text;
</tbody>
</table>
<br>
<form id="pager">
<img src="first.png" class="first" alt="First"/>
<img src="prev.png" class="prev" alt="Previous"/>
<input type="text" class="pagedisplay" title="Jump to"/>
<img src="next.png" class="next" alt="Next"/>
<img src="last.png" class="last" alt="Last"/>
<select class="pagesize" title="shRNAs per page">
  <option value="10">10</option>
  <option value="20">20</option>
  <option value="30">30</option>
  <option value="40">40</option>
  <option value="50">50</option>
  <option value="100">100</option>
  <option value="$count">ALL</option>
</select>
&nbsp;&nbsp;&nbsp; Total <span class="sense">$count</span> shRNAs
produced by Picky shRNA design suite on $now.
</form>
</body>
</html>
FOOTER

print STDERR "$count shRNA designed; use seed $seed for regeneration.\n";
