<commit_before>#!/usr/bin/env perl

=head1 NAME

activeTE_msa - evaluate MSAs for characteristics of active TEs

=head1 SYNOPSIS

activeTE_msa family.aln

=head1 DESCRIPTION

=head1 Command line arguments

 -flank   [integer] Length of flanking sequence (default 100)
 -all     Process all

=head1 AUTHORS - <NAME>

<NAME> - bradc-AT-ucr.edu

<NAME> - jason.stajich-AT-ucr.edu [contributor]

=cut

use strict;
use warnings;
use Bio::AlignIO;
use Bio::SearchIO;
use Bio::Seq;
use Bio::SeqIO;
use Bio::PrimarySeq;
use Bio::Tools::IUPAC;
use FindBin;
use File::Spec;
use Getopt::Long;
use Data::Dumper;

#set default flank length and reset if option is given

my $flank = 100;
my $trimal =
  'trimal';   # specify the full path to this if it can't be found automatically
my $all;

my $PROGRAM_NAME = "activeTE";

GetOptions(
  'flank:i'  => \$flank,
  'all'      => \$all,
  'trimal:s' => \$trimal,
);

if ( defined $all ) {
  print "Intermediate files will be kept\n";
}
else {
  print "Intermediate files will be removed\n";
}

# make sure only one valid input file is specified
if ( @ARGV != 1 ) {
  print "Please specify only one input file\n";
  help();
}
my $infile = shift @ARGV;

if ( !-f $infile ) {
  print "Supplied filepath is not valid: $!\n";
  help();
}
if ( $flank <= 25 ) {
  print "Please use alignments with flanks  >25 nucleotides long\n";
  help();
}

sub help {
  print "

usage:

activeTE_msa.pl -a -f <int> <multiple alignment file>

-a keep all intermediate files  DEFAULT = remove
-f length of sequence flanks  DEFAULT = 100

The test alignments all have f = 100.

";
  exit 1;
}

#breakdown directory and filename, create output directory
my ( $volume, $in_dir, $filename ) = File::Spec->splitpath($infile);
#my ($fname_fin) = split( '\\.', $filename );    # first bit after first '.' is the name
my @fname_fin = split( '\\.', $filename );    # first bit after first '.' is the name
pop @fname_fin;
my $fname_fin = join ('.',@fname_fin);
my $out_dir_name = "activeTE_out_" . $fname_fin;
my $out_path = File::Spec->catdir( $in_dir, $out_dir_name );
if ( !-d $out_path ) {
  mkdir($out_path)
    or die "Can't create dir:$out_path $!\n";    # use built-in Perl mkdir
}

my %element_info;    # this hash is for storing element information
my %tir_positions;

#generate trimal alignment gap summary file, using sed to remove the header infmal -in $infile -sgco

open( my $trimal_run => "$trimal -in $infile -sgc | sed -n '4,\$p' |" )
  || die "Cannot run trimal!";

# parse the trimal output file with the gap %  of each position in the full alignment
my @gap_id_array;
while ( my $line = <$trimal_run> ) {
  chomp $line;       # probably not necessary
  if ( $line =~ /\S+\s+(\S+)\s+\S+/ )
  {  
    push @gap_id_array, $1;
  }
}
close($trimal_run);

#initialize a global array to store results later

my @good_aln = ();
my ( $initial_left_TSD, $initial_right_TSD );
## replacing white space with underscores in MSA
my @in_array;
open( my $in => "<$infile" );
while ( my $line = <$in> ) {
  $line =~ s/ /_/g;

  # did you want to skip lines that don't have any whitespace?
  push @in_array, $line;
}
close($in);
open( my $fix_out => ">$infile" );
print $fix_out @in_array;
close($fix_out);
@in_array = ();

# the above could just be re-written
# `perl -i -p -e 's/ /_/g;' $infile`;

#create input object and an alignment object from it
my $in_obj = Bio::AlignIO->new( -file => $infile, -format => 'fasta' );
my $full_aln_obj;
unless ( $full_aln_obj = $in_obj->next_aln() ) {
  die "Cannot find any alignment in file $infile\n";
}

my $full_aln_len = $full_aln_obj->length();

#generate a gap column matrix from the full alignment
my $gap_cols = $full_aln_obj->gap_col_matrix();

#calculate the % nucleotide identity and the fraction of copies with sequence at each position the full alignment and print to file

my $full_id_out_path =
  File::Spec->catpath( $volume, $out_path, $filename . ".full_id" );

my @full_id_array;

if ( !-e $full_id_out_path ) {
  open( my $id_out => ">$full_id_out_path" );
  warn "Starting Full ID calculation\n";
  for ( my $i = 1 ; $i <= $full_aln_len ; $i++ ) {
    my $pos             = $full_aln_obj->slice( $i, $i, 1 );
    my $pos_id          = $pos->percentage_identity;
    my $gap_col_hashref = $gap_cols->[ $i - 1 ];
    my $total_count     = 0;
    my $base_count      = 0;<commit_msg>changed warn to a print<commit_after>#!/usr/bin/env perl

=head1 NAME

activeTE_msa - evaluate MSAs for characteristics of active TEs

=head1 SYNOPSIS

activeTE_msa family.aln

=head1 DESCRIPTION

=head1 Command line arguments

 -flank   [integer] Length of flanking sequence (default 100)
 -all     Process all

=head1 AUTHORS - <NAME>

<NAME> - bradc-AT-ucr.edu

<NAME> - jason.stajich-AT-ucr.edu [contributor]

=cut

use strict;
use warnings;
use Bio::AlignIO;
use Bio::SearchIO;
use Bio::Seq;
use Bio::SeqIO;
use Bio::PrimarySeq;
use Bio::Tools::IUPAC;
use FindBin;
use File::Spec;
use Getopt::Long;
use Data::Dumper;

#set default flank length and reset if option is given

my $flank = 100;
my $trimal =
  'trimal';   # specify the full path to this if it can't be found automatically
my $all;

my $PROGRAM_NAME = "activeTE";

GetOptions(
  'flank:i'  => \$flank,
  'all'      => \$all,
  'trimal:s' => \$trimal,
);

if ( defined $all ) {
  print "Intermediate files will be kept\n";
}
else {
  print "Intermediate files will be removed\n";
}

# make sure only one valid input file is specified
if ( @ARGV != 1 ) {
  print "Please specify only one input file\n";
  help();
}
my $infile = shift @ARGV;

if ( !-f $infile ) {
  print "Supplied filepath is not valid: $!\n";
  help();
}
if ( $flank <= 25 ) {
  print "Please use alignments with flanks  >25 nucleotides long\n";
  help();
}

sub help {
  print "

usage:

activeTE_msa.pl -a -f <int> <multiple alignment file>

-a keep all intermediate files  DEFAULT = remove
-f length of sequence flanks  DEFAULT = 100

The test alignments all have f = 100.

";
  exit 1;
}

#breakdown directory and filename, create output directory
my ( $volume, $in_dir, $filename ) = File::Spec->splitpath($infile);
#my ($fname_fin) = split( '\\.', $filename );    # first bit after first '.' is the name
my @fname_fin = split( '\\.', $filename );    # first bit after first '.' is the name
pop @fname_fin;
my $fname_fin = join ('.',@fname_fin);
my $out_dir_name = "activeTE_out_" . $fname_fin;
my $out_path = File::Spec->catdir( $in_dir, $out_dir_name );
if ( !-d $out_path ) {
  mkdir($out_path)
    or die "Can't create dir:$out_path $!\n";    # use built-in Perl mkdir
}

my %element_info;    # this hash is for storing element information
my %tir_positions;

#generate trimal alignment gap summary file, using sed to remove the header infmal -in $infile -sgco

open( my $trimal_run => "$trimal -in $infile -sgc | sed -n '4,\$p' |" )
  || die "Cannot run trimal!";

# parse the trimal output file with the gap %  of each position in the full alignment
my @gap_id_array;
while ( my $line = <$trimal_run> ) {
  chomp $line;       # probably not necessary
  if ( $line =~ /\S+\s+(\S+)\s+\S+/ )
  {  
    push @gap_id_array, $1;
  }
}
close($trimal_run);

#initialize a global array to store results later

my @good_aln = ();
my ( $initial_left_TSD, $initial_right_TSD );
## replacing white space with underscores in MSA
my @in_array;
open( my $in => "<$infile" );
while ( my $line = <$in> ) {
  $line =~ s/ /_/g;

  # did you want to skip lines that don't have any whitespace?
  push @in_array, $line;
}
close($in);
open( my $fix_out => ">$infile" );
print $fix_out @in_array;
close($fix_out);
@in_array = ();

# the above could just be re-written
# `perl -i -p -e 's/ /_/g;' $infile`;

#create input object and an alignment object from it
my $in_obj = Bio::AlignIO->new( -file => $infile, -format => 'fasta' );
my $full_aln_obj;
unless ( $full_aln_obj = $in_obj->next_aln() ) {
  die "Cannot find any alignment in file $infile\n";
}

my $full_aln_len = $full_aln_obj->length();

#generate a gap column matrix from the full alignment
my $gap_cols = $full_aln_obj->gap_col_matrix();

#calculate the % nucleotide identity and the fraction of copies with sequence at each position the full alignment and print to file

my $full_id_out_path =
  File::Spec->catpath( $volume, $out_path, $filename . ".full_id" );

my @full_id_array;

if ( !-e $full_id_out_path ) {
  open( my $id_out => ">$full_id_out_path" );
  print "Starting Full ID calculation\n";
  for ( my $i = 1 ; $i <= $full_aln_len ; $i++ ) {
    my $pos             = $full_aln_obj->slice( $i, $i, 1 );
    my $pos_id          = $pos->percentage_identity;
    my $gap_col_hashref = $gap_cols->[ $i - 1 ];
    my $total_count     = 0;
    my $base_count      = 0;
{'fraction_non_alphanumeric': 0.14313612866168868, 'fraction_numerical': 0.0068925904652498565, 'mean_word_length': 3.608787718369508, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '886439', 'n_tokens_mistral': 3319, 'n_tokens_neox': 3187, 'n_words': 1055}