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

=head1 NAME 
                                                                       
pairfq.pl - Match paired-end sequences from separate FastA/Q files

=head1 SYNOPSIS    
 
pairfq.pl -f s_1_1_trim.fq -r s_1_2_trim.fq -fp s_1_1_trim_paired.fq -rp s_1_2_trim_paired.fq -fs s_1_1_trim_unpaired.fq -rs s_1_2_trim_unpaired.fq

=head1 DESCRIPTION
     
Re-pair paired-end sequences that may have been separated by quality trimming.
This script also writes the unpaired forward and reverse sequences to separate 
files so that they may be used for assembly or mapping. The input may be FastA
or FastQ format in either Illumina 1.3+ or Illumina 1.8 format.

=head1 DEPENDENCIES

This script uses the Perl modules AnyDBM_File and AnyDBM_File::Importer. AnyDBM_File
is used to create a database with either DB_File or SQLite and AnyDBM_File::Importer
allows symbols from other packages to be imported. In this script, symbols from DB_File
are imported so that the database may be stored in memory.

=head1 LICENSE
 
   The MIT License

   Copyright (c) 2013, <NAME>.

   Permission is hereby granted, free of charge, to any person obtaining
   a copy of this software and associated documentation files (the
   "Software"), to deal in the Software without restriction, including
   without limitation the rights to use, copy, modify, merge, publish,
   distribute, sublicense, and/or sell copies of the Software, and to
   permit persons to whom the Software is furnished to do so, subject to
   the following conditions:

   The above copyright notice and this permission notice shall be
   included in all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   SOFTWARE.
 
=head1 TESTED WITH:

=over

=item *
Perl 5.14.1 (Red Hat Enterprise Linux Server release 5.7 (Tikanga))

=item *
Perl 5.14.2 (Red Hat Enterprise Linux Desktop release 6.2 (Santiago); Fedora 17)

=back

=head1 AUTHOR 

<NAME>                                                

=head1 CONTACT
 
statonse at gmail dot com

=head1 REQUIRED ARGUMENTS

=over 2

=item -f, --forward

The file of forward sequences from an Illumina paired-end sequencing run.

=item -r, --reverse                                                                                                                                                       
The file of reverse sequences from an Illumina paired-end sequencing run.

=item -fp, --forw_paired

The output file to place the paired forward reads.

=item -rp, --rev_paired                                                                                                                                                  

The output file to place the paired reverse reads. 

=item -fs, --forw_unpaired                                                                                                                                                  
The output file to place the unpaired forward reads. 

=item -rs, --rev_unpaired                                                                                                                                                  
The output file to place the unpaired reverse reads. 

=back

=head1 OPTIONS

=over 2

=item -im, --in_memory

Construct the database in memory. May be faster, but will obviously use more memory.

=item -h, --help

Print a usage statement. 

=item -m, --man

Print the full documentation.

=cut      

use v5.10;
use strict;
use warnings;
use File::Basename;
use Getopt::Long;
use Data::Dumper;
BEGIN {
  @AnyDBM_File::ISA = qw( DB_File SQLite_File )
      unless @AnyDBM_File::ISA == 1; # if loaded already, AnyDBM_File::ISA has a length of one;
}
use AnyDBM_File;
use vars qw( $DB_BTREE &R_DUP );
use AnyDBM_File::Importer qw(:bdb);
use Pod::Usage;

my ($fread, $rread, $fpread, $rpread, $fsread, $rsread, $memory, $help, $man);

GetOptions(
           'f|forward=s'        => \$fread,
           'r|reverse=s'        => \$rread,
           'fp|forw_paired=s'   => \$fpread,
           'rp|rev_paired=s'    => \$rpread,
           'fs|forw_unpaired=s' => \$fsread,
           'rs|rev_unpaired=s'  => \$rsread,
           'im|in_memory'       => \$memory,
           'h|help'             => \$help,
           'm|man'              => \$man,
          ) || pod2usage( "Try '$0 --man' for more information." );

#
# Check @ARGV
#
usage() and exit(0) if $help;

pod2usage( -verbose => 2 ) if $man;

if (!$fread  || !$rread  ||
    !$fpread || !$rpread ||
    !$fsread || !$rsread) {
    say "\nERROR: Command line not parsed correctly. Check input.\n";
    usage();
    exit(1);
}

my %rseqhash;
$DB_BTREE->{cachesize} = 100000;
$DB_BTREE->{flags} = R_DUP;
my $db_file = "pairfq.bdb";

unless (defined $memory) { 
   tie %rseqhash, 'AnyDBM_File', $db_file, O_RDWR|O_CREAT, 0666, $DB_BTREE
       or die "\nERROR: Could not open DBM file $db_file: $!\n";
}

my @raux = undef;
my ($fname, $fseq, $fqual, $fid, $rname, $rseq, $rqual, $rid);
my ($fct, $rct, $fpct, $rpct, $pct, $fsct, $rsct, $sct) = (0, 0, 0, 0, 0, 0, 0, 0);
open my $r, '<', $rread or die "\nERROR: Could not open file: $rread\n";

while (($rname, $rseq, $rqual) = readfq(\*$r, \@raux)) {
    $rct++;
    if ($rname =~ /(\/\d)$/) {
	$rname =~ s/$1//;
    } 
    elsif ($rname =~ /\s(\d\:\w\:\d\:\w+)$/) { # as below, may need to add some abstraction to this regex 
	$rid = $1; chomp $rid;
	$rid =~ s/^.//;
	$rname =~ s/\s.*//;
	$rname .= "|".$rid;
    } else {
	die "\nERROR: Could not determine Illumina encoding. Exiting.\n";
    }
    $rseqhash{$rname} = join "\t",$rseq, $rqual if defined $rqual;
    $rseqhash{$rname} = $rseq if !defined $rqual;
}
close $r;
open my $f, '<', $fread or die "\nERROR: Could not open file: $fread\n";
open my $fp, '>', $fpread or die "\nERROR: Could not open file: $fpread\n";
open my $rp, '>', $rpread or die "\nERROR: Could not open file: $rpread\n";
open my $fs, '>', $fsread or die "\nERROR: Could not open file: $fsread\n";

my ($forw_id, $rev_id);
my @faux = undef;
while (($fname, $fseq, $fqual) = readfq(\*$f, \@faux)) {
    $fct++;
    if ($fname =~ /(\/\d)$/) {
	$fname =~ s/\/\d//;
    }
    elsif ($fname =~ /\s(\d\:\w\:\d\:\w+)$/) {   # may need to add some abstraction to this regex
        $fid = $1; chomp $fid;
        $fid =~ s/^.//;
	$fname =~ s/\s.*//;
	$fname .= "|".$fid;
    } else {
	die "\nERROR: Could not determine Illumina encoding. Exiting.\n";
    }
    if ($fname =~ /\|/) {
	my ($name, $id) = split /\|/, $fname;
	$forw_id = $name." 1".$id;
	$rev_id  = $name." 2".$id;
    }
    if (exists $rseqhash{$fname}) {
	$fpct++; $rpct++;
	if (defined $fqual) {
	    my ($rread, $rqual) = split /\t/, $rseqhash{$fname};
	    if ($fname =~ /\|/) {
		say $fp join "\n","@".$forw_id, $fseq, "+", $fqual;
		say $rp join "\n","@".$rev_id, $rread, "+", $rqual;
	    } else {
		say $fp join "\n","@".$fname."/1", $fseq, "+", $fqual;
                say $rp join "\n","@".$fname."/2", $rread, "+", $rqual;
	    }
	} else {
	    if ($fname =~ /\|/) {
		say $fp join "\n",">".$forw_id, $fseq;
		say $rp join "\n",">".$rev_id, $rseqhash{$fname};
	    } else {
                say $fp join "\n",">".$fname."/1", $fseq;
                say $rp join "\n",">".$fname."/2", $rseqhash{$fname};
            }
	}
        delete $rseqhash{$fname};
    } else {
	$fsct++;
	if (defined $fqual) {
	    if ($fname =~ /\|/) {
		say $fs join "\n","@".$forw_id, $fseq, "+", $fqual;
	    } else {
		say $fs join "\n","@".$fname."/1", $fseq, "+", $fqual;
	    }
	} else {
	    if ($fname =~ /\|/) {
		say $fs join "\n",">".$forw_id, $fseq;
	    } else {
		say $fs join "\n",">".$fname."/1", $fseq;
	    }
	}
    }
}
close $f;
close $fp;
close $rp;
close $fs;
open my $rs, '>', $rsread or die "\nERROR: Could not open file: $rsread\n";

my $rev_id_up;
while (my ($rname_up, $rseq_up) = each %rseqhash) {
    $rsct++;
    if ($rname_up =~ /\|/) {
	my ($uname, $uid) = split /\|/, $rname_up;
	$rev_id_up .= $uname." "."2".$uid;
    }
    if ($rseq_up =~ /\t/) {
	my ($rread_up, $rqual_up) = split /\t/, $rseq_up;
	if ($rname_up =~ /\|/) {
	    say $rs join "\n","@".$rev_id_up, $rread_up, "+", $rqual_up;
	} else {
	    say $rs join "\n","@".$rname_up."/2", $rread_up, "+", $rqual_up;
	}
    } else {
	if ($rname_up =~ /\|/) {
	    say $rs join "\n",">".$rev_id_up, $rseq_up;
	} else {
	    say $rs join "\n",">".$rname_up."/2", $rseq_up;
	}
    }
    $rev_id_up = undef;
}
close $rs;

$pct = $fpct + $rpct;
$sct = $fsct + $rsct;
untie %rseqhash if defined $memory;

say "Total forward reads in $fread:              $fct";
say "Total reverse reads in $rread:              $rct";
say "Total paired reads in $fpread and $rpread:  $pct";
say "Total forward paired reads in $fpread:      $fpct";
say "Total reverse paired reads in $rpread:      $rpct";
say "Total upaired reads in $fsread and $rsread: $sct"; 
say "Total forward unpaired reads in $fsread:    $fsct";
say "Total reverse unpaired reads in $rsread:    $rsct";

exit;
<commit_msg>added code to remove dbm files and left comment about a potential bug<commit_after>#!/usr/bin/env perl

=head1 NAME 
                                                                       
pairfq.pl - Match paired-end sequences from separate FastA/Q files

=head1 SYNOPSIS    
 
pairfq.pl -f s_1_1_trim.fq -r s_1_2_trim.fq -fp s_1_1_trim_paired.fq -rp s_1_2_trim_paired.fq -fs s_1_1_trim_unpaired.fq -rs s_1_2_trim_unpaired.fq

=head1 DESCRIPTION
     
Re-pair paired-end sequences that may have been separated by quality trimming.
This script also writes the unpaired forward and reverse sequences to separate 
files so that they may be used for assembly or mapping. The input may be FastA
or FastQ format in either Illumina 1.3+ or Illumina 1.8 format.

=head1 DEPENDENCIES

This script uses the Perl modules AnyDBM_File and AnyDBM_File::Importer. AnyDBM_File
is used to create a database with either DB_File or SQLite and AnyDBM_File::Importer
allows symbols from other packages to be imported. In this script, symbols from DB_File
are imported so that the database may be stored in memory.

=head1 LICENSE
 
   The MIT License

   Copyright (c) 2013, <NAME>.

   Permission is hereby granted, free of charge, to any person obtaining
   a copy of this software and associated documentation files (the
   "Software"), to deal in the Software without restriction, including
   without limitation the rights to use, copy, modify, merge, publish,
   distribute, sublicense, and/or sell copies of the Software, and to
   permit persons to whom the Software is furnished to do so, subject to
   the following conditions:

   The above copyright notice and this permission notice shall be
   included in all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
   SOFTWARE.
 
=head1 TESTED WITH:

=over

=item *
Perl 5.14.1 (Red Hat Enterprise Linux Server release 5.7 (Tikanga))

=item *
Perl 5.14.2 (Red Hat Enterprise Linux Desktop release 6.2 (Santiago); Fedora 17)

=back

=head1 AUTHOR 

<NAME>                                                

=head1 CONTACT
 
statonse at gmail dot com

=head1 REQUIRED ARGUMENTS

=over 2

=item -f, --forward

The file of forward sequences from an Illumina paired-end sequencing run.

=item -r, --reverse                                                                                                                                                       
The file of reverse sequences from an Illumina paired-end sequencing run.

=item -fp, --forw_paired

The output file to place the paired forward reads.

=item -rp, --rev_paired                                                                                                                                                  
The output file to place the paired reverse reads. 

=item -fs, --forw_unpaired                                                                                                                                                  
The output file to place the unpaired forward reads. 

=item -rs, --rev_unpaired                                                                                                                                                  
The output file to place the unpaired reverse reads. 

=back

=head1 OPTIONS

=over 2

=item -im, --in_memory

Construct the database in memory. May be faster, but will obviously use more memory.

=item -h, --help

Print a usage statement. 

=item -m, --man

Print the full documentation.

=cut      

##TODO: There is a bug returning too many reads to the reverse singletons file if the dbm file is used. 7/2/13

use v5.10;
use strict;
use warnings;
use File::Basename;
use Getopt::Long;
use Data::Dumper;
BEGIN {
  @AnyDBM_File::ISA = qw( DB_File SQLite_File )
      unless @AnyDBM_File::ISA == 1; # if loaded already, AnyDBM_File::ISA has a length of one;
}
use AnyDBM_File;
use vars qw( $DB_BTREE &R_DUP );
use AnyDBM_File::Importer qw(:bdb);
use Pod::Usage;

my ($fread, $rread, $fpread, $rpread, $fsread, $rsread, $memory, $help, $man);

GetOptions(
           'f|forward=s'        => \$fread,
           'r|reverse=s'        => \$rread,
           'fp|forw_paired=s'   => \$fpread,
           'rp|rev_paired=s'    => \$rpread,
           'fs|forw_unpaired=s' => \$fsread,
           'rs|rev_unpaired=s'  => \$rsread,
           'im|in_memory'       => \$memory,
           'h|help'             => \$help,
           'm|man'              => \$man,
          ) || pod2usage( "Try '$0 --man' for more information." );

#
# Check @ARGV
#
usage() and exit(0) if $help;

pod2usage( -verbose => 2 ) if $man;

if (!$fread  || !$rread  ||
    !$fpread || !$rpread ||
    !$fsread || !$rsread) {
    say "\nERROR: Command line not parsed correctly. Check input.\n";
    usage();
    exit(1);
}

my %rseqhash;
$DB_BTREE->{cachesize} = 100000;
$DB_BTREE->{flags} = R_DUP;
my $db_file = "pairfq.bdb";
unlink $db_file if -e $db_file;

unless (defined $memory) { 
   tie %rseqhash, 'AnyDBM_File', $db_file, O_RDWR|O_CREAT, 0666, $DB_BTREE
       or die "\nERROR: Could not open DBM file $db_file: $!\n";
}

my @raux = undef;
my ($fname, $fseq, $fqual, $fid, $rname, $rseq, $rqual, $rid);
my ($fct, $rct, $fpct, $rpct, $pct, $fsct, $rsct, $sct) = (0, 0, 0, 0, 0, 0, 0, 0);
open my $r, '<', $rread or die "\nERROR: Could not open file: $rread\n";

while (($rname, $rseq, $rqual) = readfq(\*$r, \@raux)) {
    $rct++;
    if ($rname =~ /(\/\d)$/) {
	$rname =~ s/$1//;
    } 
    elsif ($rname =~ /\s(\d\:\w\:\d\:\w+)$/) { # as below, may need to add some abstraction to this regex 
	$rid = $1; chomp $rid;
	$rid =~ s/^.//;
	$rname =~ s/\s.*//;
	$rname .= "|".$rid;
    } else {
	die "\nERROR: Could not determine Illumina encoding. Exiting.\n";
    }
    $rseqhash{$rname} = join "\t",$rseq, $rqual if defined $rqual;
    $rseqhash{$rname} = $rseq if !defined $rqual;
}
close $r;
open my $f, '<', $fread or die "\nERROR: Could not open file: $fread\n";
open my $fp, '>', $fpread or die "\nERROR: Could not open file: $fpread\n";
open my $rp, '>', $rpread or die "\nERROR: Could not open file: $rpread\n";
open my $fs, '>', $fsread or die "\nERROR: Could not open file: $fsread\n";

my ($forw_id, $rev_id);
my @faux = undef;
while (($fname, $fseq, $fqual) = readfq(\*$f, \@faux)) {
    $fct++;
    if ($fname =~ /(\/\d)$/) {
	$fname =~ s/\/\d//;
    }
    elsif ($fname =~ /\s(\d\:\w\:\d\:\w+)$/) {   # may need to add some abstraction to this regex
        $fid = $1; chomp $fid;
        $fid =~ s/^.//;
	$fname =~ s/\s.*//;
	$fname .= "|".$fid;
    } else {
	die "\nERROR: Could not determine Illumina encoding. Exiting.\n";
    }
    if ($fname =~ /\|/) {
	my ($name, $id) = split /\|/, $fname;
	$forw_id = $name." 1".$id;
	$rev_id  = $name." 2".$id;
    }
    if (exists $rseqhash{$fname}) {
	$fpct++; $rpct++;
	if (defined $fqual) {
	    my ($rread, $rqual) = split /\t/, $rseqhash{$fname};
	    if ($fname =~ /\|/) {
		say $fp join "\n","@".$forw_id, $fseq, "+", $fqual;
		say $rp join "\n","@".$rev_id, $rread, "+", $rqual;
	    } else {
		say $fp join "\n","@".$fname."/1", $fseq, "+", $fqual;
                say $rp join "\n","@".$fname."/2", $rread, "+", $rqual;
	    }
	} else {
	    if ($fname =~ /\|/) {
		say $fp join "\n",">".$forw_id, $fseq;
		say $rp join "\n",">".$rev_id, $rseqhash{$fname};
	    } else {
                say $fp join "\n",">".$fname."/1", $fseq;
                say $rp join "\n",">".$fname."/2", $rseqhash{$fname};
            }
	}
        delete $rseqhash{$fname};
    } else {
	$fsct++;
	if (defined $fqual) {
	    if ($fname =~ /\|/) {
		say $fs join "\n","@".$forw_id, $fseq, "+", $fqual;
	    } else {
		say $fs join "\n","@".$fname."/1", $fseq, "+", $fqual;
	    }
	} else {
	    if ($fname =~ /\|/) {
		say $fs join "\n",">".$forw_id, $fseq;
	    } else {
		say $fs join "\n",">".$fname."/1", $fseq;
	    }
	}
    }
}
close $f;
close $fp;
close $rp;
close $fs;
open my $rs, '>', $rsread or die "\nERROR: Could not open file: $rsread\n";

my $rev_id_up;
while (my ($rname_up, $rseq_up) = each %rseqhash) {
    $rsct++;
    if ($rname_up =~ /\|/) {
	my ($uname, $uid) = split /\|/, $rname_up;
	$rev_id_up .= $uname." "."2".$uid;
    }
    if ($rseq_up =~ /\t/) {
	my ($rread_up, $rqual_up) = split /\t/, $rseq_up;
	if ($rname_up =~ /\|/) {
	    say $rs join "\n","@".$rev_id_up, $rread_up, "+", $rqual_up;
	} else {
	    say $rs join "\n","@".$rname_up."/2", $rread_up, "+", $rqual_up;
	}
    } else {
	if ($rname_up =~ /\|/) {
	    say $rs join "\n",">".$rev_id_up, $rseq_up;
	} else {
	    say $rs join "\n",">".$rname_up."/2", $rseq_up;
	}
    }
    $rev_id_up = undef;
}
close $rs;

$pct = $fpct + $rpct;
$sct = $fsct + $rsct;
untie %rseqhash if defined $memory;
unlink $db_file if -e $db_file;

say "Total forward reads in $fread:              $fct";
say "Total reverse reads in $rread:              $rct";
say "Total paired reads in $fpread and $rpread:  $pct";
say "Total forward paired reads in $fpread:      $fpct";
say "Total reverse paired reads in $rpread:      $rpct";
say "Total upaired reads in $fsread and $rsread: $sct"; 
say "Total forward unpaired reads in $fsread:    $fsct";
say "Total reverse unpaired reads in $rsread:    $rsct";

exit;

{'fraction_non_alphanumeric': 0.15590013140604467, 'fraction_numerical': 0.009461235216819974, 'mean_word_length': 2.2187447132464895, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 64, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3593853', 'n_tokens_mistral': 7027, 'n_tokens_neox': 6670, 'n_words': 2165}