We do the same thing but add a pairFixer script. Basically it goes through the pair files twice to fix the pairs. The first time to see all the headers(just the common part) that are seen twice and a second time to create these files: pair1 pair2 singles

I also added a -M12 to the clipper part to tolerate non-full length adapters but maybe 12 is a bit too short.

[HTML][HTML]

!/usr/bin/perl

use strict;

my $outName = $ARGV[0]; my $pair1 = $ARGV[1]; my $pair2 = $ARGV[2]; my $pair1FH; my $pair2FH; my $pair1outFH; my $pair2outFH; my $pairSingleFH;

open ($pair1FH, '-|', 'zcat', "$pair1") || die "Cannot Open pair 1 File"; open ($pair2FH, '-|', 'zcat', "$pair2") || die "Cannot Open pair 2 File"; my @pairedFiles = ($pair1FH,$pair2FH);

my %readCount; my $totalRecords = 0; for my $pairFH (@pairedFiles) { while (my $head = <$pairFH>) { <$pairFH>; <$pairFH>; <$pairFH>; $totalRecords++; my $key = substr($head, 0, length($head)-3); my $value = $readCount{$key}; if(defined $value) { $value++; $readCount{$key} = $value; } else { $readCount{$key} = 1; } } close($pairFH); }

open ($pair1FH, '-|', 'zcat', "$pair1") || die "Cannot Open pair 1 File"; open ($pair2FH, '-|', 'zcat', "$pair2") || die "Cannot Open pair 2 File"; @pairedFiles = ($pair1FH,$pair2FH); open ($pair1outFH, ">$outName.pair1.fastq") || die "Cannot Open pair out 1 File"; open ($pair2outFH, ">$outName.pair2.fastq") || die "Cannot Open pair out 2 File"; open ($pairSingleFH, ">$outName.single.fastq") || die "Cannot Open single File"; my @outputFiles = ($pair1outFH,$pair2outFH);

my $idx=0; my $records=0; my $pairs=0; my $single=0; print "Copying reads\n"; for my $pairFH (@pairedFiles) { my $outFH = $outputFiles[$idx]; $idx++;

while(my $head = <$pairFH>) { my $seq = <$pairFH>; my $head2 = <$pairFH>; my $qual = <$pairFH>;

my $key = substr($head, 0, length($head)-3);
my $value = $readCount{$key};
my $fh;
$records++;
if($value == 1) {
  $fh = $pairSingleFH;
  $single++;
}
elsif($value > 1) {
  $fh = $outFH;
  $pairs++;
}
else {
  die "No value found for: $key\n";
}

print $fh $head;
print $fh $seq;
print $fh $head2;
print $fh $qual;

print "\rTotal: $totalRecords, Records: $records, Pairs: $pairs, Single: $single ";

} close($outFH); close($pairFH); } close($pairSingleFH); print "Total: $totalRecords, Records: $records, Pairs: $pairs, Single: $single\n"; [HTML] [HTML]

Haibao Tang has a program that will accept a FASTA file of adaptors and trim those all in one go. It does some sort of local alignment on the adaptors so it will be more thorough tan the fastx_clipper and it can also trim low quality reads.

i think those fastx_clipper commands are not pair-friendly - i.e. if one of your sequences is "all adapter" you will orphan the mate and mess up the pairing.