LOCRC 실습 1: raw count에서 DEG 후보까지

1. 무엇을 실습하나

이 실습에서 답할 질문은 하나입니다.

같은 환자에게서 얻은 대장암 조직과 인접 정상조직을 비교하면, 어떤 유전자의 발현이 일관되게 달라지는가?

이처럼 두 조건 사이에서 발현량이 일관되게 달라진 유전자를 차등 발현 유전자(differentially expressed gene, DEG)라고 부릅니다. 이 실습에서는 Tumor와 Normal의 차이가 환자 22명에게 반복되는지 통계적으로 검사해 DEG를 찾습니다.

사용할 데이터는 NCBI GEO의 GSE251845입니다. 50세를 넘겨 진단된 후기 발병 대장암(late-onset colorectal cancer, LOCRC) 환자 22명에게서 종양 조직과 인접 정상조직을 하나씩 채취했으므로 총 44개 RNA-seq 샘플이 있습니다. FASTQ부터 다시 처리하지 않고, GEO가 제공하는 HTSeq raw count에서 시작합니다.

LOCRC는 별도의 병리 진단명이나 암 종류가 아니라, 대장암 환자를 진단 연령으로 나눈 연구·역학 용어입니다. 보통 50세 미만의 조기 발병 대장암과 구분할 때 사용하며, 정확한 연령 경계는 연구마다 확인해야 합니다. GSE251845 원 논문은 50세 미만과 50세 초과 환자군을 구분했습니다.

colorectal cancer 조직은 무엇인가

colorectal은 colon(결장)과 rectum(직장)을 함께 가리킵니다. colorectal cancer 조직은 이 부위에 생긴 종양에서 수술이나 생검으로 얻은 조직입니다. GSE251845에서는 수술로 절제한 대장 종양과 같은 환자의 인접한 비종양 조직에서 RNA를 추출했습니다.

대장은 맹장(cecum)에서 시작해 상행결장, 횡행결장, 하행결장, 구불결장(sigmoid colon)으로 이어지고, 그 아래의 직장(rectum)으로 연결됩니다. 이 colon과 rectum에 생기는 암을 함께 colorectal cancer라고 부릅니다. 출처: 미국 국립암연구소 SEER 대장 해부학 자료, Wikimedia Commons의 퍼블릭 도메인 파일 정보.

그림은 종양이 생기는 해부학적 위치를 보여 줍니다. RNA-seq가 실제로 측정하는 것은 그림 속 장기 전체가 아니라, 수술한 대장에서 떼어 낸 작은 조직 조각입니다. 따라서 count를 해석하려면 그 조각 안에 어떤 세포가 들어 있는지도 알아야 합니다.

GSE251845의 C 샘플에서는 암세포만 따로 분리하지 않았습니다. 악성 상피세포와 함께 면역세포, 섬유아세포, 혈관세포, 남아 있는 정상 상피세포도 들어 있습니다. bulk RNA-seq는 이 세포들이 낸 RNA를 합쳐 측정하므로, 이 실습의 DEG는 암세포만의 변화가 아니라 대장 종양 조직 전체의 변화입니다.

예를 들어 Tumor에서 면역 관련 유전자가 증가해도 암세포가 그 유전자를 더 많이 발현한 것인지, 조직 안의 면역세포 비율이 늘어난 것인지 bulk RNA-seq만으로는 구분할 수 없습니다. 이런 조직 구성의 한계와 인접 정상조직을 대조군으로 쓰는 이유는 종양 조직과 인접 정상조직에서 설명합니다.

분석할 데이터는 어떻게 생겼나

파일을 열면 유전자와 샘플 이름, 정수 count가 들어 있는 큰 표가 나옵니다. 아래는 실제 GSE251845_htseq_raw_counts.csv.gz에서 첫 다섯 유전자와 첫 두 환자의 열만 잘라낸 것입니다.

gene ID	24C_htseq.out	24N_htseq.out	27C_htseq.out	27N_htseq.out
ENSG00000000003	9,896	6,279	5,653	4,198
ENSG00000000005	32	49	40	43
ENSG00000000419	5,129	2,406	3,549	2,052
ENSG00000000457	929	1,288	522	649
ENSG00000000460	952	315	191	182

행 이름인 ENSG00000000003은 유전자를 구분하는 Ensembl gene ID입니다. 열 이름에서 24는 환자 번호, C는 cancer 조직, N은 같은 환자의 adjacent normal 조직을 뜻합니다. _htseq.out은 HTSeq가 만든 샘플별 결과 파일 이름에서 붙은 부분이므로 metadata를 만들 때 제거합니다.

첫 번째 숫자 9,896은 24번 환자의 종양 조직에서 ENSG00000000003 유전자에 배정된 RNA-seq fragment가 9,896개였다는 뜻입니다. 바로 옆의 6,279는 같은 환자의 인접 정상조직에서 관측된 값입니다.

위 표는 파일 맨 위의 유전자 5행만 보여 주므로 __no_feature는 보이지 않습니다. __no_feature, __ambiguous 같은 유전자가 아닌 특수 카운터는 실제 파일의 맨 아래에 있습니다. 각 행의 의미와 분석 전에 분리하는 이유는 HTSeq raw count 레퍼런스에서 확인할 수 있습니다.

raw count에서 DEG까지

9,896 ÷ 6,279를 계산하면 약 1.58이므로 이 환자에서는 Tumor count가 더 높아 보입니다. 하지만 이 비율 하나만으로 DEG라고 판단할 수는 없습니다.

첫 번째 이유는 샘플마다 읽은 RNA-seq fragment의 총량이 다르기 때문입니다. 더 깊게 시퀀싱한 샘플에서는 대부분의 유전자 count가 함께 커질 수 있으므로, 이 실습에서 사용할 PyDESeq2가 샘플별 시퀀싱 깊이를 먼저 보정합니다.

두 번째 이유는 환자마다 원래 발현 기준선이 다르기 때문입니다. 24C와 24N은 같은 사람에게서 얻은 한 쌍이므로, 다른 환자의 샘플처럼 독립적으로 취급하지 않습니다. 분석 모델에 환자 ID를 넣으면 환자별 기준선 차이를 고려한 뒤, 22명에게 반복되는 Tumor 효과를 추정할 수 있습니다.

이 paired 구조를 PyDESeq2에 전달할 때 design="~ patient + condition"을 사용합니다. 표기법은 design formula 레퍼런스에서 설명합니다.

실습에서 확인할 순서

에이전트에게 한 단계씩 요청하면서 다음 질문을 확인합니다.

1. 파일에 정말 Tumor 22개와 Normal 22개가 있으며 모든 환자의 짝이 맞는가?
2. 전체 유전자 발현 패턴만으로 Tumor와 Normal이 구분되는가?
3. 발현량이 달라진 유전자는 어느 방향으로 얼마나 움직였는가?
4. 통계적으로 확실하면서 변화량도 큰 후보는 무엇인가?
5. 상위 후보의 변화가 일부 환자에게만 생긴 것인가, 여러 환자에게 반복되는가?

이 실습의 해석 범위

이 실습이 답하는 질문은 GSE251845의 후기 발병 대장암 환자 22명에서 Tumor와 인접 Normal 조직의 발현이 어떻게 다른가입니다.

인접 정상조직은 건강한 대조군이 아니다

N은 같은 환자의 종양 옆에서 채취한, 연구진이 비종양으로 분류한 조직입니다. 종양 미세환경이나 염증의 영향을 받을 수 있으므로 완전히 건강한 사람의 대장으로 해석하지 않습니다. 자세한 비교는 종양 조직과 인접 정상조직에서 설명합니다.

이 페이지는 GSE251845 안에서 LOCRC 종양 조직과 인접 정상조직을 비교하는 데 집중합니다. 아래에는 이 데이터로 만든 PCA와 DEG 결과만 사용합니다.

데이터와 분석 방법의 출처는 GSE251845 GEO 페이지, 원 논문, PyDESeq2 레퍼런스에서 확인할 수 있습니다.

2. 에이전트에게 실습 환경 준비시키기

Codex CLI나 Claude Code를 실행한 뒤 아래 프롬프트를 복사해 전달합니다. 에이전트는 독립된 프로젝트 폴더, Python 3.11 환경, 고정된 패키지 목록, 분석 스크립트와 결과 폴더를 준비합니다. 이 단계에서는 데이터를 내려받거나 분석하지 않습니다.

▶ 실습 환경 설정 프롬프트 복사

GSE251845 RNA-seq 실습 환경을 준비해줘.

요구사항:
1. 현재 폴더가 이미 colorectal-deg 프로젝트인지 확인해. 아니라면 현재 폴더 아래에 colorectal-deg를 만들고 그 안에서 작업해.
2. uv가 설치되어 있는지 확인해. 설치되어 있지 않다면 임의로 설치하지 말고 공식 설치 방법만 알려준 뒤 멈춰.
3. uv init --bare --python 3.11로 프로젝트를 초기화해.
4. uv add "pydeseq2==0.5.4" pandas matplotlib seaborn scikit-learn gprofiler-official을 실행해.
5. analysis.py와 outputs 디렉터리를 만들어. analysis.py에는 아직 분석 코드를 넣지 마.
6. .gitignore에 .venv/, data/, __pycache__/를 추가해. outputs는 분석 결과를 확인할 수 있도록 제외하지 마.
7. uv run python으로 Python, PyDESeq2, pandas와 gprofiler-official의 버전을 출력해 설치 상태를 검증해.
8. 생성하거나 수정한 파일, 사용한 Python 버전, 설치한 주요 패키지 버전을 보고하고 멈춰.

전역 Python 환경이나 저장소의 다른 파일은 수정하지 마.

준비 직후 생기는 폴더와 파일

프롬프트가 끝나면 아래 구조가 생깁니다. 📁은 폴더, 📄는 파일입니다.

• 📁 colorectal-deg/: 이 실습만을 위한 프로젝트 루트
  • 📁 .venv/: uv가 관리하는 Python 3.11 가상환경. .gitignore에 들어가므로 Git에는 저장되지 않습니다.
  • 📁 outputs/: 이후 프롬프트에서 만드는 그래프와 결과표를 저장할 폴더. 준비 직후에는 비어 있습니다.
  • 📄 .gitignore: .venv/, data/, __pycache__/를 Git 추적 대상에서 제외합니다.
  • 📄 analysis.py: 에이전트가 분석 코드를 차례로 쌓아 갈 Python 스크립트. 준비 직후에는 비어 있습니다.
  • 📄 pyproject.toml: 이 프로젝트가 직접 사용하는 Python 버전과 패키지를 기록합니다.
  • 📄 uv.lock: 하위 의존성까지 포함한 구체적인 패키지 버전을 고정합니다.

data/는 아직 생기지 않습니다. 다음 섹션의 1단계 프롬프트가 원본 count 파일을 내려받을 때 처음 만들어집니다. 다른 컴퓨터에서는 pyproject.toml과 uv.lock이 있는 프로젝트 루트에서 uv sync를 실행해 같은 환경을 복원할 수 있습니다.

에이전트가 한 번에 모든 분석을 끝내게 하지 않는 것이 중요합니다. 아래 프롬프트를 위에서부터 하나씩 전달하고, 실제 출력과 그래프를 함께 확인한 뒤 다음 단계로 넘어가세요.

3. 프롬프트로 데이터를 한 단계씩 살펴보기

각 프롬프트는 앞 단계에서 만든 analysis.py에 코드를 이어서 작성하고 uv run python analysis.py로 전체 흐름을 다시 검증하게 합니다. 그래프는 outputs/에 저장됩니다.

예시 결과와 비교하는 방법

각 단계 아래에는 GSE251845를 분석한 예시 결과가 있습니다. 그림을 누르면 원본 크기로 열립니다. 에이전트가 Python으로 만든 결과는 패키지 버전, 필터, 그래프 설정에 따라 수치와 모양이 조금 달라질 수 있으므로, 숫자를 그대로 맞추기보다 전체 패턴과 해석이 같은지를 비교하세요.

3-1. raw count와 환자 짝 확인

첫 단계에서는 통계 검정을 하지 않습니다. 파일의 모양과 샘플 이름부터 확인합니다. 이 검사를 통과하지 못하면 이후 그래프가 그럴듯해 보여도 다른 질문을 분석한 것일 수 있습니다.

▶ 1단계 프롬프트: 데이터 다운로드와 구조 확인 복사

GSE251845 데이터 구조를 확인하는 단계만 진행해줘.

1. 프로젝트 루트에 data 디렉터리를 만들어.
2. 아래 GEO 파일을 data/GSE251845_htseq_raw_counts.csv.gz로 내려받아. 파일이 이미 있고 gzip 검증을 통과하면 다시 받지 마.
 https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE251845&file=GSE251845_htseq_raw_counts.csv.gz&format=file
3. analysis.py에 pandas로 파일을 읽는 코드를 작성해. 첫 열을 index로 사용하고, __로 시작하는 HTSeq 특수 카운터를 별도 표로 출력한 뒤 유전자 행렬에서는 제거해. 그다음 샘플이 행이 되도록 전치해.
4. 샘플 이름에서 _htseq.out을 제거하고 대문자로 통일해. 끝의 C는 Tumor, N은 Normal로 해석하고 앞의 숫자는 patient로 저장한 metadata를 만들어.
5. 다음 조건을 assert로 검증해.
 - 샘플은 총 44개
 - Tumor 22개와 Normal 22개
 - 환자는 22명
 - 모든 환자에게 C와 N이 하나씩 존재
 - count matrix와 metadata의 행 순서가 동일
6. 원본 행렬 크기, HTSeq 특수 카운터, 특수 행을 제거하고 전치한 행렬 크기, metadata 앞부분, 조건별 개수, 환자별 샘플 개수를 출력해.
7. uv run python analysis.py로 실행하고 실제 출력이 각 검증 조건과 맞는지 한국어로 설명한 뒤 멈춰.

아직 저발현 유전자 필터링, DESeq2, PCA나 다른 그래프는 실행하지 마.

정상적으로 읽었다면 샘플 축은 44가 됩니다. 원본 파일은 유전자가 행이고 샘플이 열이지만, PyDESeq2 입력은 샘플이 행이고 유전자가 열이므로 전치가 필요합니다. count 값은 TPM처럼 정규화된 값이 아니라 정수 raw count입니다.

파일 끝의 __no_feature, __ambiguous, __alignment_not_unique 등은 유전자가 아니라 어느 유전자에도 최종 배정되지 않은 read의 이유별 집계입니다. DESeq2 행렬에서는 제외하지만, 제거하기 전에 샘플별 값을 확인하면 정렬이나 annotation 문제를 찾는 단서가 됩니다. 각 행의 뜻은 HTSeq raw count 레퍼런스에서 확인하세요.

여기서 가장 중요한 출력은 환자별 샘플 수입니다. 모든 환자에게 정확히 두 샘플이 있어야 patient 효과와 condition 효과를 분리할 수 있습니다. 짝이 하나라도 빠지거나 35c가 소문자라는 이유로 Tumor에서 누락되면 paired design이 깨집니다.

3-2. 전체 발현 패턴을 PCA로 보기

PCA는 개별 유전자보다 먼저 샘플 전체 구조를 확인하는 분석입니다. raw count에는 시퀀싱 깊이에 따른 분산 차이가 있으므로, PyDESeq2의 VST로 값을 변환한 뒤 PCA를 계산합니다.

PCA를 벡터 관점에서 보면

샘플 하나는 [gene_1, gene_2, ..., gene_24116]처럼 유전자마다 VST 발현값 하나를 가진 큰 벡터입니다. 각 유전자가 벡터의 원소 하나입니다.

▶ 2단계 프롬프트: paired DESeq2 모델과 PCA 복사

이전 데이터 검증 코드를 유지하고 paired DESeq2 모델과 PCA 단계만 추가해줘.

1. 44개 중 적어도 3개 샘플에서 count가 10 이상인 유전자만 유지해.
2. PyDESeq2의 DeseqDataSet을 design="~ patient + condition"으로 만들고 n_cpus=1로 실행해.
3. condition의 기준 범주는 Normal이 되게 하고, Tumor 대 Normal contrast로 DeseqStats를 실행해.
4. 결과표를 padj 오름차순으로 정렬해 outputs/GSE251845_tumor_vs_normal_PyDESeq2.csv로 저장해.
5. VST 변환값으로 PCA를 계산해. 점 하나는 샘플 하나로 그리고 Normal은 청록색, Tumor는 주황색으로 표시해.
6. matplotlib의 비대화형 backend를 사용하고 그래프를 outputs/pca.png에 160 dpi로 저장해.
7. PC1과 PC2의 설명 분산, 각 조건의 PC1 범위, PC2에서 가장 멀리 떨어진 샘플 이름을 출력해.
8. uv run python analysis.py로 전체 스크립트를 실행해. outputs/pca.png를 직접 확인하고, 분리 방향과 이상치 후보를 한국어로 설명한 뒤 멈춰.

PCA 다음 단계의 MA plot이나 볼케이노 플롯은 아직 만들지 마.

이 실습 프롬프트를 실행해 생성한 PCA. 필터를 통과한 유전자 24,116개를 사용했으며 PC1은 34.03%, PC2는 9.46%의 분산을 설명합니다.

점 하나는 조직 샘플 하나이며, 같은 환자의 Normal과 Tumor도 서로 다른 두 점으로 표시됩니다. Normal은 PC1 오른쪽, Tumor는 왼쪽에 모이고 두 조건의 PC1 범위가 겹치지 않습니다. PCA는 조건 정보를 사용하지 않고 축을 찾는 비지도 분석이므로, 가장 큰 변화 방향이 종양 여부와 강하게 대응한다고 해석할 수 있습니다. PC1은 44개 샘플 사이에서 관찰된 전체 발현 차이의 34.03%를 한 축으로 보여 줍니다. 이 값은 종양 여부의 인과 효과나 분류 정확도가 아닙니다.

PC1의 좌우 방향에는 고정된 생물학적 의미가 없습니다. 계산 결과의 부호가 뒤집혀 Normal이 왼쪽에 나타나도 샘플 사이의 거리와 분리 정도는 같습니다. 중요한 것은 어느 쪽이 왼쪽인지가 아니라 두 조건이 얼마나 분리되는지입니다.

PC2는 PC1에 나타나지 않은 샘플 간 차이 중 다음으로 큰 방향을 보여 주며, 전체 발현 차이의 9.46%를 설명합니다. 이 그림에서는 PC1이 조건을 주로 분리하므로 PC2에서 같은 조건 안의 표본 차이가 잘 드러납니다. 50N은 PC2 방향으로 다른 Normal보다 멀리 떨어져 있어 확인이 필요한 표본입니다. 생물학적 개인차, 조직 구성, 배치 또는 품질 차이일 수 있으므로 원시 read 품질, library size와 조직 위치를 확인하기 전에는 이상치로 확정하거나 제거하지 않습니다.

이 그림에는 Tumor와 Normal의 차이뿐 아니라 환자마다 원래 다른 발현 특성도 함께 나타납니다. 환자별 개인차를 고려하는 paired design은 PCA가 아니라 뒤의 DEG 검정에 적용됩니다.

3-3. DEG 결과표를 한 행씩 읽기

같은 실행에서 paired DESeq2 검정도 완료되어 outputs/GSE251845_tumor_vs_normal_PyDESeq2.csv에 유전자별 결과표가 저장됩니다. PCA가 44개 샘플의 전체 구조를 보여 줬다면, 이 표는 각 유전자에서 Tumor와 Normal이 얼마나 다른지 수치로 보여 줍니다.

GSE251845 DEG 결과표 예시. 가로로 긴 표이므로 그림을 눌러 크게 보세요. symbol, entrez_id, regulation은 결과를 읽기 쉽게 나중에 추가한 열입니다.

표의 행 하나는 유전자 하나입니다. PyDESeq2가 직접 만드는 핵심 결과 여섯 열은 baseMean, log2FoldChange, lfcSE, stat, pvalue, padj입니다. 한꺼번에 외우지 말고 첫 행의 ETV4부터 읽어 보겠습니다.

열	ETV4의 값	먼저 읽을 뜻
baseMean	약 1,431	분석한 샘플 전체에서 어느 정도 관측됐는가
log2FoldChange	약 6.23	어느 조건에서 얼마나 높았는가
lfcSE	약 0.29	추정한 변화량이 얼마나 불확실한가
stat	약 21.67	변화량이 불확실성보다 얼마나 큰가
pvalue	약 $3.9 \times 10^{-104}$	ETV4 하나만 검사했을 때 통계적 근거가 얼마나 강한가
padj	약 $1.1 \times 10^{-99}$	모든 유전자를 함께 검사한 뒤에도 근거가 남는가

gene_id의 ENSG00000175832는 Ensembl ID이고 ETV4는 사람이 읽기 쉬운 유전자 이름입니다. symbol, entrez_id, regulation은 PyDESeq2 원본 결과가 아니라 유전자 이름과 변화 방향을 알아보기 쉽게 나중에 붙인 열입니다.

ETV4의 baseMean은 약 1,431로, 44개 샘플의 시퀀싱 깊이를 보정한 count를 평균한 관측 규모입니다. log2FoldChange는 약 6.23이며 이 분석의 비교 방향이 Tumor - Normal이므로, 모델은 ETV4가 Tumor에서 약 $2^{6.23}\approx75$배 높다고 추정했습니다.

lfcSE는 약 0.29이고 stat은 6.23 ÷ 0.29 ≈ 21.67입니다. 추정한 변화량이 표준오차의 약 21.7배이며, 모든 유전자를 함께 검사한 뒤의 padj도 약 $1.1 \times 10^{-99}$로 매우 작습니다. 따라서 이 데이터에서 ETV4의 증가 방향과 통계적 근거는 강하지만, 이 결과만으로 ETV4가 암의 원인이라고 결론 내릴 수는 없습니다.

baseMean부터 padj까지의 계산 관계, 표준오차와 표준편차의 차이, NaN과 LFC shrinkage는 DESeq2·PyDESeq2 결과표 읽기에서 자세히 설명합니다.

MA plot으로 결과표 전체 보기

결과표에는 유전자가 24,116줄 있습니다. 한 줄씩 읽기 어려우므로 유전자 하나를 점 하나로 바꾸어 한 화면에 그린 그림이 MA plot입니다. 먼저 점 전체의 모양을 확인하고, 그다음 눈에 띄는 유전자를 찾습니다.

▶ 3단계 프롬프트: MA plot 복사

이전 분석을 유지하고 DESeq2 결과로 MA plot을 만드는 단계만 추가해줘.

1. 점 하나가 유전자 하나가 되게 해.
2. x축은 log10(baseMean + 1), y축은 log2FoldChange로 사용해.
3. padj가 0.05 미만인 유전자는 주황색, 나머지는 회색으로 표시해.
4. y=0 기준선을 그리고, 극단값 때문에 중심부가 보이지 않지 않도록 표시 범위를 -5에서 5로 제한해. 범위 밖 유전자의 개수는 별도로 출력해.
5. 그래프를 outputs/ma-plot.png에 160 dpi로 저장해.
6. padj가 0.05 미만인 유전자 수, 양의 log2FC와 음의 log2FC 개수를 출력해.
7. uv run python analysis.py로 실행하고 outputs/ma-plot.png를 확인해. 저발현 구간의 퍼짐과 0선 위아래의 의미를 한국어로 설명한 뒤 멈춰.

아직 변화량 기준으로 후보를 고르거나 볼케이노 플롯을 만들지 마.

현재 실험의 outputs/ma-plot.png. 필터를 통과한 유전자 24,116개를 그렸으며 세로축은 -5에서 5까지만 표시했습니다.

점 하나가 유전자 하나입니다. 가로축은 baseMean이므로 오른쪽일수록 44개 샘플에서 전반적으로 count가 큰 유전자입니다. 세로축은 log2FoldChange입니다. 0보다 위에 있으면 Tumor에서 높고, 아래에 있으면 Normal에서 높습니다. log2FoldChange = 1은 Tumor가 약 2배, -1은 Tumor가 Normal의 약 절반이라는 뜻입니다.

먼저 전체 모양을 보는 이유

이 그림에서 가장 먼저 볼 것은 우상단이나 우하단의 특정 점이 아니라 점들이 모여 만든 전체 모양입니다.

왼쪽은 count가 작은 유전자입니다. 원래 값이 작으면 몇 count만 달라져도 비율은 크게 바뀌므로 점들이 위아래로 넓게 퍼집니다. 오른쪽은 count가 큰 유전자입니다. 몇 count의 차이에 덜 흔들리므로 점들이 상대적으로 0선 가까이에 좁게 모입니다. 현재 그림처럼 왼쪽이 넓고 오른쪽이 좁아지는 모양은 count 데이터에서 흔히 볼 수 있습니다.

또 점들이 대체로 0선 위와 아래에 모두 있는지도 봅니다. 이 과정은 특정 유전자를 고르기 전에, 발현량이 작은 유전자에서 변화량이 유난히 불안정하지 않은지와 결과 전체가 한 방향으로 심하게 치우치지 않았는지를 확인하는 단계입니다.

그다음 눈에 띄는 유전자를 찾기

전체 모양을 확인한 뒤에는 점의 위치를 다음처럼 읽습니다.

점의 위치	쉬운 해석
오른쪽 위	count가 전반적으로 크고 Tumor에서 더 높음
오른쪽 아래	count가 전반적으로 크고 Normal에서 더 높음
오른쪽의 0선 주변	count는 크지만 Tumor와 Normal의 차이는 작음
왼쪽 위·아래	변화는 커 보이지만, count가 작아서 비율이 크게 흔들렸을 수 있음

따라서 우상단과 우하단은 후속 검토할 후보를 찾기 좋은 영역입니다. 다만 그곳에 있다는 이유만으로 바로 중요한 유전자라고 결론 내리지는 않습니다. 주황색은 다중검정 보정 뒤에도 padj < 0.05인 유전자입니다. 현재 결과에서는 15,099개가 주황색이고, Tumor에서 증가한 유전자는 7,302개, 감소한 유전자는 7,797개입니다.

주황색이라고 모두 변화가 큰 것도 아닙니다. 여러 샘플에서 작은 차이가 꾸준히 나타나면 0선 가까이에 있어도 통계적으로 유의할 수 있습니다. 반대로 0선에서 멀어 보여도 소수 샘플 때문에 생긴 불안정한 값일 수 있습니다. 그래서 padj, log2FoldChange, 샘플별 count를 함께 확인해야 합니다.

실제로 세로축 표시 범위인 -5에서 5를 벗어난 유전자가 237개 있습니다. 그중 ENSG00000184811은 log2FoldChange = -27.98로 매우 커 보이지만, lfcSE = 216.88, padj = 0.918입니다. 44개 중 21개 샘플의 raw count가 0이고 값이 50N과 47N 두 샘플에 몰려 있어 추정이 매우 불안정합니다. 0선에서 멀리 떨어진 점만 골라서는 안 되는 실제 사례입니다.

정리하면 MA plot은 ① 전체 점의 모양 확인, ② 우상단과 우하단의 눈에 띄는 점 확인, ③ padj와 샘플별 count로 다시 검토하는 순서로 읽습니다. 변화량과 통계적 근거를 중심으로 후보를 좁히는 작업은 다음 볼케이노 플롯에서 이어집니다.

3-4. 볼케이노 플롯으로 후보 범위 좁히기

볼케이노 플롯은 평균 발현량 대신 변화량과 통계적 근거를 두 축에 놓습니다. MA plot에서 전체 데이터의 모양을 확인한 뒤, 후속 검토할 후보를 좁힐 때 사용합니다.

▶ 4단계 프롬프트: 볼케이노 플롯과 DEG 후보 복사

이전 분석을 유지하고 볼케이노 플롯과 이 실습의 기준을 통과한 DEG 후보를 만드는 단계만 추가해줘.

1. padj와 log2FoldChange가 결측이 아닌 유전자만 사용해.
2. g:Profiler의 identifier conversion을 사용해 결과의 Ensembl gene ID를 human gene symbol로 변환해. gprofiler-official이 없으면 uv add gprofiler-official로 설치해. 변환 결과의 Ensembl ID와 symbol을 outputs/ensembl-to-symbol.csv에 저장하고, 변환되지 않은 ID의 수도 출력해.
3. DESeq2 결과에는 Ensembl ID를 gene_id로 그대로 보존하고 symbol 열을 추가해. 하나의 Ensembl ID에 변환 결과가 여러 개면 첫 값을 조용히 고르지 말고 중복 내용을 출력해 확인해. symbol이 없는 유전자는 그림에서만 Ensembl ID를 대신 사용해.
4. x축은 log2FoldChange, y축은 -log10(padj)로 계산해. padj가 0이면 float의 가장 작은 양수로 제한해 로그 오류를 막아.
5. padj < 0.05이면서 log2FoldChange >= 1인 유전자는 빨간색, log2FoldChange <= -1인 유전자는 파란색, 나머지는 회색으로 표시해.
6. x=-1과 x=1, y=-log10(0.05)에 점선을 그어 선택 기준을 보이게 해.
7. Tumor-up과 Normal-up에서 padj가 가장 작은 유전자 3개씩, 그리고 abs(log2FoldChange)가 10보다 큰 극단값에만 라벨을 붙여. 라벨은 Ensembl ID가 아니라 symbol을 우선 사용하고, 같은 유전자를 중복 표시하지 마. symbol이 없을 때만 Ensembl ID를 표시해.
8. 그래프를 outputs/volcano-plot.png에 160 dpi로 저장해. 라벨끼리 겹치거나 그림 밖으로 잘리지 않는지 확인해.
9. 선택된 전체 후보 수, Tumor에서 증가한 후보 수, Normal에서 증가한 후보 수를 출력해. 각 방향에서 padj가 가장 작은 10개 유전자는 symbol과 Ensembl ID를 함께 출력해.
10. abs(log2FoldChange)가 10보다 큰 극단값은 symbol, Ensembl ID, baseMean, log2FoldChange, lfcSE, padj와 원래 count를 별도로 출력해 저발현 불안정성 여부를 확인해.
11. uv run python analysis.py로 실행하고 outputs/volcano-plot.png를 확인해. 네 영역과 극단값을 한국어로 설명한 뒤 멈춰.

이 기준이 보편적인 생물학적 진실이라고 표현하지 마. 이 실습에서 정한 DEG 후보 선택 기준이라고 명시해.

현재 실험의 outputs/volcano-plot.png. padj < 0.05와 |log2FoldChange| >= 1을 이 실습의 DEG 후보 선택 기준으로 표시했습니다.

점 하나는 유전자 하나입니다. 가로축의 오른쪽은 Tumor에서 높은 유전자, 왼쪽은 Normal에서 높은 유전자입니다. 0에서 좌우로 멀어질수록 두 조건의 차이가 큽니다. 세로축은 padj를 -log10으로 바꾼 값이므로, 위로 갈수록 padj가 작고 통계적 근거가 강합니다.

그림의 영역은 다음 순서로 읽을 수 있습니다.

위치	이번 그림에서의 의미
오른쪽 위의 빨간 점	Tumor에서 2배 이상 높고 padj < 0.05인 DEG 후보
왼쪽 위의 파란 점	Normal에서 2배 이상 높고 padj < 0.05인 DEG 후보
위쪽 중앙의 회색 점	차이는 2배보다 작지만 여러 샘플에서 일관될 수 있음
아래쪽의 회색 점	차이가 커 보여도 통계적 근거가 충분하지 않음

세로 점선 -1과 1은 약 2배 변화 기준이고, 가로 점선은 padj = 0.05입니다. 이 기준으로 24,116개 유전자 가운데 8,256개가 선택됐습니다. Tumor에서 증가한 후보는 3,631개, Normal에서 증가한 후보는 4,625개입니다. 빨간색과 파란색은 이 연습용 기준을 통과했다는 뜻이지, 암의 원인이나 좋은 바이오마커로 확인됐다는 뜻은 아닙니다.

오른쪽 위에서 가장 높은 빨간 점은 앞에서 한 행씩 살펴본 ETV4(ENSG00000175832)입니다. log2FoldChange = 6.23, padj는 약 $1.5 \times 10^{-99}$로, Tumor에서 크게 증가했고 통계적 근거도 강하다는 두 조건을 함께 만족합니다.

좌우로 멀리 떨어졌다고 모두 믿을 만한 후보인 것은 아닙니다. 가장 왼쪽의 TRARG1(ENSG00000184811)은 log2FoldChange = -27.98로 차이가 매우 커 보이지만, 세로축 아래쪽의 회색 점입니다. padj = 0.918이고 44개 중 21개 샘플의 raw count가 0이며 값이 두 Normal 샘플에 몰려 있어 추정이 불안정합니다.

반면 |log2FoldChange| > 10인 OTOP2, CA1, TMIGD1은 왼쪽 위의 파란 점입니다. 변화량이 크다는 공통점은 같지만 lfcSE가 약 0.63에서 0.79이고 padj도 매우 작습니다. 따라서 극단값은 무조건 버리거나 선택하지 않고, lfcSE, padj, 환자별 count를 확인해 구분해야 합니다.

볼케이노 플롯은 ① 좌우 거리로 변화량 확인, ② 높이로 통계적 근거 확인, ③ 두 기준을 통과한 후보의 환자별 count 재확인 순서로 읽습니다. 점선은 분석 목적과 후속 검증 비용에 따라 바꿀 수 있는 선택 규칙입니다.

3-5. 상위 유전자의 환자별 움직임 확인

DESeq2 결과표의 한 행은 22명에게서 나타난 차이를 하나의 log2FoldChange로 요약합니다. 이번에는 가장 작은 padj를 가진 유전자 하나를 골라, 요약되기 전의 환자별 값을 다시 펼쳐 봅니다. 이 그림으로 증가 방향이 여러 환자에게 반복되는지, 일부 환자의 큰 값에만 끌려간 결과인지를 확인합니다.

▶ 5단계 프롬프트: 상위 DEG의 paired plot 복사

이전 분석을 유지하고 상위 DEG의 환자별 paired plot을 만드는 단계만 추가해줘.

1. padj가 가장 작은 유전자 하나를 top_gene으로 선택해.
2. DESeq2의 정규화 count에서 top_gene의 값을 가져오고 metadata의 patient, condition과 결합해.
3. x축은 Normal과 Tumor, y축은 normalized count로 사용해. 같은 환자의 두 점을 회색 선으로 연결하고 Normal과 Tumor 점의 색을 다르게 표시해.
4. 그래프 제목에 gene symbol과 Ensembl ID를 모두 표시하고 outputs/top-gene-paired.png에 160 dpi로 저장해.
5. 환자마다 Tumor - Normal 차이를 계산해 증가한 환자 수, 감소한 환자 수, 차이가 가장 큰 환자와 가장 작은 환자를 출력해. 두 극단 환자와 Tumor에서 증가하지 않은 환자가 있다면 그림에 patient 번호를 표시해.
6. top_gene의 baseMean, log2FoldChange, pvalue, padj도 함께 출력해.
7. uv run python analysis.py로 실행하고 outputs/top-gene-paired.png를 확인해. 선의 방향이 환자 전반에서 반복되는지, 예외 환자가 있는지 한국어로 설명한 뒤 멈춰.

이 그림만으로 해당 유전자가 암의 원인이라고 결론 내리지 마.

현재 결과에서 padj가 가장 작은 유전자는 ETV4(ENSG00000175832)입니다. 그림은 다음처럼 읽습니다.

그림 요소	뜻
왼쪽 청록색 점	한 환자의 Normal 조직에서 측정한 ETV4 normalized count
오른쪽 주황색 점	같은 환자의 Tumor 조직에서 측정한 ETV4 normalized count
두 점을 잇는 회색 선	환자 한 명의 Normal과 Tumor 한 쌍
위로 향하는 선	그 환자에서는 Tumor의 ETV4가 더 높음
아래로 향하는 선	그 환자에서는 Normal의 ETV4가 더 높음

normalized count는 샘플마다 다른 시퀀싱 깊이를 보정한 count입니다. 회색 선 하나를 따라 왼쪽에서 오른쪽으로 읽으면 같은 환자 안에서 ETV4가 얼마나 달라졌는지 볼 수 있습니다.

이번 결과에서는 22명 모두 회색 선이 위로 향했습니다. Tumor에서 감소한 환자와 차이가 없는 환자는 각각 0명입니다. 따라서 ETV4 증가는 한두 명의 극단값 때문에 생긴 것이 아니라, 이 데이터의 모든 환자에게 같은 방향으로 나타났습니다.

다만 증가한 크기는 서로 다릅니다. 환자 32는 Tumor와 Normal의 normalized count 차이가 약 4,631로 가장 컸고, 환자 33은 약 56으로 가장 작았습니다. 그림의 32와 33 표시는 이 두 환자를 가리킵니다. 이 차이는 두 count를 뺀 절대 차이이며, 배수 변화인 log2FoldChange와는 다른 값입니다.

Normal 점들이 바닥에 뭉쳐 보이는 것은 모두 정확히 0이어서가 아닙니다. 세로축이 로그 축이 아니고 Tumor 값이 약 4,700까지 올라가므로, 상대적으로 작은 Normal 값 사이의 차이가 눌려 보입니다. 이 그림은 환자별 방향을 확인하기에는 좋지만 padj와 추정 불확실성은 보여 주지 않습니다. 그 판단에는 paired DESeq2 결과표를 함께 사용합니다.

이 그림은 연관성의 일관성을 보여 줄 뿐 인과관계를 증명하지 않습니다. Tumor에서 발현이 증가한 유전자가 암을 일으켰을 수도 있지만, 암이 생긴 결과로 증가했거나 조직의 세포 구성 차이를 반영했을 수도 있습니다. 유전자 이름으로 바꿀 때는 분석에 사용한 genome annotation 버전도 함께 고정해야 합니다.

4. 그래프를 모두 본 뒤 답해 볼 질문

1. design = ~ condition이 아니라 ~ patient + condition을 사용한 이유는 무엇인가?
2. PCA에서 PC1의 좌우가 뒤집혀도 해석이 달라지지 않는 이유는 무엇인가?
3. MA plot의 저발현 구간에서 fold change가 더 넓게 퍼지는 이유는 무엇인가?
4. padj가 작지만 log2FoldChange가 작은 유전자는 어떤 의미인가?
5. 상위 DEG가 여러 환자에게 반복되어도 암의 원인이라고 바로 말할 수 없는 이유는 무엇인가?

이제 후보 유전자 목록을 하나씩 읽는 대신, 공통 기능과 pathway로 묶어 볼 차례입니다. LOCRC 실습 2: GO와 GSEA로 기능·경로 찾기에서는 이 페이지가 만든 DEG 목록과 전체 유전자 순위를 새 입력으로 사용합니다.