#!/usr/bin/perl

$args=@ARGV;
$BIN_VIENNA = "/usr/local/bin";
$BIN_PACMIT = "/home/ray/Documents/Research/MicroRNA/Code/Pacmit_public_v2";


if ($args ne 3){
  print "######################################################################################\n";
  print "#                                    rnaplfold.pl                                    #\n";
  print "######################################################################################\n";
  print "# This script executes RNAplfold for each one of the sequences in a .fa fasta file   #\n";
  print "#                                                                                    #\n";
  print "# execution: > rnaplfold.pl [.fa 3'UTR file] [length of pieces W]                    #\n";
  print "#                           [maximum length of pairing L]                            #\n";
  print "#                                                                                    #\n";
  print "# Output:  for each 3'UTR we obtain a file with the P_u values in the current dir.   #\n";
  print "#                                                                                    #\n";
  print "# Input:                                                                             #\n";
  print "# Argument 1 => fasta 3'UTR file: e.g. example_utr.fa                                #\n";
  print "# Argument 2 => Length of window for the local folding: e.g. 80                      #\n";
  print "# Argument 3 => Maximum distance allowed between pairing nucleotides: e.g. 40        #\n";
  print "######################################################################################\n";
  print "\n";
}else{

  $fasta = $ARGV[0];
  $w_fix = $ARGV[1];
  $l_fix = $ARGV[2];

  $fasta =~ s/.fa//;
  system("$BIN_PACMIT/fasta2vienna.awk $fasta.fa > $fasta.vienna");
  open(UTRS, "$fasta.vienna");
  local $/;
  $utrs_file = <UTRS>;
  close(UTRS);
  @utrs = split(">", $utrs_file);
  shift(@utrs);

  for $utr (@utrs){
    $W=$w_fix;
    $L=$l_fix;
    @gene=split("\n", $utr);
    @title=split(" ", $gene[0]);
    $utr_name=$title[0];
    $utr_seq=$gene[1];
    $l_utr= length($utr_seq);
    if($l_utr < $w_fix) {$W=$l_utr};
    if($l_utr < $l_fix) {$L=$l_utr};
    @output_name=($utr_name, "_lunp");
    $output=join('', @output_name);
    $column=5;
    system("echo $utr_seq | $BIN_VIENNA/RNAplfold -u 4 -W $W -L $L");
    system("awk \'{if(\$1>=4) print \$1,\$$column}\' plfold_lunp > $utr_name.fa_4");
    system("rm plfold_dp.ps  plfold_lunp");
  }
}
