{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# What do quality scores mean?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Phred Quality Scores\n", "We can calculate the *phred quality score* from the probability of sequencing error (i.e. the base call is wrong) using:\n", "\n", "\n", "\\begin{equation*}\n", "Q = -10 \\log_{10} p\n", "\\end{equation*}\n", "\n", "Alternatively, we can rearrange to calculate the probability of error from the *phred quality score* using:\n", "\n", "\\begin{equation*}\n", "p = 10^{\\frac{Q}{-10}}\n", "\\end{equation*}\n", "\n", "Where $Q$ is the *phred quality score* and $p$ is the probability of error (i.e. the base call is wrong)\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Phred Prob of \n", "score Error \n", " 0 1.00000\n", " 1 0.79433\n", " 2 0.63096\n", " 3 0.50119\n", " 4 0.39811\n", " 5 0.31623\n", " 6 0.25119\n", " 7 0.19953\n", " 8 0.15849\n", " 9 0.12589\n", " 10 0.10000\n", " 11 0.07943\n", " 12 0.06310\n", " 13 0.05012\n", " 14 0.03981\n", " 15 0.03162\n", " 16 0.02512\n", " 17 0.01995\n", " 18 0.01585\n", " 19 0.01259\n", " 20 0.01000\n", " 21 0.00794\n", " 22 0.00631\n", " 23 0.00501\n", " 24 0.00398\n", " 25 0.00316\n", " 26 0.00251\n", " 27 0.00200\n", " 28 0.00158\n", " 29 0.00126\n", " 30 0.00100\n", " 31 0.00079\n", " 32 0.00063\n", " 33 0.00050\n", " 34 0.00040\n", " 35 0.00032\n", " 36 0.00025\n", " 37 0.00020\n", " 38 0.00016\n", " 39 0.00013\n", " 40 0.00010\n", " 41 0.00008\n" ] } ], "source": [ "def error_prob(quality):\n", " qval = quality\n", " return 10**(qval/-10.0)\n", "\n", "print (\"{0:^5} {1:^10}\".format(\"Phred\", \"Prob of\"))\n", "print (\"{0:^5} {1:^10}\".format(\"score\", \"Error\"))\n", "for phred in range(0,42):\n", " print (\"{0:^5} {1:03.5f}\".format(phred, error_prob(phred)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ascii Codes\n", "In FASTQ files, phred scores are represented using characters. Each character on the keyboard can be represented by a number, called an ascii code." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Character ASCII # \n", " ! 33 \n", " \" 34 \n", " # 35 \n", " $ 36 \n", " % 37 \n", " & 38 \n", " ' 39 \n", " ( 40 \n", " ) 41 \n", " * 42 \n", " + 43 \n", " , 44 \n", " - 45 \n", " . 46 \n", " / 47 \n", " 0 48 \n", " 1 49 \n", " 2 50 \n", " 3 51 \n", " 4 52 \n", " 5 53 \n", " 6 54 \n", " 7 55 \n", " 8 56 \n", " 9 57 \n", " : 58 \n", " ; 59 \n", " < 60 \n", " = 61 \n", " > 62 \n", " ? 63 \n", " @ 64 \n", " A 65 \n", " B 66 \n", " C 67 \n", " D 68 \n", " E 69 \n", " F 70 \n", " G 71 \n", " H 72 \n", " I 73 \n", " J 74 \n", " K 75 \n", " L 76 \n", " M 77 \n", " N 78 \n", " O 79 \n", " P 80 \n", " Q 81 \n", " R 82 \n", " S 83 \n", " T 84 \n", " U 85 \n", " V 86 \n", " W 87 \n", " X 88 \n", " Y 89 \n" ] } ], "source": [ "print (\"{0:^8} {1:^8}\".format(\"Character\", \"ASCII #\"))\n", "for i in range(33,90):\n", " print(\"{0:^8} {1:^8}\".format(chr(i),i))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Phred Encodings\n", "There are several different ways to encode phred scores with ascii characters. The two most common are called phred+33 and phred+64. The names are strange until you understand how then encoding works. \n", "\n", "### Phred+33\n", "To use the phred+33 encoding, take the phred quality score, add 33 to it, then use the ascii character corresponding to the sum. For example, using the phred+33 encoding, a quality score of 30 would be represented with the ascii character with the ascii code of 63 (30 + 33), which is '?'.\n", "\n", "### Phred+64\n", "The phred+64 encoding works the same as the phred+33 encoding, except you add 64 to the phred score to determine the ascii code of the quality character. You will only find phred+64 encoding on older data, which was sequenced several years ago. The tricky part is that there is no indication in the FASTQ file as to which encoding was used, you have to make an educated guess." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Phred Prob of Phred+33 Phred+64\n", "score Error Ascii Ascii \n", " 0 1.00000 ! @ \n", " 1 0.79433 \" A \n", " 2 0.63096 # B \n", " 3 0.50119 $ C \n", " 4 0.39811 % D \n", " 5 0.31623 & E \n", " 6 0.25119 ' F \n", " 7 0.19953 ( G \n", " 8 0.15849 ) H \n", " 9 0.12589 * I \n", " 10 0.10000 + J \n", " 11 0.07943 , K \n", " 12 0.06310 - L \n", " 13 0.05012 . M \n", " 14 0.03981 / N \n", " 15 0.03162 0 O \n", " 16 0.02512 1 P \n", " 17 0.01995 2 Q \n", " 18 0.01585 3 R \n", " 19 0.01259 4 S \n", " 20 0.01000 5 T \n", " 21 0.00794 6 U \n", " 22 0.00631 7 V \n", " 23 0.00501 8 W \n", " 24 0.00398 9 X \n", " 25 0.00316 : Y \n", " 26 0.00251 ; Z \n", " 27 0.00200 < [ \n", " 28 0.00158 = \\ \n", " 29 0.00126 > ] \n", " 30 0.00100 ? ^ \n", " 31 0.00079 @ _ \n", " 32 0.00063 A ` \n", " 33 0.00050 B a \n", " 34 0.00040 C b \n", " 35 0.00032 D c \n", " 36 0.00025 E d \n", " 37 0.00020 F e \n", " 38 0.00016 G f \n", " 39 0.00013 H g \n", " 40 0.00010 I h \n", " 41 0.00008 J i \n" ] } ], "source": [ "print (\"{0:^5} {1:^8} {2:^8} {3:^8}\".format(\"Phred\", \"Prob of\", \"Phred+33\", \"Phred+64\"))\n", "print (\"{0:^5} {1:^8} {2:^8} {3:^8}\".format(\"score\", \"Error\", \"Ascii\", \"Ascii\"))\n", "for phred in range(0,42):\n", " print (\"{0:^5} {1:03.5f} {2:^8} {3:^8}\".format(phred, error_prob(phred), chr(phred+33), chr(phred+64)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Why +33?\n", "ASCII 33 is the first \"normal\" ASCII character that. 1 through 32 include whitespace and non-printing characters, which cannot be identified by eye)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Character ASCII # \n", " \u0000 0 \n", " \u0001 1 \n", " \u0002 2 \n", " \u0003 3 \n", " \u0004 4 \n", " \u0005 5 \n", " \u0006 6 \n", " \u0007 7 \n", " \b 8 \n", " \t 9 \n", " \n", " 10 \n", " \u000b", " 11 \n", " \f", " 12 \n", " \r", " 13 \n", " \u000e 14 \n", " \u000f 15 \n", " \u0010 16 \n", " \u0011 17 \n", " \u0012 18 \n", " \u0013 19 \n", " \u0014 20 \n", " \u0015 21 \n", " \u0016 22 \n", " \u0017 23 \n", " \u0018 24 \n", " \u0019 25 \n", " \u001a 26 \n", " \u001b 27 \n", " \u001c", " 28 \n", " \u001d", " 29 \n", " \u001e", " 30 \n", " \u001f 31 \n", " 32 \n", " ! 33 \n", " \" 34 \n", " # 35 \n", " $ 36 \n", " % 37 \n", " & 38 \n", " ' 39 \n" ] } ], "source": [ "print (\"{0:^8} {1:^8}\".format(\"Character\", \"ASCII #\"))\n", "for i in range(0,40):\n", " print(\"{0:^8} {1:^8}\".format(chr(i),i))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.3" } }, "nbformat": 4, "nbformat_minor": 1 }