2
$\begingroup$

I have a fasta file with some DNA sequences. I would like to simulate next-generation sequencing reads from it. I'm doing it without any base error and mutation error.

wgsim -e 0 -r 0 sequence.fa seq_0_1.fq seq_0_2.fq

To my knowledge, this is a perfect simulation. Next, I give the paired alignments to bwa for alignment.

bwa mem -M sequence.fa seq_0_1.fq seq_0_2.fq > P0.sam 

Now, I check the ASCII of base quality (column 11 in the SAM format, the specification is here).

head -n 3 P0.sam | cut -f11

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

According to this page, the ASCII (in order of quality) is:

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~

Question: My experiment is supposed to be perfect, so my base-quality metrics are expected to appear near the right-end of the list (like x,y,z). However, my metric I is no-where near the top of the list. In particular, I'm unable to achieve anything from J to ~. Why?

$\endgroup$
4
$\begingroup$

When fasta reads are aligned they are by default assigned the phred score of 40, which in phred+33 encoding is represented by I. Phred+33 uses ASCII characters from ! to I:

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI

Phred+33 was originally used by sanger machines but it is now used by all popular platforms such as Illumina, Ion-torrent (also Ion-proton) and Roche 454. Since most platforms now use phred+33, most aligners also assume it by default. You can specify other quality encoding formats if you want (such as phred+64 or Solexa).

Have a look at wikipedia page on FASTQ format for the different quality encoding formats.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.