#!/usr/bin/perl
#
# Name: unscreen_shRNA.pl
# Purpose: to design shRNAs for a set of target genes based on the
#          basic shRNA design parameters. The output from this program
#          should be screened using Picky with some nontarget genes to
#          ensure gene specificity and to collect the closest
#          nontarget prediction for each shRNA candidate. After Picky
#          screening, the output from this program and *.picky files
#          can be merged together by merge_screened_shRNA.pl.
# Created: 2/18/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;

die "Usage: unscreened_shRNA.pl gene_targets.seq > unscreened_shRNAs.seq\n"
  unless @ARGV==1;

# Gibbs free energey 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"
);

# keep a set of designed shRNAs to avoid duplicates
my %set;

# subroutine to scan through each target 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 later by merge_screened_shRNA.pl, but
# note that the 3' overhang does not participate in siRNA binding with
# target, so it is disconnected from the shRNA (Picky ignores such
# short bases).
sub scan {
  my $seq = shift;
  $seq = uc($seq);
  while ($seq =~ /(?=(..G.{8}T..[^G].{5}A..))/g) {
    my ($shRNA, $overhang) = (substr($1, 0, 21), substr($1, -2));
    unless (exists($set{$shRNA}) || ($1=~tr/ACGT/ACGT/)!=length($1)) {
      print "$shRNA $overhang\n";
    }
    $set{$shRNA}++;
  }
}

# main routine to read each target sequences into $seq and call scan
my ($count, $seq) = ( 0, "");
while (<>) {
  if (/^>/) { # FASTA header
    $count++;
    &scan($seq) if length($seq);
    $seq = "";
  } else {
    chomp;
    s/\r//g; # remove DOS or MacOS end-of-line characters just in case
    $seq .= $_;
  }
}
scan($seq) if length($seq);

print STDERR "$count target sequences seen; ",
    scalar keys(%set), " shRNA candidates generated.\n";
