[Ntobe]

@GenoMax, we suspected that our reads might be somehow corrupted during the ftp transer. We used fastQC to check for quality and trimming didn't seem necessary. I'm still exploring trying different tophat threds and options to speed up mapping. Thanks for your suggestions

[Brian Bushnell]

Thanks so much Brian Bushnell. The reason we trying to explore cleaning the reads is because our tophat jobs were running out of time in the server before completion. Our raw read files are too big (~13GB per .gz file) and we were thinking that some of the reads might not be of good quality. If that the case, why not remove them and map only good quality reads? Again, this might not be a good approach and that why I'm seeking help.

You can quality-filter, but it will increase bias if you do it aggressively. Quality-trimming is a much better approach for quantitative analysis. Incidentally, BBMap is substantially faster than Tophat, and STAR is (from what I have heard) substantially faster than BBMap. Switching to a faster aligner would be a better approach than discarding enough of the data to finish in a fixed time window.

[cement_head]

If it's useful to anyone this is a small script I knocked up when we had to process some fastq files which were corrupted during an FTP transfer. You can pipe data through it and it does some basic sanity checks to ensure that the file looks like valid fastq data. It will remove any entries which look broken and leave you just the good stuff.

#!/usr/bin/perl
use warnings;
use strict;

while (<>) {

  unless (/^\@/) {
    warn "$_ should have had an \@ at the start and it didn't\n";
    next;
  }
  my $id1 = $_;
  my $seq = <>;
  my $id2 = <>;
  my $qual = <>;

  if ($seq =~/^[@+]/) {
    warn "Sequence '$seq' looked like an id";
    next;
  }
  if ($qual =~/^[@+]/) {
    warn "Quality '$qual' looked like an id";
    next;
  }
  if ($id2 !~ /^\+/) {
    warn "Midline '$id2' didn't start with a +";
    next;
  }

  if ($qual =~ /[GATCN]{20,}/) {
    warn "Quality '$qual' looked like sequence";
    next;
  }

  if (length($seq) != length($qual)) {
    warn "Seq $seq and Qual $qual weren't the same length";
    next;
  }

  print $id1,$seq,$id2,$qual;


}

Hi,

This is my lack of PERL shining through, but is there a way to get a two output files, one with the "good" reads and one with the "bad" reads?

Thanks,
Andor

[GenoMax]

Save code below in a file (check.pl in this example) and then run as

$ perl check.pl input_seq.fq good.fq bad.fq

Note: This code assumes that every fastq entry has 4 lines in the file. It does not check for gross problems (i.e. missing an entire line in a record). Use at your own risk, not extensively tested :-)

-----------------------------------

#!/usr/bin/perl
use warnings;
use strict;

my $infile = $ARGV[0];
my $outfile1 = $ARGV[1];
my $outfile2 = $ARGV[2];

open (IN, "$infile") or die "can't open the outputfile: $infile\n";
open (OUT1, ">$outfile1") or die "can't open the outputfile: $outfile1\n";
open (OUT2, ">$outfile2") or die "can't open the outputfile: $outfile2\n";


while (<>) {

  my $id1 = $_;
  my $seq = <>;
  my $id2 = <>;
  my $qual = <>;
  
  if ($id1 !~ /^[\@]/) {
        print OUT2 $id1,$seq,$id2,$qual;
    next;
  }

  if ($seq !~ /[GATCN]+/g) {
        print OUT2 $id1,$seq,$id2,$qual;
    next;
  }
  if ($id2 !~ /^\+/) {
        print OUT2 $id1,$seq,$id2,$qual;
    next;
  }

  if (length($seq) != length($qual)) {
    warn "Seq $seq and Qual $qual weren't the same length";
        print OUT2 $id1,$seq,$id2,$qual;
    next;
  }

  print OUT1 $id1,$seq,$id2,$qual;

}
close IN;
close OUT1;
close OUT2;

Edit : The regex for base checking is not working so a record with sequence characters other than ACGTN will get through.