RSEM

FASTQ to Paper 시리즈 2부 마감

이 편으로 "FASTQ to Paper" 시리즈 2부의 정량 단계가 마무리됩니다. 다음 편(S19)의 DESeq2 볼케이노 플롯이 이 시리즈의 결승 장면입니다. 여기까지 오시느라 수고하셨습니다.

지난 편 Salmon/kallisto의 유사 정렬이 표준이 된 현재, 전통 카운팅 도구(RSEM · featureCounts · HTSeq)는 왜 아직 쓰일까요. 세 가지 이유가 있습니다.

논문 재현: 옛 논문의 파이프라인을 정확히 재현해야 할 때.
정렬이 이미 있음: STAR/HISAT2로 정렬한 BAM이 있다면 정렬을 버리기 아까움.
특수 요구: 아이소폼 수준 EM(RSEM) 또는 엄격한 카운팅(HTSeq).

이번 편에서 세 도구를 비교하고 실무 결정 규칙을 잡습니다.

featureCounts — 속도의 왕

Subread 패키지의 일부. Yang Liao와 Wei Shi가 개발.

원리

가장 단순한 접근입니다.

정렬된 BAM에서 각 리드를 스캔.
리드가 특정 유전자(GTF의 exon 유니온)와 겹치면 그 유전자에 카운트 +1.
리드가 여러 유전자에 걸치면 옵션에 따라 배분 또는 폐기.

인터벌 트리로 구현되어 있어 극도로 빠릅니다.

실행 명령

bash

featureCounts \
    -a gencode.v44.annotation.gtf.gz \
    -o counts.tsv \
    -p --countReadPairs \
    -s 2 \
    -T 16 \
    -Q 30 \
    sample1.bam sample2.bam sample3.bam

-a: GTF 주석 파일.
-p --countReadPairs: paired-end, 페어 하나 = 카운트 1.
-s 2: strandedness. 0=unstranded, 1=forward, 2=reverse.
-Q 30: 최소 MAPQ.
-T: 스레드.

결과 TSV의 각 행이 유전자, 각 열이 샘플.

Strandedness — 자주 실수하는 옵션

Illumina TruSeq Stranded 프로토콜은 리드가 mRNA의 역방향(antisense) 을 읽습니다. featureCounts -s 2가 이 반전을 반영. 잘못 잡으면 유전자 카운트가 절반으로 떨어지고, 안티센스 유전자(예: 마우스 XIST)만 유전자로 잡히는 이상 결과가 나옵니다.

Strandedness 판정 방법:

bash

# infer_experiment.py (RSeQC 패키지)
infer_experiment.py -i sample.bam -r hg38_gencode.bed

이 결과가 "1++,1--,2+-,2-+"가 우세면 reverse (=-s 2).

강점

속도: 다른 도구의 5～10배.
다중 샘플 동시 처리: 여러 BAM을 한 번에.
가벼움: C로 짜여 있어 메모리 부담 적음.

약점

아이소폼 미지원: 유전자 수준만.
애매한 다중 매핑 리드 처리: 옵션이 많아 신중 필요.

HTSeq — 엄격한 카운팅

Simon Anders와 Wolfgang Huber가 개발. DESeq2 저자와 같은 팀.

원리

featureCounts와 유사하지만 옵션이 엄격합니다.

union 모드 (기본): 리드가 한 유전자에만 걸치면 카운트, 여러 유전자에 걸치면 폐기.
intersection-strict 모드: 리드의 모든 부분이 한 유전자에 있어야 카운트.
intersection-nonempty 모드: 중간.

기본 union 모드가 대부분 실무에 잘 맞습니다.

실행 명령

bash

htseq-count \
    --format=bam \
    --order=name \
    --stranded=reverse \
    --mode=union \
    sample.name_sorted.bam \
    gencode.v44.annotation.gtf.gz > sample.counts.tsv

--order=name: 이름순 정렬된 BAM 요구 (samtools sort -n으로 준비).
--stranded: yes/no/reverse.
--mode: union/intersection-strict/intersection-nonempty.

강점

엄격: DESeq2 저자가 만든 만큼 통계 가정을 엄격히.
표준 참조 구현: featureCounts의 결과와 비교되는 gold standard.

약점

속도: 다른 도구의 1/5 이하.
이름순 정렬 요구: 파이프라인 유연성 저하.

RSEM — 아이소폼 수준 EM

Bo Li가 위스콘신에서 개발. 유사 정렬 이전 시대의 아이소폼 정량 표준.

원리

Salmon과 유사한 EM이지만 정렬을 먼저 합니다. Salmon과 다른 점.

RSEM: STAR나 Bowtie2로 트랜스크립트에 정렬 → EM으로 아이소폼 배분.
Salmon: 정렬 없이 k-mer 해시 → EM으로 아이소폼 배분.

정렬 단계가 있어 훨씬 느립니다.

실행 명령

bash

# 참조 인덱스
rsem-prepare-reference \
    --gtf gencode.v44.annotation.gtf.gz \
    --star \
    --star-path $(which STAR) \
    -p 16 \
    hg38.fa hg38_rsem

# 정량
rsem-calculate-expression \
    --paired-end --star --star-path $(which STAR) \
    -p 16 \
    --output-genome-bam \
    clean_R1.fq.gz clean_R2.fq.gz \
    hg38_rsem \
    SAMPLE01

결과에 유전자 수준(.genes.results)과 아이소폼 수준(.isoforms.results) 파일.

강점

아이소폼 수준 정량: 각 트랜스크립트의 IsoPct(isoform percent) 산출.
참조 구현: 많은 논문이 RSEM 결과를 인용.

약점

속도: 매우 느림. Salmon의 20～50배.
정렬 요구: 정렬 인프라 필수.

세 도구 실무 매트릭스

상황	1순위	2순위
신규 프로젝트 유전자 수준 정량	Salmon (S17)	featureCounts
STAR/HISAT2 BAM이 이미 있음	featureCounts	HTSeq
아이소폼 수준 정량 (알렙터너티브 스플라이싱)	Salmon	RSEM
논문 재현 (원 논문 도구)	원 도구 그대로	—
DESeq2 원 논문 스타일 재현	HTSeq	featureCounts
속도 우선	Salmon	featureCounts

핵심 결론: 새 프로젝트라면 Salmon(S17), 정렬이 이미 있으면 featureCounts, 참조 재현이면 원 도구.

Salmon vs featureCounts — 결과 비교

두 도구 모두 실행해 결과 비교하면 유전자 수준 카운트가 대략 90～95% 일치합니다. 차이의 주 원인:

다중 매핑 리드 처리: Salmon은 EM으로 배분, featureCounts는 폐기.
아이소폼 편향: 특정 아이소폼이 우세할 때 Salmon이 정확.

DGEA 결과에도 영향이 있지만, 대부분 강력한 신호는 두 도구 결과가 일치합니다.

손 계산 — Union 모드가 실제로 하는 일

가상 시나리오. 리드 R이 GTF 상 두 유전자 g1과 g2에 걸침 (겹친 엑손 등).

HTSeq union 모드: R을 폐기. 카운트 0.
featureCounts 기본: R을 폐기 (multi-overlap).
featureCounts --allowMultiOverlap: g1과 g2 모두 +1 (또는 소수 배분 옵션).

기본 관례는 폐기입니다. 실무 노트 → 한 유전자가 다른 유전자 안에 있는 경우 (예: 미토콘드리아 tRNA in rRNA) 카운트 유실이 발생할 수 있음. 이런 유전자 특별 취급이 필요합니다.

R로 DESeq2 입력 만들기

세 도구 모두 R에서 DESeq2로 넘길 수 있습니다.

featureCounts

library(DESeq2)

counts <- read.table("counts.tsv", header=TRUE, skip=1, row.names=1)
counts <- counts[, 6:ncol(counts)]  # meta 열 제거
colnames(counts) <- gsub(".bam", "", colnames(counts))

coldata <- data.frame(row.names=colnames(counts), 
                      condition=c("normal","normal","tumor","tumor"))

dds <- DESeqDataSetFromMatrix(countData=counts, colData=coldata, design=~condition)

HTSeq

library(DESeq2)

files <- list.files("counts_dir", pattern="*.tsv", full.names=TRUE)
sampleTable <- data.frame(
    sampleName = gsub(".counts.tsv", "", basename(files)),
    fileName = files,
    condition = c("normal","normal","tumor","tumor"))

dds <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable,
                                  directory=".",
                                  design=~condition)

library(tximport)
files <- list.files("rsem", pattern="*.genes.results", full.names=TRUE)
txi <- tximport(files, type="rsem", txIn=FALSE, txOut=FALSE)

dds <- DESeqDataSetFromTximport(txi, colData=coldata, design=~condition)

세 접근 모두 결국 DESeq2의 DESeqDataSet 객체로 수렴.

Colab에서 소형 실습

bash

!apt-get install -y subread >/dev/null

# STAR 정렬 BAM이 있다고 가정 (S04에서 준비)
!featureCounts \
    -a gencode.v44.annotation.gtf.gz \
    -o counts.tsv \
    -p --countReadPairs \
    -s 2 \
    -T 2 \
    -Q 30 \
    SAMPLE01_Aligned.sortedByCoord.out.bam

!head -5 counts.tsv

결과 열: Geneid, Chr, Start, End, Strand, Length, sample1.bam. 유전자별 카운트가 마지막 열.

CS 매핑

인터벌 트리: featureCounts의 GTF exon 조회. DryBench "인터벌 트리 기초" 편(interval-tree-basics) 참고.
다중 매핑 처리: 리드가 여러 인터벌에 걸칠 때의 결정 규칙.
EM 재사용: RSEM은 Salmon과 같은 EM을 정렬 후에 적용.

FASTQ to Paper 시리즈 여정 정리

여기까지의 13편(S06～S18)이 "FASTQ to Paper" 시리즈의 정량·변이 파이프라인을 완주했습니다.

S06～S11 (GATK4 6편): DNA 변이의 표준 파이프라인.
S17～S18 (RNA-seq 정량 2편): 발현량 계산.
앞으로 S19～S21 (DGEA 3편): 실제 볼케이노 플롯 완성.

여러분은 이제 원시 FASTQ에서 논문 그림까지의 전 여정을 하나의 흐름으로 이해할 수 있습니다.

마무리

RSEM · featureCounts · HTSeq은 유사 정렬 이전의 표준이었고 지금도 특정 상황에서 유효합니다. 도구가 여러 개인 이유가 각각의 강점 때문임을 감으로 잡아두시면 됩니다. 다음 편(S19)이 시리즈의 하이라이트 — DESeq2로 볼케이노 플롯을 완성합니다. 이미 파일럿 편으로 배포되어 있으니 오늘 배운 정량 결과를 그대로 이어받아 확인해보세요.

더 깊게 파고 싶다면

Liao Y. et al. (2014), featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30:923. featureCounts 원 논문.
Anders S. et al. (2015), HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31:166. HTSeq 원 논문.
Li B. & Dewey C.N. (2011), RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. RSEM 원 논문.
Xiaole Shirley Liu Harvard STAT115 W4 — RNA-seq 정량 강의. 자막 완비. featureCounts와 Salmon의 결과 차이 실무 사례.
Bioconductor tximport 튜토리얼: 세 도구 결과를 DESeq2로 통합하는 표준 방법.

한 샘플로 featureCounts · Salmon 둘 다 돌려 카운트를 비교해보시면 실무 감이 옵니다. 다음 편에서 파일럿 S19 DESeq2 볼케이노 플롯으로 시리즈의 결승을 확인해봅시다.

BioPlayground

RSEM · featureCounts · HTSeq: 전통 정렬 후 카운팅의 세 갈래

FASTQ to Paper 시리즈 2부 마감

featureCounts — 속도의 왕

원리

실행 명령

Strandedness — 자주 실수하는 옵션

강점

약점

HTSeq — 엄격한 카운팅

원리

실행 명령

강점

약점

RSEM — 아이소폼 수준 EM

원리

실행 명령

강점

약점

세 도구 실무 매트릭스

Salmon vs featureCounts — 결과 비교

손 계산 — Union 모드가 실제로 하는 일

R로 DESeq2 입력 만들기

featureCounts

HTSeq

RSEM

Colab에서 소형 실습

CS 매핑

FASTQ to Paper 시리즈 여정 정리

마무리

더 깊게 파고 싶다면