bb.tassel

#!/usr/bin/perl
#----------------------------------------------------------#
#        Author: Douglas Senalik dsenalik@wisc.edu         #
#             https://github.com/dsenalik/bb               #
#----------------------------------------------------------#
# "Black Box" program series
# version 1 2016-03-08 initial commit
# version 2 2019-02-04 Add --barcodefile, --column, and --gzlevel
# version 3 2019-07-25 Add --force
# version 4 2022-04-09 Add --directory and --parallel
# version 5 2022-09-04 Add --outdirectory
# version 6 2023-04-25 Support overhang length other than 4, support longer barcodes
# version 7 2024-06-03 Automatically generate an MD5 checksum file
=bb
Add barcodes to reverse reads for Tassel GBS pipeline
=cut bb
use strict;
use warnings;
use Getopt::Long;      # for getting command line parameters
use File::Basename;
use POSIX ":sys_wait_h";  # to define WNOHANG constant



############################################################
# configuration variables
############################################################
my @defaultoverhangs = ( "CAGC", "CTGC" );  # ApeKI



############################################################
# global variables
############################################################
my $ansiup    = "\033[1A";  # terminal control
(my $prognopath = $0) =~ s/^.*[\/\\]//;
my @bch;  # array of hashes for looking up barcode hits
my $minbarcodelength;    # length of shortest barcode
my $maxbarcodelength;    # length of longest barcode
my $overhanglength;      # length of strings in overhangs



############################################################
# command line parameters
############################################################
my $directory;           # process entire directory
my $outdirectory;        # destination of process entire directory
my $parallel     = 4;    # process this many in parallel
my $finfilename  = "";   # input file name
my $rinfilename  = "";   # input file name
my $outfilename  = "";   # output file name
my $barcodefilename;     # optional file with barcodes, one per line
my $barcodecolumn = 1;   # if using --barcodefile
my $bcset        = 1;    # predefined barcode set number
my @overhangs;
my $gzlevel      = 9;    # gzip compression level, default to best compression
my $force;               # used to overwrite existing output files
my $nochecksum   = 0;    # suppress generation of am MD5 checksum file when finished
my $help         = 0;    # print help and exit
my $quiet        = 0;    # only show errors
my $debug        = 0;    # print extra debugging information
GetOptions (
            "directory=s"    => \$directory,         # string
            "outdirectory=s" => \$outdirectory,      # string
            "finfile=s"      => \$finfilename,       # string
            "rinfile=s"      => \$rinfilename,       # string
            "outfile=s"      => \$outfilename,       # string
            "overhang=s"     => \@overhangs,         # string
            "barcodefile=s"  => \$barcodefilename,   # string
            "bcset=i"        => \$bcset,             # integer
            "column=i"       => \$barcodecolumn,     # integer
            "force"          => \$force,             # flag
            "gzlevel=i"      => \$gzlevel,           # integer
            "parallel=i"     => \$parallel,          # integer
            "nochecksum"     => \$nochecksum,        # flag
            "help"           => \$help,              # flag
            "quiet"          => \$quiet,             # flag
            "debug"          => \$debug);            # flag
# debug implies not quiet
if ( $debug ) { $quiet = 0; }
if ( ( $finfilename ) or ( $rinfilename ) or ( $outfilename ) )
  {
    unless ( $finfilename ) { print "Error, --finfile not specified\n"; $help = 1; }
    unless ( $rinfilename ) { print "Error, --rinfile not specified\n"; $help = 1; }
    unless ( $outfilename ) { print "Error, --outfile not specified\n"; $help = 1; }
  }
elsif ( !$directory )
  { $help = 1; }
unless ( @overhangs ) { @overhangs = @defaultoverhangs; }
if ( $outfilename eq "-" ) { $quiet = 1 }



############################################################
# print help screen
############################################################
if ( $help )
  {
    print "$prognopath
This program will preprocess paired-end Illumina reads from
GBS (Genotyping By Sequencing) experiments to make them
compatible with TASSEL. This involves copying the barcode
from the forward reads to the beginning of the reverse reads,
since TASSEL cannot otherwise identify the reverse reads which
do not have a barcode.

Required parameters:
  either
  --finfile=xxx     input FASTQ file name of forward reads,
                    can be gzip compressed if .gz extension
  --rinfile=xxx     input FASTQ file name of reverse reads,
  --outfile=xxx     output file name, compressed if .gz extension
                    or use - for STDOUT
  or
  --directory=xxx   process all read pairs in this directory
Optional parameters:
  --outdirectory=xxx  when using --directory, destination for
                    processed output files. Default is same directory
  --parallel=xxx    when using --directory, process up to this
                    many pairs of reads files in parallel,
                    default is 4
  --overhang=xxx    sequence(s) of restriction enzyme
                    overhang, e.g. for ApeKI (this is the
                    default), use --overhang twice:
                    --overhang=CAGC --overhang=CTGC
                    they must all be the same length
  --bcset=xxx       predefined barcode set, 1 (default) or 2
  --barcodefile=xxx specify a custom set of barcodes
  --column=xxx      which column within file --barcodefile
                    contains the barcode, default=1 (1-based)
  --force           overwrite output files if they exist
  --gzlevel=xxx     gzip compression level, 1(fast) to 9(best),
                    default=9
  --nochecksum      suppress generation of an MD5 checksum file
  --help            print this screen
  --quiet           only print error messages
  --debug           print extra debugging information
";
    exit 1;
  } # if ( $help )



############################################################
# initialization
############################################################
init();

my $r;
if ( $directory )
  {
    $r = processdirectory( $directory, $outdirectory, $parallel, $force );
  }
else
  {
    $r = processreads( $finfilename, $rinfilename, $outfilename, $force );
  }



############################################################
# cleanup and end program
############################################################
exit $r;



############################################################
sub processreads { my ( $finfilename, $rinfilename, $outfilename, $force ) = @_;
############################################################
if ( ( -s $outfilename ) and ( ! $force ) )
  {
    unless ( $quiet )
      { print "Warning, output file \"$outfilename\" exists, not processing. Use --force to override\n"; }
    return 1;
  }
my $FINF = stdopen ( "<", $finfilename );
my $RINF = stdopen ( "<", $rinfilename );
my $FOUTF = stdopen ( ">", $outfilename );
my $ROUTF = $FOUTF;  # for future, in case want non-interleaved file
my $nfound = 0;
my $nnotfound = 0;
my $nlines = 0;
while ( ! eof $FINF )
  {
    my @flines;
    my @rlines;
    for my $i ( 0..3 )
      {
        # if either file ends prematurely, it is probably corrupted, so stop with error
        if ( eof $FINF )
          { die "Error, premature ending of forward file line ".($nlines+$i)." \"$finfilename\"\n"; }
        if ( eof $RINF )
          { die "Error, premature ending of reverse file line ".($nlines+$i)." \"$rinfilename\"\n"; }
        ( $flines[$i] = <$FINF> ) =~ s/[\r\n]//g;
        ( $rlines[$i] = <$RINF> ) =~ s/[\r\n]//g;
      }

    # find barcode
    my $bclength;
    for (my $i=$maxbarcodelength; $i>=$minbarcodelength; $i--)
      {
        if ( $bch[$i]->{substr($flines[1],0,$i+$overhanglength)} )
          {
            $bclength = $i;
            last;
          }
      }

    # add barcode to reverse read. Reverse read should already have enzyme site
    if ( $bclength )
      {
        $nfound++;
        substr($rlines[1],0,0) = substr($flines[1],0,$bclength);
        substr($rlines[3],0,0) = substr($flines[3],0,$bclength);

        # save to interleaved file
        for my $i ( 0..3 )
          { print $FOUTF $flines[$i], "\n"; }
        for my $i ( 0..3 )
          { print $ROUTF $rlines[$i], "\n"; }
      }
    else # forward read lacking any known barcode+enzyme, discard both F and R reads
      { $nnotfound++; }

    $nlines += 4;

  } # while <$INF>
stdclose ( $FOUTF );
stdclose ( $RINF );
stdclose ( $FINF );
unless ($nochecksum)
  { generatemd5($outfilename); }
unless ( $quiet ) { print "$outfilename: Found " . commify($nfound) . " read pairs with barcodes, "
                    . commify($nnotfound) . " pairs without barcodes discarded\n"; }
} # sub processreads



############################################################
sub processdirectory { my ( $directory, $outdirectory, $parallel, $force ) = @_;
############################################################
# the assumption is that reverse reads files are identical to forward except
# that _R1_ is replaced with _R2_, and for output _R1_001 is replaced with _tassel
# or _R1_ is replaced with _tassel_
my @fwdreads = glob( $directory . '/*_R1*.fastq*' );
push( @fwdreads, glob( $directory . '/*_R1*.fq*' ) );
unless ( @fwdreads )
  {
    print "Error, no files found in \"$directory\" matching *_R1*.fastq* or *_R1*.fq*\n";
    return 1;
  }
# append slash if not present
if ( $outdirectory )
  { $outdirectory =~ s|/*$|/|; }

my @pids;
my $i = 0;  # for debug message
foreach my $finfilename ( @fwdreads )
  {
    # exclude md5 checksum files that may have been picked up by the glob
    next if ($finfilename =~ m|\.md5$|);

    $i++;
    # this regex pattern is guaranteed to be in the name based on the glob
    (my $rinfilename = $finfilename ) =~ s|_R1_|_R2_|;
    (my $outfilename = $finfilename ) =~ s|_R1_\d+|_tassel|;
    $outfilename =~ s|_R1_|_tassel_|;
    if ( $outdirectory )
      { $outfilename = $outdirectory . basename( $outfilename ); }

    # limit number of concurrent processes
    while ( countchildren(@pids) >= $parallel )
      { sleep(60); }

    # split into threads
    my $pid = fork();

    if ( not defined $pid )
      { die "Fork failed thread $i $!\n"; }
    elsif ( $pid == 0 )  # i.e., child
      {
        debugmsg ( "Child $i starting \"$outfilename\" ".timestr() );
        processreads( $finfilename, $rinfilename, $outfilename, $force );
        debugmsg ( "Child $i finished \"$outfilename\" ".timestr() );
        exit 0;
      }
    else  # parent
      {
        push ( @pids, $pid );
        debugmsg ( "Launched child $i pid $pid \"$outfilename\" ".timestr() );
        # short sleep so that if child fails immediately, we
        # can avoid the longer sleep in the loop
        sleep( 1 );
      }
  }

# wait for all threads to finish
debugmsg ( "Waiting for threads to finish ".timestr() );
foreach ( @pids )
  { waitpid ( $_, 0 ); }
debugmsg ( "All threads finished ".timestr() );

} # sub processdirectory



###############################################################
sub countchildren { my ( @pids ) = @_;
###############################################################
my $nrunning = 0;
foreach my $pid ( @pids )
  {
    my $returncode = waitpid($pid, WNOHANG);
    unless ( $returncode ) { $nrunning++; }
  }
debugmsg( "countchildren returns $nrunning" );
return( $nrunning );
} # sub countchildren



############################################################
sub generatemd5 { my ($filename) = @_;
############################################################
  my $md5 = $filename . '.md5';
  unless (-s $md5)
    {
      my $cmd = "md5sum \"$filename\" > \"$md5\"";
      run($cmd);
    }
} # sub generatemd5



############################################################
sub init {
############################################################
# initialize barcode array of hashes @bch from list of barcodes
# the order of barcodes within @bc is not important;
# initializes global variables $overhanglength and @bch

  # define the default barcode sets
  my @bc1 = qw/ AACCGAGA ACAGGGAA ACGACTAC ACGTGGTA CCATGGGT CCGGATAT CGCCTTAT CGCGGAGA CGTGTGGT
                GCTGTGGA GGATTGGT GTGAGGGT TAGCATGC TAGGCCAT TATCGGGA TCTCAGTC TGCAAGGA TGGTACGT TTCCTGGA
                AAAAGTT AACGCCT AATATGC ACGTGTT ATGAAAC ATTAATT ATTGGAT CATAAGT CGCTGAT CGGTAGA CTACGGA
                CTTGCTT GAACTTC GAATTCA GCGGAAT GGACCTA GTCGATT TAGCGGA TATTTTT TCGAAGA TCTGTGA TGCTGGA
                ACCTAA AGTGGA ATATGT ATCGTA ATGCCT CATCGT CCACAA CCAGCT CGCGGT CTATTA CTTCCA
                GAGATA GCCAGT GCTCTA GGAAGA GGTTGT GTACTT GTTGAA TAACGA TAGGAA TGGCTA TTCAGA
                ACAAA ACCGT AGCCC AGGAT ATTGA CATCT CCTAC CGCTT CTAGC CTGTA GAGGA
                GCTTA GGAAC GGTGT GTATT GTCAA TAATA TACAT TCACC TCGTT TGCGA TTCTC
                AACT ACTA AGGC CAGA CGAT CTCC GATC GCGT GTAA TCAC TGCA /;
  my @bc2 = qw/ TGCT TTGGCA CGCCAT CGCTCA CGTACC CTATGGA TACTGAT GTCTGAA TCAGTAT ACCGAGT GCAACGT
                GCCAGAT CAGCT GATGTC GTCACT ACGCTC TTCAAGT TACGGTA CTTGAGA ATCAGTT GTTACGA ATGTCAA
                ATAGTCA CCAGGTA GACTC TGTTAC CATCGC CTGATC GAGCAGT ACTATGT GATCATA TGTGCAA TGCATAT
                GGTGGCA TAGCCAT TGGCAAC ACGTC TGACCT AGGTCT CCAGTC ACACGGT CGTGAAT GCATTGA CGACAGT
                AGTTCGA TCGTTAA AGGAGTC ATGATCT CACGTT TTGCAC CCTGCA TACGCC GACGTGA AACTTGT CAGGTAT
                ACGTGTA ACGCTGA TGATCAT GTAGAGC CATGTTA GTTAGC CTAGCT ACCTGC CCGTAC TCATAGT TGACGTA
                TGCAATA GATGCAT TAAGCTT CGTCGGA GACCTAT TGTAAGC AGCATT TAGTGC CTTGAC TCGACC TTACGAT
                GCTATAA ATATCGT CTAATGT CTTAGAT GACTATT ATTAGCA GGACCTT CTCCGA TGAATC GCTCAC GCATCC
                GGCTAGA ATCGTAT AGTCTAT GCGTAAT TGATACA TAGGTCA TGTCATT GTGTACA /;

  # optional override of default barcodes with a file of barcodes
  my @bc;
  if ( $barcodefilename )
    { @bc = loadbarcodefile( $barcodefilename, $barcodecolumn ); }
  elsif ( $bcset == 1 )  # default
    { @bc = @bc1; }
  elsif ( $bcset == 2 )
    { @bc = @bc2; }
  else
    { die "No barcode file, or set is invalid\n"; }

  # determine range in lengths of barcodes
  foreach my $barcode (@bc)
    {
      my $l = length($barcode);
      if (!defined($minbarcodelength) or $l < $minbarcodelength) { $minbarcodelength = $l; }
      if (!defined($maxbarcodelength) or $l > $maxbarcodelength) { $maxbarcodelength = $l; }
    }

  # check that all --overhang parameters are the same length
  $overhanglength = length($overhangs[0]);
  foreach my $overhang ( @overhangs )
    {
      if ( length($overhang) != $overhanglength )
        { die "Error, all overhang sequences must be the same length\n"; }
    }

  # initialize global @bch barcode+enzyme lookup hashes
  foreach my $abarcode ( @bc )
    {
      my $l = length($abarcode);
      foreach my $overhang ( @overhangs )
        { $bch[$l]->{ $abarcode . uc($overhang) } = 1; }
    }

} # sub init



############################################################
sub loadbarcodefile { my ( $infilename, $column ) = @_;
############################################################
# read file of barcode sequences and return as array. $column is 1-based
  my $INF = stdopen( '<', $infilename );
  my @bc;
  while ( my $aline = <$INF> )
    {
      $aline =~ s/[\r\n]//g;
      next if ( $aline =~ m/^#/ );  # skip comment lines
      my @cols = split( /\t/, $aline );
      if ( $cols[$column-1] )
        { push( @bc, $cols[$column-1] ); }
    }
  stdclose( $INF );
  unless( $quiet ) { print "Read ".scalar(@bc)." barcodes from file \"$infilename\"\n"; }
  return( @bc );
} # sub loadbarcodefile



############################################################
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 -$gzlevel -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
      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



############################################################
sub run { my ( $command, $errorokay ) = @_;
############################################################
# Run a system command, but use bash. Perl by default uses sh.
# Command may be piped multiple commands or redirected from or to files.
  debugmsg( "Running command \"$command\"" );
  my $result = system( 'bash', '-o', 'pipefail', '-c', $command );
  if ( ( $result ) and ( ! $errorokay ) )
    {
      # exitvalue 141 can occur if piping to head or grep -q
      my $exitvalue = $result >> 8;
      my $signal = $result & 255;
      die( "Error $exitvalue:$signal running command \"$command\"\n" );
    }
  return( $result );
} # sub run



# eof