#!/usr/bin/perl
#
# Name: merge_screened_shRNA.pl
# Purpose: to merge shRNA candidates screened by Picky against various
#          nontarget gene sets. There is no limit on the number of
#          nontarget gene sets but duplicative nontarget gene sets
#          should be avoided. The program rigorously check for
#          consistency among all input data. Thus, if two genes were
#          found to have the same name but different DNA sequences
#          this program may consider that a fatal error even if the
#          two genes are from two different screening output. The
#          shRNA candidates were previously selected by
#          unscreened_shRNA.pl. The output from this program includes
#          many parameters that can be used to prioritize the shRNAs.
# Created: 2/19/2008
# Copyright: (c) 2008, 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: merge_screened_shRNA.pl {opt. arguments} unscreened_shRNAs.seq *.picky > output.html
  -text: generate text output instead; not fancy HTML table
  -loopsize: the shRNA loop length beyond the two 3'-overhang bases
  -hygro: a DNA fragment unique to your host genomes 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 both human and mouse gene sets as
# nontargets and found to be unique to them; make sure you use Picky
# to find a DNA fragment that's unique to your host genome
my ($text, $loopsize, $hygro, $seed) =
    (0, 7, 'GGCGTGGATATGTCCTGCGGGTAAATAGCTGCGCCGATGGTTTCTACAAAGATCGT', time^$$);

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

# sanity check of arguments
die $usage if @ARGV<2;
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 Ensemble links.
sub outlink {
  if ($_[0] =~ /^ENST\d/) {
    return '<a href="http://may2009.archive.ensembl.org/Homo_sapiens/geneview?transcript=', $_[0],
           '" target="ensembl">', $_[0], "</a>";
  } elsif ($_[0] =~ /^ENSMUST\d/) {
    return '<a href="http://jul2009.archive.ensembl.org/Mus_musculus/geneview?transcript=', $_[0],
           '" target="ensembl">', $_[0], "</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, $multi_loc, $min_range, $tail_eng,
    $bind_eng, $preferred_features);
sub check_temp_diff {
  return $temp_diff > 10.0;
}
sub check_AT_count {
  return $at_count > 3;
}
sub check_multi {
  return $multi_loc > 1;
}
sub check_min_range {
  return $min_range < 16;
}
sub check_tail {
  return $tail_eng > -10.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_min_range, \$min_range],
                  [\&check_multi, \$multi_loc], [\&check_bind, \$bind_eng]);

# read in the original unscreened shRNAs generated by unscreened_shRNA.pl
my %list;
my @ordered_list;
open FILE, $ARGV[0] or die "can't open file $ARGV[0]";
while (<FILE>) {
  chomp;
  my ($shRNA, $overhang) = split;
  push @ordered_list, $list{$shRNA} = [ $shRNA, $overhang ];
}
close FILE;
shift @ARGV;

# read in all Picky screening output against different nontarget gene sets;
# the many 'die' statements perform rigorous checks of data consistency.
for (@ARGV) {
  my $header_skipped = 0;
  open FILE, $_ or die "can't open file $_";
  while (<FILE>) {
    unless ($header_skipped) {
      $header_skipped=1;
      next;
    }
    chomp;
    my @fields = split/\t/;
    die "unknown shRNA candidate $fields[0]" unless exists($list{$fields[0]});
    die "sequence target of shRNA candidate $fields[0] is unknown" if
	$fields[7] eq "unknown!";
    my $ptr = $list{$fields[0]};
    if (defined($ptr->[2])) { # previously we have seen this shRNA candidate
      my $hsh = $ptr->[2]; # so load its target information hash
      # target melting temperatures can differ slightly but not too
      # much between different gene sets due to dangling end base difference
      die "unmatched target temperatures $fields[2]!=",
                   $hsh->{"melting\ttemperature"}, " for shRNA $fields[0]"
	  if abs($fields[2]-$hsh->{"melting\ttemperature"})>2.0;
      # any new shRNA targets will be stored in %new_targets so their
      # locations will NOT be checked until they are added to the hash
      my %new_targets;
      # now check the target list to see if ALL target locations match
      for (my $i=7; $i<@fields; $i+=4) {
	if (exists($hsh->{$fields[$i]})) { # check existing target locations
	  die "unmatched target locations $fields[$i]:$fields[$i+2] for shRNA $fields[0]"
	      unless exists($hsh->{$fields[$i]}{$fields[$i+2]});
	} else {
	  # a new target is found; its locations will just be added
	  unless (exists($new_targets{$fields[$i]})) {
	    $new_targets{$fields[$i]} = { };
	  }
	  $new_targets{$fields[$i]}{$fields[$i+2]} = 1;
	}
      }
      # finally copy new target information into the hash for this shRNA
      for my $target (keys %new_targets) {
	$hsh->{$target} = $new_targets{$target};
      }
      $hsh = $ptr->[3]; # now load its nontarget information hash
      if (exists($hsh->{$fields[3]})) {
	die "unmatched nontarget temperature $fields[6]!=",
	                       $hsh->{$fields[3]}[1], " for shRNA $fields[0]"
	    if abs($fields[6]-$hsh->{$fields[3]}[1]) > 0.01;
	die "unmatched nontarget location $fields[3]:$fields[5] for shRNA $fields[0]"
	    if $fields[5] != $hsh->{$fields[3]}[0];
      } elsif ($fields[3] ne "not seen!") {
	$hsh->{$fields[3]} = [ @fields[5, 6] ]; # add the new nontarget info
      }
    } else {
      # this is the first time we run into this shRNA candidate;
      # create a new anonymous hash to store its target information
      my $hsh = $ptr->[2] = { "melting\ttemperature" => $fields[2] };
      # store all targets and their locations using hash to hashes
      for (my $i=7; $i<@fields; $i+=4) {
	unless (exists($hsh->{$fields[$i]})) {
	  $hsh->{$fields[$i]} = { }; # each new target has a location hash
	}
	$hsh->{$fields[$i]}{$fields[$i+2]} = 1; # add to location hash
      }
      # create a new anonymous hash to store its nontarget information
      $ptr->[3] = $fields[3] eq "not seen!" ? { } : { $fields[3] => [ @fields[5, 6] ] };
    }
  }
  close FILE;
}

# 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}, 14: {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(s)<br>(Click for details)</th>
<th>Nontarget<br>Loc(s)</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>Min<br>Range</th>
<th>Multi<br>Loc</th>
<th>Binding<br>Energy</th>
<th>Score<br>(0-6)</th>
<th>shRNA&nbsp;Schematic&nbsp;Drawings&nbsp;(Oligos&nbsp;hidden&nbsp;under&nbsp;URL&nbsp;links)</th>
</tr>
</thead>
<tbody>
HEADER

my $count=1;
for (@ordered_list) {

  die "unscreened shRNA candidate $_->[0]" unless defined($_->[2]);

  # The following six parameter values are computed for each shRNA
  # according to its 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 minimum difference between the Picky computed melting
  #      temperatures of the shRNA target(s) and their closest
  #      nontarget(s). This value tells how specific a shRNA will be.
  #   5) The minimum range within which a shRNA candidate satisfies
  #      the monotonic energy ladder requirement. This value tells how
  #      easily it will unwind and be incorporated into RISC complex.
  #   6) The maximum count of target sites on each target gene. If the
  #      same shRNA targets a gene multiple times, it's considered
  #      more potent. Note that multiple targeted genes is not
  #      considered the same as multiple target sites.
  # These values are output and can be used to rate the shRNAs.

  $temp_diff = 100.0;  # find minimal temperature difference
  for my $ary (values %{$_->[3]}) {
    my $t_d = $_->[2]{"melting\ttemperature"}-$ary->[1];
    $temp_diff = $t_d if $t_d<$temp_diff;
  }
  $temp_diff = int($temp_diff*100)/100; # we don't want extra digits
  my $tail = substr($_->[0], -6);  # last six bases on the target site;
  $at_count = ($tail =~ tr/AT/AT/);     # count AT bases
  # calculate tail unwinding energy
  $tail_eng = 0.0;
  for (my $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 (my $i=0; $i<20; $i++) {
    $bind_eng += $energies{substr($_->[0], $i, 2)};
  }
  # discover the minimum range that allows the shRNA candidate to
  # satisfy the monotonic energy ladder requirement such that it is
  # always increasing toward 3'-end when averaged in the $min_range
  $min_range=18;
  while ($min_range>0) {
    $min_range--;
    my ($i, $range_sum) = (2, 0.0);
    for (my $j=0; $j<$min_range; $i++, $j++) {
      $range_sum += $energies{substr($_->[0], $i, 2)};
    }
    for ( ; $i<20; $i++) {
      my $eng = $energies{substr($_->[0], $i, 2)};
      last if $eng<$range_sum/$min_range;
      $range_sum -= $energies{substr($_->[0], $i-$min_range, 2)};
      $range_sum += $eng;
    }
    last if $i<20; # if not, skip this shRNA candidate
  }
  $min_range++;
  # determine the maximum number of multiple target sites on each gene
  $multi_loc = 0;
  for my $hsh (values %{$_->[2]}) {
    next unless ref $hsh;
    my $c = keys %{$hsh};
    $multi_loc = $c if $c>$multi_loc;
  }

  # construct the shRNA parts
  my $sens = substr($_->[0], 2, 19);
  my $loop = &find_loop($_->[1]);
  my $anti = &reverse_complement($_->[0]);
  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 my $target (sort keys %{$_->[2]}) {
      next if $target eq "melting\ttemperature";
      print " $target:";
      my $c = 0;
      for my $loc (sort {$a <=> $b} keys %{$_->[2]{$target}}) {
	print "," if $c;
	print "$loc";
	$c++;
      }
    }
    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 my $target (sort keys %{$_->[2]}) {
      next if $target eq "melting\ttemperature";
      for (my $i=0; $i<keys %{$_->[2]{$target}}; $i++) {
	print outlink($target), "<br>\n";
      }
    }
    # target locations
    print "</td><td>\n";
    for my $target (sort keys %{$_->[2]}) {
      next if $target eq "melting\ttemperature";
      for my $loc (sort {$a <=> $b} keys %{$_->[2]{$target}}) {
	print "$loc<br>\n";
      }
    }
    # closest nontarget genes (there is only one for each screened
    # gene set and duplicative ones will not be shown)
    print "</td><td>\n";
    for my $nontarget (sort keys %{$_->[3]}) {
      print outlink($nontarget), "<br>\n";
    }
    # approximate nontarget locations
    print "</td><td>\n";
    for my $nontarget (sort keys %{$_->[3]}) {
      print $_->[3]{$nontarget}[0], "<br>\n";
    }
    # target melting temperature
    print "</td><td>", $_->[2]{"melting\ttemperature"}, "\n";
    # nontarget melting temperatures
    print "</td><td>\n";
    for my $nontarget (sort keys %{$_->[3]}) {
      print $_->[3]{$nontarget}[1], "<br>\n";
    }
    print "</td>\n";
    # check and print the following shRNA features
    #   melting temperature difference
    #   3'-end AT count
    #   3'-end tail unwind energy
    #   minimum range for satisfying energy ladder
    #   multiple target sites on the same gene
    #   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";
