Lesson 15 of 20

읽기를 유전자에 붙이고 세기

정렬할 것인가, 세기만 할 것인가

이제 FASTQ 를 유전자별 숫자로 바꿉니다. 길이 두 갈래인데, 나오는 표는 비슷해요.

① 정렬 — STAR

읽기를 유전체의 정확한 자리에 붙인다
BAM 파일이 나온다
메모리 32 GB · 시료당 20분
스플라이싱·변이·커버리지를 볼 수 있다

② 세기 — Salmon · kallisto

자리는 안 찾고 어느 전사체인지만 맞힌다
카운트 표가 바로 나온다
메모리 4 GB · 시료당 3분
발현량만 볼 거라면 충분하다

차등발현만 볼 거라면 Salmon 이 답입니다. 그림에 커버리지를 그리거나 변이를 볼 거라면 STAR 가 필요해요.

High school정렬할 것인가, 세기만 할 것인가

	정렬 (STAR · HISAT2)	유사정렬 (salmon · kallisto)
하는 일	유전체 어디에 붙는지 정확히	어느 전사체에서 왔을지 추정
속도	느림	훨씬 빠름
메모리	수십 GB	적음
결과물	BAM — 위치 정보	카운트만
변이·스플라이싱	가능	불가

발현만 필요하면 유사정렬이 훨씬 효율적입니다.

High school스플라이싱 인식 정렬

RNA read는 엑손 경계를 넘습니다. 유전체에 붙이려면 중간에 인트론을 건너뛰어야 해요.

일반 정렬 도구	스플라이싱 인식 도구
연속된 구간만 찾음	read를 쪼개서 붙임
경계 read를 놓침	접합부를 검출

STAR·HISAT2가 그 일을 합니다 — DNA용 도구(BWA 등)는 RNA에 부적합해요.

읽기 길이를 늘리면 오히려 DNA 정렬기로는 더 못 씁니다.

Python Runs in your browser

# 엑손 경계에 걸치는 읽기가 얼마나 되나
exon_len = 145.0          # 사람 엑손 중앙 길이 (bp)
read_len = 100.0

# 읽기 시작 자리가 엑손 안에 고르게 흩어진다면
inside = max(exon_len - read_len, 0)
p_span = 1 - inside / exon_len
print("엑손 길이 %.0f bp, 읽기 길이 %.0f bp" % (exon_len, read_len))
print("경계에 걸치는 읽기의 비율: %.0f%%" % (p_span * 100))
print()
for rl in [50, 75, 100, 150]:
    inside = max(exon_len - rl, 0)
    p = 1 - inside / exon_len
    print("읽기 %3d bp -> 경계에 걸침 %.0f%%" % (rl, p * 100))
print()
print("읽기가 길수록 더 많이 걸칩니다.")
print("DNA 정렬기를 쓰면 이만큼을 통째로 버립니다.")

실제 명령 두 줄

# STAR — 유전체에 붙인다
STAR --genomeDir /ref/GRCh38_star \
     --readFilesIn R1.trim.gz R2.trim.gz --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode GeneCounts \
     --outFileNamePrefix Ctrl_1_

# Salmon — 전사체를 맞힌다
salmon quant -i /ref/GRCh38_salmon -l A \
     -1 R1.trim.gz -2 R2.trim.gz \
     --validateMappings --gcBias \
     -o Ctrl_1_quant
    

깃발	무엇
`--quantMode GeneCounts`	STAR 가 카운트까지 함께 낸다
`-l A`	라이브러리 방향을 알아서 알아내라 (10강의 dUTP)
`--gcBias`	GC 함량에 따른 치우침을 보정한다

High school가닥 정보 쓰기

10강에서 dUTP로 가닥을 기억시켰습니다. 이제 그 정보를 씁니다.

설정	내용
unstranded	가닥 무시 — 정보 낭비
stranded (정방향)	read가 mRNA와 같은 방향
reverse (역방향)	read가 반대 방향 — dUTP 방식은 이것

설정을 틀리면 정량이 크게 틀립니다 — 겹치는 유전자에서 특히요.

확신이 없으면 도구로 자동 판별할 수 있어요(infer_experiment 등). 반드시 확인하세요.

여기서 나온 표가 17강부터 계속 쓰입니다.

Python Runs in your browser

# 카운트 표가 만들어지는 마지막 단계
reads = [
    ("read1", "chr1", 1000, "+"), ("read2", "chr1", 1040, "+"),
    ("read3", "chr1", 5200, "-"), ("read4", "chr2",  300, "+"),
    ("read5", "chr1", 1010, "+"), ("read6", "chr1", 9999, "+"),
]
genes = [
    ("ACTB",  "chr1",  900, 1200, "+"),
    ("CD3D",  "chr1", 5100, 5400, "-"),
    ("GAPDH", "chr2",  200,  500, "+"),
]

counts = {g[0]: 0 for g in genes}
unassigned = 0
for rid, chrom, pos, strand in reads:
    hit = None
    for name, gc, s, e, gs in genes:
        if chrom == gc and s <= pos <= e and strand == gs:
            hit = name
    if hit:
        counts[hit] += 1
    else:
        unassigned += 1
    print("%-6s %s:%-6d %s  ->  %s" % (rid, chrom, pos, strand, hit or "배정 안 됨"))
print()
print("카운트 표:")
for g, c in counts.items():
    print("  %-8s %d" % (g, c))
print("  배정 안 됨: %d" % unassigned)
print()
print("이것이 featureCounts 가 하는 일의 전부입니다 — 자리와 방향을 보고 세기.")

High school다중 정렬 read

어떤 read는 여러 곳에 똑같이 잘 붙습니다. 어떻게 셀까요?

처리	내용	영향
버림	고유 정렬만	패밀리 유전자가 과소평가
나눠 줌	확률적으로 배분 (EM)	유사정렬의 기본 방식
무작위로 하나	—	재현성 문제

16강에서 이 문제가 돌아옵니다 — 유사 유전자가 많은 집안은 정량이 어렵습니다.

🤔 Check yourself

RNA-seq 에 DNA 정렬기(BWA)를 쓰면 안 되는 이유는?

🤔 Check yourself

차등발현만 볼 계획이다. 무엇을 쓰는 것이 나은가?