3강 / 20강

품질 점수 — 글자 하나가 확률 하나

Phred 점수와 ASCII · 그리고 이 자료의 품질이 여섯 가지뿐인 이유

FASTQ 의 네 번째 줄은 서열과 글자 수가 똑같습니다. 서열의 i번째 염기에 대해, 품질 줄의 i번째 글자가 "이 염기가 틀렸을 확률"을 담고 있기 때문입니다.

Phred 점수 Q = −10 · log₁₀ P
Q 20 → 100번에 1번 틀림 · Q 30 → 1,000번에 1번 · Q 40 → 10,000번에 1번

그 숫자를 파일에 넣을 때는 ASCII 문자 한 글자로 바꿉니다. Illumina 는 Q + 33 번째 문자를 씁니다 (Phred+33).

문자 ! " # $ … 5 … A … E … I Q 0 1 2 3 … 20 … 32 … 36 … 40

그래서 # 은 Q2, E 는 Q36 입니다.

고등학생글자 하나가 확률 하나

FASTQ 4번 줄의 이상한 문자들은 암호가 아니라 숫자 다. 문자를 ASCII 코드로 바꾸고 33을 빼면 Phred 점수가 나온다.

Phred 점수가 뜻하는 것

Q	틀릴 확률	맞을 확률
10	1/10	90%
20	1/100	99%
30	1/1,000	99.9%
40	1/10,000	99.99%

Q 가 10 오를 때마다 오류 확률이 10분의 1 이 된다. 지진의 규모나 소리의 데시벨과 같은 로그 눈금이다. Q30 과 Q40 의 차이는 Q10 과 Q20 의 차이와 같은 "한 칸" 이지만 절대적인 오류 개수로는 천 배 차이다.

Q = -10 log₁₀(오류 확률). 기계가 자기 확신을 숫자로 적어 준다

문자 → 숫자 → 확률. ord(c) - 33 한 줄이 전부입니다.

Python 브라우저에서 실행

import bu

h, s, q = bu.reads("WT_1", 1)[0]
scores = bu.phred(q)

print("서열 :", s[:40])
print("품질 :", q[:40])
print()
print(f"{'위치':>4}{'염기':>5}{'문자':>5}{'Q':>5}{'오류 확률':>12}")
for i in range(10):
    print(f"{i+1:>4}{s[i]:>5}{q[i]:>5}{scores[i]:>5}{bu.err_rate(scores[i]):>12.5f}")
print()
print(f"이 read 의 평균 Q : {sum(scores)/len(scores):.1f}")

고등학생왜 하필 ASCII 로 적나

Q 점수를 그냥 숫자로 적으면 되지 않나? "30 28 35 40 …" 이렇게. 안 된다. 자리를 너무 많이 먹는다.

방식	76글자 read 하나당	2,500만 read
숫자 + 공백	약 228 바이트	약 5.7 GB
ASCII 한 글자	76 바이트	약 1.9 GB

세 배 차이다. 그리고 서열과 품질의 길이가 정확히 같아진다 는 이점이 더 크다. i 번째 글자의 품질은 i 번째 문자. 찾을 필요 없이 위치가 곧 대응이다.

33을 빼는 이유

ASCII 0~32 는 줄바꿈 같은 보이지 않는 제어 문자 다. 파일에 넣으면 파일이 깨진다. 그래서 33부터 시작한다. Q0 = 문자 '!', Q40 = 문자 'I'. 오프셋 33 이라고 부른다.

형식의 설계에는 언제나 이유가 있다 — 대개는 용량과 편의

직접 세어 보면 이 자료의 정체가 드러납니다.

Python 브라우저에서 실행

import bu
import numpy as np
from collections import Counter

# 이 자료에 실제로 나오는 품질 점수는 몇 가지일까?
c = Counter()
for _, _, q in bu.reads("WT_1"):
    c.update(bu.phred(q))

tot = sum(c.values())
print(f"품질 점수의 종류 : {len(c)} 가지")
print()
print(f"{'Q':>4}{'문자':>5}{'염기 수':>12}{'비율':>9}{'오류 확률':>12}")
for qv in sorted(c):
    print(f"{qv:>4}{chr(qv+33):>5}{c[qv]:>12,}{c[qv]/tot*100:>8.2f}%{bu.err_rate(qv):>12.4f}")

여섯 가지뿐입니다

Q2 · Q14 · Q21 · Q27 · Q32 · Q36. 그 사이 값은 하나도 없습니다.

이것은 고장이 아니라 설계입니다. Illumina 는 2014년 무렵부터 품질 점수 뭉치기(quality score binning)를 씁니다 — 40가지 점수를 8개 칸으로 몰아넣는 것이죠. 파일 크기가 크게 줄고, downstream 결과는 거의 안 바뀐다는 것이 그들의 주장입니다.

얻는 것

품질 줄이 잘 압축됩니다. 같은 글자가 반복되니까요. 이 자료도 78%가 Q36 한 값입니다.

잃는 것

Q36 이라고 적힌 염기들의 실제 품질은 Q34~Q40 으로 제각각이었을 텐데, 그 구분이 영영 사라졌습니다.

그리고 Q2 는 특별합니다. Illumina 는 "이 염기는 아예 믿지 말라"는 뜻으로 Q2(#)를 씁니다. 1강에서 본 read 끝의 ##### 가 그것이었죠.

고등학생품질이 여섯 가지밖에 없다

이 자료에서 실제로 나타나는 품질 점수를 세어 보면 여섯 가지뿐 이다. 원래는 0부터 40까지 마흔한 가지가 나와야 한다.

기계가 반올림을 했다

일루미나는 어느 시점부터 품질 점수를 몇 개 구간으로 뭉뚱그리는 방식을 쓴다. 정확히 Q34 인지 Q35 인지는 하류 분석에서 거의 쓸모가 없는데, 저장 공간은 크게 잡아먹기 때문이다. 파일 크기를 크게 줄이려고 일부러 정보를 버린 것이다.

사진을 JPEG 으로 저장하는 것과 같다. 눈에 안 띄는 정보를 버려서 용량을 줄인다. 되돌릴 수 없지만, 대부분의 용도에는 지장이 없다.

언제 문제가 되나

품질 점수를 정밀하게 쓰는 분석에서는 문제가 된다. 변이 검출(variant calling)에서 염기 품질을 재보정하는 절차가 그 예다. RNA-seq 카운트에는 거의 영향이 없다.

"값이 여섯 개뿐" 은 자료가 나쁜 게 아니라 압축된 것이다

오른쪽 그림은 자료가 아니라 예시입니다 — 뭉치기 전 모양을 짐작해 본 것뿐이에요.

Python 브라우저에서 실행

import bu
import matplotlib.pyplot as plt
import numpy as np
from collections import Counter

c = Counter()
for _, _, q in bu.reads("WT_1"):
    c.update(bu.phred(q))

qs = sorted(c)
fig, ax = plt.subplots(1, 2, figsize=(10, 3.8))
ax[0].bar([str(q) for q in qs], [c[q] for q in qs], color="#4c6ef5")
ax[0].set_title("이 자료 (NextSeq · 뭉친 점수)"); ax[0].set_xlabel("Q")

# 뭉치기 전이라면 이렇게 이어져 있었을 것이다 (모양만 보이는 예시)
x = np.arange(2, 41)
y = np.exp(-((x - 36) ** 2) / 18) * c[36]
ax[1].bar(x, y, color="#adb5bd")
ax[1].set_title("뭉치지 않았다면 (예시)"); ax[1].set_xlabel("Q")
plt.tight_layout()

print(f"실제 종류 {len(qs)}가지 — 오른쪽은 '이랬을 것이다' 를 그려 본 예시입니다.")

🎮 품질 변환기

문자 · Q · 오류 확률 사이를 오가며 감을 잡으세요.

고등학생평균 Q 를 그냥 평균 내면 안 되는 이유

read 하나의 품질을 한 숫자로 요약하고 싶다. Q 값들을 산술평균하면 될까? 안 된다.

Q 는 로그 눈금이다

Q40 과 Q20 의 평균을 산술적으로 내면 Q30 이다. 그런데 실제 오류 확률로 보면 0.0001 과 0.01 의 평균은 0.00505, 이것을 Q 로 바꾸면 약 Q23 이다. Q30 이 아니다.

	계산	결과
틀린 방법	(40+20)/2	Q30 — 실제보다 좋게 보임
맞는 방법	확률로 바꿔 평균 → 다시 Q	Q23

일반화하면

로그 눈금 위의 값은 더하고 나눌 수 없다. pH, 데시벨, 지진 규모, 그리고 Phred 점수가 모두 그렇다. 평균을 내려면 원래 눈금(확률, 농도, 에너지)으로 되돌아가야 한다. 나쁜 쪽이 평균을 지배한다는 점도 기억할 만하다 — Q40 이 아무리 많아도 Q2 몇 개가 read 전체를 망친다.

단위를 확인하지 않은 평균은 계산이 아니라 사고다

🤔 확인 문제

품질 문자 # 의 Phred 점수와 뜻은?

해설 # 은 ASCII 35 이므로 Q = 35 − 33 = 2 입니다. 오류 확률 63% — 사실상 무작위죠. Illumina 는 이것을 read segment quality control indicator 로 씁니다.

🤔 확인 문제

이 자료의 품질 점수가 여섯 가지뿐인 이유는?

해설 quality score binning 입니다. NextSeq·NovaSeq 세대의 기본 동작이고, 2강에서 본 NS500424(NextSeq) 와 아귀가 맞습니다 — 메타데이터의 HiSeq 2000 이었다면 이런 모양이 안 나옵니다.

🤔 확인 문제

Phred 점수를 산술 평균 내는 것이 엄밀히 틀린 이유는?