GEOquery로 공개 발현 데이터 받아서 분석하기 — GEO에서 내 PC로, 재분석까지 (R 실전)
TL;DR — 논문에 "데이터는 GEO에 GSE12345로 올렸다"는 한 줄이 보이면, 그 발현 데이터를 내 R 세션으로 그대로 가져와 재분석할 수 있습니다. GEOquery는 그 다리입니다 —
getGEO("GSExxxxx")한 줄이면 발현행렬(expression matrix)·샘플 메타데이터(metadata)·프로브 주석(annotation)이 한 객체로 들어옵니다. 워크플로는 짧습니다 — GEO 검색 →getGEO→exprs/pData/fData→ 그룹 라벨 정리 → 다운스트림(limma·DESeq2).단, 마이크로어레이(microarray)와 RNA-seq는 GEO에서 오는 형태가 다릅니다 — 어레이는 정규화된 발현행렬이 바로 오지만, RNA-seq는 보통 원시 카운트(raw count)가 보충 파일(supplementary file)로 따로 옵니다. 이 글은 그 차이까지 포함해, 복붙하면 그대로 도는 GEOquery 코드를 따라갑니다. (코드는 Bioconductor GEOquery 비네트 기준입니다.)
🔗 관련 글 · 받은 데이터를 바로 분석에 투입: DESeq2로 차등발현 분석 · edgeR 차등발현 분석 · RNA-seq 정규화란?
0. 왜 GEOquery인가 — 1분 직관
블로그의 RNA-seq 분석 글들(DESeq2·edgeR·정규화·PCA)은 모두 "발현 데이터가 이미 손에 있다"는 전제로 시작합니다. 그런데 실전에서 가장 먼저 막히는 건 "그 데이터를 대체 어디서 가져오나"입니다. 답의 8할은 GEO (Gene Expression Omnibus) — NCBI가 운영하는 세계 최대 공개 발현 데이터 저장소입니다.
- 논문 Methods나 Data Availability에 적힌 GSE 번호 하나면, 그 연구의 발현행렬과 샘플 정보를 통째로 받을 수 있습니다. 내 실험 데이터가 없어도 남의 공개 데이터로 재분석·검증·연습이 가능합니다.
- 문제는 GEO 웹페이지에서 손으로 파일을 받아 R로 읽으면 포맷이 제각각이라 번거롭다는 점입니다. GEOquery는 이 과정을 함수 한 줄로 줄입니다 — 다운로드·파싱·객체화를 한 번에.
- 받은 결과는 R의 표준 발현 객체인 ExpressionSet (S4 클래스) 또는 SummarizedExperiment라, limma·DESeq2·edgeR 같은 분석 패키지에 바로 연결됩니다.
1. 준비물 — 설치와 GEO 개념
GEOquery는 Bioconductor 패키지입니다. 어레이 데이터를 다룰 limma, 카운트를 다룰 DESeq2도 함께 깔아 둡니다.
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("GEOquery", "limma", "DESeq2")) # 최초 1회
library(GEOquery)먼저 GEO의 세 가지 식별자를 구분해야 코드가 헷갈리지 않습니다. 셋은 포함 관계입니다 — GSE 하나가 여러 GSM을 묶고, 각 GSM은 하나의 GPL 위에서 측정됩니다.
| 식별자 | 뜻 | 비유 | 받는 함수 |
|---|---|---|---|
| GSE (Series) | 한 연구 전체(샘플 묶음) | 논문 한 편의 데이터셋 | getGEO("GSE...") |
| GSM (Sample) | 개별 샘플 하나 | 환자/조직 한 개 | getGEO("GSM...") |
| GPL (Platform) | 측정 플랫폼(칩·시퀀서) | 어떤 마이크로어레이 기종인가 | getGEO("GPL...") |
2. getGEO — 시리즈 데이터 받기
핵심 함수는 단 하나, getGEO()입니다. 재현 가능한 공개 예시로 GSE10072 (폐 선암종 종양 vs 정상, Affymetrix HG-U133A 어레이, 107 샘플)를 받아 보겠습니다. 마이크로어레이라 정규화된 발현행렬이 series matrix에 바로 들어 있습니다.
# GSE 시리즈 다운로드 — 처음엔 시간이 걸림(파일을 받아 캐시에 저장)
gse <- getGEO("GSE10072", GSEMatrix = TRUE) # GSEMatrix=TRUE가 기본(빠른 series matrix)
# ★ getGEO는 '리스트'를 반환 — 한 GSE가 여러 플랫폼을 가질 수 있기 때문
length(gse) # 플랫폼 개수(여기선 1)
eset <- gse[[1]] # 첫 번째 플랫폼의 ExpressionSet을 꺼냄
eset # 객체 요약: assayData(발현), phenoData(메타), featureData(주석)ExpressionSet (storageMode: lockedEnvironment)
assayData: 22283 features, 107 samples
element names: exprs
phenoData: sampleNames: GSM254625 GSM254626 ... (107 total)
varLabels: title geo_accession ... characteristics_ch1
featureData: featureNames: 1007_s_at 1053_at ... (22283 total)요약을 보면 22283 프로브 × 107 샘플의 발현행렬, 샘플별 표현형(phenotype) 표, 프로브 주석이 한 객체에 들어 있는 게 보입니다. 이게 ExpressionSet의 세 칸입니다.
getGEO()는 단일 ExpressionSet이 아니라 리스트를 돌려줍니다(gse[[1]]로 꺼내기). 시리즈가 두 플랫폼(예: 어레이 + miRNA 칩)에 걸쳐 있으면 원소가 2개가 되니, length(gse)를 먼저 확인하는 습관을 들이세요.3. exprs / pData / fData — 세 칸 꺼내기
ExpressionSet에서 실제로 쓸 세 덩어리를 뽑는 접근자(accessor)는 외우기 쉽습니다 — 발현은 exprs, 표현형은 pData, 피처(프로브) 주석은 fData입니다.
# ① 발현행렬 — 행=프로브, 열=샘플 (어레이는 보통 log2 정규화된 강도)
expr <- exprs(eset)
dim(expr) # 22283 × 107
expr[1:4, 1:3] # 좌상단만 살짝
# ② 샘플 메타데이터(표현형) — 행=샘플, 열=변수
meta <- pData(phenoData(eset)) # pData(eset)으로 줄여도 동일
dim(meta) # 107 × (수십 개 변수)
colnames(meta)[1:8] # title, geo_accession, source_name_ch1, characteristics_ch1 ...
# ③ 피처(프로브) 주석 — 프로브 ID ↔ 유전자 심볼 등
feat <- fData(eset)
colnames(feat) # ID, "Gene Symbol", "Gene Title", ENTREZ_GENE_ID ...> expr[1:4, 1:3]
GSM254625 GSM254626 GSM254627
1007_s_at 9.84 9.71 10.02
1053_at 6.12 6.33 6.05
117_at 7.45 7.20 7.51
121_at 9.03 8.88 9.14세 줄이 끝입니다. expr은 분석에 넣을 발현행렬, meta는 그룹을 가를 표현형 표, feat은 프로브를 유전자 이름으로 바꿀 사전입니다. 이제 "정상 vs 종양" 같은 비교 라벨만 정리하면 분석 준비가 끝납니다.
eset@assayData 처럼 슬롯(slot)에 직접 손대지 말고 exprs()·pData()·fData() 함수를 쓰는 게 안전합니다. S4 객체의 내부 구조가 바뀌어도 접근자는 동작을 보장합니다.4. 그룹 라벨 정리 — 메타데이터에서 비교군 만들기
다운스트림 분석은 "어느 샘플이 어느 그룹인가"를 알아야 합니다. 그 정보는 meta의 characteristics_ch1나 source_name_ch1 같은 열에 자유 텍스트(free text)로 들어 있어, 깔끔한 요인(factor)으로 바꿔 줘야 합니다.
# 어떤 열에 그룹 정보가 있는지 눈으로 확인
meta$source_name_ch1[1:6]
# [1] "Adenocarcinoma of the Lung" "Adenocarcinoma of the Lung"
# "Normal Lung Tissue" ...
# 텍스트에서 종양/정상 라벨을 뽑아 요인으로 — grepl로 패턴 매칭
group <- ifelse(grepl("Normal", meta$source_name_ch1), "Normal", "Tumor")
group <- factor(group, levels = c("Normal", "Tumor")) # 대조군을 앞에
table(group) # 그룹별 샘플 수 확인group
Normal Tumor
49 58characteristics_ch1에 "disease state: tumor" 식으로 들어 있기도 하고, 여러 열에 흩어져 있기도 합니다. 항상 직접 열어 보고(View(meta) 또는 head) 패턴을 확인한 뒤 라벨링하세요.5. 탐색 — 차원·결측·로그 변환 점검
분석에 넣기 전, 발현행렬이 멀쩡한지 빠르게 점검합니다. 핵심은 세 가지 — 차원이 맞나, 결측(NA)이 있나, 이미 로그 변환됐나입니다.
dim(expr) # 프로브 × 샘플 — 라벨 길이(107)와 열 수가 같아야
sum(is.na(expr)) # 결측값 개수 — 어레이는 종종 NA가 섞임
range(expr, na.rm = TRUE) # 값의 범위 — log2면 대략 2~16, 원시강도면 수천~수만
# 이미 log 스케일인지 빠른 판정(값의 최댓값이 크면 원시 강도일 가능성)
qx <- quantile(expr, c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
qx 0% 25% 50% 75% 100%
2.30 5.41 6.58 8.12 14.92 <- 범위가 2~15 → 이미 log2 정규화됨위처럼 값이 대략 2~15 사이면 이미 log2로 정규화된 series matrix라 추가 변환이 필요 없습니다. 만약 최댓값이 수만 단위로 크면 원시 강도이니 log2(expr + 1)로 변환합니다. 결측이 많은 프로브는 expr <- expr[complete.cases(expr), ]로 솎아 냅니다.
# 원시 강도일 때만(범위가 매우 클 때) 로그 변환
if (max(expr, na.rm = TRUE) > 100) {
expr <- log2(expr + 1)
}
# 빠른 시각 점검 — 샘플 간 분포가 비슷한지(boxplot) · 이상 샘플 탐지
boxplot(expr[, 1:10], las = 2, main = "Sample distributions (log2)")이제 expr(발현행렬) + group(비교 라벨)이 준비됐습니다. 이 둘이 곧 limma(어레이) 또는 DESeq2/edgeR(카운트)의 입력입니다.
6. 마이크로어레이 vs RNA-seq — GEO에서 뭐가 다른가
여기가 GEOquery 실전의 가장 중요한 분기점입니다. getGEO()로 받았을 때 exprs()에 값이 바로 차 있으면 보통 마이크로어레이입니다. 반대로 exprs()가 비거나(0행) 이상하면, 십중팔구 RNA-seq라 원시 카운트가 보충 파일에 따로 있습니다.
| 항목 | 마이크로어레이 (예: Affymetrix GSE) | RNA-seq (카운트 기반 GSE) |
|---|---|---|
| series matrix의 발현행렬 | 정규화된 강도가 바로 옴 | 비어 있거나 불완전한 경우가 많음 |
exprs(eset) 결과 | log2 정규화 행렬(바로 사용) | 빈 행렬 / NA → 추가 처리 필요 |
| 원시 데이터 위치 | series matrix 안에 포함 | 보충 파일(supplementary, .txt.gz/.csv.gz) |
| 받는 함수 | getGEO() 하나로 끝 | getGEO() + getGEOSuppFiles() |
| 값의 형태 | 연속형 정규화 강도 | 정수 원시 카운트(raw count) |
| 다운스트림 | limma | DESeq2 / edgeR |
| 전처리 | 대개 정규화 완료 | 직접 정규화·필터링 필요 |
exprs()가 비는가 — GEO의 series matrix는 본래 마이크로어레이 시대에 설계된 "행렬" 포맷입니다. RNA-seq 연구자들은 정규화 방식이 제각각이라(TPM·FPKM·raw count) 표준 행렬 대신 원시 카운트를 보충 파일로 첨부하는 게 관행이 됐습니다. 그래서 RNA-seq GSE는 exprs()가 비어 있고, 카운트 파일을 따로 받아야 합니다.RNA-seq — getGEOSuppFiles로 원시 카운트 받기
RNA-seq 시리즈(예: GSE164073)라면 보충 파일을 받아 직접 읽습니다.
# ① 먼저 어떤 보충 파일이 있는지 '목록만' 확인(fetch_files=FALSE)
supp <- getGEOSuppFiles("GSE164073", fetch_files = FALSE)
supp[, c("fname")] # 예: GSE164073_raw_counts.txt.gz
# ② 실제 다운로드 — 현재 작업폴더 아래 GSE164073/ 폴더로 받아짐
getGEOSuppFiles("GSE164073")
# ③ 받은 gzip 카운트 파일을 행렬로 읽기 — 첫 열을 유전자 ID(행이름)로
counts <- read.delim(gzfile("GSE164073/GSE164073_raw_counts.txt.gz"),
row.names = 1, check.names = FALSE)
counts <- as.matrix(counts)
dim(counts) # 유전자 × 샘플 정수 행렬 → DESeq2/edgeR 입력getRNASeqData("GSE...")는 이를 SummarizedExperiment 객체로 바로 받아 줍니다(hasRNASeqQuantifications("GSE...")로 가능 여부 확인). 보충 파일 포맷이 제각각이라 파싱이 번거로울 때 유용합니다.7. 자주 발생하는 에러 & 해결
exprs(eset)가 0행이거나 비어 있음 — 거의 항상 RNA-seq 시리즈입니다. series matrix엔 카운트가 없으니getGEOSuppFiles()로 보충 파일을 받아 읽으세요(6절).- 다운로드 타임아웃 / 중간에 끊김 — 큰 series matrix는 기본 타임아웃(60초)을 넘길 수 있습니다.
options(timeout = 600)으로 늘리고 다시 시도하세요. 캐시 덕에 받다 만 파일은 이어받습니다. Error: ... VROOM_CONNECTION_SIZE— 시리즈가 아주 클 때 파서 버퍼가 부족하다는 신호입니다.Sys.setenv(VROOM_CONNECTION_SIZE = 500000)후getGEO()를 다시 부르면 해결됩니다.length(gse)가 2 이상(여러 플랫폼) — 한 GSE가 두 칩에 걸친 경우입니다. 각 원소가 독립 ExpressionSet이니gse[[1]],gse[[2]]를 따로 다루거나, 분석할 플랫폼(GPL)만 고르세요.- 프로브 ID는 있는데 유전자 심볼이 없음 —
fData(eset)에"Gene Symbol"열이 없으면, 플랫폼 주석을 따로 받습니다 —getGEO("GPL96")의 Table에 프로브↔유전자 매핑이 있고, 어레이 기종이면 Bioconductor 주석 패키지(예:hgu133a.db)로도 매핑됩니다. - 같은 프로브가 여러 개(중복 유전자) — 한 유전자에 프로브가 여러 개일 수 있습니다. 분석 전 발현 평균/최댓값으로 대표 프로브 하나만 남기거나(
limma::avereps), 프로브 수준 그대로 분석 후 마지막에 매핑하세요.
8. 다운스트림 연결 — limma / DESeq2 · edgeR
받은 데이터의 형태에 따라 분석 도구가 갈립니다. 정규화된 어레이 강도는 limma, 정수 카운트는 DESeq2/edgeR입니다.
# ── (A) 마이크로어레이 → limma ── expr=log2 정규화 행렬, group=요인
library(limma)
design <- model.matrix(~ group) # Normal(기준) vs Tumor
fit <- lmFit(expr, design) # 선형모형 적합
fit <- eBayes(fit) # 경험적 베이즈 보정
topTable(fit, coef = 2, number = 10) # 상위 차등발현 프로브(logFC·adj.P)# ── (B) RNA-seq 카운트 → DESeq2 ── counts=정수 행렬, coldata=샘플 정보
library(DESeq2)
coldata <- data.frame(group = group, row.names = colnames(counts))
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = coldata,
design = ~ group)
dds <- DESeq(dds)
results(dds) # logFC·padj 차등발현 결과여기까지 오면 GEO의 공개 데이터가 완전한 분석 파이프라인으로 들어왔습니다 — GEO 검색 → GEOquery 다운로드 → 발현행렬·메타 → 그룹 라벨 → limma/DESeq2. 이 블로그의 차등발현·정규화·농축분석 글들이 모두 이 입력에서 출발합니다.
• (어레이 분석) RNA-seq 정규화란? — 카운트를 다룰 때 정규화의 이유
• (카운트 분석) DESeq2로 차등발현 분석 · edgeR 차등발현 분석 — 받은 카운트의 다음 단계
• (QC·해석) PCA로 샘플 QC — 받은 데이터의 배치·이상치 점검 · clusterProfiler 농축분석 · p값과 FDR
보너스 — 데이터셋 고르기: GEO에서 무엇을 받을까
재분석할 GSE를 고를 땐 세 가지를 봅니다. ① 플랫폼 — 같은 종·같은 기술(어레이 vs RNA-seq)을 골라야 비교가 깔끔합니다. ② 샘플 수 — 그룹당 최소 3개, 통계적 힘을 원하면 그룹당 6개 이상이 안전합니다. ③ 실험 설계 — 비교군이 명확한지(종양 vs 정상, 처리 vs 대조), 메타데이터에 그룹 정보가 잘 적혀 있는지 GEO 페이지에서 미리 확인하세요. GEO 검색창에서 "lung cancer"[Title] AND "Homo sapiens"[Organism] AND "expression profiling by array"[DataSet Type] 같은 필터로 좁히면 원하는 시리즈를 빠르게 찾습니다.
References
- Davis, S., Meltzer, P. S. (2007). GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 23(14), 1846–1847.
- GEOquery Vignette (Davis, S., Bioconductor). getGEO·GSEMatrix·exprs·pData·fData·getGEOSuppFiles·getRNASeqData·hasRNASeqQuantifications.
- Barrett, T., et al. (2013). NCBI GEO: archive for functional genomics data sets — update. Nucleic Acids Research, 41(D1), D991–D995. (GSE·GSM·GPL 구조)
- Edgar, R., Domrachev, M., Lash, A. E. (2002). Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Research, 30(1), 207–210.
- Ritchie, M. E., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47.
- 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.
Pipette & Pipeline · A bio portfolio journal
'Code Bench > R Snippets' 카테고리의 다른 글
| ggplot2 기초 — 바이오 데이터 그래프 그리기 (산점도·막대·박스플롯, R 실전) (0) | 2026.06.17 |
|---|---|
| limma-voom으로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2·edgeR과 무엇이 다른가 (R 실전) (0) | 2026.06.17 |
| WGCNA로 유전자 공발현 네트워크 분석하기 — 모듈부터 허브 유전자까지 (R 실전) (0) | 2026.06.16 |
| edgeR로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2와 무엇이 다른가 (R 실전) (0) | 2026.06.15 |
| 주성분분석(PCA)으로 RNA-seq 샘플 QC하기 — vst·plotPCA·prcomp·배치효과 (R 실전) (0) | 2026.06.11 |
