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$

1 Answer 1

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$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .