base quality encoding changed after "bwa samse" command

tomjan · 02-25-2013, 02:45 PM

hello,

Please look at base quality string in my sample2New.fq file

@EBRI093151_0051:4:55:2998:9540#0/1
ACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGAT
+
FFFFFFFFFFDFFBFDEAEEEFFFFFFFFCFFEFFCEEEDDFFEEEFEADDFDFDEEDFFE@FCDDD>ACDFADD?CCECDB<?@047:9@?BB+B@@@]]

after commands

opt/bwa-0.6.2/bwa index -a bwtsw -p ref reference.fa
/opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai -I ref sample2New.fq
/opt/bwa-0.6.2/bwa samse -f sample2New.sam -r "@RG\tID:sample2\tPL:ILLUMINA\tPU

u1\tLB:sample2\tSM:sample2" ref sample2New.sai sample2New.fq

I can see changed base quality string in the sample2New.sam file

EBRI093151_0051:4:55:2998:9540#0 0 Chr10 377653 0 101M * 0 0 ACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGAT ''''''''''%''#'%&"&&&''''''''$''&''$&&&%%''&&&'&"%%'%'%&&%''&!'$%%%^_"$%'"%% $$&$%#^] !^Q^U^XESC^Z! ##^L#!!!>> RG:Z:sample2 XT:A:R NM:i:0 X0:i:3 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101 XA:Z:Chr10,+33,101M,0;Chr10,+242847,101M,0;

and ofcourse the command

java -Xmx8g -jar /opt/picard-tools-1.85/SortSam.jar SO=coordinate INPUT=sample2New.sam OUTPUT=sample2New.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true

fails with error

Exception in thread "main" java.lang.IllegalArgumentException: Invalid fastq character:

Why "bwa samse" is changing quality encoding??
Do you have an idea what Im doing wrong?

thanks

maubp · 02-26-2013, 08:35 AM

Could you edit your post to use the [ code ] and [ /code ] tags? This is easily done via the advanced editor view where there is a button for this in the tool bar (not shown in the quick reply edit box).

maubp · 02-26-2013, 08:41 AM

I've not checked all the bases (due to the forum formatting), however, it would appear to be down to a FASTQ encoding problem. It appears bwa defaulted to assuming the obsolete Illumina specific ASCII encoding of PHRED+64, while your data was actually the original standard Sanger ASCII encoding of PHRED+33 (now adopted by Illumina). For background, see:
http://dx.doi.org/10.1093/nar/gkp1137

In your FASTQ file, the first base has quality code 'F', ASCII character 70. Under the Sanger FASTQ scheme that means 70-33 = quality 37. However, if read in as the obsolete Illumina scheme it would be 70-64 = 6 quality, which when output again in SAM format (which uses the Sanger FASTQ scheme) becomes 6+33 = ASCII 39 = ' (single quote).

Solution - there is a command line option to tell bwa you have a Sanger style FASTQ file. Use it, otherwise you get a bad SAM/BAM file.

tomjan · 02-26-2013, 01:12 PM

Thank you for your help Maubp,

Your explanations helped me to find the solutions

The problem was in "bwa aln" cmmand
/opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai -I ref sample2New.fq

from the documentation we can se "-I The input is in the Illumina 1.3+ read format (quality equals ASCII-64). ". So, everything is OK when I ommit the -I option.

/opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai ref sample2New.fq

Once again, Thank You for your help.

maubp · 02-26-2013, 01:23 PM

Well done

- and thank you for posting back with the details for anyone searching about this again in the future.

(I couldn't remember the details about the switch, and wasn't at a machine where I could quickly check - but this way you'll probably remember the problem and solution

)

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
What are "Base-Space" reads? [BWA]	dp05yk	Bioinformatics	3	03-21-2012 04:44 PM
"allele balance ratio" and "quality by depth" in VCF files	efoss	Bioinformatics	2	10-25-2011 12:13 PM
FastQC "Per Base Sequence Content": systematic deviation at 3' end of reads	d f	Illumina/Solexa	4	09-28-2010 10:46 AM
"Systems biology and administration" & "Genome generation: no engineering allowed"	seb567	Bioinformatics	0	05-25-2010 01:19 PM
Help: what a "*" indicates in the reads base column of pileup format?	rainbow	Illumina/Solexa	1	01-29-2010 04:59 AM

02-25-2013, 02:45 PM	#1
tomjan Junior Member Location: pl Join Date: Apr 2012 Posts: 8	base quality encoding changed after "bwa samse" command hello, Please look at base quality string in my sample2New.fq file @EBRI093151_0051:4:55:2998:9540#0/1 ACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGAT + FFFFFFFFFFDFFBFDEAEEEFFFFFFFFCFFEFFCEEEDDFFEEEFEADDFDFDEEDFFE@FCDDD>ACDFADD?CCECDB<?@047:9@?BB+B@@@]] after commands opt/bwa-0.6.2/bwa index -a bwtsw -p ref reference.fa /opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai -I ref sample2New.fq /opt/bwa-0.6.2/bwa samse -f sample2New.sam -r "@RG\tID:sample2\tPL:ILLUMINA\tPUu1\tLB:sample2\tSM:sample2" ref sample2New.sai sample2New.fq I can see changed base quality string in the sample2New.sam file EBRI093151_0051:4:55:2998:9540#0 0 Chr10 377653 0 101M * 0 0 ACAACACAGTGGGTTGGAGTAGAGCATCTCCAAAGGCCCTTTCCAATCCAACATGAGTAACTCAAGCTCTGCACCAGCCACGAAAAGGCAAGGCTTTGGAT ''''''''''%''#'%&"&&&''''''''$''&''$&&&%%''&&&'&"%%'%'%&&%''&!'$%%%^_"$%'"%% $$&$%#^] !^Q^U^XESC^Z! ##^L#!!!>> RG:Z:sample2 XT:A:R NM:i:0 X0:i:3 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101 XA:Z:Chr10,+33,101M,0;Chr10,+242847,101M,0; and ofcourse the command java -Xmx8g -jar /opt/picard-tools-1.85/SortSam.jar SO=coordinate INPUT=sample2New.sam OUTPUT=sample2New.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true fails with error Exception in thread "main" java.lang.IllegalArgumentException: Invalid fastq character: Why "bwa samse" is changing quality encoding?? Do you have an idea what Im doing wrong? thanks

02-26-2013, 08:35 AM	#2
maubp Peter (Biopython etc) Location: Dundee, Scotland, UK Join Date: Jul 2009 Posts: 1,543	Could you edit your post to use the [ code ] and [ /code ] tags? This is easily done via the advanced editor view where there is a button for this in the tool bar (not shown in the quick reply edit box).

02-26-2013, 08:41 AM	#3
maubp Peter (Biopython etc) Location: Dundee, Scotland, UK Join Date: Jul 2009 Posts: 1,543	I've not checked all the bases (due to the forum formatting), however, it would appear to be down to a FASTQ encoding problem. It appears bwa defaulted to assuming the obsolete Illumina specific ASCII encoding of PHRED+64, while your data was actually the original standard Sanger ASCII encoding of PHRED+33 (now adopted by Illumina). For background, see: http://dx.doi.org/10.1093/nar/gkp1137 In your FASTQ file, the first base has quality code 'F', ASCII character 70. Under the Sanger FASTQ scheme that means 70-33 = quality 37. However, if read in as the obsolete Illumina scheme it would be 70-64 = 6 quality, which when output again in SAM format (which uses the Sanger FASTQ scheme) becomes 6+33 = ASCII 39 = ' (single quote). Solution - there is a command line option to tell bwa you have a Sanger style FASTQ file. Use it, otherwise you get a bad SAM/BAM file.

02-26-2013, 01:12 PM	#4
tomjan Junior Member Location: pl Join Date: Apr 2012 Posts: 8	Thank you for your help Maubp, Your explanations helped me to find the solutions The problem was in "bwa aln" cmmand /opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai -I ref sample2New.fq from the documentation we can se "-I The input is in the Illumina 1.3+ read format (quality equals ASCII-64). ". So, everything is OK when I ommit the -I option. /opt/bwa-0.6.2/bwa aln -t 10 -f sample2New.sai ref sample2New.fq Once again, Thank You for your help.

02-26-2013, 01:23 PM	#5
maubp Peter (Biopython etc) Location: Dundee, Scotland, UK Join Date: Jul 2009 Posts: 1,543	Well done - and thank you for posting back with the details for anyone searching about this again in the future. (I couldn't remember the details about the switch, and wasn't at a machine where I could quickly check - but this way you'll probably remember the problem and solution )