code snip to decide Phred encoding of FASTQ file

It's a common task to check the Phred quality score encoding of fastq file. Sure, it's fairly easy to check this by eyes (according to the rul below), and also there might be already some tools for that purpose.  The task is essentially a problem of converting ASCII character into its decimal format.

  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS.....................................................
  ..........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX......................
  ...............................IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII......................
  .................................JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ......................
  LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL....................................................
  !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~
  |                         |    |        |                              |                     |
 33                        59   64       73                            104                   126

 S - Sanger        Phred+33,  raw reads typically (0, 40)
 X - Solexa        Solexa+64, raw reads typically (-5, 40)
 I - Illumina 1.3+ Phred+64,  raw reads typically (0, 40)
 J - Illumina 1.5+ Phred+64,  raw reads typically (3, 40)
    with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator (bold) 
    (Note: See discussion above).
 L - Illumina 1.8+ Phred+33,  raw reads typically (0, 41)

To be simple, I just wrote a line of code for that, by only looking the first lines of reads. Here it is:

head input.fastq | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}'

To use it, you can also save it into a bash, like

#!/bin/sh

inputfile=$1

# Update: The following part for checking the file extension can be simplified (Thanks to the comment from Unknown) 
[[ "$inputfile" =~ ".*\.(fq\.gz$|fastq\.gz$)" ]] && zcat $inputfile | head -n40 ...
[[ "$inputfile" =~ ".*\.(fq$|fastq$)" ]] && head -n40 $inputfile | ...

less $inputfile | head -n40 | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}'

The only trick is the command od, which is to display/convert file content in a specific format. Here is the more detail for od (http://publib.boulder.ibm.com/infocenter/pseries/v5r3/index.jsp?topic=/com.ibm.aix.cmds/doc/aixcmds4/od.htm).

simple way to get reads length distribution of FASTQ files

• Using perl: 

cat input.fq | perl -ne '$s=<>;<>;<>;chomp($s);print length($s)."\n";' > input.readslength.txt

• Using awk:

cat input.fq | awk '{if(NR%4==2) print length($1)}' > input.readslength.txt

• if zipped file, using:

 zcat input.fq.gz | ...

• get length statistics:

sort input.readslength.txt | uniq -c

• or using fancy code from Jim Kent's utility:

textHistogram


So, one line code for all input fastq files would be:

find *.fq.gz -not -name \*raw\* -printf "zcat %p | awk '{if(NR%%4==2) print length(\$1)}' | textHistogram -maxBinCount=59 stdin \n" | sh


Note that you have to use double % to escape the % character for printf formatting control, just like in C. (Thanks for zencuke's answer here)

You will get something like this:

RNAseq.20E_library.result_primary.clean.fa

large values truncated: need 35 bins or larger binSize than 1

Maximum value 53.000000

 18 ********* 27730

 19 ******************* 58997

 20 ************************************************************ 186919

 21 ************************ 74536

 22 ************** 45171

 23 ************* 39107

 24 *************** 45560

 25 *********************** 70452

 26 ************************************************* 154030

 27 ************************************************************ 187704

 28 ************************** 81198

 29 ***** 17016

 30 ** 5341

 31 * 2439

 32  0

 33  1

 34  0

 35  0

 36  0

 37  0

 38  0

 39  1

<minVal or >= 40  173

RNAseq.Day_13_library.result_primary.clean.fa

large values truncated: need 35 bins or larger binSize than 1

Maximum value 53.000000

 18 ***** 22570

 19 ********* 40335

 20 ***************************** 127999

 21 ********************* 95179

 22 ****************** 79808

 23 ********* 39596

 24 ******* 29423

 25 ************* 55438

 26 ************************************* 164868

 27 ************************************************************ 265722

 28 *************************** 120353

 29 ********* 38625

 30 *** 14684

 31 * 5214

 32  0

 33  0

 34  1

<minVal or >= 35  140