<commit_before>
=head1 NAME

Bio::SecreTary::TMpred 

=head1 DESCRIPTION

An object to run the trans-membrane helix prediction program tmpred.

=head1 AUTHOR

<NAME> (<EMAIL>)

=cut

package Bio::SecreTary::TMpred;

use strict;
use warnings;

use IO::File;
use File::Temp;
use CGI ();
use Bio::SecreTary::Table;
use Bio::SecreTary::Helix;

use Readonly;
Readonly my $FALSE => 0;
Readonly my $TRUE  => 1;
Readonly my $PASCAL_STYLE =>
  $TRUE;    # if this is true, does the same as the old pascal code.
Readonly my $TMHLOFFSET => ($PASCAL_STYLE)
  ? 1
  : 0;      
# TMHLOFFSET gets added to the max_tmh_length,
# and (1 - TMHLOFFSET) gets subtracted from min_tmh_length
# TMHLOFFSET => 1  makes it agree with pascal code.
# finds helices with length as large as next bigger odd number, e.g. if
# $max_tmh_length is 33, will find some helices of lengths 34 and 35.
# with TMHLOFFSET => 0, 33->33, 32->33, 31->31, etc. i.e. goes up to
# next greater OR EQUAL odd number, rather than to next STRICTLY greater odd.
# similarly the min length is affected. with TMHLOFFSET => 1, 17->17, 16->17, 15->15
# i.e. if you specify min length of 16 you will never see helices shorter than 17
# but with TMHLOFFSET => 0, 17->17, 16->15, 15->15, etc. now you find the length 16 ones, (as well as length 15 ones which are discarded in good_solutions).



=head2 function new()

Synopsis:
my $tmpred_obj = Bio::SecreTary::TMpred->new(); # using defaults
    or
    my $tmpred_obj = Bio::SecreTary::TMpred->new({min_score => 600});
	Arguments: $arg_hash_ref holds some parameters describing which solutions will be found by tmpred and kept by good_solutions: min_score, min_tm_length, max_tm_length, min_beg, max_beg.
	Returns:	an instance of a TMpred object
	Side effects:	Creates the object.
  Description: Creates a TMpred object with certain parameters which determine which trans-membrane helices to keep.

=cut

sub new {
    my $class        = shift;
    my $self         = bless {}, $class;
    my $arg_hash_ref = shift;
    $self->{version} = (exists $arg_hash_ref->{version})? $arg_hash_ref->{version} : 'perl';
    $self->{min_score} =
      ( exists $arg_hash_ref->{min_score} ) ? $arg_hash_ref->{min_score} : 500;
    $self->{min_tm_length} =
      ( exists $arg_hash_ref->{min_tm_length} )
      ? $arg_hash_ref->{min_tm_length}
      : 17;
    $self->{max_tm_length} =
      ( exists $arg_hash_ref->{max_tm_length} )
      ? $arg_hash_ref->{max_tm_length}
      : 33;
    $self->{min_beg} =
      ( exists $arg_hash_ref->{min_beg} ) ? $arg_hash_ref->{min_beg} : 0;
    $self->{max_beg} =
      ( exists $arg_hash_ref->{max_beg} ) ? $arg_hash_ref->{max_beg} : 35;

    $self->{min_halfw} =
      int( ( ( $self->{min_tm_length} ) - ( 1 - $TMHLOFFSET ) ) / 2 );
    $self->{max_halfw} =
      int( ( ( $self->{max_tm_length} ) + $TMHLOFFSET ) / 2 );

    $self->setup_tables();
    my @aa_array = split( '', 'ACDEFGHIKLMNPQRSTVWYX' ); # first 20 are std aa's; X means unknown aa
    my %aa_row_hash = ()
      ; # keys: letters representing aa's; values: numbers 0-20 (20 for unknown)
    for ( my $i = 0 ; $i < scalar @aa_array ; $i++ ) {
        $aa_row_hash{ $aa_array[$i] } = $i;
    }
    $self->set_aa_row_hash( \%aa_row_hash );

    return $self;
}

sub run_tmpred {
    my $self     = shift;
    my ($sequence_id, $sequence, $version, $do_long_output) =
	('Anon_prot_seq', undef, 'perl', $FALSE);
    my $arg_hash_ref = undef;
    my $arg1 = shift; # can be either sequence, or ref to  hash holding arguments
    if(ref($arg1)){
	$arg_hash_ref = $arg1;
    }
    else{ # arg1 is not a ref; should be sequence
	$sequence = $arg1; 
	my $arg2 = shift;
	if(ref($arg2)){
	    $arg_hash_ref = $arg2;
	}else{ # arg2 is not a ref; should be sequence_id
	    $sequence_id = (defined($arg2))? $arg2: 'A_prot_seq';
	}
    }
    if(defined($arg_hash_ref)){
if(ref($arg_hash_ref) eq 'HASH'){
	$sequence = $arg_hash_ref->{sequence} if(exists $arg_hash_ref->{sequence});
	$sequence_id = $arg_hash_ref->{sequence_id} if(exists $arg_hash_ref->{sequence_id});
	$version = $arg_hash_ref->{version} if(exists $arg_hash_ref->{version});
	$do_long_output = $arg_hash_ref->{do_long_output} if(exists $arg_hash_ref->{do_long_output});
    }
    else{
	die "Argument to run_tmpred is ref but not hash ref.\n";
}
    }
    # check and fix up the sequence.
    $sequence =~ s/\s+//g;    # remove whitespace
    $sequence =~ s/\*$//;     # remove * at end if present
    $sequence = uc $sequence;

    #  die  "sequence contains non A-Z characters:[$sequence]\n"
    if ( $sequence =~ /[^A-Z]/ ) {
        warn
          "substituting 'X' for non A-Z character. \n$sequence_id  $sequence\n";
        $sequence =~ s/[^A-Z]/X/g;
    }

    $sequence =~ s/[BJOUZ]/X/g;    # to get just 20 aa's plus X for unknown.
    if ($PASCAL_STYLE) {
	# the pascal version of tmpred simply substitutes A for X.
        $sequence =~ s/X/A/g;
    }
    my $good_solns;
    my $tmpred_raw_out = undef;
    if ( $version eq 'pascal' ) {
        $tmpred_raw_out =
          $self->run_tmpred_pascal( $sequence_id, $sequence );

        # process tmpred output to get summary,
        # i.e. score, begin and end positions for tmh's satisfying limits.
        $good_solns =
          ( $self->good_solutions_pascal($tmpred_raw_out) );    # , $tmpred_raw_out);
    }
    else {                                               # use perl code
        $good_solns = $self->run_tmpred_perl( $sequence_id, $sequence );
	if($do_long_output){
	    $tmpred_raw_out = 'Long output not implemented in perl version.';
	}
    }
  #  return $good_solns;
    return ($good_solns, $tmpred_raw_out);
}


sub run_tmpred_pascal {                     # this uses pascal code.
    my $self        = shift;
    my $sequence_id = shift;
    my $sequence    = shift;

    my $temp_file_dir    = '/home/tomfy/tempfiles';
    my $temp_file_handle = File::Temp->new(
        TEMPLATE => 'tmpred_input_XXXXXX',
        DIR      => $temp_file_dir,
        UNLINK   => 0
    );

    my $temp_file = $temp_file_handle->filename;
    print $temp_file_handle ">$sequence_id\n$sequence\n";
    $temp_file_handle->close();

    my $max_tm_length = $self->{max_tm_length};
    my $min_tm_length = $self->{min_tm_length};
    my $tmpred_dir    = "/home/tomfy/tmpred";
    my $tmpred_pascal_out =
`$tmpred_dir/tmpred  -def -in=$temp_file  -out=-  -par=$tmpred_dir/matrix.tab -max=$max_tm_length  -min=$min_tm_length`;

    return $tmpred_pascal_out;
}


=head2 function good_solutions_pascal()
	Synopsis: my $solution_string = $tmpred->good_solutions_pascal(); 
Arguments: none.
Returns: String of "good" solutions' (score,beg,end) info.
Description: Parses tmpred output; finds the tmhs which
are consistent with limits, and store them in a string,
    e.g.: "(2000,17,34)  (1500,12,36)".

=cut

sub good_solutions_pascal {
    my $self       = shift;
    my $tmpred_out = shift;

    my $solutions = "";
    my $ok        = 0;
    while ( $tmpred_out =~ /(.*?\n)/ ) {
        my $line = $1;
        $tmpred_out = substr( $tmpred_out, length $line );    #
        last  if ( $line =~ /Table of correspondences/ );
        $ok++ if ( $line =~ /Possible transmembrane helices/ );
        $ok++ if ( $line =~ /Inside to outside helices/ );
        if ( $ok == 2 ) {

            if ( $line =~
                /(\d+)\s*\(\s*(\d+)\)\s*(\d+)\s*\(\s*(\d+)\)\s*(\d+)\s*(\d+)/ )
            {
                my ( $begin, $begcore, $end, $endcore, $score, $center ) =
                  ( $1, $2, $3, $4, $5, $6 );

                my $length = $end + 1 - $begin;

                if (   $length > $self->{max_tm_length}
                    or $length < $self->{min_tm_length} )
                {
                    warn "length out of allowed range: $length \n";
                }
                if (    $score >= $self->{min_score}
                    and $begin <= $self->{max_beg}
                    and $begin >= $self->{min_beg}
                    and $length <= $self->{max_tm_length}
                    and $length >= $self->{min_tm_length} )
                {
                    $solutions .= "($score,$begin,$end)  ";
                }
		else{
		    warn "Solution rejected in good_solutions_pascal.\n($score,$begin,$length)";
		}
            }
        }
    }

    if ( $solutions eq "" ) { $solutions = "(-10000,0,0)"; }
    else { $solutions =~ s/(\s+)$//; }    # eliminate whitespace at end.
    return $solutions;
}

sub run_tmpred_perl {
    my $self        = shift;
    my $sequence_id = shift;
    my $sequence    = shift;

    my $aa_row_hash             = $self->get_aa_row_hash();
    my @seq_array               = split( "", $sequence );
    my @sequence_aanumber_array = ();
    foreach (@seq_array) {
        if ( exists $aa_row_hash->{$_} ) {
            push @sequence_aanumber_array, $aa_row_hash->{$_};
        }
        else {
            die "no row in aa_row_hash for aa: [$_] \n";
        }
    }

    my $length = scalar @sequence_aanumber_array;

    # get the 6 profiles:
    # in-to-out
    my $io_center_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_center_table() );
    my $io_nterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_nterm_table() );
    my $io_cterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_cterm_table() );

    # out-to-in
    my $oi_center_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_center_table() );
    my $oi_nterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_nterm_table() );
    my $oi_cterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_cterm_table() );

    ( $length == scalar @$io_center_prof )
      || die "seq length: ", $length, "  io_center_profile size: ",
      scalar @$io_center_prof, "\n";
    if (0) {
        for ( my $i = 0 ; $i < scalar @$io_center_prof ; $i++ ) {
            print $sequence_aanumber_array[$i], "  ",
              $io_center_prof->[$i], "  ",
              $io_nterm_prof->[$i],  "  ",
              $io_cterm_prof->[$i],  "  ",
              $oi_center_prof->[$i], "  ",
              $oi_nterm_prof->[$i],  "  ",
              $oi_cterm_prof->[$i],  "\n";
        }
    }

    # get score array for in-to-out using 3 profiles (center, nterm, cterm)
    my $io_score =
      $self->make_curve( $io_center_prof, $io_nterm_prof, $io_cterm_prof );

    # same for out-to-in
    my $oi_score =
      $self->make_curve( $oi_center_prof, $oi_nterm_prof, $oi_cterm_prof );

    my @io_helices = ();
    {
        my $fhresult;
        my $helix;
        my $start = 0;
        while ( $start < $length ) {
            ( $fhresult, $start, $helix ) =
              $self->find_helix( $length, $start, $io_score, $io_center_prof,
                $io_nterm_prof, $io_cterm_prof );
            $helix->set_nt_in('TRUE');    # io is inside-to-outside
            if ($fhresult) {
                push @io_helices, $helix;
            }
            last if ( $start >= $length );
        }
    }

    my @oi_helices = ();
    {
        my $fhresult;
        my $start = 0;
        my $helix;
        while ( $start < $length ) {
            ( $fhresult, $start, $helix ) =
              $self->find_helix( $length, $start, $oi_score, $oi_center_prof,
                $oi_nterm_prof, $oi_cterm_prof );
            $helix->set_nt_in('FALSE');
            if ($fhresult) {
                push @oi_helices, $helix;
            }
            last if ( $start >= $length );
        }
    }

    my @helices = ( @io_helices, @oi_helices );
    return $self->good_solutions_perl( \@helices );

    #	return (\@io_helices, \@oi_helices);
    #   return ( scalar @io_helices, scalar @oi_helices );

}

sub good_solutions_perl {    # get the good solutions starting from the
                             # output of perl tmpred
    my $self           = shift;
    my $tmpred_helices = shift;    # ref to array of predicted helices

    my $good_string = "";
    foreach (@$tmpred_helices) {
        my $beg    = $_->get_nterm()->[0] + 1;
        my $end    = $_->get_cterm()->[0] + 1;
        my $score  = $_->get_score();
        my $length = $end + 1 - $beg;

        if (    $score >= $self->{min_score}
            and $length >= $self->{min_tm_length}
            and $length <= $self->{max_tm_length}
            and $beg >= $self->{min_beg}
            and $beg <= $self->{max_beg} )
        {
            $good_string .= "($score,$beg,$end)  ";
        }
    }
    if ( $good_string eq "" ) { $good_string = "(-10000,0,0)"; }
    else { $good_string =~ s/(\s+)$//; }    # eliminate whitespace at end.
    return $good_string;
}

sub make_profile {                        # makes a profile, i.e. an array
                                          # containing ...
                                          # my $sequence     = shift;
    my $seq_aanumber_array = shift;       # ref to array of numbers
    my $table              = shift;
    my $ref_position = $table->get_marked_position();
    my $matrix       = $table->get_table();
    my $ncols   = scalar @{ $matrix->[0] };   # ncols is # elements in first row
    my @profile = ();
    my $length = scalar @$seq_aanumber_array;    #length $sequence;

# need to be careful here to make it consistent with pascal code,
# which has 1-based arrays. ref_position is 1 less than in pascal code for this reason.
    for ( my $k = 0 ; $k < $length ; $k++ ) {
        my $m    = 0;
        my $kmrf = $k - $ref_position;
        my $plo  = $kmrf;                        # $k - $ref_position;
        $plo = 0 if ( $plo < 0 );
        my $pup = $ncols + $kmrf;                # $k - $ref_position;
        $pup = $length if ( $pup > $length );

        for (
            my ( $p, $i ) = ( $plo, $plo - $kmrf ) ;  # $k - $ref_position ) ) ;
            $p < $pup ; $p++, $i++
          )
        {
            $m += $matrix->[ $seq_aanumber_array->[$p] ]->[$i];
        }

        my $round_m =
          ( $m < 0 ) ? int( $m * 100 - 0.5 ) : int( $m * 100 + 0.5 );
        push @profile, $round_m;<commit_msg>Adding long output like that produced by original pascal code for tmpred<commit_after>
=head1 NAME

Bio::SecreTary::TMpred 

=head1 DESCRIPTION

An object to run the trans-membrane helix prediction program tmpred.

=head1 AUTHOR

<NAME> (<EMAIL>)

=cut

package Bio::SecreTary::TMpred;

use strict;
use warnings;

use Carp;
use IO::File;
use File::Temp;
use CGI ();
use Bio::SecreTary::Table;
use Bio::SecreTary::Helix;

use Readonly;
Readonly my $FALSE    => 0;
Readonly my $TRUE     => 1;
Readonly my %defaults => (
    'version'       => 'perl',
    'min_score'     => 500,
    'min_tm_length' => 17,
    'max_tm_length' => 33,
    'min_beg'       => 0,
    'max_beg'       => 35,
'lo_orientational_threshold' => 80,
  'hi_orientational_threshold' => 200,
  'avg_orientational_threshold' => 80);

Readonly my $PASCAL_STYLE =>
  $TRUE;    # if this is true, does the same as the old pascal code.
Readonly my $TMHLOFFSET => ($PASCAL_STYLE)
  ? 1
  : 0;

# TMHLOFFSET gets added to the max_tmh_length,
# and (1 - TMHLOFFSET) gets subtracted from min_tmh_length
# TMHLOFFSET => 1  makes it agree with pascal code.
# finds helices with length as large as next bigger odd number, e.g. if
# $max_tmh_length is 33, will find some helices of lengths 34 and 35.
# with TMHLOFFSET => 0, 33->33, 32->33, 31->31, etc. i.e. goes up to
# next greater OR EQUAL odd number, rather than to next STRICTLY greater odd.
# similarly the min length is affected. with TMHLOFFSET => 1, 17->17, 16->17, 15->15
# i.e. if you specify min length of 16 you will never see helices shorter than 17
# but with TMHLOFFSET => 0, 17->17, 16->15, 15->15, etc. now you find
# the length 16 ones, (as well as length 15 ones which are discarded in good_solutions).
# set up defaults for tmpred parameters:

=head2 function new

  Synopsis : my $tmpred_obj = Bio::SecreTary::TMpred->new();    # using defaults
  or my $tmpred_obj =
  Bio::SecreTary::TMpred->new( { min_score => 600 } );
  Arguments: $arg_hash_ref holds some parameters describing which 
      solutions will be found by tmpred :
      min_score, min_tm_length, max_tm_length, min_beg, max_beg . 
  Returns: an instance of a TMpred object 
  Side effects: Creates the object . 
  Description: Creates a TMpred object with certain parameters which 
      determine which trans-membrane helices to find.

=cut

sub new {
    my $class        = shift;
    my $self         = bless {}, $class;
    my $arg_hash_ref = shift;

    foreach my $param ( keys %defaults ) {    # set params to defaults
        $self->{$param} = $defaults{$param};
    }

    if ( defined $arg_hash_ref ) {            # otherwise just use defaults
        if ( ref($arg_hash_ref) eq 'HASH' ) {    # check it's really a hash ref
            foreach my $param ( keys %$arg_hash_ref ) {
                if ( exists $self->{$param} ) {
                    if ( defined $arg_hash_ref->{$param} ) {
                        $self->{$param} = $arg_hash_ref->{$param};
                    }
                }
                else {    # a key of the arg hash is an unknown parameter.
                    carp
"In TMpred constructor. Attempt to set unknown parameter $param.\n"
                      . "Known params are:  "
                      . join( " ", keys %$arg_hash_ref ), "\n";
                }
            }
        }
        else {
            carp(
"Argument to TMpred constructor is not hash ref as it should be.\nUsing defaults.\n"
            );
        }
    }
    $self->{min_halfw} =
      int( ( ( $self->{min_tm_length} ) - ( 1 - $TMHLOFFSET ) ) / 2 );
    $self->{max_halfw} =
      int( ( ( $self->{max_tm_length} ) + $TMHLOFFSET ) / 2 );

    $self->setup_tables();
    my @aa_array = split( '', 'ACDEFGHIKLMNPQRSTVWYX' )
      ;    # first 20 are std aa's; X means unknown aa
    my %aa_row_hash = ()
      ; # keys: letters representing aa's; values: numbers 0-20 (20 for unknown)
    for ( my $i = 0 ; $i < scalar @aa_array ; $i++ ) {
        $aa_row_hash{ $aa_array[$i] } = $i;
    }
    $self->set_aa_row_hash( \%aa_row_hash );

    return $self;
}

sub run_tmpred {
    my $self = shift;
    my ( $sequence_id, $sequence, $version, $do_long_output ) =
      ( 'Anon_prot_seq', undef, 'perl', $FALSE );
    my $arg_hash_ref = undef;
    my $arg1 =
      shift;    # can be either sequence, or ref to  hash holding arguments
    if ( ref($arg1) ) {
        $arg_hash_ref = $arg1;
    }
    else {      # arg1 is not a ref; should be sequence
        $sequence = $arg1;
        my $arg2 = shift;
        if ( ref($arg2) ) {
            $arg_hash_ref = $arg2;
        }
        else {    # arg2 is not a ref; should be sequence_id
            $sequence_id = ( defined($arg2) ) ? $arg2 : 'A_prot_seq';
        }

        #	print "id, sequence: $sequence_id, $sequence \n";
    }
    if ( defined($arg_hash_ref) ) {
        if ( ref($arg_hash_ref) eq 'HASH' ) {
            $sequence = $arg_hash_ref->{sequence}
              if ( exists $arg_hash_ref->{sequence} );
            $sequence_id = $arg_hash_ref->{sequence_id}
              if ( exists $arg_hash_ref->{sequence_id} );
            $version = $arg_hash_ref->{version}
              if ( exists $arg_hash_ref->{version} );
            $do_long_output = $arg_hash_ref->{do_long_output}
              if ( exists $arg_hash_ref->{do_long_output} );
        }
        else {
            die "Argument to run_tmpred is ref but not hash ref.\n";
        }
    }

    # check and fix up the sequence.
    $sequence =~ s/\s+//g;    # remove whitespace
    $sequence =~ s/\*$//;     # remove * at end if present
    $sequence = uc $sequence;

    #  die  "sequence contains non A-Z characters:[$sequence]\n"
    if ( $sequence =~ /[^A-Z]/ ) {
        warn
          "substituting 'X' for non A-Z character. \n$sequence_id  $sequence\n";
        $sequence =~ s/[^A-Z]/X/g;
    }

    $sequence =~ s/[BJOUZ]/X/g;    # to get just 20 aa's plus X for unknown.
    if ($PASCAL_STYLE) {

        # the pascal version of tmpred simply substitutes A for X.
        $sequence =~ s/X/A/g;
    }
    my $good_solns;
    my $tmpred_raw_out = undef;
    if ( $version eq 'pascal' ) {
        $tmpred_raw_out = $self->run_tmpred_pascal( $sequence_id, $sequence );

        # process tmpred output to get summary,
        # i.e. score, begin and end positions for tmh's satisfying limits.
        $good_solns = ( $self->good_solutions_pascal($tmpred_raw_out) );
    }
    else {    # use perl code
        my ( $io_helices_ref, $oi_helices_ref ) =
          $self->run_tmpred_perl( $sequence_id, $sequence );
        $good_solns = $self->good_solutions_perl($io_helices_ref);
        $good_solns .= $self->good_solutions_perl($oi_helices_ref);
        if ( 1 or $do_long_output ) {
            print "io/oi helices: ", scalar @$io_helices_ref, "  ",
              scalar @$oi_helices_ref, "\n";
            print $self->short_output(
                $sequence_id, $sequence, $io_helices_ref, $oi_helices_ref
              ),
              "\n";
            $tmpred_raw_out = $self->_long_output($sequence_id, $sequence, $io_helices_ref, $oi_helices_ref);
        }
    }

    #  print "good_solutions: $good_solns \n";
    #  return $good_solns;
    return ( $good_solns, $tmpred_raw_out );
}

sub run_tmpred_pascal {    # this uses pascal code.
    my $self        = shift;
    my $sequence_id = shift;
    my $sequence    = shift;

    my $temp_file_dir    = '/home/tomfy/tempfiles';
    my $temp_file_handle = File::Temp->new(
        TEMPLATE => 'tmpred_input_XXXXXX',
        DIR      => $temp_file_dir,
        UNLINK   => 0
    );

    my $temp_file = $temp_file_handle->filename;
    print $temp_file_handle ">$sequence_id\n$sequence\n";
    $temp_file_handle->close();

    my $max_tm_length = $self->{max_tm_length};
    my $min_tm_length = $self->{min_tm_length};
    my $tmpred_dir    = "/home/tomfy/tmpred";
    my $tmpred_pascal_out =
`$tmpred_dir/tmpred  -def -in=$temp_file  -out=-  -par=$tmpred_dir/matrix.tab -max=$max_tm_length  -min=$min_tm_length`;

    return $tmpred_pascal_out;
}

=head2 function good_solutions_pascal()
	Synopsis: my $solution_string = $tmpred->good_solutions_pascal(); 
Arguments: none.
Returns: String of "good" solutions' (score,beg,end) info.
Description: Parses tmpred output; finds the tmhs which
are consistent with limits, and store them in a string,
    e.g.: "(2000,17,34)  (1500,12,36)".

=cut

sub good_solutions_pascal {
    my $self       = shift;
    my $tmpred_out = shift;

    my $solutions = "";
    my $ok        = 0;
    while ( $tmpred_out =~ /(.*?\n)/ ) {
        my $line = $1;
        $tmpred_out = substr( $tmpred_out, length $line );    #
        last  if ( $line =~ /Table of correspondences/ );
        $ok++ if ( $line =~ /Possible transmembrane helices/ );
        $ok++ if ( $line =~ /Inside to outside helices/ );
        if ( $ok == 2 ) {

            if ( $line =~
                /(\d+)\s*\(\s*(\d+)\)\s*(\d+)\s*\(\s*(\d+)\)\s*(\d+)\s*(\d+)/ )
            {
                my ( $begin, $begcore, $end, $endcore, $score, $center ) =
                  ( $1, $2, $3, $4, $5, $6 );

                my $length = $end + 1 - $begin;

                if (   $length > $self->{max_tm_length}
                    or $length < $self->{min_tm_length} )
                {
                    warn "length out of allowed range: $length \n";
                }
                if (    $score >= $self->{min_score}
                    and $begin <= $self->{max_beg}
                    and $begin >= $self->{min_beg}
                    and $length <= $self->{max_tm_length}
                    and $length >= $self->{min_tm_length} )
                {
                    $solutions .= "($score,$begin,$end)  ";
                }
                else {
                    warn
"Solution rejected in good_solutions_pascal.\n($score,$begin,$length)";
                }
            }
        }
    }

    if ( $solutions eq "" ) { $solutions = "(-10000,0,0)"; }
    else { $solutions =~ s/(\s+)$//; }    # eliminate whitespace at end.
    return $solutions;
}

sub run_tmpred_perl {
    my $self        = shift;
    my $sequence_id = shift;
    my $sequence    = shift;

    my $aa_row_hash             = $self->get_aa_row_hash();
    my @seq_array               = split( "", $sequence );
    my @sequence_aanumber_array = ();
    foreach (@seq_array) {
        if ( exists $aa_row_hash->{$_} ) {
            push @sequence_aanumber_array, $aa_row_hash->{$_};
        }
        else {
            die "no row in aa_row_hash for aa: [$_] \n";
        }
    }

    my $length = scalar @sequence_aanumber_array;

    # get the 6 profiles:
    # in-to-out
    my $io_center_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_center_table() );
    my $io_nterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_nterm_table() );
    my $io_cterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_io_cterm_table() );

    # out-to-in
    my $oi_center_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_center_table() );
    my $oi_nterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_nterm_table() );
    my $oi_cterm_prof =
      make_profile( \@sequence_aanumber_array, $self->get_oi_cterm_table() );

    ( $length == scalar @$io_center_prof )
      || croak "seq length: ", $length, "  io_center_profile size: ",
      scalar @$io_center_prof, " Should be the same but aren't.\n";

    # get score array for in-to-out using 3 profiles (center, nterm, cterm)
    my $io_score =
      $self->make_curve( $io_center_prof, $io_nterm_prof, $io_cterm_prof );

    # same for out-to-in
    my $oi_score =
      $self->make_curve( $oi_center_prof, $oi_nterm_prof, $oi_cterm_prof );

    my @io_helices = ();
    {
        my $fhresult;
        my $helix;
        my $start = 0;
        while ( $start < $length ) {
            ( $fhresult, $start, $helix ) =
              $self->find_helix( $length, $start, $io_score, $io_center_prof,
                $io_nterm_prof, $io_cterm_prof );
            $helix->set_nt_in('TRUE');    # io is inside-to-outside
            if ($fhresult) {
                push @io_helices, $helix;
            }
            last if ( $start >= $length );
        }
    }

    my @oi_helices = ();
    {
        my $fhresult;
        my $start = 0;
        my $helix;
        while ( $start < $length ) {
            ( $fhresult, $start, $helix ) =
              $self->find_helix( $length, $start, $oi_score, $oi_center_prof,
                $oi_nterm_prof, $oi_cterm_prof );
            $helix->set_nt_in('FALSE');
            if ($fhresult) {
                push @oi_helices, $helix;
            }
            last if ( $start >= $length );
        }
    }

#  print "number of io, oi helices: ", scalar @io_helices, "  ", scalar @oi_helices, "\n";
# $self->raw_output(\@io_helices, \@oi_helices);

    #  my @helices = ( @io_helices, @oi_helices );
    #   return $self->good_solutions_perl( \@helices );

    return ( \@io_helices, \@oi_helices );

    #   return ( scalar @io_helices, scalar @oi_helices );

}

sub short_output {
    my $self           = shift;
    my $id             = shift;
    my $sequence       = shift;
    my $io_helices_ref = shift;
    my $oi_helices_ref = shift;
    my $out            = "$id";
    $out .= substr( $sequence, 0, 3 ) . '...' . substr( $sequence, -3 ) . "\n";
    $out .=
        "Inside to outside helices : "
      . scalar @$oi_helices_ref . "\n"
      . "    from     to   score  center \n";

    for (@$io_helices_ref) {
        $out .= $_->get_descriptor_string() . "\n";
    }
    $out .=
        "Outside to inside helices : "
      . scalar @$oi_helices_ref . "\n"
      . "    from     to   score  center \n";
    for (@$oi_helices_ref) {
        $out .= $_->get_descriptor_string() . "\n";
    }

    return $out;
}

sub _long_output {
    my $self           = shift;
    my $id             = shift;
    my $sequence       = shift;
    my $io_helices_ref = shift;
    my $oi_helices_ref = shift;
    my $long_out       = "TMpred prediction output for : $id \n\n";
    $long_out .= "Sequence: " . substr( $sequence, 0, 3 ),
      '...' . substr( $sequence, -3 ) . "   length: " . length $sequence . "\n";
    $long_out .=
        "Prediction parameters: TM-helix length between "
      . $self->{min_tm_length} . " and "
      . $self->{max_tm_length}
      . "\n\n\n";

    $long_out .=
        "1.) Possible transmembrane helices\n"
      . "==================================\n"
      . "The sequence positions in brackets denominate the core region.\n"
      . "Only scores above "
      . $self->{min_score}
      . " are considered significant.\n\n";
    $long_out .=
        "Inside to outside helices : "
      . scalar @$io_helices_ref . "\n"
      . "    from     to   score  center \n";

    for (@$io_helices_ref) {
        $long_out .= $_->get_descriptor_string() . "\n";
    }
    $long_out .=
        "Outside to inside helices : "
      . scalar @$oi_helices_ref . "\n"
      . "    from     to   score  center \n";
    for (@$oi_helices_ref) {
        $long_out .= $_->get_descriptor_string() . "\n";
    }

    $long_out .=
        "2.) Table of correspondences\n"
      . "============================\n"
      . "Here is shown, which of the inside->outside helices correspond\n"
      . "to which of the outside->inside helices.\n";
    $long_out .=
        "  Helices shown in brackets are considered insignificant.\n"
      . 'A "+"  symbol indicates a preference of this orientation.' . "\n"
      . 'A "++" symbol indicates a strong preference of this orientation.'
      . "\n\n";

    $long_out .= "           inside->outside | outside->inside \n";

  
    my ( $n_io, $n_oi ) = ( scalar @$io_helices_ref, scalar @$oi_helices_ref );
    my ( $io_count, $oi_count ) = 0;
    my $padding = '         ';
    while (1) {
        if ( $io_count >= $n_io ) {
            $long_out .= "$padding | "
              . $oi_helices_ref->[$oi_count]->get_descriptor_string() . "\n";
            $oi_count++;
        }
        elsif ( $oi_count >= $n_oi ) {
            $long_out .= $io_helices_ref->[$io_count]->get_descriptor_string() . " | $padding" . "\n";
            $io_count++;
        }
        elsif ( $io_helices_ref->[$io_count]->get_center()->[0] <=
            $oi_helices_ref->[$oi_count]->get_center()->[0] )
        {    # io center is first (or same pos)
            if ( $oi_helices_ref->[$oi_count]->get_sh_nterm()->[0] <=
                $io_helices_ref->[$io_count]->get_sh_cterm()->[0] )
            {    # cores overlap
 $long_out .= $self->_io_oi_correspondence_string($io_helices_ref->[$io_count], $oi_helices_ref->[$oi_count]);
                # $long_out .=
                #     $io_helices_ref->[$io_count]->get_descriptor_string() . '|'
                #   . $oi_helices_ref->[$oi_count]->get_descriptor_string() . "\n";
                $io_count++;
                $oi_count++;
            }
            else {
 $long_out .= $self->_io_oi_correspondence_string($io_helices_ref->[$io_count], undef);
                # $long_out .= $io_helices_ref->[$io_count]->get_descriptor_string() . " | $padding " . "\n";
                $io_count++;
            }
        }
        else {    # oi center pos < io center pos

            if ( $io_helices_ref->[$io_count]->get_sh_nterm()->[0] <=
                $oi_helices_ref->[$oi_count]->get_sh_cterm()->[0] )
            {     # cores overlap
                $long_out .= $self->_io_oi_correspondence_string($io_helices_ref->[$io_count], $oi_helices_ref->[$oi_count]);
      #              $io_helices_ref->[$io_count]->get_descriptor_string() . '|'
      #            . $oi_helices_ref->[$oi_count]->get_descriptor_string() . "\n";
                $io_count++;
                $oi_count++;
            }
            else {
                $long_out .= $self->_io_oi_correspondence_string(undef, $oi_helices_ref->[$oi_count]);
		#  "$padding | " . $oi_helices_ref->[$oi_count]->get_descriptor_string() . "\n";
                $oi_count++;
            }

	  }
	last if($io_count >= $n_io  and $oi_count >= $n_oi);
      }

  $long_out .=
        "3.) Suggested models for transmembrane topology \n"
      . "============================================== \n"
      . "These suggestions are purely speculative and should be used with\n"
      . "EXTREME CAUTION since they are based on the assumption that\n"
      . "all transmembrane helices have been found.\n";
    $long_out .=
        "In most cases, the Correspondence Table shown above or the\n"
      . "prediction plot that is also created should be used for the\n"
      . "topology assignment of unknown proteins.\n\n";
    $long_out .=
      "2 possible models considered, only significant TM-segments used\n\n";

    return $long_out;
}

sub _io_oi_correspondence_string{
  my $self = shift;
  my $io_helix = shift;
my $oi_helix = shift;

#print "XXX: ", $self->{lo_orientational_threshold}, "  ", $self->{hi_orientational_threshold}, "\n";

my $io_string = '                       ';
my ($io_beg, $io_end, $io_score);
if(defined $io_helix){
 ($io_beg, $io_end, $io_score) = (
$io_helix->get_nterm()->[0] + 1, 
$io_helix->get_cterm()->[0] + 1, 
$io_helix->get_score()) ;

$io_string = sprintf("%6i-%4i (%2i) %5i ", 
		     $io_beg, $io_end, $io_end-$io_beg+1, $io_score);
}
my $oi_string = '                       ';
  my ($oi_beg, $oi_end, $oi_score);
if(defined $oi_helix){
  ($oi_beg, $oi_end, $oi_score) = (
$oi_helix->get_nterm()->[0] + 1, 
$oi_helix->get_cterm()->[0] + 1, 
$oi_helix->get_score()) ;

$oi_string = sprintf("%6i-%4i (%2i) %5i ", 
		     $oi_beg, $oi_end, $oi_end-$oi_beg+1, $oi_score);
}
my $score_diff;
my ($io_pref, $oi_pref) = '  ';
  if(!defined $io_helix){
    $oi_pref = '++';
  }elsif(!defined $oi_helix){
    $io_pref = '++';
  }else{ # both defined
    $score_diff = $io_score - $oi_score;
    if($score_diff > $self->{hi_orientational_threshold}){
      $io_pref = '++';
    }elsif($score_diff > $self->{lo_orientational_threshold}){
      $io_pref = ' +';
    }elsif($score_diff < -1*$self->{hi_orientational_threshold}){
      $oi_pref = '++';
    }elsif($score_diff < -1*$self->{lo_orientational_threshold}){
      $oi_pref = ' +';
    }else{
      # leave them both as '  '
    }
  }
return $io_string . $io_pref . '  |' . $oi_string . $oi_pref . "\n";
}

sub good_solutions_perl {    # get the good solutions starting from the
                             # output of perl tmpred
    my $self               = shift;
    my $tmpred_helices_ref = shift;    # ref to array of predicted helices

    my $good_string = "";
    foreach (@$tmpred_helices_ref) {
        my $beg    = $_->get_nterm()->[0] + 1;
        my $end    = $_->get_cterm()->[0] + 1;
        my $score  = $_->get_score();
        my $length = $end + 1 - $beg;

        if (    $score >= $self->{min_score}
            and $length >= $self->{min_tm_length}
            and $length <= $self->{max_tm_length}
            and $beg >= $self->{min_beg}
            and $beg <= $self->{max_beg} )
        {
            $good_string .= "($score,$beg,$end)  ";
        }
    }
    if ( $good_string eq "" ) { $good_string = "(-10000,0,0)"; }
    else { $good_string =~ s/(\s+)$//; }    # eliminate whitespace at end.
    return $good_string;
}

sub make_profile {                          # makes a profile, i.e. an array
                                            # containing ...
                                            # my $sequence     = shift;
    my $seq_aanumber_array = shift;         # ref to array of numbers
    my $table              = shift;
    my $ref_position = $table->get_marked_position();
    my $matrix       = $table->get_table();
    my $ncols   = scalar @{ $matrix->[0] };   # ncols is # elements in first row
    my @profile = ();
    my $length = scalar @$seq_aanumber_array;    #length $sequence;

# need to be careful here to make it consistent with pascal code,
# which has 1-based arrays. ref_position is 1 less than in pascal code for this reason.
    for ( my $k = 0 ; $k < $length ; $k++ ) {
        my $m    = 0;
        my $kmrf = $k - $ref_position;
        my $plo  = $kmrf;                        # $k - $ref_position;
        $plo = 0 if ( $plo < 0 );
        my $pup = $ncols + $kmrf;                # $k - $ref_position;
        $pup = $length if ( $pup > $length );

        for (
            my ( $p, $i ) = ( $plo, $plo - $kmrf ) ;  # $k - $ref_position ) ) ;
            $p < $pup ; $p++, $i++
          )
        {
            $m += $matrix->[ $seq_aanumber_array->[$p] ]->[$i];
        }

        my $round_m =
          ( $m < 0 ) ? int( $m * 100 - 0.5 ) : int( $m * 100 + 0.5 );
        push @profile, $round_m;
{'fraction_non_alphanumeric': 0.15114859917007947, 'fraction_numerical': 0.01046893222315641, 'mean_word_length': 2.0987477939322634, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 320, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2238164', 'n_tokens_mistral': 13523, 'n_tokens_neox': 12831, 'n_words': 3041}