tximport로 salmon·kallisto 정량 불러오기 — RNA-seq 카운트의 출발점 (R 실전)

TL;DR — salmon·kallisto는 빠른 pseudo/selective-alignment로 transcript 수준의 정량⁠(quantification)을 냅니다⁠(quant.sf / abundance.h5). 그런데 DESeq2·edgeR·limma 같은 차등발현⁠(differential expression, DE) 도구는 보통 gene 수준 카운트 행렬에서 출발하죠. 그 둘 사이를 잇는 게 tximport입니다 — transcript→gene 집계와 평균 transcript 길이 보정을 한 번에 처리하고, DESeqDataSetFromTximport() 한 줄로 DE 입력까지 연결합니다.

이 글은 salmon 정량 폴더에서 시작해 tx2gene 매핑 → tximport() → 반환 구조 확인 → DESeq2 입력까지, 복붙하면 도는 코드로 따라갑니다. kallisto·edgeR·limma 경로와 countsFromAbundance 선택 기준, 자주 터지는 에러도 정리합니다. (코드는 tximport·DESeq2 공식 vignette 기준입니다.)

🔗 관련 글  ·  카운트 행렬이 준비된 다음 단계: DESeq2로 차등발현 분석  ·  edgeR 차등발현  ·  GEOquery로 공개데이터 받기

그림 1. salmon·kallisto의 transcript 정량을 tximport가 gene 카운트로 모아 DESeq2·edgeR 입력으로 잇는 업스트림 파이프라인.

이 글에서 할 것

이 블로그의 DE 글들⁠(DESeq2·edgeR·limma-voom)은 전부 카운트 행렬이 이미 있다는 데서 출발합니다. 그런데 그 카운트는 어디서 오나요? 요즘 RNA-seq 정량은 salmon·kallisto로 transcript 수준 추정치를 먼저 뽑고, 그걸 gene 수준으로 모으는 흐름이 표준입니다. 이 글은 바로 그 'DESeq2 글의 한 단계 앞'을 채웁니다.

데모는 단순합니다 — salmon이 샘플마다 만들어 둔 quant.sf 폴더가 있다고 가정하고, tximport로 읽어 gene 수준 카운트로 모은 뒤, DESeqDataSetFromTximport()에 그대로 넣습니다. 최종 결과물은 DESeq2가 바로 받을 수 있는 dds 객체입니다.

txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ condition)
# 이 한 흐름이 transcript 정량을 gene 수준 DE 입력으로 바꾼다

0. 왜 tximport인가 — 1분 직관

salmon (selective alignment)·kallisto (pseudo-alignment)는 리드를 게놈에 일일이 정렬하지 않고, transcriptome을 기준으로 어느 transcript에서 나왔을 확률이 높은지를 빠르게 추정합니다. 출력은 transcript마다의 추정 카운트와 TPM (transcripts per million)이죠. 문제는 두 가지입니다.

• 수준이 다르다. salmon·kallisto는 transcript 수준인데, DE 분석은 보통 gene 수준에서 합니다. 한 유전자의 여러 isoform을 합쳐야 합니다.
• 길이 보정이 필요하다. 같은 유전자라도 샘플마다 우세하게 발현된 isoform이 다르면 평균 transcript 길이가 달라집니다. 단순히 카운트만 더하면 이 차이가 묻혀 편향이 생깁니다.

tximport는 이 둘을 함께 처리합니다 — transcript 추정치를 gene 단위로 합산하고, 샘플별 평균 transcript 길이를 offset으로 넘겨 DESeq2·edgeR이 정규화에 반영하도록 합니다⁠(Soneson et al. 2015). 덕분에 단순 합산보다 정확한 gene 수준 추론이 됩니다.

📌 한 줄 요약 — salmon·kallisto는 transcript를, DE는 gene을 본다. tximport가 그 사이에서 집계 + 길이 보정을 맡아 카운트 행렬과 offset을 만든다.

1. 준비물 — 설치

Bioconductor 패키지 tximport가 핵심이고, 연결할 DE 패키지로 DESeq2⁠(edgeR·limma도 동일하게 가능)를 깝니다. 메타데이터를 자동으로 붙이고 싶으면 tximeta도 유용합니다. kallisto의 .h5를 읽으려면 rhdf5가 필요합니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("tximport", "DESeq2", "tximeta", "rhdf5"))   # 최초 1회

library(tximport)   # 정량 불러오기 + transcript→gene 집계
library(DESeq2)     # DESeqDataSetFromTximport로 바로 연결

여기서 salmon·kallisto 출력이 무엇인지 먼저 짚어 둡시다. salmon은 샘플 폴더마다 quant.sf⁠(transcript별 추정 카운트·TPM·유효 길이)를, kallisto는 abundance.h5·abundance.tsv·run_info.json을 남깁니다. tximport는 이 파일들을 가리키는 경로 벡터만 받으면 됩니다.

2. 파일 경로 벡터 만들기

tximport의 첫 입력은 각 샘플의 정량 파일 경로를 담은 벡터입니다. 핵심은 벡터에 이름⁠(names)을 붙이는 것 — 그 이름이 곧 카운트 행렬의 열 이름이 됩니다.

# 샘플 메타데이터⁠(보통 별도 .csv/.tsv로 관리)
dir     <- "/path/to/project"
coldata <- read.csv(file.path(dir, "samples.csv"))   # run, condition 등 컬럼
coldata$condition <- factor(coldata$condition, levels = c("control", "treated"))

# salmon 출력 구조: <dir>/salmon/<run>/quant.sf
files <- file.path(dir, "salmon", coldata$run, "quant.sf")
names(files) <- coldata$run     # 열 이름이 될 샘플 ID

all(file.exists(files))         # 모두 TRUE여야 진행

[1] TRUE

경로가 하나라도 FALSE면 거기서 멈추세요 — 폴더 구조나 압축 여부⁠(quant.sf vs quant.sf.gz)부터 확인하는 게 빠릅니다. salmon을 --gzip으로 돌렸다면 파일명이 quant.sf.gz입니다.

3. tx2gene — transcript ID를 gene ID로 잇는 지도

tximport가 transcript를 gene으로 모으려면 어떤 transcript가 어떤 gene에 속하는지 알려 주는 2열 표가 필요합니다 — 1열 transcript ID, 2열 gene ID. 이게 tx2gene입니다. GENCODE·Ensembl GTF에서 만들거나, TxDb·EnsDb 같은 주석 패키지에서 뽑습니다.

# 방법 A — TxDb 주석 패키지에서 (예: 사람 hg38)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v44.annotation.gtf")   # 또는 기존 TxDb 사용
k    <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")          # 1열 TXNAME, 2열 GENEID

head(tx2gene, 3)

              TXNAME            GENEID
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000488147.1 ENSG00000227232.5

# 방법 B — 이미 만들어 둔 매핑 표를 읽기 (GENCODE 동봉 파일 등)
tx2gene <- read.csv(file.path(dir, "tx2gene.gencode.v44.csv"))

💡 버전 접미사 주의 — GENCODE transcript ID는 ENST00000456328.2처럼 점 뒤에 버전이 붙습니다. salmon 인덱스의 ID와 tx2gene의 ID에서 이 버전이 어긋나면 매핑이 통째로 비어 버립니다. tximport(..., ignoreTxVersion = TRUE)로 점 이하를 무시하거나, 양쪽 ID 표기를 처음부터 맞추세요.

4. tximport 실행 — salmon

이제 본론입니다. type = "salmon"으로 지정하고 tx2gene을 넘기면, tximport가 모든 샘플의 quant.sf를 읽어 gene 수준으로 집계합니다.

txi <- tximport(files,
                type      = "salmon",
                tx2gene   = tx2gene,
                ignoreTxVersion = TRUE)   # ID 버전 불일치 방지⁠(필요 시)

names(txi)

[1] "abundance"           "counts"              "length"
[4] "countsFromAbundance"

tximport 한 줄이 진행 상황을 찍으며⁠(reading in files 1 2 3 ...) transcript 추정치를 gene 단위로 모읍니다. 반환값은 리스트입니다.

kallisto일 때는?

도구만 바꾸면 됩니다. kallisto는 abundance.h5⁠(권장, rhdf5 필요) 또는 abundance.tsv를 읽습니다.

# kallisto 출력: <dir>/kallisto/<run>/abundance.h5
files_k <- file.path(dir, "kallisto", coldata$run, "abundance.h5")
names(files_k) <- coldata$run

txi_k <- tximport(files_k, type = "kallisto", tx2gene = tx2gene)
# .tsv를 쓰면 GENCODE의 'ENST|ENSG|...' 헤더 분리를 위해 ignoreAfterBar = TRUE

salmon·kallisto·RSEM 모두 같은 인터페이스로 읽힙니다 — type만 바뀌고 나머지는 동일합니다.

5. 반환 구조 — counts · abundance · length

tximport가 돌려주는 리스트의 세 핵심 원소를 이해하면, 뒤 단계가 전부 쉬워집니다.

그림 2. tximport 리스트의 세 행렬과, 여러 transcript가 한 gene으로 모이는 집계 개념.

dim(txi$counts)              # gene × sample 카운트 행렬
txi$counts[1:3, 1:2]         # DE 도구가 받을 추정 카운트
txi$length[1:3, 1:2]         # gene별 평균 transcript 길이⁠(샘플마다 다름)

[1] 38194     6
                  sample1 sample2
ENSG00000000003   742.30  689.11
ENSG00000000419   531.07  604.88
ENSG00000000457   198.44  221.50

• $counts — gene × sample 추정 카운트 행렬. DESeq2·edgeR·limma가 받는 그 행렬입니다. salmon의 확률적 할당 때문에 정수가 아닌 실수일 수 있습니다⁠(정상).
• $abundance — 같은 모양의 TPM 행렬. 라이브러리 크기·길이가 보정된 발현량이라, 빠른 탐색·시각화나 샘플 간 단순 비교에 씁니다⁠(DE 검정 입력으로 직접 쓰진 않습니다).
• $length — gene별 평균 transcript 길이 행렬. 샘플마다 우세 isoform이 달라 값이 다릅니다. 이 행렬이 길이 offset으로 DE 도구에 넘어가, isoform 전환에 따른 편향을 보정합니다.

📈 왜 $length가 따로 있나 — 같은 유전자라도 샘플 A는 긴 isoform, 샘플 B는 짧은 isoform이 우세할 수 있습니다. 그러면 카운트가 같아도 실제 발현은 다르죠. tximport는 이 평균 길이를 offset으로 넘겨, DESeq2·edgeR이 "길이 차이를 감안한 뒤" 비교하게 합니다. 이게 transcript 정량을 gene DE에 쓸 때 단순 합산보다 나은 핵심 이유입니다.

6. DESeqDataSetFromTximport — DESeq2로 연결

가장 깔끔한 길입니다. txi 리스트를 통째로 넘기면, DESeq2가 $counts를 카운트로, $length를 길이 offset으로 자동 등록합니다. 별도 정규화 코드가 필요 없습니다.

# coldata의 행 순서가 txi$counts의 열 순서와 같아야 함
rownames(coldata) <- coldata$run
stopifnot(all(rownames(coldata) == colnames(txi$counts)))

dds <- DESeqDataSetFromTximport(txi,
                                colData = coldata,
                                design  = ~ condition)
dds <- DESeq(dds)
res <- results(dds)                  # 이후는 일반 DESeq2 워크플로와 동일

여기서부터는 DESeq2 글과 완전히 같습니다 — results()로 logFC·padj를 뽑고, 시각화하면 됩니다. tximport는 그 출발점인 dds를 만들어 주는 역할까지입니다.

edgeR·limma로 연결하려면

edgeR·limma도 길이 offset만 제대로 넘기면 됩니다. 최신 edgeR⁠(≥ 4.10)은 전용 함수가 있고, 그 이전 버전이라면 countsFromAbundance로 길이 보정된 카운트를 만들어 DGEList에 넣는 방식이 이식성이 좋습니다.

# (a) 최신 edgeR — 길이 offset을 알아서 등록
library(edgeR)
y <- DGEListFromTximport(txi)        # edgeR >= 4.10
y <- normLibSizes(y)

# (b) 이식성 있는 방식 — 길이 보정 카운트로 DGEList 구성
txi_ls <- tximport(files, type = "salmon", tx2gene = tx2gene,
                   countsFromAbundance = "lengthScaledTPM")
y <- DGEList(txi_ls$counts, group = coldata$condition)
y <- normLibSizes(y)                 # 이후 voom 또는 glmQLFit

lengthScaledTPM으로 만든 카운트는 길이·라이브러리 크기가 이미 반영돼 있어, 별도 offset 없이 edgeR·limma-voom에 바로 넣을 수 있습니다.

7. countsFromAbundance — 언제 무엇을 쓰나

countsFromAbundance는 "카운트를 TPM⁠(abundance)에서 다시 만들지" 정하는 옵션입니다. DE 도구마다 권장값이 다릅니다.

값	무엇을 하나	언제
"no" (기본)	원래 추정 카운트 + 별도 길이 offset	DESeqDataSetFromTximport · 최신 edgeR⁠(offset 자동)
"lengthScaledTPM"	TPM을 평균 길이·라이브러리 크기로 스케일한 카운트	offset을 따로 못 넘길 때 — edgeR·limma 직접 입력
"scaledTPM"	TPM을 라이브러리 크기로만 스케일⁠(길이 보정 X)	길이 보정이 불필요하거나 별도로 할 때
"dtuScaledTPM"	transcript 수준 DTU 분석용⁠(txOut = TRUE)	isoform/DTU 분석⁠(DRIMSeq·DEXSeq)

규칙은 단순합니다 — DESeq2로 갈 거면 "no" (기본)로 두고 offset을 자동으로 쓰게 하고, edgeR·limma에 카운트만 직접 넣을 거면 "lengthScaledTPM"으로 길이 보정을 카운트에 녹이세요. transcript 수준 차등사용⁠(DTU)을 볼 때만 "dtuScaledTPM" + txOut = TRUE입니다.

🧭 tximeta 한 줄 팁 — salmon 출력에 메타데이터가 함께 있으면, tximport 대신 tximeta()를 쓰면 게놈·주석 버전을 자동 인식해 tx2gene을 직접 만들 필요가 줄어듭니다. 재현성이 중요한 프로젝트에서 편리합니다⁠(Love et al. 2020).

8. salmon vs 정렬 기반⁠(featureCounts) — 무엇이 다른가

카운트를 만드는 길은 크게 둘입니다. salmon·kallisto 같은 pseudo/selective-alignment 방식과, STAR·HISAT2로 게놈에 정렬한 뒤 featureCounts·htseq-count로 세는 정렬 기반 방식입니다. tximport는 전자의 출력을 받습니다.

항목	salmon / kallisto (+ tximport)	정렬 기반⁠(STAR → featureCounts)
정렬	transcriptome에 pseudo/selective	genome에 full alignment
출력 수준	transcript → tximport로 gene 집계	gene⁠(또는 exon) 카운트 직접
속도·자원	빠름·메모리 적음	느림·BAM 저장 큼
isoform 처리	다중매핑 리드를 확률적으로 분배	보통 다중매핑 제외⁠(편향 가능)
길이 보정	tximport가 평균 길이 offset 제공	없음⁠(별도 처리 필요)
강점	빠른 대규모 정량·isoform 인지	게놈 좌표·신규 transcript 탐색

그림 3. transcriptome pseudo-alignment (salmon) vs genome 정렬로 세는 카운트⁠(featureCounts) — 같은 카운트, 다른 경로.

둘 다 유효하지만, 요즘 표준 RNA-seq 정량은 속도와 isoform 인지 때문에 salmon·kallisto로 기우는 추세입니다. 그리고 그 출력을 DE로 잇는 다리가 정확히 tximport입니다. 게놈 좌표가 필요한 분석⁠(신규 transcript 발굴, 변이 연계 등)에서는 정렬 기반이 여전히 제 몫을 합니다.

9. 자주 발생하는 에러 & 해결

• tximport 후 카운트가 거의 0이거나 gene이 몇 개뿐 — tx2gene의 transcript ID와 salmon 인덱스 ID의 버전 불일치가 1순위입니다. head(rownames(read.delim(files[1])))로 실제 ID를 보고, ignoreTxVersion = TRUE를 주거나 tx2gene ID를 맞추세요.
• Error: all(file.exists(files)) is not TRUE — 경로가 틀렸거나 압축⁠(quant.sf.gz)입니다. files를 직접 출력해 폴더 구조와 파일명을 눈으로 확인하세요.
• tx2gene 관련 경고로 transcript가 누락 — GTF·TxDb 버전이 salmon 인덱스를 만든 주석과 다른 경우입니다. 정량과 같은 버전의 주석으로 tx2gene을 만드는 게 원칙입니다.
• kallisto .h5를 못 읽음 — rhdf5 미설치입니다. BiocManager::install("rhdf5") 후 다시 시도하거나, abundance.tsv로 읽으세요⁠(GENCODE면 ignoreAfterBar = TRUE).
• DESeqDataSetFromTximport에서 행/열 불일치 — colData의 행 순서가 txi$counts의 열 순서와 달라서입니다. rownames(coldata)와 colnames(txi$counts)가 같은지 stopifnot으로 먼저 확인하세요.
• 카운트가 정수가 아니라고 경고 — salmon의 확률적 할당 때문에 실수 카운트는 정상입니다. DESeqDataSetFromTximport는 이를 알아서 처리합니다⁠(직접 반올림하지 마세요).

🎓 다음 학습
• (바로 다음) DESeq2로 차등발현 분석 — DESeqDataSetFromTximport로 만든 dds를 그대로 이어서 · edgeR 차등발현 · limma-voom 차등발현
• (개념 보강) RNA-seq 정규화란? — TPM·CPM·길이 보정이 왜 필요한지
• (데이터 확보) GEOquery로 공개데이터 받기 — 분석할 RNA-seq 데이터를 어디서 구하나

References

1. Soneson, C., Love, M. I., Robinson, M. D. (2015). Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research, 4, 1521. (tximport 방법론·길이 offset)
2. Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods, 14(4), 417–419.
3. Bray, N. L., Pimentel, H., Melsted, P., Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification (kallisto). Nature Biotechnology, 34(5), 525–527.
4. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. (DESeqDataSetFromTximport)
5. Love, M. I., Soneson, C., et al. (2020). Tximeta: Reference sequence checksums for provenance identification in RNA-seq. PLoS Computational Biology, 16(2), e1007664.
6. tximport vignette (Bioconductor). tximport·tx2gene·countsFromAbundance·ignoreTxVersion·DGEListFromTximport.

Pipette & Pipeline · A bio portfolio journal