################################################################################
#
#  kPerl Sequence Laboratory
#  Object Library for PWM-Coded Sequence Motif, Interface to PWMsuite/PWMscan*
#
#  copyright (c)
#    Fritz Lipmann Institute Jena, Genome Analysis Group, 2010
#    Institute of Molecular Biotechnology Jena, Dept. Genome Analysis, 2004
#  author
#    Karol Szafranski, karol.szafranski@leibniz-fli.de
#
################################################################################
#
#  DESCRIPTION
#
# - purpose:
#   this class works as an interface for the positional weight matrix suite
#   "PWMsuite".
#
# - dependencies:
#   - The search engine binaries, written in C, have to be installed
#     in $CorePath{call}{PWMscanNt} and $CorePath{call}{PWMscanAa}
#   - The motif definition format is described in the help output of the
#     motif search engines.
#   - The inherited class "DataFileSync" manages the synchronization between
#     definition file and the internal data structure representation of the
#     motif definition.
#
# - individual description of functions can be found at the beginning of the
#   code blocks
#
################################################################################
#
#  OBJECT OPERATORS
#
#   *** NONE ***
#
#
#  OBJECT METHODS  for external access
#
# - housekeeping
#   new           create object, see MainLib::DefaultObjHash.pm
#   ini           initialize object with motif definition from file
#   AddSwitch     modify object options
#   Load          load (overwrite) motif definition from file
#                 Better use method Sync() than this function to trigger
#                 economic synchronization to file.
#   Save          save (overwrite) motif definition to file
#                 Better use method Sync() than this function to trigger
#                 economic synchronization to file.
#   Sync          test for need of synchronization action, possibly
#                 flush data between poles of synchronization.
#                 See class "MainLib::DataFileSync" for details
#   SyncConnect   --"--  auto-set on successful initialization
#   SyncMaster    --"--  auto-set according to the data type of the
#                        initialization argument
#   SyncValid     --"--
#   Clone         return copy of object
#
# - data management
#   Def           read motif definition => return copy of motif definition
#                 data structure
#   DefFile       read/write name of definition file
#                 alias to MainLib::DataFileSync::SyncFile()
#   DefType       return type of motif definition (last name of package)
#   DefSeqType    read sequence type of motif definition
#   ID            read/write motif identifier
#   Thresh        read/write scoring threshold for motif match search
#   Valid         check validity of motif definition
#                 *** fake implementation ***
#   Width         read/write matrix width
#   RevCompl      return object for reverse-complement matrix
#                 *** not implemented ***
#
# - functionality
#   Search        return hits for search on sequence using PWM
#
#
#  OBJECT DATA STRUCTURE
#  (hash)
#
#   def           reference on motif definition data
#                 On initialization, a hash data structure has to be provided,
#                 which contains at least the fields:
#                   id, matrix, thresh.
#                 Additional fields are optional:
#                   offset (default: 1), orient (+1), ...[more?]
#                 substructure of motif definition data:
#     alph          sequence alphabet covered by the matrix definition
#     AlphType      sequence type: nucleotide/protein
#     backgr        background distribution of sequence symbols
#     descr         matrix description
#     id            motif identifier
#     matrix        reference on matrix array, access hierarchy:
#                     1st dimension (array): sequence position
#                     2nd dimension (hash): symbols
#     offset        motif position which corresponds to matrix position 1
#                   NOTE: matrix position counting starts at 1
#                   for all calulations, value "1" is used
#     orient        motif orientation: +1 / -1
#                   always implicitly +1 in file format
#     thresh        scoring threshold
#   width         matrix width
#   SeqHits       count of hits in scanned sequences
#   SeqScanned    count of scanned sequences
#   switch        hash reference for object switches, cf. method AddSwitch()
#     -debug        print debug protocol to STDERR
#     -SyncOnce     lazy synchronization
#     -SyncOnDestr  on destruction, automatically save definition changes back
#                   to file.
#     -TmpPreserve  do not unlink temporary files, propagate option -preserve
#                   to global manager of temporary files
#   Sync*         see class "MainLib::DataFileSync" for details
#
################################################################################
#
#  FUNCTIONS, DATA
#
#   @ISA
#   %_LibParam
#
# - housekeeping
#   $_LibParam{TmpManag}
#   &new  see MainLib::DefaultObjHash.pm
#   &ini
#   &AddSwitch
#   &_LocalSwitch  see MainLib::DefaultObjHash.pm
#   $_LibParam{*}  file format resources
#   &Load
#   &Save
#   &Clone
#   &DESTROY
#
# - data management
#   &Def
#   &DefFile
#   &DefSeqType
#   &DefType
#   &ID
#   &Width
#   &Valid
#   &RevCompl  *** not implemented ***
#
# - search
#   motif search result data struture  see MotifLib.pm
#   &Search
#   &SearchStrand
#
#
#  STD OPTIONS
#
#  -debug      print debug protocol to STDERR
#
################################################################################
#
#  DEBUG, CHANGES, ADDITIONS
#
# - look also for notes in the header of each function block
#
################################################################################

package SeqLab::MotifPWMscan;

# includes
use strict; #use warnings;  # OK 20040820
use FileHandle;
use MainLib::DefaultObjHash;
use MainLib::Data qw(&DataClone &DataPrint);
use MainLib::Path;
use MainLib::File;
use MainLib::FileTmp;
use MainLib::DataFileSync;
use MainLib::Misc qw(&MySub);
use SeqLab::SeqBench qw (%SeqSmb &SeqType &SeqStrPure);

# inheritance
our @ISA;
push @ISA, qw(MainLib::DefaultObjHash MainLib::DataFileSync);

# package-wide constants and variables
my %_LibParam;


################################################################################
# housekeeping
################################################################################


# handle temporary files globally
# - encapsulation in a sub{} allows MainLib::FileTmp to identify the calling
#   package
my $pcFT = sub{ $_LibParam{TmpManag} = MainLib::FileTmp->new(); };
&$pcFT;


# parametric initialization
#
# INTERFACE
# - argument 1: source argument:
#               - path of sequence input file
#               - reference to motif definition data structure
#                 see explanation on top passage about $this->{def}
#               - reference to an object of __PACKAGE__
#                 This will work like $this->Clone()
#
# - options:    handled by $this->AddSwitch()
#
# - return val: - object reference
#               - undef if an error occured
#
sub ini {
  my ($this, $arg, %opt) = @_;

  # initialize object, keep it blessed
  %$this = ();
  # enter object switches
  $this->AddSwitch(%opt);
  my $debug = $this->{switch}{-debug};

  if (0) { }

  # initialize object with definition file
  elsif (! ref($arg)) {
    if (! -r $arg) {
      die sprintf "%s. ERROR: unable to read motif source file %s\n", &MySub, $arg||"''";
    }
    $this->{SyncFile} = $arg;
    $this->{SyncMaster} = 'file';
    $this->SyncConnect($this->SyncValid());
  }

  # initialize object with definition data structure
  elsif (ref($arg) eq 'HASH') {
    $this->Def($arg);
    $this->{SyncMaster} = 'memory';
  }

  # initialize object from existing object
  # WARNING: there are two conflicting sync bindings now!
  elsif (ref($arg) eq __PACKAGE__) {
    $debug and printf STDERR "%s. initialising with existing object (ref: %s, ID: %s)\n", &MySub,
      ref($arg)||"''", $arg->ID()||"''";
    return $arg->Clone();
  }

  # initialization error
  else {
    printf STDERR "%s. ERROR: unknown initialization argument\n", &MySub;
    return undef;
  }

  # return
  return $this;
}


# enter object switches
#
# INTERFACE
# - argument 1*: hash of switches
# - return val:  success status (boolean)
#
# DESCRIPTION
# - for description of object switches see top
#
sub AddSwitch {
  my ($this,%oopt) = @_;
  my $bErr=0;

  # loop over switches
  while (my($key,$val) = each(%oopt)) {
    if (0) {}

    # delegate -debug to "MainLib::FileTmp"
    elsif ($key eq '-debug' and ($val||0)>1) {
      $bErr ||= ! $_LibParam{TmpManag}->AddSwitch($key=>$val-1,-preserve=>1);
    }
    # delegate -TmpPreserve to "MainLib::FileTmp"
    elsif ($key eq '-TmpPreserve') {
      $bErr ||= ! $_LibParam{TmpManag}->AddSwitch(-preserve=>$val);
    }

   # from here on, each case block exclusively defines the action of object
   # re-shaping associated with that switch
    if (0) {}

    # delegate these to "MainLib::DataFileSync"
    elsif ($key =~ m/^-Sync/) {
      $bErr ||= ! &MainLib::DataFileSync::AddSwitch ($this, $key, $val);
    }

    #options that we just have to enter
    else {
      if (defined $val) { $this->{switch}{$key} = $val; }
      else       { delete $this->{switch}{$key}; }
    }
  }

  # return success status
  return !$bErr;
}


# translation of PWM file header fields to object fields
# - the "default" field specifies a value that is written by method
#   $this->Save() if no certain value was specified.
#
$_LibParam{HeadLabel} = [
  { file=>'ALPH',   code=>'AlphType', default=>'nucleotide' },
  { file=>'BACKGR', code=>'backgr',   default=>undef },
  { file=>'DESCR',  code=>'descr',    default=>undef },
  { file=>'NAME',   code=>'id',       default=>undef },
  { file=>'THRESH', code=>'thresh',   default=>undef },
  { file=>'WIDTH',  code=>'width',    default=>undef },
  ];

# PWMsuite's syntax of sequence type
#
$_LibParam{SeqType2PWMscan} = {
  nucleotide => 'nt',
  protein    => 'aa',
  };
$_LibParam{SeqTypeFromPWMscan} = {
  protein => 'protein',
  aa      => 'protein',
  nucleotide => 'nucleotide',
  nt         => 'nucleotide',
  ACGT       => 'nucleotide',
  };


# read motif definition file
#
# INTERFACE
# - options:
#   -debug      [STD]
#   -file       source file, default $this->{SyncFile}
#
# - return val: success status (boolean)
#
# DESCRIPTION
# - Note that this method is free of arguments!
# - see PWMsuite programmes for specification of the matrix definition format.
#
# DEVELOPER'S NOTES
# - Loading the motif file is necessary for: validity check, preparation of
#   reverse-complement matrix (essential for reverse-strand search), etc.
# - Beware of calling $this->AnyAttribute() methods! They will probably result
#   in an endless recursion of this function via $this->Sync().
#
sub Load {
  my ($this,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};

  # open source file
  my $ArgIn = $lopt{-file} || $this->{SyncFile};
  my $hIn;
  if (ref($ArgIn) =~ m/\b(FileHandle|GLOB)/) {
    $hIn = $ArgIn;
  } else {
    $hIn = FileHandle->new($ArgIn) or return 0;
  }

  # read definition header
  # - translate file format field descriptors to internal keys
  my $line;
  $this->{def} = {};
  while (defined ($line=<$hIn>)) {
    if ($line =~ m/^#/) { next }
    if ($line =~ m/^PWM\b/) { last }
    chomp ($line);
    if ($line =~ m/^(\w+)\s+(.+)/) {
      my ($pLbl) = grep{ $_->{file} eq uc($1) } @{$_LibParam{HeadLabel}};
      unless($pLbl){ next }
      $pLbl->{code} or next;
      if ($pLbl->{file} eq 'ALPH') {
        $this->{def}{$pLbl->{code}} = $_LibParam{SeqTypeFromPWMscan}{$2};
        if ($this->{def}{$pLbl->{code}} !~ m/^(nucleotide|protein)$/) {
          printf STDERR "%s. WARNING: cannot handle alphabet %s\n  %s", &MySub,
            $2||"''";
          return 0;
        }
      } else {
        $this->{def}{$pLbl->{code}} = $2;
      }
    }
  }
  if ($debug) {
    printf STDERR "%s. motif definition header parsed:\n", &MySub;
    print  STDERR map { "  $_->{file} -> $_->{code}: $this->{def}{$_->{code}}\n" }
      @{$_LibParam{HeadLabel}};
  }
  $this->{def}{id} ||= ($ArgIn=~m/^(.+\/)?([^ \t.]+)/)[1];
  $this->{def}{orient} ||= +1;
  $this->{def}{thresh} ||= 0.0;
  unless ($this->{def}{width}) {
    printf STDERR "%s %s. ERROR: NULL-length motif definition in file %s\n", &MySub,
      $this->{def}{id}||"''", $ArgIn||"''";
    return 0;
  }

  # data lines
  # - field "width" is ensured to be defined
  my $pMtf = [];
  for (my $i=0; $i<$this->{def}{width}; ++$i) {

    # read next line
    unless (defined ($line=<$hIn>)) {
      printf STDERR "%s %s. ERROR: premature end of motif definition file, "
        . "matrix length %d/%d of %d\n", &MySub, $this->{def}{id}||"''",
        $i, int(@$pMtf), $this->{def}{width};
      return 0;
    }
    if ($line =~ m/^#/) { redo }
    chomp ($line);

    # split line into array of positional weights
    # definition: line = matrix column = single seq position
    $$pMtf[$i] = [ grep{ length($_) } split(/\s+/,$line) ];
  }
  $debug and printf STDERR "%s. read motif definition, file %s, width %d\n", &MySub,
    $ArgIn||"''", int(@$pMtf);

  # determine sequence type and reformat motif defintion data structure
  my @smb;
  unless ($this->{def}{AlphType}) {
    $debug and printf STDERR "%s %s. assigning default seq. alphabet, file %s\n", &MySub,
      $this->{def}{id}||"''", $ArgIn||"''";
    if (@{$$pMtf[0]} > 5) {
      $this->{def}{AlphType} = 'protein';
      @smb = sort split (//, $SeqSmb{prot});
    } else {
      $this->{def}{AlphType} = 'nucleotide';
      @smb = sort split (//, $SeqSmb{DNA});
    }
  }
  if ($this->{def}{AlphType} eq 'protein') {
    @smb = sort split (//, $SeqSmb{prot});
  } else {
    @smb = sort split (//, $SeqSmb{DNA});
  }
  $debug and printf STDERR "%s. motif %s: type %s, width %d\n",
    &MySub, $this->{def}{id}, $this->{def}{AlphType}, $this->{def}{width};
  $this->{def}{alph} = [ @smb ];
  $debug and printf STDERR "%s %s. sequence symbols: %s\n",
    &MySub, $this->{def}{id}||"''", join(' ',@smb);
  for (my $i=0; $i<$this->{def}{width}; ++$i) {
    @{$this->{def}{matrix}[$i]}{@smb} = @{$$pMtf[$i]};
  }

  # return data
  return $this->Def();
}


# save object's definition to file
#
# INTERFACE
# - options:
#   -debug      [STD]
#   -file       target file, default $this->{SyncFile}
#
# - return val: success status (boolean)
#
sub Save {
  my ($this,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $ArgOut = $lopt{-file} || $this->{SyncFile};
  $debug and printf STDERR "%s %s. saving motif definition to file %s\n", &MySub,
    $this->ID()||"''", $ArgOut;

  # anything to do?
  unless ($ArgOut) { return undef }
  unless ($this->Valid()) {
    printf STDERR "%s. ERROR: cannot save invalid motif definition %s\n", &MySub,
      $this->ID()||"''";
    return 0;
  }

  # open file, write definition header
  my $hOutMtf;
  if (ref($ArgOut) =~ m/^(FileHandle|GLOB)$/) {
    $hOutMtf = $ArgOut;
  } else {
    $hOutMtf = FileHandle->new($ArgOut,'w') or return 0;
  }
  printf $hOutMtf "PWMsuite v2.0\n";
  # write header fields which have a value defined or which have a default value
  foreach my $pLbl (grep{ (defined($_->{code}) and length($this->{def}{$_->{code}})) or $_->{default} }
    @{$_LibParam{HeadLabel}}
  ) {
    if ($pLbl->{file} eq 'ALPH') {
      printf $hOutMtf "%-8s%s\n", $$pLbl{file}, $_LibParam{SeqType2PWMscan}{$this->{def}{$$pLbl{code}}};
    } else {
      printf $hOutMtf "%-8s%s\n", $$pLbl{file},
        $$pLbl{code} ? $this->{def}{$$pLbl{code}}:$$pLbl{default};
    }
  }
  printf $hOutMtf "PWM\n";

  # data lines
  printf $hOutMtf "#\n#%s\n", join('',map { sprintf('%6s',$_) } @{$this->{def}{alph}});
  for (my $i=0; $i<$this->{def}{width}; ++$i) {
    print  $hOutMtf ' ', map { sprintf('%6.2f',$_) }
      @{$this->{def}{matrix}[$i]}{@{$this->{def}{alph}}};
    print  $hOutMtf "\n";
  }

  # finish matrix definition
  print  $hOutMtf "//\n";
  if ($hOutMtf != $ArgOut) { close($hOutMtf) }

  # return successfully
  return 1;
}


# return copy of object
#
# INTERFACE
# - return val: object reference
#
# DEBUG, CHANGES, ADDITIONS
# - this function possibly does not make sense since synchronization is still
#   done to the same file as was done by the template object.
#
sub Clone {
  my ($this) = @_;

  # prepare copy
  my $pCopy = &DataClone({%$this}) or return undef;
  bless $pCopy;

  # refine copy
  delete $pCopy->{TmpPath};

  return $pCopy;
}


# destroy object
#
sub DESTROY {
  my ($this) = @_;
  &MainLib::DataFileSync::DESTROY ($this);
}


################################################################################
# data management
################################################################################


# enter motif definition via data structure
#
# INTERFACE
# - argument 1*: reference to motif definition data structure (write mode)
#
# - options beside object options:
#   -add         allow partial assignments, add argument values that are
#                defined, default: replace complete definition
#
# - return val:  reference to motif definition data structure
#                The motif definition data structure is returned as a copy.
#                Use re-assignment via Def() to modify the object's data.
#
sub Def {
  my ($this,$pMtf,%opt) = @_;

  # synchronize object with file
  # synchronization does not need to be active
  $this->SyncValid() and $this->Sync();

  # change motif definition
  if ($pMtf) {
    if ($opt{-add}) {
      foreach my $ItField (%$pMtf) {
        $this->{def}{$ItField} = $pMtf->{$ItField};
      }
    } else {
      $this->{def} = &DataClone ($pMtf);
    }
    $this->{def}{orient} ||= +1;

    # check definition
    $this->Valid() or printf STDERR "%s. ERROR: in motif definition assignment\n", &MySub;
    $this->SyncMemTouch(-warn=>1);
  }

  return &DataClone ($this->{def});
}


# read/write name of motif defintion file
#
sub DefFile { return SyncFile(@_); }


# read defintion type
#
sub DefType { return (split('::',__PACKAGE__))[-1] }


# sequence type of motif definition
#
# INTERFACE
# - return val: sequence type
#
sub DefSeqType {
  my ($this) = @_;

  # synchronize object with file
  # synchronization does not need to be active
  $this->SyncValid() and $this->Sync();

  # determine sequence type
  # this is needed if the object data structure is the master representation
  #   and this data structure is being gradually constructed
  unless ($this->{def}{AlphType}) {
    $this->{def}{AlphType} = (keys(%{$this->{def}{matrix}[0]}) > 5) ? 'protein':'nucleotide';
  }

  return $this->{def}{AlphType};
}


# defintion identifier
#
# INTERFACE
# - argument 1*: motif identifier (write mode)
# - return val:  motif identifier
#
sub ID {
  my ($this,$arg) = @_;

  # synchronize object with file
  # synchronization does not need to be active
  $this->SyncValid() and $this->Sync();

  # set motif identifier
  if (defined $arg) {
    $this->{def}{id} = $arg;
    $this->SyncMemTouch(-warn=>1);
  }

  # determine motif identifier from motif definition file
  else {
    $this->{def}{id} ||=
      &PathSplit ($this->SyncFile())->{nameroot} || undef;
  }

  return $this->{def}{id};
}


# score threshold
#
# INTERFACE
# - argument 1*: score threshold (write mode)
# - return val:  score threshold
#
sub Thresh {
  my ($this,$arg) = @_;

  # synchronize object with file
  # synchronization does not need to be active
  $this->SyncValid() and $this->Sync();

  # set name of file
  if (defined $arg) {
    $this->{def}{thresh} = $arg;
    $this->SyncMemTouch(-warn=>1);
  }

  return $this->{def}{thresh};
}


# check validity of motif definition
#
# INTERFACE
# - return val:  validity (boolean)
#
# DESCRIPTION
# - do not confuse this function with SyncValid()
#
# DEBUG, CHANGES, ADDITIONS
# - currently, this is a fake implementation
#
sub Valid {
  my ($this) = @_;
  return 1;
}


# width of motif definition
#
# INTERFACE
# - argument 1*: motif width (write mode)
# - return val:  motif width
#
sub Width {
  my ($this,$arg) = @_;

  # synchronize object with file
  # synchronization does not need to be active
  $this->SyncValid() and $this->Sync();

  # set matrix width
  if (defined $arg) {
    $this->{def}{width} = $arg;
    $#{$this->{def}{matrix}} = $arg - 1;
    $this->SyncMemTouch(-warn=>1);
  }

  # default value
  $this->{def}{width} ||= int (@{$this->{def}{matrix}});

  return $this->{def}{width};
}


# reverse-complement PWM
#
# INTERFACE
# - return val: - object reference
#               - undef if an error occurred
#
# DESCRIPTION
# - For motif search, this function is not necessary, since "PWMscanNt" supports
#   reverse-complement search mode.
#
# DEBUG, CHANGES, ADDITIONS
# - implement me! This might be achieved by calling the PWM suite conversion
#   program "PWMconvert".
#
sub RevCompl {
  my ($this,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  if ($this->DefSeqType() ne 'nucleotide') { return undef }
  # synchronization done

  # reverse-complement motif definition
  my $pMtfRev = $this->Clone();
  # ...

  return $this->{RevCompl} = $pMtfRev;
}



################################################################################
# search
################################################################################


# identify and locate instances of search pattern
#
# INTERFACE
# - argument 1: query sequence string
#
# - options:
#   -best       return n best scoring hit sites only
#   -debug      [STD]
#   -HitSurrd   lengths of stored surrounding sequences, default: 0
#               This takes effect in &SearchStrand
#   -isPure     no need to purify sequence string
#   -strands    which strand to search
#               0   both strands (default)
#               -1  minus strand only
#               1   plus strand only
#
# - return val: - reference to result data structure
#               - undef if an error occurred
#
# DESCRIPTION
# - the query sequence string will be purified prior to search
#
sub Search {
  my ($this,$sSeq,%opt) = @_;

  # function parameters
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  $opt{-strands} ||= 0;

  # process / test search sequence
  unless ($lopt{isPure}) {
    $sSeq = &SeqStrPure ($sSeq, -upper=>1);
  }
  unless ($sSeq) { return undef }

  # sync with motif defintion file
  unless ($this->SyncConnect() and $this->Sync()) {
    die sprintf "%s %s. ERROR: missing sync connection\n", &MySub, $this->ID();
  }

  # test for fitting sequence alphabets
  if (&SeqType($sSeq,-basic=>1) ne $this->DefSeqType()) {
    printf STDERR "%s %s. ERROR: seq type %s doesn't fit to motif seq type\n",
      &MySub, $this->ID(), &SeqType($sSeq,-type=>'basic');
    return undef;
  }

  # locate positively oriented instances
  my @MtfHit;
  if ($opt{-strands} >= 0) {
    push @MtfHit, @{ $this->SearchStrand($sSeq,%opt)||[] };
    $debug and printf STDERR "%s. motif %s, seq #%d, strand +1: %d hits\n", &MySub,
      $this->ID(), $this->{SeqScanned}, int @MtfHit;
  }
  # locate negatively oriented instances
  if ($opt{-strands} <= 0) {
    push @MtfHit, @{ $this->SearchStrand($sSeq,-strand=>-1,%opt)||[] };
    $debug and printf STDERR "%s. motif %s, seq #%d, both strands: %d hits\n", &MySub,
      $this->ID(), $this->{SeqScanned}, int @MtfHit;
  }

  # select best scoring hits
  if ($lopt{-best}) {
    # *** implement me ***
  }

  # sort complete list of matches, return result
  $this->{SeqScanned} ++;
  $this->{SeqHits} += int @MtfHit;
  return [ sort{ $a->{offset}<=>$b->{offset} } @MtfHit ];
}


# identify and locate instances of search pattern in one strand
#
# INTERFACE
# - argument 1: sequence
#
# - options:
#   -debug      [STD]
#   -HitSurrd   lengths of stored surrounding sequences, default: 0
#   -strand     which sequence strand to search, default: +1
#
# - return val: - reference to instance data structure
#               - undef if an error occurred
#
# DESCRIPTION
# - match offset value is stated in biological counting logics - always "1"
#   here, though the motif definition file format allows other values.
# ? match offset value refers to the beginning of the match in the
#   sense-directed template sequence. The offset of the motif may be at
#   the end of the match if the motif matches as an inverted instance.
#
sub SearchStrand {
  my ($this,$sSeq,%opt) = @_;
  my %lopt = $this->_LocalSwitch(%opt);
  my $debug = $lopt{-debug};
  my $strand = $lopt{-strand} || 1;
  my $HitSurrd = $lopt{-HitSurrd};

  # ensure sync with motif defintion file
  # assuming memory->file sync done in Search()
  if ($strand==-1 and $this->SyncMaster() eq 'file') {
    unless ($this->SyncConnect()) {
      die sprintf "%s %s. ERROR: missing sync connection\n", &MySub, $this->ID();
    }
    $this->Sync() or return undef;
  }

  # prepare sequence file, reserve error log, prepare call string, start call
  my $PathSeq = $_LibParam{TmpManag}->Create();
  &WriteFile ($PathSeq, sprintf(">%s\n%s\n",$this->DefType().'-arg',uc($sSeq)) );
  my $PathErr = $_LibParam{TmpManag}->Create();
  my $CallLbl = ($this->DefSeqType() eq 'protein') ? 'PWMscanAa' : 'PWMscanNt';
  my $CallProg = $CorePath{call}{$CallLbl};
  if (! -x $CallProg) {
    printf STDERR "%s. missing executable \"%s\" at %s\n", &MySub,
      $CallLbl, $CallProg;
    exit 1;
  }
  my $call = join (' ', $CallProg, $HitSurrd?('-f',$HitSurrd):(),
    ($strand==-1)?('-r'):(), $this->{SyncFile}, "< $PathSeq");
  my ($CallIn);
  if (not $CallIn = FileHandle->new("$call 2>$PathErr |")) {
    printf STDERR "%s. ERROR: motif search call failed, call:\n  %s\n", &MySub, $call||"''";
    printf STDERR "  error log: %s\n%s",
      -s($PathErr)?'':'NONE', -s($PathErr)?scalar &ReadFile($PathErr):'';
    return [];
  }

  # parse search result on the fly
  my (@hit);
  while (defined (my $line=<$CallIn>)) {
    my %HitInst = (
      MotifID => $this->ID(),
      length  => $this->Width(),
      orient  => $strand,
    );
    @HitInst{qw(NULL NULL offset score ScoreExp instance)} = split(/\s+/,$line);
    # flanks were reported by search engine, central match is upper case
    if ($HitSurrd and ($HitInst{instance}||'') =~ m/[A-Z]+/) {
      @HitInst{qw(ante post)} = ($`, $');
      $HitInst{instance} = $&;
    }
    push @hit, { %HitInst };
  }
  if (-s $PathErr) {
    printf STDERR "%s. ERROR in search function, system call:\n  %s\n", &MySub, $call||"''";
    print  STDERR "  error log:\n", &ReadFile($PathErr);
  } else {
    $debug and printf STDERR "%s. seq #%d, strand %d: %d result lines from search machine\n", &MySub,
      $this->{SeqScanned}, $strand, $.;
  }

  # tidy up, return hits
  if (!$debug and !$lopt{-TmpPreserve}) { unlink $PathSeq, $PathErr; }
  return \@hit;
}


1;
# $Id: MotifPWMscan.pm,v 1.22 2018/06/05 18:02:56 szafrans Exp $
