SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 02-25-2013, 01:45 PM   #1
tomjan
Junior Member
 
Location: pl

Join Date: Apr 2012
Posts: 8
Default 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
tomjan is offline   Reply With Quote
Old 02-26-2013, 07:35 AM   #2
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

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 is offline   Reply With Quote
Old 02-26-2013, 07:41 AM   #3
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

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.
maubp is offline   Reply With Quote
Old 02-26-2013, 12:12 PM   #4
tomjan
Junior Member
 
Location: pl

Join Date: Apr 2012
Posts: 8
Default

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.
tomjan is offline   Reply With Quote
Old 02-26-2013, 12:23 PM   #5
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

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 )
maubp is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:01 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO