bb.scaftigs

#!/usr/bin/perl
#----------------------------------------------------------#
#        Author: Douglas Senalik dsenalik@wisc.edu         #
#----------------------------------------------------------#
# "Black Box" program series
=bb
Split FASTA sequences into scaftigs
=cut bb
use strict;
use warnings;
use Getopt::Long;      # for getting command line parameters



############################################################
# configuration variables
############################################################
my $defaultminlen = 1;
my $defaultprefix = ".";
my $defaultwrap = 100;



############################################################
# global variables
############################################################
my $ansiup    = "\033[1A";  # terminal control
(my $prognopath = $0) =~ s/^.*[\/\\]//;
my $ntotallines = 0;
my $nseq = 0;
my $nbases = 0;
my $nscaftigs = 0;
my $nscafbases = 0;



############################################################
# command line parameters
############################################################
my @infilenames;         # input file names
my $outfilename;         # output file name
my $minlen        = $defaultminlen;  # number of Ns needed to split
my $prefix        = $defaultprefix;  # scaftig number prefix
my $wrap          = $defaultwrap;
my $notrim;              # trim header after first white space
my $ucase;               # translate all bases to upper case
my $help;                # print help and exit
my $quiet;               # only show errors
my $debug;               # print extra debugging information
GetOptions (
            "infiles=s"      => \@infilenames,       # string
            "outfile=s"      => \$outfilename,       # string
            "minlen=i"       => \$minlen,            # integer
            "prefix=s"       => \$prefix,            # string
            "notrim"         => \$notrim,            # flag
            "ucase"          => \$ucase,             # flag
            "wrap=i"         => \$wrap,              # integer
            "help"           => \$help,              # flag
            "quiet"          => \$quiet,             # flag
            "debug"          => \$debug);            # flag
# debug implies not quiet
if ( $debug ) { $quiet = 0; }
unless ( @infilenames ) { $help = 1 }
unless ( $outfilename ) { $help = 1 }
# quiet mode if printing to STDOUT
if ( $outfilename eq "-" ) { $quiet = 1 }



############################################################
# print help screen
############################################################
if ( $help )
  {
    print "$prognopath
This program will return scaftigs, all portions of a final assembly consisting
of contiguous sequence, with sequences split at every occurrence of gaps of
unknown bases (Ns)

Required parameters:
  --infile=xxx      input assembly FASTA file name (all scaffolds, contigs, etc.),
                    .gz compressed okay, or use - for STDIN,
                    multiple uses of --infile are allowed
  --outfile=xxx     output FASTA file name (all scaftigs),
                    .gz compressed okay, or use - for STDOUT
Optional parameters:
  --notrim          default is to trim original header after first
                    white space, this disables that behavior
  --prefix=xxx      prefix scaftig number with this text,
                    default prefix is \"$defaultprefix\"
  --minlen=xxx      gap must be at least this many Ns long,
                    default is split if $defaultminlen or more Ns
  --wrap=xxx        wrap sequence lines to this length, default
                    is $defaultwrap, or use --wrap=0 for no wrap
  --help            print this screen
  --quiet           only print error messages
  --debug           print extra debugging information
";
    exit 1;
  } # if ( $help )



############################################################
# main processing loop
############################################################
{
my $OUTF = stdopen ( ">", $outfilename );
for my $infilename ( @infilenames )
  {
    my $INF = stdopen ( "<", $infilename );
    debugmsg ( "Processing file \"$infilename\"" );
    my $nlines = 0;
    my $currhdr;
    my $currseq;
    while ( my $aline = <$INF> )
      {
        $nlines++;
        $aline =~ s/[\r\n]//g;
        if ( ( $nlines == 1 ) and ( $aline !~ m/^>/ ) )
          { die "Error, input file does not appear to be a FASTA file, it should start with \">\"\n"; }
        if ( $aline =~ m/^>/ )
          {
            writeone ( $OUTF, $currhdr, $currseq );
            $currhdr = $aline;
            undef $currseq;
          }
        else  # sequence line
          {
            if ( $ucase ) { $aline =~ tr/a-z/A-Z/; }
            $currseq .= $aline;
          }
      } # while <$INF>
    writeone ( $OUTF, $currhdr, $currseq );
    stdclose ( $INF );
    $ntotallines += $nlines;
  } # for $infilename
stdclose ( $OUTF );
}



############################################################
# cleanup and end program
############################################################
unless ( $quiet )
  {
    print commify($ntotallines), " input lines processed from ",
          commify(scalar @infilenames), " input file", ( ( scalar @infilenames ) > 1 )?"s":"", "\n";
    print commify($nseq), " input sequences, ",
          commify($nbases), " input bases\n";
    print commify($nscaftigs), " output scaftigs, ",
          commify($nscafbases), " output bases\n";
  }
exit 0;



############################################################
sub writeone { my ( $OUTF, $hdr, $seq ) = @_;
############################################################
  # uses global variables $nseq, $nbases, $nscaftigs, $nscafbases, $prefix, $wrap
  my $scaftigidx = 0;

  # first pass will be undefined $hdr and $seq
  return if ( ( ! defined $hdr ) and ( ! defined $seq ) );

  # update global input sequence counter
  $nseq++;
  $nbases += length ( $seq );

  # trim header after first white space, if desired
  unless ( $notrim ) { $hdr =~ s/\s.*$//; }

  # while match non-gap loop
  while ( $seq =~ m/([^Nn]{$minlen,})/g )
    {
      my $tmpseq = $1;
      $nscafbases += length ( $tmpseq );
      $scaftigidx++;
      $nscaftigs++;
      my $tmphdr = $hdr . $prefix . $scaftigidx;

      # wrap sequence lines, if desired, $wrap=0 disables wrapping
      if ( $wrap )
        { $tmpseq =~ s/(.{$wrap})/$1\n/g; }
      unless ( $tmpseq =~ m/\n$/ ) { $tmpseq .= "\n"; }

      print $OUTF $tmphdr, "\n", $tmpseq;
    } # while match non-gap

} # sub writeone



############################################################
sub debugmsg { my ( $text, $noreturn, $nolinenum ) = @_;
############################################################
  if ( $debug )
    {
      my ($package, $filename, $line, $sub) = caller(0);
      unless ( $nolinenum ) { $text = "Line $line: " . $text; }
      if ( ! ( $noreturn ) ) { $text .= "\n"; }
      print $text;
    } # if ( $debug )
} # sub debugmsg



###############################################################
sub timestr {
###############################################################
  @_ = localtime(shift || time);
  return(sprintf("%04d/%02d/%02d %02d:%02d", $_[5]+1900, $_[4]+1, $_[3], @_[2,1]));
} # sub timestr



###############################################################
sub commify {
###############################################################
# http://perldoc.perl.org/perlfaq5.html#How-can-I-output-my-numbers-with-commas
  local $_ = shift;
  1 while s/^([-+]?\d+)(\d{3})/$1,$2/;
  return $_;
} # commify



###############################################################
sub stdopen { my ( $mode, $filename, $extratext ) = @_;
###############################################################
# a replacement for the three-parameter open which also allows
# the use of "-" as the file name to mean STDIN or STDOUT
  my $fh;  # the file handle
  if ( $filename eq "-" )  # only exact match to "-" has special meaning
    {
      if ( $mode =~ m/>/ )
        { $fh = *STDOUT }
      else
        { $fh = *STDIN }
    }
  else
    {
      # supplemental passed text for error messages, need one more space
      if ( defined $extratext )
        { $extratext .= " " }
      else
        { $extratext = "" }

      my $text;  # this is only used for error message
      if ( $mode =~ m/^\+?>>/ )  # ">>" or "+>>"
        { $text = "append" }
      elsif ( $mode =~ m/^\+?>/ )  # ">" or "+>"
        { $text = "output" }
      elsif ( $mode =~ m/^\+?</ )  # "<" or "+<"
        { $text = "input" }
      elsif ( $mode eq "-|" )
        { $text = "piped input" }
      elsif ( $mode eq "|-" )
        { $text = "piped output" }
      else
        { die "Error, unsupported file mode \"$mode\" specified to stdopen( $mode, $filename, $extratext )\n"; }

      # if file name ends in ".gz", gzip compression is assumed, and handle it transparently
      if ( $filename =~ m/\.gz$/ )
        {
          if ( $mode =~ m/^>$/ ) # output mode
            { $mode = "|-"; $filename = "gzip -c > \"$filename\""; }
          elsif ( $mode =~ m/^<$/ ) # input mode
            { $mode = "-|"; $filename = "gunzip -c \"$filename\""; }
          else
            { die "Error, can't handle gzip compression with mode \"$mode\" for file \"filename\"\n"; }
        } # if gzip compressed file
      elsif ( $filename	=~ m/\.bz2$/ )
       	{
          if ( $mode =~ m/^>$/ ) # output mode
            { $mode = "|-"; $filename = "bzip2 -c > \"$filename\""; }
          elsif ( $mode =~ m/^<$/ ) # input mode
            { $mode = "-|"; $filename = "bunzip2 -c \"$filename\""; }
          else
            { die "Error, can't handle bzip2 compression with mode \"$mode\" for file \"filename\"\n"; }
       	}
      open ( $fh, $mode, $filename ) or die ( "Error opening ${extratext}file \"$filename\" for $text: $!\n" );
    }
  # return the opened file handle to the caller
  return $fh;
} # sub stdopen



###############################################################
sub stdclose { my ( $fh ) = @_;
###############################################################
# same as built-in close, except in case of STDIN or STDOUT,
# and in this case the file handle is not closed

  unless ( fileno($fh) <= 2 )  # if file number is this low, is stdin or stdout or stderr
    { close ( $fh ) or die ( "Error closing file handle: $!\n" ); }

} # sub stdclose



# eof