부록 E13: RNA-seq 분석 파이프라인

7 분 소요

이 장에서 배울 것

이번 장에서는 RNA-seq 분석 파이프라인(RNA-seq analysis pipeline)을 배웁니다. RNA-seq은 세포 안에 어떤 RNA가 얼마나 있는지 측정하는 실험입니다. 프로그래밍 관점에서 보면 RNA-seq 분석은 “실험 기계가 만든 원시 파일을 읽고, 품질을 확인하고, 유전자별 숫자 표로 바꾸고, 두 조건을 비교하는 과정”입니다.

핵심 용어를 먼저 정리하겠습니다.

RNA-seq: RNA의 양을 대량으로 측정하는 시퀀싱 실험입니다.
read: 시퀀싱 기계가 읽어 낸 짧은 서열 조각입니다.
FASTQ: read 서열과 품질 점수를 함께 저장하는 파일 형식입니다.
품질관리(QC, quality control): 데이터가 분석할 만한 상태인지 확인하는 단계입니다. 앞으로는 필요한 곳에서 QC라고도 부르겠습니다.
trimming: 품질이 낮은 끝부분이나 어댑터(adapter, 실험 과정에서 붙는 짧은 인공 서열)를 잘라내는 작업입니다.
정렬(alignment): read가 기준 유전체나 전사체의 어느 위치와 맞는지 찾는 작업입니다.
정량(quantification): 각 유전자나 전사체에 read가 얼마나 배정되는지 세는 작업입니다.
카운트 행렬(count matrix): 행은 유전자, 열은 샘플, 값은 read 수인 숫자 표입니다.
정규화(normalization): 샘플마다 읽힌 전체 양이 다른 문제를 보정하는 작업입니다.
차등발현분석(differential expression analysis): 두 조건 사이에서 발현량이 달라진 유전자를 찾는 분석입니다.

RNA-seq 분석 파이프라인

가장 쉬운 비유: 도서관 책 조각 맞추기

RNA-seq 데이터를 도서관에서 찢어진 책 조각을 많이 받은 상황으로 생각해 봅시다. 각 조각에는 짧은 문장이 적혀 있습니다. 먼저 조각이 너무 흐리거나 찢어졌는지 확인합니다. 그다음 조각을 원래 책의 어느 페이지에 붙일 수 있는지 찾습니다. 마지막으로 “어떤 책에서 조각이 많이 나왔는지”를 세면, 사람들이 어떤 책을 많이 읽었는지 짐작할 수 있습니다.

RNA-seq도 비슷합니다. read는 RNA에서 나온 짧은 조각입니다. read를 유전자에 연결하고, 유전자별 read 수를 세면 어떤 유전자가 많이 발현되었는지 알 수 있습니다.

전체 흐름 한눈에 보기

RNA-seq 파이프라인은 보통 다음 흐름으로 진행됩니다.

FASTQ 파일
→ 품질관리(QC)
→ 필요하면 trimming
→ 정렬 또는 준정렬
→ 유전자별 read 수 계산
→ count matrix 만들기
→ 정규화
→ 차등발현분석
→ 결과 해석과 그림 만들기

여기서 준정렬(pseudoalignment)은 read를 기준 서열에 정확한 위치까지 붙이는 대신, 어떤 전사체에서 왔는지를 빠르게 추정하는 방식입니다. 입문 단계에서는 “정렬보다 빠른 정량 방식도 있다” 정도로 이해하면 충분합니다.

1단계: FASTQ 파일 확인

FASTQ 파일은 대략 다음처럼 생겼습니다.

@read001
ACGUACGUACGU
+
IIIIIIIIIIII

첫 줄은 read 이름, 둘째 줄은 서열, 셋째 줄은 구분자, 넷째 줄은 품질 점수입니다. 실제 파일은 DNA 알파벳 A, T, G, C로 저장되는 경우가 많습니다. RNA에서 온 정보라도 실험과 파일 처리 과정에서는 DNA 문자처럼 다루는 일이 흔합니다.

2단계: 품질관리

품질관리는 데이터를 분석하기 전에 “이 파일이 괜찮은가?”를 보는 단계입니다. 예를 들어 read 끝부분의 품질이 낮거나, 어댑터 서열이 너무 많이 섞여 있으면 이후 분석 결과가 흔들릴 수 있습니다.

터미널에서는 보통 이런 식으로 도구를 실행합니다.

fastqc sample1.fastq.gz

이 명령은 sample1.fastq.gz 파일의 품질 보고서를 만듭니다. 여기서 중요한 것은 명령어를 외우는 것이 아니라, “원시 데이터를 바로 믿지 말고 먼저 검사한다”는 태도입니다.

3단계: trimming

품질이 낮은 끝부분이나 어댑터가 있으면 잘라낼 수 있습니다. trimming은 데이터를 예쁘게 꾸미는 작업이 아니라, 잘못된 서열 조각이 정렬이나 정량을 방해하지 않도록 줄이는 작업입니다.

다만 무조건 많이 자르는 것이 좋은 것은 아닙니다. 너무 많이 자르면 read가 너무 짧아져서 어디서 온 조각인지 알기 어려워집니다. 따라서 trimming은 품질 보고서를 보고 필요한 만큼만 해야 합니다.

4단계: 정렬 또는 정량

read를 기준 유전체(reference genome)나 기준 전사체(reference transcriptome)에 연결해야 합니다. 정렬 도구로는 STAR, HISAT2 같은 도구가 자주 등장합니다.

아주 단순화한 명령어 모양은 다음과 같습니다.

hisat2 -x genome_index -U sample1.fastq.gz -S sample1.sam

뜻은 “sample1.fastq.gz의 read를 genome_index라는 기준에 맞추고, 결과를 sample1.sam에 저장하라”입니다. 실제 연구에서는 paired-end read, 스레드 수, annotation 파일 등 더 많은 옵션을 씁니다.

정렬을 거치지 않고 빠르게 정량하는 도구도 있습니다. Salmon, kallisto 같은 도구가 대표적입니다. 이런 도구는 전사체별 양을 빠르게 추정할 때 많이 쓰입니다.

5단계: count matrix 만들기

분석의 핵심 결과 중 하나는 유전자별 read 수 표입니다.

gene     sample1  sample2  sample3
BRCA1    120      140      30
TP53     80       75       200
EGFR     10       12       90

이 표를 카운트 행렬이라고 부릅니다. 이제부터는 생물학 문제가 숫자 표 문제로 바뀝니다. 어떤 유전자가 어느 조건에서 많이 나왔는지, 차이가 우연인지 아닌지 계산할 수 있습니다.

6단계: 정규화

샘플마다 전체 read 수가 다르면 단순 count만 비교하면 위험합니다. 예를 들어 sample1은 전체 1,000만 read가 읽혔고 sample2는 전체 2,000만 read가 읽혔다면, sample2의 count가 전체적으로 커 보일 수 있습니다. 이 차이는 생물학적 차이가 아니라 실험에서 읽힌 양의 차이일 수 있습니다.

정규화는 이런 차이를 보정하는 작업입니다. 입문 단계에서는 다음 정도만 기억하면 됩니다.

정규화 전: read 수 자체
정규화 후: 샘플 간 비교가 조금 더 공정해진 값

7단계: 차등발현분석

차등발현분석은 두 조건 사이에서 발현량이 달라진 유전자를 찾는 작업입니다. 예를 들어 정상 세포와 암세포를 비교해 어떤 유전자가 더 많이 발현되는지 볼 수 있습니다.

R에서는 DESeq2 같은 패키지가 자주 쓰입니다.

library(DESeq2)
# count matrix와 샘플 정보를 넣어 조건별 차이를 분석합니다.

여기서 중요한 것은 “p-value가 작다”만 보는 것이 아닙니다. 여러 유전자를 동시에 검사하므로 다중검정 보정이 필요하고, 발현 차이가 통계적으로 유의한지뿐 아니라 생물학적으로도 의미 있는지 봐야 합니다.

작은 코드 감각: count matrix 읽기

파이썬에서는 결과 표를 이렇게 읽을 수 있습니다.

import pandas as pd

counts = pd.read_csv("counts.csv")
print(counts.head())

counts.csv가 유전자별 read 수 표라면, 판다스로 표를 읽고 이후 필터링, 평균 계산, 그림 그리기를 할 수 있습니다.

실전 보강: count 숫자를 해석하는 최소 감각

RNA-seq에서 가장 위험한 초보자 실수는 raw count를 그대로 비교하는 것입니다. 예를 들어 다음 표를 봅시다.

gene    control_total_reads  treatment_total_reads  control_count  treatment_count
GeneA   1,000,000            2,000,000              100            160

처리군의 GeneA count가 160이고 대조군은 100이므로 처리군에서 더 높아 보입니다. 하지만 처리군은 전체 read 수도 2배입니다. 단순 count만 보면 전체적으로 더 많이 읽힌 샘플을 “발현이 높다”고 착각할 수 있습니다. 그래서 정규화가 필요합니다.

아주 거칠게 백만 read당 count를 계산하면 다음과 같습니다.

control:   100 / 1,000,000 × 1,000,000 = 100
treatment: 160 / 2,000,000 × 1,000,000 = 80

raw count는 처리군이 높지만, 전체 read 수를 고려하면 오히려 처리군의 상대적 발현은 낮아집니다. 실제 분석에서는 DESeq2 같은 도구가 더 정교한 방식으로 library size와 분산을 보정합니다. 입문 단계에서는 “샘플마다 읽힌 총량이 다르면 raw count 직접 비교는 위험하다”는 감각을 잡으면 됩니다.

실전 보강: log2 fold change와 FDR

차등발현 결과에서 자주 보는 값은 log2FoldChange와 padj입니다.

log2FoldChange = 1: 처리군 발현이 대조군보다 약 2배 높다는 뜻입니다.
log2FoldChange = -1: 처리군 발현이 대조군보다 약 1/2배 낮다는 뜻입니다.
padj: 여러 유전자를 동시에 검사할 때 생기는 우연한 발견을 보정한 p-value입니다. 흔히 FDR 조정값이라고 생각하면 됩니다.

유전자 20,000개를 동시에 검사하면 우연히 p-value가 작게 나오는 유전자가 생길 수 있습니다. 그래서 RNA-seq에서는 p-value만 보고 결론내리지 않고, padj와 fold change를 함께 봅니다.

초보자가 자주 하는 오해

오해 1: count가 크면 무조건 생물학적으로 중요하다. count가 큰 유전자는 원래 많이 발현되는 housekeeping gene일 수 있습니다. 조건 간 차이를 봐야 합니다.
오해 2: p-value가 작으면 원인 유전자다. 통계적으로 차이가 있다는 뜻이지, 질병의 원인이라는 뜻은 아닙니다.
오해 3: biological replicate 없이 차등발현을 믿을 수 있다. 반복 샘플이 없으면 조건 차이와 샘플 우연 차이를 구분하기 어렵습니다.
오해 4: QC는 형식적인 단계다. QC에서 어댑터 오염이나 품질 저하를 놓치면 뒤의 통계가 깔끔해 보여도 결론이 흔들립니다.

이전 개념과 다음 개념의 연결

RNA-seq 파이프라인은 이 부록의 뒤쪽 항목들과 계속 연결됩니다. count matrix는 E17의 SQL/데이터베이스 사고와 이어지고, 많은 샘플을 반복 처리할 때는 E21 워크플로우가 필요합니다. 데이터가 커지면 E22 HPC가 필요하고, 결과를 다시 만들 수 있게 하려면 E20 환경관리와 E24 재현성이 필요합니다.

생물정보학에서 왜 중요한가

RNA-seq은 생물정보학에서 가장 널리 쓰이는 분석 중 하나입니다. 유전자 발현, 질병 상태, 약물 반응, 세포 유형, 조직 차이를 연구할 때 계속 등장합니다. RNA-seq 파이프라인을 이해한다는 것은 “실험 데이터가 어떻게 숫자 표로 바뀌는지”를 이해한다는 뜻입니다.

어려운 개념 보강: count matrix, metadata, 정규화를 한 줄로 연결하기

RNA-seq에서 가장 중요한 표는 count matrix와 metadata table입니다. count matrix는 “각 유전자가 각 샘플에서 몇 번 읽혔는가”를 담고, metadata는 “각 샘플이 어떤 조건인가”를 담습니다. 둘은 sample_id로 정확히 연결되어야 합니다.

count matrix
            S1   S2   S3
GeneA      100  120  500
GeneB       20   30   25

metadata
sample_id   condition
S1          normal
S2          normal
S3          cancer

여기서 S3의 GeneA count가 높다는 사실만으로는 부족합니다. S3가 cancer 샘플이라는 metadata가 연결되어야 “암 조건에서 GeneA가 높아 보인다”는 해석이 가능합니다.

정규화가 필요한 이유는 샘플마다 전체 read 수가 다를 수 있기 때문입니다. 예를 들어 두 샘플의 GeneA count가 같아도 전체 read 수가 다르면 의미가 달라질 수 있습니다.

S1: GeneA = 100, 전체 read = 10,000
S2: GeneA = 100, 전체 read = 100,000

S1에서는 GeneA가 전체의 1%이고, S2에서는 0.1%입니다. raw count만 보면 둘 다 100이지만, 샘플 전체 규모를 고려하면 다르게 해석됩니다. 이것이 정규화의 직관입니다.

아주 단순한 비율 계산은 다음과 같습니다.

비율 = 특정 유전자 count / 샘플 전체 count

실제 차등발현분석 도구는 더 정교한 size factor나 분포 모델을 사용합니다. 하지만 입문 단계에서는 “샘플마다 읽힌 총량이 달라서 raw count를 그대로 비교하면 위험하다”는 점을 잡으면 됩니다.

주의할 점도 있습니다. DESeq2 같은 도구는 보통 raw count를 입력으로 받고 내부에서 정규화를 수행합니다. 그래서 사용자가 임의로 정규화한 값을 다시 넣으면 분석 가정이 깨질 수 있습니다. 즉, 정규화가 필요하다는 말과 “아무 정규화 값을 아무 도구에 넣어도 된다”는 말은 다릅니다.

미니 실습 블록: RNA-seq 분석 순서와 입력·출력 연결하기

이 실습은 RNA-seq 분석 순서와 입력·출력 연결하기를 직접 손으로 확인하는 연습입니다. 왜 필요한가 하면, RNA-seq은 여러 도구가 이어지는 파이프라인이므로 각 단계의 입력과 출력이 맞아야 다음 단계가 정상 동작하기 때문입니다.

FASTQ 확인
품질관리(QC)
trimming
reference/genome 또는 transcriptome 정렬
gene별 count 생성
count matrix와 metadata 연결
정규화와 차등발현 분석
그림과 보고서 작성

각 코드 요소의 의미를 풀어보면 다음과 같습니다. FASTQ는 원본 read이고, 정렬 후에는 BAM 또는 transcript quantification 결과가 만들어집니다. gene별 count가 모이면 count matrix가 되고, metadata의 condition 정보와 연결해 차등발현을 수행합니다.

생물정보학/계산생물학에서 쓰이는 장면은 분명합니다. 처리군-대조군 bulk RNA-seq에서 DEG를 찾는 기본 흐름입니다.

흔한 오해 또는 주의점도 있습니다. biological replicate 없이 차등발현을 과신하면 조건 차이와 샘플 우연 차이를 구분하기 어렵습니다.

핵심 정리

RNA-seq 분석은 FASTQ 파일에서 시작해 품질관리, trimming, 정렬 또는 정량, count matrix 생성, 정규화, 차등발현분석으로 이어집니다. 핵심은 read를 유전자나 전사체에 연결하고, 그 수를 세어 조건 간 차이를 해석하는 것입니다.