sqp_seq.pm

#!/usr/bin/perl
#
# sqp_seq.pm
# Eric Nawrocki
# EPN, Tue Mar 19 13:29:58 2019 [incept, in vadr]
# EPN, Tue Jul  2 11:47:41 2019 [migrated from vadr's epn-seq.pm (as of commit 69b003d)]]
# version: 0.10
#
use strict;
use warnings;
use Time::HiRes qw(gettimeofday);

# NOTE: do not add any 'require' statements here, e.g. 'require
# sqp_utils.pm' because the program that uses sequip must handle that
# so each program can specify sequip from a specific directory defined
# by a specific environment variable. This is how, for example,
# ribovore can require a specific version of sequip on the same file
# system that has vadr installed with a potentially different version
# of sequip.

#################################################################
# Subroutine:  seq_SqstringAddNewlines()
# Incept:      EPN, Thu Mar 14 06:12:11 2019
#
# Purpose:     Add newlines to $sqstring after every $linelen
#              characters and return result.
#
# Arguments: 
#   $sqstring: the sequence string
#   $linelen:  interval for newlines
# 
# Returns:     $sqstring with newlines inserted every $linelen 
#              characters and at end of string.
# 
# Dies:        Never.
#
################################################################# 
sub seq_SqstringAddNewlines { 
  my $nargs_expected = 2;
  my $sub_name = "seq_SqstringAddNewlines";
  if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); } 
  my ($sqstring, $linelen) = @_;

  my $retstr = "";
  my $sqpos = 0;
  my $sqlen = length($sqstring);

  while($sqpos < $sqlen) { 
    $retstr .= substr($sqstring, $sqpos, $linelen) . "\n";
    $sqpos += $linelen;
  }

  return $retstr;
}

#################################################################
# Subroutine: seq_SqstringCapitalize
# Incept:     EPN, Fri Mar 15 13:32:36 2019
# 
# Purpose:    Capitalize a string in place.
# 
# Arguments:
#   $sqstring_R: REF to sequence string to capitalize
#
# Returns:    void
# 
# Dies:       never
#
#################################################################
sub seq_SqstringCapitalize {
  my $sub_name = "seq_SqstringCapitalize";
  my $nargs_expected = 1;
  if(scalar(@_) != $nargs_expected) { die "ERROR $sub_name entered with wrong number of input args" }
 
  my ($sqstring_R) = @_;
  
  $$sqstring_R =~ tr/a-z/A-Z/;
  return;
}

#################################################################
# Subroutine: seq_SqstringDnaize
# Incept:     EPN, Fri Mar 15 13:33:39 2019
# 
# Purpose:    Convert a RNA/DNA sqstring to DNA in place.
# 
# Arguments:
#   $sqstring_R: REF to sequence string to capitalize
#
# Returns:    void
# 
# Dies:       never
#
#################################################################
sub seq_SqstringDnaize {
  my $sub_name = "seq_SqstringDnaize";
  my $nargs_expected = 1;
  if(scalar(@_) != $nargs_expected) { die "ERROR $sub_name entered with wrong number of input args" }
 
  my ($sqstring_R) = @_;
  
  $$sqstring_R =~ tr/Uu/Tt/;
  return;
}


#################################################################
# Subroutine: seq_SqstringReverseComplement
# Incept:     EPN, Fri Mar 15 15:35:10 2019
# 
# Purpose:    Reverse complement a RNA/DNA sqstring in place.
# 
# Arguments:
#   $sqstring_R: REF to reverse complement
#
# Returns:    void
# 
# Dies:       never
#
#################################################################
sub seq_SqstringReverseComplement {
  my $sub_name = "seq_SqstringReverseComplement";
  my $nargs_expected = 1;
  if(scalar(@_) != $nargs_expected) { die "ERROR $sub_name entered with wrong number of input args" }
 
  my ($sqstring_R) = @_;

  # DNA-ize it
  seq_SqstringDnaize($sqstring_R);
  # reverse it 
  $$sqstring_R = reverse $$sqstring_R;
  # complement it
  $$sqstring_R =~ tr/ACGTRYMKHBVDacgtrymkhbvd/TGCAYRKMDVBHtgcayrkmdvbh/;
  # see esl_alphabet.c::set_complementarity()
  # note that S, W, N are omitted as they are their own complements

  return;
}


#################################################################
# Subroutine: seq_SqstringDiffSummary
# Incept:     EPN, Fri Mar 15 13:35:28 2019
# 
# Purpose:    Return a string summarizes the differences between
#             two sqstrings.
# 
# Arguments:
#   $sqstring1: sqstring 1
#   $sqstring2: sqstring 2
#
# Returns:    String with N newlines for N differences.
#             "" if sqstrings are identical.
# 
# Dies:       never
#
#################################################################
sub seq_SqstringDiffSummary {
  my $sub_name = "seq_SqstringDiffSummary";
  my $nargs_expected = 2;
  if(scalar(@_) != $nargs_expected) { die "ERROR $sub_name entered with wrong number of input args" }
 
  my ($sqstring1, $sqstring2) = @_;

  if(! defined $sqstring1) { 
    return "sequence 1 is undefined\n"; 
  }
  if(! defined $sqstring2) { 
    return "sequence 2 is undefined\n"; 
  }
  seq_SqstringCapitalize(\$sqstring1);
  seq_SqstringCapitalize(\$sqstring2);
  seq_SqstringDnaize(\$sqstring1);
  seq_SqstringDnaize(\$sqstring2);
  if($sqstring1 eq $sqstring2) { 
    return "";
  }
  my $n1 = length($sqstring1); 
  my $n2 = length($sqstring2); 
  if($n1 != $n2) {
    return "sequence lengths mismatch ($n1 != $n2)\n"; 
  }
  my $ret_str = "";
  my @sqstring1_A = split("", $sqstring1); 
  my @sqstring2_A = split("", $sqstring2); 
  my $n = ($n1 > $n2) ? $n1 : $n2;
  for(my $i = 0; $i < $n; $i++) { 
    if($i >= $n1) { 
      $ret_str .= " char " . ($i+1) . " seq1: off-end seq2: " . $sqstring2_A[$i] . "\n";
    }
    elsif($i >= $n2) { 
      $ret_str .= " char " . ($i+1) . " seq1: " . $sqstring1_A[$i] . " seq2: off-end\n";
    }
    elsif($sqstring1_A[$i] ne $sqstring2_A[$i]) { 
      $ret_str .= " char " . ($i+1) . " seq1: " . $sqstring1_A[$i] . " seq2: " . $sqstring2_A[$i] . "\n";
    }
  }

  return $ret_str;
}

#################################################################
# Subroutine: seq_CodonValidateStartCapDna()
# Incept:     EPN, Sat Feb 23 10:01:55 2019
# 
# Purpose:    Given an already capitalized DNA codon, return '1' 
#             if it's a valid start codon, else return 0.
#
# Args:
#  $codon:    the codon
#  $tt:       the translation table ('1' for standard)
#  $atg_only: return '1' only for 'ATG'
#
# Returns:    '1' if $codon is valid start codon for translation table $tt
#             else '0'
#
#################################################################
sub seq_CodonValidateStartCapDna {
  my $sub_name = "seq_CodonValidateStartCapDna";
  my $nargs_exp = 3;
  if(scalar(@_) != $nargs_exp) { die "ERROR $sub_name entered with wrong number of input args"; }

  my ($codon, $tt, $atg_only) = @_;

  if($atg_only) { 
    if($codon eq "ATG") {
      return 1; 
    }
    else { 
      return 0;
    }
  }
  else { 
    if($tt == 1) { # TTG CTG ATG
      if($codon =~ m/^[TCAYWMH]TG$/) { 
        return 1;
      }
    }
    elsif($tt == 2) { # ATT ATC ATA ATG GTG
      if(($codon =~ m/^AT[ACGTRYSWKMBDHVN]$/) || # ATN
         ($codon =~ m/^[AGR]TG$/)) { 
        return 1;
      }
    }
    elsif($tt == 3) { # ATA ATG GTG
      if(($codon =~ m/^AT[AGR]$/) || 
         ($codon =~ m/^[AGR]TG$/)) { 
        return 1;
      }
    }
    elsif($tt == 4) { # TTA TTG CTG ATT ATC ATA ATG GTG
      if(($codon =~ m/^AT[ACGTRYSWKMBDHVN]$/) || # ATN
         ($codon =~ m/^TT[AGR]$/) || 
         ($codon =~ m/^[CGS]TG$/)) { 
        return 1;
      }
    }
    elsif($tt == 5) { # TTG ATT ATC ATA ATG GTG
      if(($codon =~ m/^AT[ACGTRYSWKMBDHVN]$/) || # ATN
         ($codon =~ m/^[TGK]TG$/)) { 
        return 1;
      }
    }
    elsif(($tt == 6) || ($tt == 10) || ($tt == 14) || ($tt == 16) || ($tt == 22) || ($tt == 27) || ($tt == 28) || ($tt == 29) || ($tt == 30) || ($tt == 31)) { #ATG
      if(($codon eq "ATG")) {
        return 1;
      }
    }
    elsif(($tt == 9) || ($tt == 21)) { # ATG GTG
      if($codon =~ m/^[AGR]TG$/) { 
        return 1;
      }
    }
    elsif($tt == 11) { # TTG CTG ATT ATC ATA ATG GTG
      if(($codon =~ m/^AT[ACGTRYSWKMBDHVN]$/) || # ATN
         ($codon =~ m/^[TCGYKSB]TG$/)) { 
        return 1;
      }
    }
    elsif(($tt == 12) || ($tt == 26)) { # CTG ATG
      if($codon =~ m/^[CAM]TG$/) { 
        return 1;
      }
    }
    elsif($tt == 13) { # TTG ATA ATG GTG
      if(($codon =~ m/^[TAGWKRD]TG$/) || 
         ($codon =~ m/^AT[AGR]$/)) { 
        return 1;
      }
    }
    elsif($tt == 23) { # ATT ATG GTG
      if(($codon =~ m/^AT[TGK]$/) || 
         ($codon =~ m/^[AGR]TG$/)) { 
        return 1;
      }
    }
    elsif(($tt == 24) || ($tt == 33)) { # TTG CTG ATG GTG
      if($codon =~ m/^[ACGTRYSWKMBDHVN]TG$/) { # NTG
        return 1;
      }
    }
    elsif($tt == 25) { # TTG ATG GTG
      if($codon =~ m/^[TAGWKRD]TG$/) { 
        return 1;
      }
    }
  }

  return 0;
}

#################################################################
# Subroutine: seq_CodonValidateStopCapDna()
# Incept:     EPN, Mon Mar 14 13:47:57 2016
# 
# Purpose:    Given an already capitalized DNA codon, return '1' 
#             if it's a valid stop codon, else return 0.
#             Any codon with an ambiguous nt will return 0.
#
# Args:
#  $codon:  the codon
#  $tt:     the translation table ('1' for standard)
#
# Returns:    '1' if codon is valid stop, else '0'
#
# Ref: https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes
#################################################################
sub seq_CodonValidateStopCapDna {
  my $sub_name = "seq_CodonValidateStopCapDna";
  my $nargs_exp = 2;
  if(scalar(@_) != $nargs_exp) { die "ERROR $sub_name entered with wrong number of input args"; }

  my ($codon, $tt) = @_;
  
  if(($tt == 1) || ($tt == 11) || ($tt == 12) || ($tt == 26) || ($tt == 28)) { # TAA TAG TGA
    if(($codon =~ m/^TA[AGR]$/) || 
       ($codon =~ m/^T[AGR]A$/)) { 
      return 1;
    }
  }
  elsif($tt == 2) { # TAA TAG AGA AGG
    if(($codon =~ m/^TA[AGR]$/) || 
       ($codon =~ m/^AG[AGR]$/)) { 
      return 1;
    }
  }
  elsif(($tt == 3) || ($tt == 4) || ($tt == 5) || ($tt == 9) || ($tt == 10) || ($tt == 13) || ($tt == 21) || ($tt == 24) || ($tt == 25) || ($tt == 31)) { # TAA TAG
    if($codon =~ m/^TA[AGR]$/) {
      return 1;
    }
  }
  elsif(($tt == 6) || ($tt == 27) || ($tt == 29) || ($tt == 30)) { # TGA
    if($codon eq "TGA") { 
      return 1;
    }
  }
  elsif(($tt == 14) || ($tt == 33)) { # TAG
    if($codon eq "TAG") { 
      return 1;
    }
  }
  elsif($tt == 16) { # TAA TGA
    if($codon =~ m/^T[AGR]A$/) { 
      return 1;
    }
  }
  elsif($tt == 22) { # TCA TAA TGA
    if($codon =~ m/^T[ACGMSRV]A$/) { 
      return 1;
    }
  }
  elsif($tt == 23) { # TTA TAA TAG TGA
    if(($codon =~ m/^T[TAGWKRD]A$/) || 
       ($codon =~ m/^TA[AGR]$/)) { 
      return 1;
    }
  }
  
  return 0;
}

#################################################################
# Subroutine: seq_Overlap()
# Incept:     EPN, Mon Mar 14 13:47:57 2016 [dnaorg.pm]
#
# Purpose:    Calculate number of nucleotides of overlap between
#             two regions.
#
# Args:
#  $start1: start position of hit 1 (must be <= $end1)
#  $end1:   end   position of hit 1 (must be >= $end1)
#  $start2: start position of hit 2 (must be <= $end2)
#  $end2:   end   position of hit 2 (must be >= $end2)
#  $FH_HR:  REF to hash of file handles, including "log" and "cmd"
#
# Returns:  Two values:
#           $noverlap:    Number of nucleotides of overlap between hit1 and hit2, 
#                         0 if none
#           $overlap_reg: region of overlap, "" if none
#
# Dies:     if $end1 < $start1 or $end2 < $start2.
#
#################################################################
sub seq_Overlap {
  my $sub_name = "seq_Overlap";
  my $nargs_exp = 5;
  if(scalar(@_) != 5) { die "ERROR $sub_name entered with wrong number of input args"; }

  my ($start1, $end1, $start2, $end2, $FH_HR) = @_; 

  # printf("in $sub_name $start1..$end1 $start2..$end2\n");

  if($start1 > $end1) { ofile_FAIL("ERROR in $sub_name start1 > end1 ($start1 > $end1)", 1, $FH_HR); }
  if($start2 > $end2) { ofile_FAIL("ERROR in $sub_name start2 > end2 ($start2 > $end2)", 1, $FH_HR); }

  # Given: $start1 <= $end1 and $start2 <= $end2.
  
  # Swap if nec so that $start1 <= $start2.
  if($start1 > $start2) { 
    my $tmp;
    $tmp   = $start1; $start1 = $start2; $start2 = $tmp;
    $tmp   =   $end1;   $end1 =   $end2;   $end2 = $tmp;
  }
  
  # 3 possible cases:
  # Case 1. $start1 <=   $end1 <  $start2 <=   $end2  Overlap is 0
  # Case 2. $start1 <= $start2 <=   $end1 <    $end2  
  # Case 3. $start1 <= $start2 <=   $end2 <=   $end1
  if($end1 < $start2) { return (0, ""); }                                           # case 1
  if($end1 <   $end2) { return (($end1 - $start2 + 1), ($start2 . "-" . $end1)); }  # case 2
  if($end2 <=  $end1) { return (($end2 - $start2 + 1), ($start2 . "-" . $end2)); }  # case 3
  die "ERROR in $sub_name, unforeseen case in $sub_name, $start1..$end1 and $start2..$end2";

  return; # NOT REACHED
}

#################################################################
# Subroutine: seq_StripVersion()
# Incept:     EPN, Thu Feb 11 14:25:52 2016
#
# Purpose:    Given a ref to an accession.version string, remove the version.
#
# Arguments: 
#   $accver_R: ref to accession version string to remove version from
#
# Returns:    Nothing, $$accver_R has version removed
#
# Dies:       never
#################################################################
sub seq_StripVersion {
  my $sub_name  = "seq_StripVersion()";
  my $nargs_expected = 1;
  if(scalar(@_) != $nargs_expected) { printf STDERR ("ERROR, $sub_name entered with %d != %d input arguments.\n", scalar(@_), $nargs_expected); exit(1); } 
  
  my ($accver_R) = (@_);

  $$accver_R =~ s/\.[0-9]*$//; # strip version

  return;
}

####################################################################
# the next line is critical, a perl module must return a true value
return 1;
####################################################################