DNA에도 기억이 새겨집니다

S38(ChIP-seq), S39(ATAC-seq), S40(Hi-C)으로 히스톤 변형, 크로마틴 접근성, 3D 구조를 읽었습니다. 후성유전체의 마지막 층위가 하나 남았습니다 — **DNA 자체에 새겨진 화학적 표식, 메틸화(methylation)**입니다.

DNA 메틸화는 시토신(C) 염기의 5번 탄소에 메틸기(-CH₃)가 공유결합으로 추가되는 것입니다. 주로 CpG 디뉴클레오타이드(5'-CG-3')에서 일어나며, 프로모터 CpG 섬(CpG island)의 메틸화는 일반적으로 유전자 침묵과 연관됩니다. 발달 과정에서 세포 유형별 유전자 ON/OFF 패턴이 메틸화로 "잠금"되고, 이 패턴은 세포 분열을 거쳐도 유지됩니다 — 그래서 "화학적 기억"이라 부릅니다.

이 편에서는 메틸화를 단일 염기 해상도로 읽는 바이설파이트 시퀀싱의 화학적 원리, Bismark 정렬 도구의 독특한 전략, 그리고 차등 메틸화 영역(DMR) 검출까지를 다룹니다.

바이설파이트 전환의 화학

바이설파이트(sodium bisulfite, NaHSO₃) 처리의 핵심은 비메틸화 시토신(C)을 우라실(U)로 전환하는 것입니다. 메틸화된 시토신(mC)은 메틸기의 보호로 전환되지 않고 C로 남습니다.

시퀀싱하면 우라실은 티민(T)으로 읽히므로:

• C → T: 이 위치의 시토신은 메틸화되지 않았음
• C → C: 이 위치의 시토신은 메틸화되어 있음

이것이 단일 염기 해상도로 메틸화를 읽는 원리입니다.

전환 효율과 불완전 전환

이론적으로는 비메틸화 C가 100% U로 전환되어야 하지만, 실제 전환 효율은 99～99.5% 정도입니다. 불완전 전환은 메틸화를 과대 추정하게 합니다. 람다 파지 DNA(완전 비메틸화)를 스파이크인으로 첨가해 전환 효율을 측정하는 것이 표준 관행입니다.

WGBS vs RRBS: 커버리지와 비용의 트레이드오프

• WGBS (Whole-Genome Bisulfite Sequencing): 전체 게놈의 모든 C를 읽습니다. 가장 포괄적이지만, 30× 커버리지를 위해 ～10억 리드가 필요합니다.
• RRBS (Reduced Representation Bisulfite Sequencing): MspI 같은 제한효소로 CpG 밀집 영역만 잘라내어 시퀀싱합니다. 게놈의 ～1～3%만 읽지만, CpG 섬 대부분을 커버하며 비용이 훨씬 적습니다.

Bismark: 3문자 알파벳 정렬

바이설파이트 처리된 리드를 일반 정렬기(BWA, Bowtie2)로 정렬하면 문제가 생깁니다 — 원래 C였던 염기가 T로 바뀌어 있어, 참조 게놈과 불일치(mismatch)가 대거 발생합니다. Bismark은 이 문제를 in silico 전환으로 해결합니다.

1. 참조 게놈을 두 가지로 준비: (a) 모든 C→T 전환, (b) 모든 G→A 전환 (역상보 가닥 처리)
2. 리드도 in silico 전환: C→T 또는 G→A
3. 전환된 리드를 전환된 참조에 정렬 (Bowtie2 백엔드)
4. 원래 서열로 되돌려 각 C 위치의 **전환 여부(메틸화 여부)**를 판정

이 전략은 사실상 DNA 알파벳을 **4문자(ACGT)에서 3문자(A/T/G 또는 A/C/T)**로 축소한 뒤 정렬하는 것입니다. 정보 손실이 있으므로 매핑률이 일반 정렬보다 낮은 것이 정상입니다(60～80%).

# 1단계: Bismark 게놈 인덱스 준비
bismark_genome_preparation --bowtie2 genome/

# 2단계: 정렬 (페어드엔드)
bismark --bowtie2 --genome genome/ \
  -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz \
  -o aligned/

# 3단계: 메틸화 추출
bismark_methylation_extractor --paired-end --comprehensive \
  --bedGraph --CX \
  aligned/reads_R1_bismark_bt2_pe.bam

# 4단계: 커버리지 보고서
bismark2report --alignment_report aligned/reads_R1_bismark_bt2_PE_report.txt

출력은 각 CpG 위치별 메틸화 리드 수 / 총 리드 수 = 메틸화 비율(β-value, 0～1)입니다.

차등 메틸화 영역(DMR) 검출

두 조건(예: 종양 vs 정상, 유형 A 세포 vs 유형 B 세포) 사이에서 메틸화 수준이 유의하게 다른 영역을 **DMR(Differentially Methylated Region)**이라 합니다.

methylKit (R)

library(methylKit)

# 메틸화 데이터 로드 (Bismark 출력 또는 bedGraph)
myobj <- methRead(
  list("tumor.txt", "normal.txt"),
  sample.id = list("tumor", "normal"),
  treatment = c(1, 0),
  assembly = "hg38",
  context = "CpG"
)

# 필터링: 최소 10× 커버리지
filtered <- filterByCoverage(myobj, lo.count = 10, hi.perc = 99.9)

# 타일링: 1000bp 창으로 집계 (개별 CpG 대신 영역 단위)
tiles <- tileMethylCounts(filtered, win.size = 1000, step.size = 1000)

# 차등 메틸화 검정
meth_diff <- calculateDiffMeth(tiles)

# DMR 추출 (25% 이상 차이, q < 0.01)
dmr <- getMethylDiff(meth_diff, difference = 25, qvalue = 0.01)

DSS/bsseq (Bioconductor)

DSS(Dispersion Shrinkage for Sequencing data)는 베이즈 접근으로, CpG 단위의 메틸화 비율을 베타-이항(beta-binomial) 모델로 추정합니다. 이웃 CpG 사이의 공간적 상관을 이용해 커버리지가 낮은 위치의 추정을 안정화합니다.

library(DSS)
library(bsseq)

# BSseq 객체 생성
bs <- BSseq(chr = chr_vec, pos = pos_vec,
            M = cbind(tumor_M, normal_M),
            Cov = cbind(tumor_cov, normal_cov),
            sampleNames = c("tumor", "normal"))

# 차등 메틸화 검정
dml_test <- DMLtest(bs, group1 = "tumor", group2 = "normal")

# DMR 검출 (연속된 유의 CpG 병합)
dmr <- callDMR(dml_test, p.threshold = 0.001)

메틸화의 생물학적 의미

• 프로모터 CpG 섬 메틸화 → 유전자 침묵: 메틸화된 CpG에 MBD(메틸-CpG 결합 도메인) 단백질이 결합해 전사 기구의 접근을 차단합니다.
• 유전자 본체(gene body) 메틸화 → 전사 활성: 역설적이지만, 유전자 본체의 메틸화는 오히려 활발히 전사되는 유전자에서 높습니다.
• X 불활성화: 여성의 두 X 염색체 중 하나가 메틸화로 침묵됩니다.
• 게놈 각인(imprinting): 부모 유래에 따라 메틸화 패턴이 달라, 한쪽 부모의 대립유전자만 발현됩니다.
• 암: 프로모터 과메틸화(종양 억제 유전자 침묵) + 전 게놈 저메틸화(게놈 불안정성)가 암의 후성유전적 특성입니다.

CS 매핑

• 알파벳 확장/축소: 바이설파이트 전환은 DNA의 4문자 알파벳에서 C/T 구분을 제거해 3문자로 축소합니다. 이것은 문자열 알고리즘에서 알파벳 크기가 정렬 복잡도와 정확도에 미치는 영향의 실증입니다.
• 비대칭 정렬: Bismark의 in silico 전환 + 이중 참조 전략은 일반 정렬의 변형으로, 정보 이론적으로 매핑 불확실성이 증가하는 대가를 측정 가능한 형태로 수용합니다.
• 베타 분포 / 베이즈 추정: 메틸화 비율(성공/시행)의 불확실성을 베타 분포로 모델링하고, 이웃 정보로 사전 분포를 업데이트하는 것은 베이즈 추론의 교과서적 응용입니다.

자주 만나는 결함

• 불완전 전환 미측정 — 람다 파지 스파이크인 없이 전환 효율을 가정하면 메틸화를 체계적으로 과대 추정합니다. 반드시 측정합시다.
• 낮은 커버리지에서 메틸화 비율 보고 — CpG 하나에 리드 2개(1개 메틸화)면 β = 0.5이지만, 신뢰도는 극히 낮습니다. 최소 5～10× 커버리지 필터를 적용합시다.
• CpG만 분석 — 포유류에서는 CpG 메틸화가 주류이지만, 식물과 초기 발달에서는 CHG, CHH 맥락(H = A/T/C)의 메틸화도 중요합니다. --CX 옵션으로 모든 맥락을 추출합시다.
• 산화적 메틸시토신(5hmC) 미구분 — 표준 바이설파이트는 5mC와 5hmC를 구분하지 못합니다(둘 다 전환 안 됨). 5hmC가 중요한 경우 산화적 바이설파이트(oxBS-seq) 또는 TAB-seq를 사용합시다.
• PCR 편향: 바이설파이트 처리로 DNA가 손상되어 GC 편향이 심해집니다. 중복 제거와 GC 편향 보정이 필요합니다.

더 깊게 파고 싶다면

본문은 BPD가 직접 재구성한 서술입니다. 심화는 아래 정통 자료를 활용합시다.

• 바이설파이트 시퀀싱 리뷰: Krueger et al. (2012), DNA methylome analysis using short bisulfite sequencing data, Nature Methods 9:145.
• Bismark 공식 문서: Krueger & Andrews (2011), Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications, Bioinformatics 27:1571.
• methylKit: Akalin et al. (2012), methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles, Genome Biology 13:R87.
• DSS: Park & Wu (2016), Differential methylation analysis for BS-seq data under general experimental design, Bioinformatics 32:1446.
• Harvard STAT115 후성유전학 강의: Xiaole Shirley Liu 교수 — DNA 메틸화의 통계적 분석을 다룹니다.
• ENCODE WGBS 데이터: encodeproject.org — 무료 공개 WGBS 데이터셋.

이로써 후성유전체 시리즈 — ChIP-seq(S38), ATAC-seq(S39), Hi-C(S40), 바이설파이트 시퀀싱(S41) — 을 마쳤습니다. 유전자 발현을 조절하는 네 가지 층위를 모두 읽을 수 있게 되었습니다. 다음 섹션에서는 시선을 개별 유전자에서 집단으로 돌려, GWAS와 집단 유전학으로 넘어갑니다.