CODE BENCH · NGS PIPELINES

deeptools — 정렬된 리드를 신호 트랙으로: bamCoverage·computeMatrix·plotHeatmap

TL;DR — deeptools는 정렬된 리드(BAM)를 유전체 위의 신호 트랙(bigWig)으로 바꾸고, 그 신호를 특정 구간 주변에서 메타플롯·히트맵으로 시각화하는 파이썬 도구 모음입니다. ChIP-seq·ATAC-seq·CUT&RUN처럼 '어디에 신호가 얼마나 쌓였나'를 봐야 하는 분석의 표준 후처리죠.

뼈대는 네 갈래입니다. bamCoverage로 BAM을 정규화된 bigWig 커버리지 트랙으로 만들고, multiBamSummary + plotCorrelation·plotPCA로 반복 실험(replicate)의 재현성을 확인하고, bamCompare로 ChIP 대 input을 log2 비율로 견주고, computeMatrix → plotHeatmap·plotProfile로 TSS 같은 기준점 둘레의 신호를 히트맵과 평균 곡선으로 그립니다.

작은 BAM 예제로 명령을 한 줄씩 따라간 뒤, 정렬된 BAM에서 bamCoverage로 트랙을 만들고 computeMatrix로 행렬을 뽑아 plotHeatmap까지 잇는 전체 워크플로를 완성합니다. 0-based 좌표와 bigWig 포맷도 함께 짚죠. 명령은 deeptools 3.5 기준입니다.

🔗 관련 글 · deeptools가 신호로 바꾸는 리드는 대개 앞선 도구의 산출물입니다 — 정렬 결과의 포맷은 SAM/BAM 파일이란?, 그 BAM을 정렬·인덱싱하는 samtools·bcftools 실전, 신호 트랙의 대표 사용처인 ChIP-seq, 그리고 같은 방식으로 트랙을 읽는 또 다른 후성유전체 분석 Hi-C·3D 유전체에서 이어집니다.

이 글에서 만들 것

ChIP-seq나 ATAC-seq를 돌리면 결국 손에 남는 건 정렬된 리드가 담긴 BAM 파일입니다. 그런데 BAM을 IGV로 열어 봐도 리드 하나하나가 늘어서 있을 뿐, "이 프로모터에 신호가 얼마나 몰렸나", "샘플 두 개가 서로 닮았나" 같은 질문에는 바로 답이 안 나오죠. 리드를 세어 트랙으로 요약하고, 그 트랙을 관심 구간에 맞춰 정렬해 그림으로 만들어야 비로소 읽힙니다.

deeptools는 그 요약과 시각화를 전담합니다. BAM을 유전체 커버리지 트랙(bigWig)으로 바꾸고(bamCoverage), 여러 샘플의 커버리지를 모아 상관·주성분으로 재현성을 보고(multiBamSummary·plotCorrelation·plotPCA), ChIP과 input을 비율로 견주고(bamCompare), 유전자 시작점(TSS) 같은 기준점 주변의 신호를 히트맵·메타플롯으로 펼칩니다(computeMatrix·plotHeatmap·plotProfile). 입력은 BAM과 bigWig, 출력은 bigWig·PNG·npz라 앞뒤 도구와 자연스럽게 이어지죠.

이 글은 그 deeptools를 처음부터 익힙니다. 구체적으로 ① 설치와 예제 BAM을 준비하고, ② bamCoverage로 정규화된 신호 트랙을 만들고, ③ multiBamSummary로 샘플 간 재현성을 확인하고, ④ bamCompare로 ChIP 대 input을 비교하고, ⑤ computeMatrix로 TSS 주변 신호를 히트맵·메타플롯으로 그린 뒤, ⑥ 이 조각들을 하나의 워크플로로 묶습니다.

# 핵심 흐름 한눈에
bamCoverage -b aln.sorted.bam -o signal.bw --normalizeUsing RPKM   # ① BAM → bigWig 트랙
computeMatrix reference-point -S signal.bw -R genes.bed -o mat.gz  # ② TSS 주변 신호 행렬
plotHeatmap -m mat.gz -o heatmap.png                              # ③ 히트맵·메타플롯

1. 사전 준비 — 설치와 예제 데이터

deeptools는 파이썬으로 짠 명령줄 도구 모음입니다. conda로 까는 게 가장 간편한데, 의존성 해결을 위해 conda-forge와 bioconda 채널을 함께 지정합니다.

# conda — conda-forge·bioconda 두 채널을 함께(권장)
conda install -c conda-forge -c bioconda deeptools

# 또는 pip
pip install deeptools

# 설치 확인(버전 출력)
deeptools --version

deeptools 3.5.6

deeptools의 입력은 좌표순으로 정렬된 BAM이고, 대부분의 명령은 그 옆에 인덱스 파일(.bai)이 있어야 돕니다. 정렬·인덱싱은 samtools가 맡습니다.

# 좌표순 정렬 → 인덱싱(.bai 생성) — deeptools가 요구
samtools sort aln.bam -o aln.sorted.bam
samtools index aln.sorted.bam

예제 데이터는 두 갈래로 구할 수 있습니다. 하나는 deeptools 저장소가 함께 배포하는 테스트 파일(deeptools/test/test_data/의 test1.bam·testA.bw·test.bed3 등)이고, 다른 하나는 ENCODE 포털에서 내려받는 공개 ChIP-seq BAM·bigWig입니다. 아래 예제는 이런 작은 BAM 한두 개가 있다고 치고 진행합니다.

출력으로 자주 나오는 bigWig는 유전체 전체에 걸친 연속적인 수치(커버리지·점수)를 담는 색인된 이진 포맷입니다. 텍스트인 BedGraph와 달리 압축·색인이 되어 있어 IGV·UCSC 같은 브라우저가 특정 구간만 빠르게 꺼내 그릴 수 있죠. 좌표는 BED와 마찬가지로 0-based 반열림이라, 다른 포맷과 맞댈 때는 좌표계를 먼저 확인해야 합니다.

2. bamCoverage — BAM을 bigWig 신호 트랙으로

가장 먼저 쓰는 명령이 bamCoverage입니다. 유전체를 일정 크기의 구간(bin)으로 잘라 각 bin에 겹치는 리드 수를 세고, 이를 정규화해 bigWig 트랙으로 내보냅니다.

# BAM → bigWig, 10bp bin, RPKM 정규화, 리드를 fragment 길이로 확장
bamCoverage -b aln.sorted.bam -o signal.bw \
  --binSize 10 --normalizeUsing RPKM --extendReads

--binSize는 해상도를 정합니다(기본 50bp). 작게 잡으면 신호가 더 촘촘해지지만 파일이 커지죠. --extendReads는 단일-끝(single-end) 리드를 실제 fragment 길이만큼 늘려, 리드가 짧아 생기는 신호 끊김을 메웁니다. ChIP처럼 fragment 전체가 신호인 실험에서 특히 중요합니다.

결과가 이진 포맷인 bigWig라 눈으로 바로 보긴 어려우니, 확인용으로 같은 계산을 텍스트인 BedGraph로 한 번 뽑아 봅니다.

# 확인용 — BedGraph로 출력해 앞 몇 줄 들여다보기
bamCoverage -b aln.sorted.bam -o signal.bedgraph -of bedgraph \
  --binSize 50 --normalizeUsing CPM
head -3 signal.bedgraph

chr1	0	50	0.00
chr1	50	100	1.83
chr1	100	150	4.57

각 줄은 '이 50bp 구간의 정규화된 커버리지가 이 값'이라는 뜻입니다. 정규화 방식은 --normalizeUsing으로 고르는데, 실험·목적에 따라 갈립니다.

RPGC (1x 정규화)만은 매핑 가능한 유전체 크기를 알려 줘야 합니다. --effectiveGenomeSize에 값을 넘기는데, 이 값은 리드 길이에 따라 달라집니다. 예컨대 사람 GRCh38·리드 100bp 기준은 2,805,636,231이고, 정확한 표는 deeptools 공식 문서에 정리돼 있습니다.

# RPGC 1x 정규화 — 유효 유전체 크기 필요(GRCh38·100bp)
bamCoverage -b aln.sorted.bam -o signal_1x.bw \
  --normalizeUsing RPGC --effectiveGenomeSize 2805636231 --binSize 10

3. multiBamSummary — 샘플 간 재현성(상관·PCA)

트랙을 만들었다면, 그림을 그리기 전에 반복 실험이 서로 닮았는지부터 확인하는 게 순서입니다. multiBamSummary는 여러 BAM의 커버리지를 한꺼번에 요약해 압축 배열(.npz)로 저장합니다. 유전체를 일정 bin으로 훑는 bins 모드가 기본이고, 특정 구간만 볼 땐 BED-file 모드를 씁니다.

# 여러 BAM의 커버리지를 bin 단위로 요약 → .npz
multiBamSummary bins \
  -b chip_rep1.bam chip_rep2.bam input.bam \
  --labels ChIP_rep1 ChIP_rep2 input \
  -o readCounts.npz

이 .npz를 plotCorrelation에 넘기면 샘플 간 상관 히트맵이, plotPCA에 넘기면 주성분 그림이 나옵니다.

# 스피어만 상관을 히트맵으로(숫자도 표시)
plotCorrelation -in readCounts.npz \
  --corMethod spearman --whatToPlot heatmap \
  --plotNumbers -o correlation.png

# 주성분 분석(PCA)
plotPCA -in readCounts.npz -o pca.png

잘 된 실험이라면 같은 조건의 반복끼리(ChIP_rep1·ChIP_rep2) 상관이 높게 뭉치고, 배경인 input과는 뚜렷이 갈립니다. PCA에서도 반복은 가까이, 조건이 다른 샘플은 멀찍이 떨어지죠. 여기서 한 반복만 엉뚱하게 튄다면, 히트맵을 그리기 전에 그 샘플부터 들여다보는 게 낫습니다. --corMethod는 pearson으로, --whatToPlot은 scatterplot으로도 바꿀 수 있습니다.

4. bamCompare — ChIP 대 input을 log2 비율로

ChIP-seq 신호에는 진짜 결합 신호뿐 아니라 열린 염색질·매핑 편향 같은 배경이 섞입니다. 이 배경을 걷어 내려면 대조군(input)과 견줘야 하는데, 그 일을 bamCompare가 합니다. 두 BAM을 같은 bin에서 비교해 비율 트랙 하나로 내보내죠.

# ChIP 대 input의 log2 비율 트랙
bamCompare -b1 chip.bam -b2 input.bam \
  --operation log2 --binSize 20 \
  -o log2_chip_vs_input.bw

--operation은 기본이 log2(두 값의 log2 비율)이고, ratio·subtract·mean 등으로도 바꿀 수 있습니다. 깊이가 다른 두 BAM을 공정하게 맞추기 위해 bamCompare는 먼저 스케일 인자를 계산하는데, 그 방식은 --scaleFactorsMethod로 고릅니다(기본 readCount, ChIP 농축이 강할 땐 SES). 출력 트랙에서 값이 0보다 크면 ChIP이 input보다 높은 구간, 곧 신호가 실제로 농축된 자리입니다.

5. computeMatrix → plotHeatmap·plotProfile — 기준점 둘레의 신호

이제 핵심입니다. 만든 트랙(bigWig)을 관심 구간 집합(BED/GTF) 둘레에서 정렬해, 히트맵과 평균 곡선으로 펼치는 단계죠. 먼저 computeMatrix가 신호 행렬을 만들고, 그 행렬을 plotHeatmap·plotProfile이 그림으로 바꿉니다. computeMatrix는 두 가지 모드로 씁니다.

reference-point 모드는 각 구간의 한 지점(TSS·TES·center)을 기준으로 위아래 일정 구간을 잘라 냅니다. 유전자 시작점(TSS) 둘레의 신호를 볼 때 씁니다.

# TSS 기준 ±3kb 구간의 신호 행렬
computeMatrix reference-point --referencePoint TSS \
  -S signal.bw -R genes.bed \
  -b 3000 -a 3000 --skipZeros \
  -o matrix_TSS.gz

-S는 점수가 담긴 bigWig, -R는 구간이 담긴 BED/GTF, -o는 gzip으로 압축된 행렬입니다. -b(before)와 -a(after)는 기준점 위·아래로 몇 bp를 볼지 정하고, --skipZeros는 신호가 아예 없는 구간을 빼 그림을 깔끔하게 합니다.

scale-regions 모드는 길이가 제각각인 구간(유전자 몸통)을 같은 길이로 늘이거나 줄여 맞춥니다. 유전자 전체(TSS→TES)에 걸친 신호 모양을 볼 때 유용하죠.

# 유전자 몸통을 5kb로 스케일 + 양옆 3kb
computeMatrix scale-regions \
  -S signal.bw -R genes.bed \
  --regionBodyLength 5000 -b 3000 -a 3000 \
  -o matrix_body.gz

행렬이 준비되면 두 가지로 그립니다. plotProfile은 모든 구간의 신호를 위치별로 평균 낸 메타플롯(평균 곡선)을, plotHeatmap은 구간을 행으로, 위치를 열로 놓은 히트맵을 그립니다.

# 메타플롯(평균 신호 곡선)
plotProfile -m matrix_TSS.gz -o profile_TSS.png --perGroup

# 히트맵(구간 × 위치) + 위쪽에 메타플롯
plotHeatmap -m matrix_TSS.gz -o heatmap_TSS.png \
  --colorMap Blues --sortUsing sum

plotProfile에서 곡선이 TSS 자리에서 봉우리를 이루면, 그 신호가 전사 시작점에 몰려 있다는 뜻입니다. plotHeatmap은 같은 정보를 구간별로 펼쳐, 어떤 유전자에서 신호가 세고 약한지를 한눈에 보여 주죠. 히트맵 행은 --sortUsing으로 신호 세기순으로 정렬하고, --kmeans를 주면 신호 패턴이 비슷한 구간끼리 자동으로 묶어 줍니다.

6. 실전 — BAM에서 히트맵까지 한 번에

이제 조각을 이어 붙입니다. 정렬된 ChIP BAM에서 시작해 트랙을 만들고, TSS 주변 신호를 뽑아 히트맵까지 가는 전체 흐름입니다.

# 1) 정렬·인덱싱(아직 안 했다면)
samtools sort chip.bam -o chip.sorted.bam
samtools index chip.sorted.bam

# 2) BAM → 정규화된 bigWig 트랙
bamCoverage -b chip.sorted.bam -o chip.bw \
  --binSize 10 --normalizeUsing RPKM --extendReads

# 3) TSS ±3kb 신호 행렬
computeMatrix reference-point --referencePoint TSS \
  -S chip.bw -R genes.bed -b 3000 -a 3000 --skipZeros \
  -o matrix.gz

# 4) 히트맵 + 메타플롯
plotHeatmap -m matrix.gz -o chip_TSS_heatmap.png --colorMap Blues

ChIP과 input을 함께 보고 싶다면, 2번에서 bamCompare로 log2 비율 트랙을 만들어 -S에 넣으면 됩니다. 트랙을 여러 개 넣으면(-S a.bw b.bw) computeMatrix가 나란히 처리해, plotHeatmap이 조건별 패널을 한 그림에 그려 주죠. bigWig만 있으면 원본 BAM 없이도 3~4단계를 반복할 수 있어, 파라미터를 바꿔 가며 그림을 다듬기 좋습니다.

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

• The file ... does not have an index — deeptools는 BAM에 인덱스(.bai)가 함께 있어야 돕니다. samtools index aln.sorted.bam으로 인덱스를 먼저 만드세요. 인덱싱은 좌표순 정렬이 끝난 뒤라야 되므로, 정렬이 안 됐으면 samtools sort부터 돌립니다.
• --effectiveGenomeSize 관련 오류(RPGC) — RPGC (1x) 정규화는 유효 유전체 크기가 반드시 필요합니다. 값은 유전체·리드 길이마다 다르니 공식 문서 표에서 맞는 값을 찾아 넣으세요(GRCh38·100bp면 2,805,636,231).
• 히트맵이 텅 비거나 신호가 0 — bigWig와 BED의 염색체 이름 규칙이 다른 경우가 흔합니다. 한쪽은 chr1, 다른 쪽은 1이면 겹침이 0으로 나오죠. 이름을 한쪽으로 통일하고, computeMatrix에 --skipZeros가 과하게 걸리지 않았는지도 확인합니다.
• computeMatrix가 너무 느리거나 메모리 부족 — 구간·트랙이 많으면 무거워집니다. -p(또는 --numberOfProcessors)로 코어를 늘리고, bamCoverage 단계의 --binSize를 키워 트랙을 가볍게 만드세요.
• 신호에 잡음·중복이 섞임 — 낮은 품질 매핑이나 PCR 중복이 트랙을 흐립니다. bamCoverage에 --minMappingQuality 30으로 낮은 MAPQ 리드를, --ignoreDuplicates로 중복 리드를 걸러 냅니다.

8. 확장 — 더 해 볼 것

• plotFingerprint로 ChIP 품질 보기 — plotFingerprint는 리드가 소수 구간에 얼마나 몰렸는지를 누적 곡선으로 그려, ChIP의 신호 대 잡음(signal-to-noise)을 한눈에 판단하게 해 줍니다. 히트맵 전 QC 단계로 좋죠.
• computeGCBias·correctGCBias로 GC 편향 보정 — 라이브러리 준비에서 생기는 GC 편향을 진단(computeGCBias)하고 보정(correctGCBias)합니다. 커버리지가 GC 함량에 휘둘릴 때 씁니다.
• multiBigwigSummary로 트랙끼리 비교 — BAM이 아니라 이미 만든 bigWig 여러 개의 상관을 볼 땐 multiBigwigSummary가 짝입니다. 공개 ENCODE bigWig와 내 트랙을 견줄 때 편하죠.
• alignmentSieve로 ATAC-seq 전처리 — alignmentSieve는 리드를 필터링하거나 ATAC-seq용 Tn5 이동 보정(shift)을 넣어 줍니다. ATAC 신호 트랙을 만들기 전 단계로 자주 씁니다.
• pyGenomeTracks로 출판용 트랙 그림 — deeptools 계열인 pyGenomeTracks는 bigWig·BED·유전자 모델을 한 창에 겹쳐 그려, 특정 유전자좌위(locus)의 트랙 그림을 논문 그림 수준으로 만듭니다.

다음 학습

References

1. Ramírez, F., Ryan, D. P., Grüning, B., Bhardwaj, V., Kilpert, F., Richter, A. S., Heyne, S., Dündar, F., & Manke, T. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Research, 44(W1), W160–W165. doi:10.1093/nar/gkw257 (deeptools 원논문)
2. Ramírez, F., Dündar, F., Diehl, S., Grüning, B. A., & Manke, T. (2014). deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Research, 42(W1), W187–W191. doi:10.1093/nar/gku365 (초판)
3. Bhardwaj, V., et al. (2024). deepTools documentation (v3.5) — tools to process and analyze deep sequencing data. deeptools.readthedocs.io
4. deepTools. Effective genome size — mappable genome sizes per read length. deeptools.readthedocs.io/effectiveGenomeSize

Pipette & Pipeline · A bio portfolio journal

CODE BENCH · NGS PIPELINES

bedtools — 유전체 구간 연산: intersect·merge·coverage와 BED 포맷

TL;DR — bedtools는 유전체 위의 좌표 구간(genomic interval)을 다루는 명령줄 도구 모음으로, '유전체 산술(genome arithmetic)의 스위스 아미 나이프'라 불립니다. BED·BAM·VCF·GFF에 담긴 구간을 교집합(intersect)·병합(merge)·커버리지(coverage) 같은 연산으로 한 줄에 처리하죠.

핵심은 BED 포맷과 그 좌표계입니다. BED는 chrom·start·end 세 칸이 기본인데, start가 0-based이고 end는 배제되는 반열림(half-open) 좌표라 chr1 0 100은 첫 100개 염기를 뜻합니다. 이 0-based가 1-based인 VCF·GFF와 자주 헷갈리는 지점이라 이 글에서 특히 강조합니다.

작은 예제 BED 파일로 intersect의 -wa·-wb·-v·-c, merge·sort, coverage·genomecov, subtract·closest·slop을 한 명령씩 따라간 뒤, ChIP-seq 피크와 프로모터의 교집합, 변이(VCF)와 엑손의 교집합을 samtools와 파이프(|)로 잇는 실전까지 만듭니다. 명령은 bedtools 2.31 기준입니다.

🔗 관련 글 · bedtools가 다루는 구간은 대개 다른 도구의 산출물입니다 — 정렬 결과는 SAM/BAM 파일이란?과 서열 정렬이란?, 변이 정보는 VCF 파일이란?, 그 BAM·VCF를 직접 다루는 samtools·bcftools 실전, 그리고 구간의 대표 사례인 ChIP-seq 피크에서 이어집니다.

이 글에서 만들 것

NGS 분석을 하다 보면 결국 좌표를 가진 '구간'을 비교하는 일이 반복됩니다. ChIP-seq 피크가 프로모터 안에 있나, 변이가 엑손에 걸리나, 이 유전자 구간의 커버리지는 얼마인가 같은 질문이죠. 스프레드시트로 좌표를 일일이 맞춰 보긴 어렵고, 직접 파이썬으로 겹침을 계산하자니 경계 조건에서 실수하기 쉽습니다.

bedtools는 이런 '구간 대 구간' 질문을 명령 한 줄로 풉니다. 두 파일의 구간이 겹치는지(intersect), 한 파일 안에서 인접한 구간을 하나로 합칠지(merge), 어떤 구간이 리드로 얼마나 덮이는지(coverage)를 각각 전용 서브명령으로 처리하죠. 입력은 BED뿐 아니라 BAM·VCF·GFF도 되고, 출력은 다시 BED라 파이프(|)로 samtools 같은 도구와 자연스럽게 이어집니다.

이 글은 그 bedtools를 처음부터 익힙니다. 구체적으로 ① BED 포맷과 0-based 좌표를 짚고, ② intersect로 교집합·겹침을 뽑고, ③ merge·sort로 구간을 병합하고, ④ coverage·genomecov로 커버리지를 계산하고, ⑤ subtract·closest·slop으로 빼기·최근접·확장을 다룬 뒤, ⑥ ChIP-seq 피크 ∩ 프로모터와 변이 ∩ 엑손을 실전으로 묶습니다.

# 핵심 흐름 한눈에
bedtools intersect -a peaks.bed -b promoters.bed -u   # ① 겹치는 구간(교집합)
bedtools merge -i sorted.bed                          # ② 인접·겹치는 구간 병합
bedtools coverage -a windows.bed -b reads.bed         # ③ 구간별 커버리지

1. 사전 준비 — 설치와 BED 포맷

bedtools는 리눅스·macOS용 명령줄 도구입니다. conda로 까는 게 가장 간편합니다.

# conda (bioconda) — 권장
conda install -c bioconda bedtools

# 설치 확인(버전 출력)
bedtools --version

bedtools v2.31.0

이제 포맷을 봅시다. BED (Browser Extensible Data)는 구간을 적는 표준 텍스트 포맷으로, 탭으로 나뉜 세 칸이 기본입니다.

• chrom — 염색체 이름(chr1)
• chromStart — 시작 좌표, 0-based
• chromEnd — 끝 좌표, 배제됨(half-open)

여기서 좌표계가 가장 중요합니다. BED는 0-based 반열림(half-open) 좌표라, 시작은 0부터 세고 끝은 포함하지 않습니다. 그래서 chr1 0 100은 1번째부터 100번째까지 100개 염기를 가리키고, 단일 염기 하나(1-based로 5번 위치)는 chr1 4 5로 적습니다. 반면 VCF·GFF·SAM은 1-based라 같은 위치를 5로 씁니다. bedtools를 쓸 때 좌표가 1씩 어긋나 보이면 대개 이 0-based를 잊은 경우죠.

칸을 더 붙이면 정보가 늘어납니다. BED3는 위 세 칸, BED6은 여기에 name·score·strand를 더한 여섯 칸, BED12는 exon 블록 구조(blockCount·blockSizes·blockStarts)까지 담아 유전자 하나를 통째로 표현합니다.

예제로 쓸 작은 파일 두 개를 만듭니다. peaks.bed는 ChIP-seq 피크, genes.bed는 유전자 구간이라 치죠.

# peaks.bed — 피크 3개
printf "chr1\t100\t200\tpeak1\nchr1\t400\t500\tpeak2\nchr1\t700\t900\tpeak3\n" > peaks.bed

# genes.bed — 유전자 2개
printf "chr1\t150\t350\tgeneA\nchr1\t800\t1000\tgeneB\n" > genes.bed

2. intersect — 구간 교집합과 겹침

가장 많이 쓰는 서브명령이 intersect입니다. -a와 -b 두 파일을 받아, A의 각 구간이 B의 구간과 겹치는 부분을 찾습니다. 옵션 없이 돌리면 겹치는 구간(교집합) 자체를 출력합니다.

# A(peaks)와 B(genes)가 겹치는 부분만 출력
bedtools intersect -a peaks.bed -b genes.bed

chr1	150	200	peak1
chr1	800	900	peak3

peak1(100–200)은 geneA(150–350)와 150–200에서 겹치고, peak3(700–900)은 geneB(800–1000)와 800–900에서 겹칩니다. peak2(400–500)는 어느 유전자와도 안 겹쳐 빠졌죠. 출력의 좌표는 겹친 구간으로 잘려 나오고, peak 이름 같은 A의 나머지 칸은 그대로 따라옵니다.

원본 구간을 통째로 보고 싶으면 옵션을 붙입니다. 자주 쓰는 조합은 이렇습니다.

• -wa — 겹친 A 구간을 원본 그대로 출력(잘리지 않음)
• -wb — 짝이 된 B 구간도 함께 출력
• -v — B와 안 겹치는 A만(여집합)
• -c — A마다 겹친 B의 개수를 뒤에 붙임
• -u — 겹침이 하나라도 있으면 A를 한 번만 출력(중복 제거)

# -wa -wb: 원본 A와 짝이 된 B를 나란히
bedtools intersect -a peaks.bed -b genes.bed -wa -wb

chr1	100	200	peak1	chr1	150	350	geneA
chr1	700	900	peak3	chr1	800	1000	geneB

# -v: 어느 유전자와도 안 겹치는 피크(여집합)
bedtools intersect -a peaks.bed -b genes.bed -v

chr1	400	500	peak2

# -c: 피크마다 겹친 유전자 수를 마지막 칸에
bedtools intersect -a peaks.bed -b genes.bed -c

chr1	100	200	peak1	1
chr1	400	500	peak2	0
chr1	700	900	peak3	1

-a에는 파일 여럿을, 겹침 기준에는 -f(A의 최소 겹침 비율)·-r(양쪽 상호 비율)을 더할 수 있습니다. 예컨대 -f 0.5는 A가 절반 이상 겹칠 때만 셈에 넣죠.

3. merge·sort — 인접·겹치는 구간 병합

merge는 한 파일 안에서 겹치거나 맞닿은 구간을 하나로 합칩니다. 단, merge는 입력이 정렬돼 있어야 하므로 먼저 sort로 좌표순(염색체→시작)으로 맞춥니다.

# regions.bed — 겹치거나 인접한 구간들
printf "chr1\t100\t200\nchr1\t180\t300\nchr1\t400\t500\nchr1\t480\t600\nchr1\t900\t1000\n" > regions.bed

# 좌표순 정렬(염색체별, 시작 위치 오름차순)
bedtools sort -i regions.bed > regions.sorted.bed

# 겹치거나 맞닿은 구간을 하나로 병합
bedtools merge -i regions.sorted.bed

chr1	100	300
chr1	400	600
chr1	900	1000

100–200과 180–300은 겹쳐 100–300으로, 400–500과 480–600은 겹쳐 400–600으로 합쳐지고, 멀리 떨어진 900–1000은 그대로 남습니다. 조금 떨어진 구간까지 합치려면 -d로 허용 간격을 줍니다. -d 350이면 600과 900 사이(300 간격)도 이어 붙죠.

병합하면서 값을 요약할 수도 있습니다. -c로 대상 칸을, -o로 연산(count·sum·mean·collapse 등)을 지정합니다.

# 병합된 각 구간이 원래 몇 개 구간에서 왔는지 세기
bedtools merge -i regions.sorted.bed -c 1 -o count

chr1	100	300	2
chr1	400	600	2
chr1	900	1000	1

4. coverage·genomecov — 커버리지 계산

coverage는 A의 각 구간이 B의 구간(리드 등)으로 얼마나 덮이는지를 셉니다. 기본 출력은 각 A 구간 뒤에 겹친 B의 개수, 덮인 염기 수, 구간 길이, 덮인 비율을 덧붙입니다.

# windows.bed의 각 구간을 reads.bed가 얼마나 덮는지
bedtools coverage -a windows.bed -b reads.bed

chr1	0	500	12	480	500	0.9600000
#                ↑개수 ↑덮인bp ↑길이 ↑비율

마지막 열의 비율(0.96)은 이 창(window)의 96%가 리드로 덮였다는 뜻입니다. 창 단위 커버리지를 훑을 때 유용하죠. 구간별이 아니라 염기별(per-base) 커버리지가 필요하면 genomecov를 씁니다. BAM을 바로 받아 BedGraph로 뽑을 수 있습니다.

# BAM에서 염기별 커버리지를 BedGraph로(-bg: 0인 구간은 생략)
bedtools genomecov -ibam aln.sorted.bam -bg > coverage.bedgraph

chr1	0	120	3
chr1	120	340	7
chr1	340	500	2

각 줄은 '이 구간의 모든 염기가 N개 리드로 덮였다'는 뜻입니다. genomecov는 BAM이 좌표순으로 정렬돼 있어야 하고, 리드가 아닌 BED 입력이면 -i와 함께 염색체 크기 파일(-g)을 줍니다.

5. subtract·closest·slop — 빼기·최근접·확장

나머지 자주 쓰는 서브명령 셋을 묶어 봅니다.

subtract는 A에서 B와 겹치는 부분을 잘라 냅니다. 블랙리스트 구간이나 반복서열을 제거할 때 씁니다.

# peaks에서 genes와 겹치는 부분을 제거
bedtools subtract -a peaks.bed -b genes.bed

chr1	100	150	peak1
chr1	400	500	peak2
chr1	700	800	peak3

peak1은 150부터 geneA와 겹쳐 100–150만 남고, peak3은 800부터 겹쳐 700–800만 남습니다. peak2는 겹침이 없어 통째로 유지되죠.

closest는 A의 각 구간에 가장 가까운 B 구간을 찾습니다. 겹치는 게 없어도 제일 가까운 것을 돌려주죠. sort된 입력이 필요하고, -d를 붙이면 거리(bp)를 마지막 칸에 적습니다.

# 각 피크에서 가장 가까운 유전자와 그 거리
bedtools closest -a peaks.sorted.bed -b genes.sorted.bed -d

chr1	100	200	peak1	chr1	150	350	geneA	0
chr1	400	500	peak2	chr1	150	350	geneA	51
chr1	700	900	peak3	chr1	800	1000	geneB	0

거리 0은 겹친다는 뜻이고, peak2는 geneA에서 51bp 떨어져 있습니다(0-based라 400과 349의 차이).

slop은 각 구간을 양옆으로 넓힙니다. TSS 주변 프로모터를 만들 때처럼 자주 쓰죠. 염색체 밖으로 삐져나가지 않게 크기 파일(-g)이 필요합니다.

# genome.txt: 염색체 크기(chrom <TAB> size)
printf "chr1\t248956422\n" > genome.txt

# 각 구간을 양쪽으로 50bp씩 확장(-b), 한쪽만 넓히려면 -l/-r
bedtools slop -i genes.bed -g genome.txt -b 50

chr1	100	400	geneA
chr1	750	1050	geneB

한쪽으로만 크기 0인 flanking 구간을 만들려면 flank를 씁니다. slop이 원본을 포함해 넓힌다면, flank는 원본 양옆의 새 구간만 돌려주죠.

6. 실전 — ChIP-seq 피크 ∩ 프로모터, 변이 ∩ 엑손

이제 조각을 모아 실전 질문에 답합니다. 첫째, 어떤 ChIP-seq 피크가 프로모터에 있나? 전사 시작점(TSS) 좌표에서 slop으로 프로모터(TSS ±2kb)를 만든 뒤, intersect의 -u로 프로모터에 걸친 피크만 추립니다.

# 1) TSS ±2kb로 프로모터 구간 만들기
bedtools slop -i tss.bed -g genome.txt -b 2000 > promoters.bed

# 2) 프로모터에 겹치는 피크만(한 번씩)
bedtools intersect -a peaks.bed -b promoters.bed -u > peaks_in_promoters.bed

둘째, 어떤 변이가 엑손에 있나? intersect는 VCF도 그대로 받습니다. -header를 주면 VCF 헤더를 살려 결과가 다시 유효한 VCF가 되죠.

# 엑손에 걸친 변이만(VCF 헤더 유지)
bedtools intersect -a variants.vcf -b exons.bed -header > exonic_variants.vcf

마지막으로, bedtools의 출력은 다시 BED라 파이프로 이어 붙일 수 있습니다. 예컨대 특정 구간에 정렬된 리드만 BAM으로 뽑아 바로 커버리지를 볼 수 있죠.

# 표적 구간(-L)의 리드만 뽑아 그대로 커버리지 계산
samtools view -b -L targets.bed aln.bam \
  | bedtools coverage -a targets.bed -b - > target_coverage.txt

-b -의 -는 '앞 명령의 출력을 표준입력으로 받으라'는 뜻입니다. 중간 파일 없이 samtools와 bedtools가 곧장 연결되죠.

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

• ***** ERROR: Requested column ... , but the number of columns ... — 입력이 탭이 아니라 공백으로 나뉘었을 때 흔합니다. BED는 탭 구분이 원칙이라, 스프레드시트에서 붙여넣은 공백을 sed 's/ /\t/g'나 printf "...\t..."로 탭으로 바꿔야 합니다.
• Error: Sorted input specified, but the file ... has ... out of order record — merge·closest처럼 정렬을 요구하는 명령에 정렬 안 된 파일을 넣은 경우입니다. 먼저 bedtools sort -i in.bed(또는 sort -k1,1 -k2,2n)로 정렬하세요.
• 좌표가 1씩 어긋나 보임 — 0-based를 잊은 경우입니다. BED start는 0-based, VCF·GFF는 1-based라, 같은 위치라도 BED가 1 작습니다. 다른 포맷과 맞대볼 때는 좌표계를 먼저 확인하세요.
• Error: Unable to open file ... genome file — slop·genomecov에 크기 파일(-g)을 안 줬거나 경로가 틀린 경우입니다. samtools faidx ref.fa로 만든 ref.fa.fai의 앞 두 칸을 쓰거나, chrom<TAB>size 형식으로 직접 만듭니다.
• 염색체 이름이 안 맞음 — 한쪽은 chr1, 다른 쪽은 1이면 겹침이 0으로 나옵니다. 이름 규칙(UCSC식 chr1 vs Ensembl식 1)을 한쪽으로 통일해야 합니다.

8. 확장 — 더 해 볼 것

• complement로 빈 구간 찾기 — bedtools complement는 크기 파일 기준으로 구간이 없는 영역을 돌려줍니다. 캡처 패널 밖 영역이나 갭을 찾을 때 유용하죠.
• makewindows로 유전체 창 나누기 — bedtools makewindows -g genome.txt -w 10000은 유전체를 10kb 창으로 잘라, coverage와 짝지어 창 단위 커버리지·밀도 프로파일을 만듭니다.
• getfasta로 서열 뽑기 — bedtools getfasta -fi ref.fa -bed peaks.bed는 구간에 해당하는 실제 염기서열을 FASTA로 추출합니다. 모티프 분석의 입력으로 바로 이어지죠.
• pybedtools로 파이썬에서 — 같은 연산을 파이썬 스크립트 안에서 하려면 pybedtools가 bedtools를 감싸 줍니다. pandas와 엮어 대량 구간을 프로그램적으로 다루기 좋습니다.
• 대안 도구와 비교 — 아주 큰 데이터에서 속도가 급하면 bedops가 정렬 기반으로 더 빠를 때가 있고, R에서는 GenomicRanges가 같은 구간 연산을 객체로 제공합니다. 셸에서 빠르게 훑기엔 bedtools가 여전히 표준이죠.

다음 학습

References

1. Quinlan, A. R., & Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 26(6), 841–842. doi:10.1093/bioinformatics/btq033 (bedtools 원논문)
2. Quinlan, A. R. (2014). BEDTools: The Swiss-Army Tool for Genome Feature Analysis. Current Protocols in Bioinformatics, 47, 11.12.1–11.12.34. doi:10.1002/0471250953.bi1112s47
3. Quinlan lab. (2023). bedtools documentation (v2.31.0) — a powerful toolset for genome arithmetic. bedtools.readthedocs.io
4. UCSC Genome Browser. Frequently Asked Questions: Data File Formats (BED · 0-based half-open 좌표). genome.ucsc.edu/FAQ/FAQformat.html

Pipette & Pipeline · A bio portfolio journal

TL;DR — Snakemake는 파이썬 기반 워크플로 관리자로, FastQC·fastp·정렬·samtools처럼 손으로 순서대로 돌리던 NGS 도구들을 하나의 Snakefile에 규칙(rule)으로 적어 자동·재현 가능하게 엮습니다. 규칙은 입력(input)→출력(output)→명령(shell)의 세 칸으로 정의하고, 파일 이름의 {sample} 같은 와일드카드(wildcard)로 여러 샘플에 같은 규칙을 일반화합니다.

Snakemake는 규칙들의 입출력을 훑어 의존성 그래프(DAG)를 스스로 그리고, 이 그래프대로 실행합니다. 그래서 snakemake -n으로 무엇이 돌지 미리 보고(dry-run), snakemake --dag | dot -Tsvg로 그래프를 그림으로 뽑을 수 있죠. 파일 타임스탬프를 비교해 바뀐 것만 다시 돌리는 증분 실행, --cores N으로 독립 작업을 동시에 돌리는 병렬화, conda·컨테이너로 도구 버전까지 고정하는 재현성이 핵심 강점입니다.

이 글은 애기장대(Arabidopsis thaliana) 공개 데이터를 예로, conda 설치부터 첫 규칙 작성·와일드카드·expand()·config.yaml·DAG 시각화·병렬 실행까지 한 줄씩 따라간 뒤, FASTQ → FastQC → fastp → 정렬 → samtools 정리를 하나의 Snakefile로 묶습니다. 명령은 Snakemake 8/9 기준입니다(현재 안정 버전 9.x).

🔗 관련 글  ·  이 글이 자동화하는 개별 도구들은 이미 따로 다뤘습니다 — 리드 품질검사·트리밍은 FastQC·fastp 편, 참조 정렬의 원리는 서열 정렬이란?, 정렬 결과 정리는 samtools·bcftools 편, 변이 호출은 GATK 편에서 이어집니다  ·  정렬을 우회하는 유사정렬 정량은 salmon으로 전사체 정량에서 다뤘습니다

이 글에서 만들 것

NGS 분석을 처음 손으로 돌려 보면 대체로 이런 식입니다. FastQC로 품질을 보고, fastp로 어댑터를 자르고, BWA나 STAR로 정렬하고, samtools로 BAM을 정렬·인덱싱하고… 샘플이 하나면 견딜 만한데, 12개·96개로 늘면 명령을 복사·붙여넣기 하다 어느 샘플에서 한 단계를 빠뜨렸는지도 헷갈립니다. 중간에 참조를 바꾸면 처음부터 다 다시 돌려야 하나 고민되고, 몇 달 뒤 "그때 그 결과 어떻게 냈더라"를 재현하기도 어렵죠.

Snakemake는 이 문제를 정면으로 풉니다. 각 단계를 규칙(rule)으로 한 번만 적어 두면, 어떤 샘플이든 같은 규칙이 알아서 적용되고(와일드카드), 도구 사이의 의존성은 Snakemake가 파일 이름을 보고 스스로 잇습니다. 바뀐 입력만 골라 다시 돌리고, 독립적인 샘플들은 병렬로 처리하며, Snakefile 하나가 곧 분석의 명세이자 재현 레시피가 됩니다.

이 글은 그 Snakemake를 처음부터 만듭니다. 구체적으로 ① conda로 Snakemake를 깔고, ② input→output→shell로 첫 규칙을 쓰고, ③ {sample} 와일드카드와 expand()로 여러 샘플에 일반화하고, ④ config.yaml로 경로·샘플 목록을 밖으로 빼고, ⑤ snakemake -n(dry-run)·--dag로 실행 계획을 확인하고, ⑥ --cores로 병렬 실행한 뒤, ⑦ FASTQ → FastQC → fastp → 정렬 → samtools를 하나로 묶은 실전 Snakefile을 완성합니다.

손으로 도구를 이어 붙이는 방식과 뭐가 다른지 한 줄로 보면 이렇습니다. 셸 스크립트는 "명령을 이 순서로 실행하라"는 명령형이지만, Snakemake는 "이 출력 파일을 만들려면 이 입력과 이 명령이 필요하다"는 선언형입니다. 무엇을 만들지만 적으면 순서·의존성·재실행 판단은 Snakemake가 알아서 하죠.

# 핵심 흐름 한눈에
snakemake -n                              # ① dry-run — 무엇이 실행될지 미리 확인
snakemake --dag | dot -Tsvg > dag.svg     # ② 의존성 그래프(DAG)를 그림으로
snakemake --cores 8                       # ③ 실제 실행 (8코어 병렬, 바뀐 것만)

1. 사전 준비 — 설치와 예제 데이터

Snakemake는 파이썬 패키지라 conda (권장) 또는 pip로 설치합니다. bioconda 채널의 conda를 쓰면 뒤에서 규칙별 conda 환경을 자동으로 만들어 주는 기능까지 매끄럽게 이어져, 여기서는 conda로 통일합니다.

# conda (bioconda) — 권장, 전용 환경에 격리 설치
conda create -n smk -c conda-forge -c bioconda snakemake
conda activate smk

# 설치 확인 (버전 출력)
snakemake --version

9.5.1

실습용 도구와 데이터도 같은 환경에 준비합니다. 이 글의 실전 예제(FastQC·fastp·BWA·samtools)를 그대로 돌리려면 도구들이 있어야 하지만, 개념을 익히는 앞 절들은 도구 없이도 따라올 수 있게 구성했습니다.

# 실전 파이프라인용 도구 (개념 학습만 할 거면 건너뛰어도 됨)
conda install -c bioconda fastqc fastp bwa samtools

작업 폴더 구조는 이렇게 잡습니다. 원시 리드는 data/에, 결과는 results/에 쌓고, 파이프라인 명세인 Snakefile과 설정 config.yaml을 폴더 맨 위에 둡니다.

project/
├── Snakefile          # 파이프라인 명세 (규칙 모음)
├── config.yaml        # 경로·샘플 목록 등 설정
├── data/
│   ├── sampleA.fastq.gz
│   └── sampleB.fastq.gz
└── results/           # 실행하면 여기에 결과가 쌓인다

2. 첫 규칙 — input → output → shell

Snakemake의 기본 단위는 규칙(rule)입니다. 규칙 하나가 "어떤 입력으로, 어떤 출력을, 어떤 명령으로 만드는가"를 적습니다. 가장 단순한 예로, FASTQ에 FastQC를 돌려 품질 리포트를 만드는 규칙을 봅시다. Snakefile이라는 이름의 파일에 다음을 적습니다.

# Snakefile — 규칙 하나: FASTQ → FastQC 리포트
rule fastqc:
    input:
        "data/sampleA.fastq.gz"          # 입력 파일
    output:
        "results/qc/sampleA_fastqc.html" # 만들 출력 파일
    shell:
        "fastqc {input} -o results/qc/"  # 실제 실행할 명령

세 칸의 뜻이 곧 Snakemake의 문법 전부에 가깝습니다. input:은 이 규칙이 필요로 하는 파일, output:은 이 규칙이 만들어 낼 파일, shell:은 실제로 돌릴 셸 명령입니다. 명령 안의 {input}·{output}은 위에 적은 경로로 자동 치환됩니다 — 경로를 두 번 적지 않아도 되죠. 규칙 이름(fastqc)은 사람이 알아보기 위한 라벨입니다.

이제 실행해 봅니다. Snakemake는 인자로 만들고 싶은 출력 파일을 받고, 그 파일을 만들 수 있는 규칙을 거꾸로 찾아 실행합니다.

# 원하는 출력 파일을 목표(target)로 지정해 실행
snakemake --cores 1 results/qc/sampleA_fastqc.html

Building DAG of jobs...
Job stats:
job       count
------    -----
fastqc    1
total     1

rule fastqc:
    input: data/sampleA.fastq.gz
    output: results/qc/sampleA_fastqc.html
Finished job 0.
1 of 1 steps (100%) done

여기서 두 가지가 중요합니다. 첫째, --cores(또는 -c)는 Snakemake 8부터 반드시 지정해야 하는 옵션입니다. 몇 개의 코어를 써도 되는지 알려 주지 않으면 실행을 거부하죠(예전 5~7.x에서는 생략 시 1코어였지만, 8부터는 명시가 필수입니다). 둘째, 목표 파일이 이미 최신 상태면 Snakemake는 "Nothing to be done"이라며 다시 돌리지 않습니다 — 이 판단이 뒤에서 볼 증분 실행의 핵심입니다.

shell 대신 파이썬 코드로 처리하고 싶으면 run:을, 파이썬 스크립트나 R 스크립트를 부르려면 script:를 쓸 수 있습니다. 명령이 길거나 옵션이 많으면 params:로 값을 빼두고, 스레드 수는 threads:, 로그는 log:로 관리합니다.

# 옵션을 갖춘 규칙: fastp 트리밍 (params·threads·log 활용)
rule fastp:
    input:
        "data/{sample}.fastq.gz"
    output:
        trimmed="results/trimmed/{sample}.fastq.gz",
        report="results/trimmed/{sample}.fastp.html"
    params:
        quality=20                       # 품질 컷오프 (params로 분리)
    threads: 4                           # 이 규칙이 쓸 스레드 수
    log:
        "logs/fastp/{sample}.log"        # 로그 파일
    shell:
        "fastp -i {input} -o {output.trimmed} "
        "-h {output.report} -q {params.quality} "
        "--thread {threads} 2> {log}"

출력이 여럿이면 이름을 붙여(trimmed=·report=) {output.trimmed}처럼 골라 쓸 수 있습니다. {params.quality}·{threads}·{log}도 같은 방식으로 명령 안에서 참조됩니다. 문자열을 여러 줄로 나눠 적으면 파이썬이 이어 붙이니, 긴 명령을 읽기 좋게 쪼갤 수 있죠.

3. 와일드카드와 expand() — 여러 샘플로 일반화

앞 규칙에 이미 답이 나와 있습니다. data/{sample}.fastq.gz의 {sample} — 이게 와일드카드(wildcard)입니다. 특정 샘플 이름을 박아 두는 대신 {sample}이라는 자리표시자를 쓰면, 하나의 규칙이 모든 샘플에 적용됩니다.

작동 방식이 꽤 영리한데, Snakemake는 출력 파일 이름에서 와일드카드 값을 거꾸로 알아냅니다. 예컨대 results/trimmed/sampleA.fastq.gz를 만들라는 요청이 오면, 출력 패턴 results/trimmed/{sample}.fastq.gz와 맞춰 {sample}=sampleA로 정하고, 그 값을 입력 패턴 data/{sample}.fastq.gz에도 그대로 채워 data/sampleA.fastq.gz를 입력으로 잡습니다. 즉 출력에서 입력으로 와일드카드가 전파되죠.

# 와일드카드 규칙: {sample} 하나로 모든 샘플 처리
rule fastqc:
    input:
        "data/{sample}.fastq.gz"
    output:
        "results/qc/{sample}_fastqc.html"
    shell:
        "fastqc {input} -o results/qc/"

이제 snakemake --cores 1 results/qc/sampleB_fastqc.html을 하면 같은 규칙이 sampleB에 적용됩니다. 그런데 샘플이 열 개면 목표 파일을 열 번 타이핑할까요? 그래서 expand()가 있습니다. 샘플 목록을 주면 파일 이름 리스트를 한 번에 펼쳐 줍니다.

# Snakefile 맨 위 — 샘플 목록과 최종 목표 정의
SAMPLES = ["sampleA", "sampleB", "sampleC"]

# 규칙 all: 최종적으로 원하는 파일들을 모아 두는 '목표 규칙'
rule all:
    input:
        expand("results/qc/{sample}_fastqc.html", sample=SAMPLES)

두 가지 관례가 여기 담겨 있습니다. 첫째, expand("....{sample}....", sample=SAMPLES)는 리스트의 각 값을 끼워 넣어 ["results/qc/sampleA_fastqc.html", ...]처럼 펼칩니다. 둘째, Snakefile의 맨 첫 규칙이 기본 목표(default target)가 됩니다. 관례상 아무것도 만들지 않고 input에 "최종 산출물"만 모아 두는 rule all을 맨 위에 두면, 그냥 snakemake --cores 8만 쳐도 all이 요구하는 모든 파일이 만들어집니다.

# 목표 파일 없이 실행 → 첫 규칙(all)이 기본 목표가 되어 전체 샘플 처리
snakemake --cores 8

4. config.yaml — 설정을 코드 밖으로

샘플 목록·참조 경로·파라미터를 Snakefile에 박아 두면, 프로젝트가 바뀔 때마다 코드를 고쳐야 합니다. 이런 값은 config.yaml로 빼는 게 관례입니다. YAML 파일에 설정을 적고, Snakefile에서 configfile:로 불러오면 config 딕셔너리로 접근할 수 있죠.

# config.yaml — 경로·샘플·파라미터 설정
samples:
  - sampleA
  - sampleB
  - sampleC
genome: "ref/TAIR10.fa"
fastp_quality: 20

# Snakefile — config.yaml 불러와 값 참조
configfile: "config.yaml"

rule all:
    input:
        expand("results/qc/{sample}_fastqc.html",
               sample=config["samples"])       # config에서 샘플 목록 가져오기

이제 샘플이 늘거나 참조가 바뀌어도 config.yaml만 고치면 됩니다. config["samples"]·config["genome"]·config["fastp_quality"]처럼 키로 값을 꺼내 규칙 곳곳에서 재사용하죠. 명령줄에서 --configfile other.yaml로 다른 설정을 바꿔 끼우거나, --config genome=ref/hg38.fa로 개별 값을 덮어쓸 수도 있습니다.

5. DAG·dry-run — 실행 계획을 미리 본다

규칙을 여럿 적으면 Snakemake는 어느 규칙의 출력이 어느 규칙의 입력인지를 파일 이름으로 이어, 방향성 비순환 그래프(DAG, directed acyclic graph)를 스스로 만듭니다. 이 그래프가 곧 실행 계획입니다 — 노드는 작업(job), 화살표는 의존 관계죠.

실제로 돌리기 전에 계획부터 확인하는 게 안전합니다. snakemake -n(또는 --dry-run)은 아무것도 실행하지 않고 무슨 작업이 어떤 순서로 돌지만 보여 줍니다. 오타·경로 실수·빠진 입력을 실행 전에 잡는 필수 습관이죠.

# dry-run — 실행하지 않고 계획만 출력
snakemake -n --cores 8

Building DAG of jobs...
Job stats:
job       count
------    -----
all       1
fastqc    3
fastp     3
total     7

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

그래프를 눈으로 보고 싶으면 --dag로 DOT 형식을 뽑아 Graphviz의 dot으로 이미지를 만듭니다. 작업이 많아 복잡하면 규칙 수준으로만 요약하는 --rulegraph가 더 읽기 좋습니다.

# 작업 수준 DAG를 SVG 이미지로 (Graphviz 필요: conda install -c conda-forge graphviz)
snakemake --dag | dot -Tsvg > dag.svg

# 규칙 수준으로 간단히 (샘플이 많을 때 권장)
snakemake --rulegraph | dot -Tsvg > rulegraph.svg

6. 증분 실행·병렬화 — 바뀐 것만, 동시에

Snakemake가 손으로 돌리기보다 나은 결정적 이유 두 가지가 증분 실행과 병렬화입니다.

증분 실행은 기본적으로 Make식 타임스탬프 비교에서 출발합니다. 출력 파일이 이미 있고 그 입력보다 새것(타임스탬프가 더 최근)이면, 그 작업은 이미 끝난 것으로 보고 건너뜁니다. 반대로 입력을 고치거나(더 최근이 됨) 출력을 지우면, 그 작업과 그 뒤에 딸린 작업들만 다시 돌죠. 그래서 100개 샘플 중 3개 리드만 교체하면, 그 3개 계보만 다시 흐릅니다 — 나머지 97개는 그대로 두고요.

다만 요즘 Snakemake는 타임스탬프 하나만 보지 않는다는 점이 Make와 다릅니다. 기본값으로 다섯 가지 방아쇠(trigger) — 파일 수정 시각(mtime)·입력 집합(input)·params 값·규칙 코드(code)·소프트웨어 환경(software-env) — 중 무엇이라도 바뀌면 재실행합니다. 그래서 입력 타임스탬프가 그대로여도 규칙의 shell 명령이나 params를 고치면 다시 돕니다. Make만 써 본 사람이 흔히 놀라는 지점인데, 결과가 항상 코드·설정과 일치하도록 보장하려는 설계죠. 순수하게 옛 Make처럼 타임스탬프만 보게 하려면 --rerun-triggers mtime을 줍니다.

# sampleB의 입력만 새로 바뀐 상황 → 그 계보만 재실행
snakemake --cores 8

Building DAG of jobs...
Job stats:
job       count
------    -----
fastp     1
total     1

# sampleB만 다시 돌고, sampleA·sampleC는 최신이라 건너뜀

병렬화는 --cores N(또는 -c N, -j N) 한 방입니다. Snakemake는 DAG에서 서로 의존하지 않는 작업을 골라 동시에 돌립니다. 코어 8개를 주면 독립적인 샘플들의 FastQC를 한꺼번에 처리하죠. 규칙에 threads:를 준 경우엔 그 값과 --cores를 함께 저울질해 배분합니다.

# 8코어로 병렬 실행 — 독립 작업을 동시에
snakemake --cores 8

강제로 전부 다시 돌리고 싶으면 --forceall(또는 -F), 특정 규칙만 다시 돌리려면 --forcerun 규칙이름을 씁니다.

7. 실전 — FASTQ→QC→트리밍→정렬→정리 한 파일로

이제 앞 조각들을 모아 실전 파이프라인을 만듭니다. FASTQ에서 시작해 FastQC·fastp 품질검사·트리밍, BWA-MEM 정렬, samtools 정렬·인덱싱까지 다섯 단계를 하나의 Snakefile로 묶습니다. 각 단계는 규칙 하나이고, 앞 규칙의 출력이 뒤 규칙의 입력으로 자연히 이어집니다. 각 도구의 옵션·개념은 앞선 글들에서 다뤘으니, 여기서는 그것들을 어떻게 자동으로 엮는지에 집중합니다.

# Snakefile — FASTQ → QC → 트리밍 → 정렬 → BAM 정렬·인덱싱
configfile: "config.yaml"

# 최종 목표: 모든 샘플의 정렬·인덱싱된 BAM과 QC 리포트
rule all:
    input:
        expand("results/sorted/{sample}.bam.bai", sample=config["samples"]),
        expand("results/qc/{sample}_fastqc.html", sample=config["samples"])

# 1) FastQC — 원시 리드 품질검사
rule fastqc:
    input:
        "data/{sample}.fastq.gz"
    output:
        "results/qc/{sample}_fastqc.html"
    shell:
        "fastqc {input} -o results/qc/"

# 2) fastp — 어댑터·저품질 염기 트리밍
rule fastp:
    input:
        "data/{sample}.fastq.gz"
    output:
        "results/trimmed/{sample}.fastq.gz"
    threads: 4
    log:
        "logs/fastp/{sample}.log"
    shell:
        "fastp -i {input} -o {output} --thread {threads} 2> {log}"

# 3) BWA-MEM — 트리밍된 리드를 참조에 정렬 → BAM
#    (미리 bwa index {genome} 로 참조 색인이 만들어져 있어야 함)
rule bwa_map:
    input:
        reads="results/trimmed/{sample}.fastq.gz",
        ref=config["genome"]
    output:
        "results/mapped/{sample}.bam"
    threads: 8
    shell:
        "bwa mem -t {threads} {input.ref} {input.reads} "
        "| samtools view -b - > {output}"

# 4) samtools sort — 좌표순 정렬
rule samtools_sort:
    input:
        "results/mapped/{sample}.bam"
    output:
        "results/sorted/{sample}.bam"
    threads: 4
    shell:
        # -T: 임시파일 접두어 (파일 목록이 아니라 접두어라 {wildcards.sample} 사용)
        "samtools sort -@ {threads} -T results/sorted/{wildcards.sample} "
        "{input} -o {output}"

# 5) samtools index — BAM 인덱싱(.bai)
rule samtools_index:
    input:
        "results/sorted/{sample}.bam"
    output:
        "results/sorted/{sample}.bam.bai"
    shell:
        "samtools index {input}"

이 한 파일이 파이프라인 전부입니다. snakemake -n --cores 8로 계획을 확인하고, 문제없으면 snakemake --cores 8로 돌리면 됩니다. Snakemake는 rule all이 요구하는 .bam.bai를 만들려면 정렬된 BAM이, 그러려면 매핑된 BAM이, 그러려면 트리밍된 리드가 필요하다는 걸 파일 이름만으로 역추적해 fastp→bwa_map→samtools_sort→samtools_index 순서를 스스로 세웁니다. 순서를 어디에도 명시하지 않았는데도요.

재현성을 한 단계 더 올리려면 규칙마다 conda:로 도구 환경을 고정합니다. 규칙에 환경 YAML을 물리고 --software-deployment-method conda(짧게 --sdm conda)로 실행하면, Snakemake가 규칙별로 그 환경을 만들어 그 안에서 명령을 돌립니다. 도구 버전까지 파이프라인에 박히니, 몇 달 뒤 다른 컴퓨터에서도 같은 결과가 재현되죠.

# 규칙에 conda 환경 지정 (envs/fastp.yaml에 fastp 버전 명시)
rule fastp:
    input:  "data/{sample}.fastq.gz"
    output: "results/trimmed/{sample}.fastq.gz"
    conda:  "envs/fastp.yaml"
    shell:  "fastp -i {input} -o {output}"

# 규칙별 conda 환경을 자동 생성·사용해 실행 (재현성↑)
snakemake --cores 8 --software-deployment-method conda

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

• --cores를 안 줘서 실행 거부 — Snakemake 8부터는 코어 수를 반드시 지정해야 합니다. 그냥 snakemake만 치면 실행을 멈추고 코어를 명시하라고 요구하죠. snakemake --cores 8(또는 전부 쓰려면 --cores all)처럼 항상 코어 수를 붙이세요.
• MissingInputException — 입력 파일이 없음 — Snakemake가 목표를 만들려면 어떤 입력이 필요한데 그 파일도, 그걸 만들 규칙도 없을 때 납니다. 대개 경로 오타이거나(data/ vs data), 파일 이름 규칙이 실제 파일과 안 맞는 경우입니다. 메시지에 어떤 파일을 찾다 실패했는지 나오니, 그 경로가 실제로 존재하는지·와일드카드 패턴과 맞는지 확인합니다(경로의 ~는 확장되지 않으니 작업 폴더 기준 상대경로를 쓰세요).
• AmbiguousRuleException — 규칙이 겹침 — 같은 출력 파일을 만들 수 있는 규칙이 둘 이상이면, Snakemake가 어느 것을 쓸지 몰라 멈춥니다. 출력 패턴이 서로 겹치지 않게 하거나, ruleorder: 규칙A > 규칙B로 우선순위를 정해 해결합니다.
• 와일드카드가 예상과 다르게 잡힘 — {sample}은 기본적으로 /를 포함해 욕심껏(greedy) 매칭하므로, 경로가 복잡하면 엉뚱하게 잡힐 수 있습니다. wildcard_constraints로 정규식을 걸어(예 sample="[A-Za-z0-9]+") 매칭 범위를 좁히면 안전합니다.
• 중간 파일이 자꾸 지워지거나 남음 — 용량 큰 중간 산출물은 temp("...")로 감싸면, 그걸 쓰는 뒷 작업이 끝난 뒤 Snakemake가 자동으로 지웁니다. 반대로 실수로 덮어쓰면 안 되는 최종 결과는 protected("...")로 감싸 쓰기 방지를 걸 수 있죠.

9. 확장 — 더 해 볼 것

• 모듈화(include·rules/) — 규칙이 많아지면 Snakefile 하나가 길어집니다. 규칙을 주제별 .smk 파일로 나눠 include: "rules/mapping.smk"로 불러오면 관리가 쉬워집니다.
• --use-conda에서 --software-deployment-method로 (버전 주의) — 재현 환경을 붙이는 플래그가 Snakemake 8에서 이름이 바뀌었습니다. 예전 --use-conda·--use-singularity는 8부터 --software-deployment-method conda·apptainer(짧게 --sdm)로 통합됐죠(옛 --use-conda도 당분간 별칭으로 동작하지만 권장되지 않습니다). 인터넷의 오래된 튜토리얼은 --use-conda를 쓰니, 버전에 맞는 플래그를 확인하세요.
• 컨테이너로 완전 재현(container:) — conda보다 더 강하게 고정하려면 규칙에 container: "docker://..."로 도커·Apptainer 이미지를 물려, OS 라이브러리까지 통째로 담습니다. --sdm apptainer로 실행하면 그 이미지 안에서 명령이 돕니다.
• 클러스터·클라우드 실행 — Snakemake 8부터 실행 백엔드를 플러그인으로 분리했습니다. 예전 --cluster "sbatch ..."·--kubernetes 대신, snakemake-executor-plugin-slurm 같은 플러그인을 깔고 --executor slurm으로 같은 Snakefile을 SLURM 클러스터나 클라우드에 그대로 던져 대규모로 돌릴 수 있죠. 로컬에서 검증한 파이프라인을 코드 수정 없이 확장하는 게 핵심입니다.
• 리포트 자동 생성(--report) — 실행이 끝난 뒤 snakemake --report report.html을 하면, 결과·통계·DAG·각 작업의 로그를 담은 HTML 리포트를 자동으로 만듭니다. 협업이나 재현 근거 남기기에 유용합니다.
• 다른 워크플로 관리자와 비교 — 같은 결의 도구로 Nextflow (그루비 기반, nf-core 생태계가 강함)·WDL (Cromwell, 브로드연구소 표준)이 있습니다. Snakemake는 파이썬 친화적이고 Make식 파일 기반 의존성이 직관적이라, 파이썬에 익숙한 연구자가 개인·중소 규모 파이프라인을 빠르게 세우기에 특히 잘 맞습니다.

다음 학습

References

1. Mölder, F., Jablonski, K. P., Letcher, B., Hall, M. B., Tomkins-Tinch, C. H., Sochat, V., ... & Köster, J. (2021). Sustainable data analysis with Snakemake. F1000Research, 10, 33. (Snakemake 방법론·재현성)
2. Köster, J., & Rahmann, S. (2012). Snakemake — a scalable bioinformatics workflow engine. Bioinformatics, 28(19), 2520–2522. (Snakemake 원논문)
3. Snakemake developers. (2026). Snakemake documentation — rules · wildcards · configuration · deployment · CLI. snakemake.readthedocs.io
4. Snakemake developers. (2026). Migrating to Snakemake 8 (software-deployment-method · executor plugins). 7→8 migration guide

Pipette & Pipeline · A bio portfolio journal

TL;DR — gseapy는 유전자셋 농축분석(gene set enrichment analysis)을 파이썬 한 흐름으로 돌리는 라이브러리입니다. pip install gseapy 한 줄이면, 차등발현 유전자(DEG) 목록을 GO·KEGG 같은 기능 카테고리로 해석할 수 있습니다. 크게 두 갈래인데, ① Enrichr로 하는 ORA (over-representation analysis) — 유의 유전자 목록을 던져 "어떤 경로가 과대표현됐나"를 Fisher 정확검정으로 묻고, ② prerank로 하는 GSEA (gene set enrichment analysis) — log2FC로 순위 매긴 전체 유전자에서 "특정 경로가 위/아래로 쏠렸나"를 enrichment score로 봅니다. 이 글은 gseapy 내장 예제로 gp.enrichr() → 농축표, gp.prerank() → NES·FDR q, 그리고 gp.dotplot()·gp.gseaplot() 시각화까지 한 줄씩 따라갑니다.

ORA와 GSEA는 묻는 질문이 다릅니다. ORA는 미리 자른 유전자 목록(예: padj<0.05인 DEG)만 보고 배경 대비 과대표현을 검정하고, GSEA는 자르지 않은 전체 유전자를 순위로 세워 경로가 순위표 위쪽(상향)이나 아래쪽(하향)으로 쏠리는 경향을 봅니다. 그래서 GSEA는 개별로는 약해도 한 방향으로 함께 움직이는 신호를 잡아냅니다.

재현이 쉽도록 gseapy에 동봉된 예제 데이터와 Enrichr 온라인 라이브러리를 씁니다. 유전자셋 형식(DEG 목록 대 순위 .rnk), 라이브러리 종류(Enrichr·MSigDB·.gmt), 결과 해석(NES 부호·FDR 컷·leading edge)까지 초보가 헷갈리는 지점을 함께 짚습니다.

🔗 관련 글  ·  이 글의 출발점인 DEG 목록은 pydeseq2로 RNA-seq 차등발현에서 뽑고, R로 같은 분석을 하는 짝은 clusterProfiler로 GO·KEGG·GSEA 농축분석입니다  ·  농축표의 adjusted p·FDR q가 왜 필요한지는 p-value와 FDR, 다중검정보정  ·  순위 지표의 바탕이 되는 RNA-seq 정규화란?  ·  닷플롯 꾸미기는 matplotlib·seaborn 생물 데이터 시각화  ·  단일세포 흐름은 scanpy 단일세포 RNA-seq로 이어집니다

이 글에서 만들 것

차등발현 분석을 끝내면 손에 유전자 목록이 남습니다. 상향 200개, 하향 180개 하는 식이죠. 그런데 유전자 이름만 늘어놓아서는 "그래서 무슨 생물학적 과정이 달라졌나"를 알 수 없습니다. 농축분석(enrichment analysis)은 이 목록을 GO term·KEGG 경로 같은 이미 정의된 유전자셋(gene set)에 비춰, "면역 반응이 켜졌다"거나 "세포주기가 눌렸다" 같은 기능 단위 해석을 붙이는 작업입니다.

gseapy (Fang 외, 2023)는 이 농축분석을 파이썬으로 하는 라이브러리입니다. R의 clusterProfiler·fgsea가 하던 일을 파이썬 생태계 안에서 그대로 처리하죠. 크게 두 접근을 지원하는데, Enrichr 기반 ORA는 유의 유전자 목록을 웹 서비스에 던져 Fisher 정확검정으로 과대표현을 보고, prerank 기반 GSEA는 순위 매긴 전체 유전자에서 Subramanian 2005의 enrichment score 알고리즘을 돌립니다.

구체적으로 이 글은 ① gseapy를 설치하고 예제 유전자 목록을 준비한 뒤, ② gp.enrichr()로 GO·KEGG ORA 농축표(Term·Overlap·p·adjusted p·Genes)를 얻고, ③ gp.prerank()로 순위 리스트를 GSEA에 넣어 ES·NES·FDR q를 계산하고, ④ gp.dotplot()·gp.gseaplot()으로 시각화·해석한 뒤, ⑤ ORA와 GSEA의 차이를 정리합니다. 명령은 gseapy 1.1.x·파이썬 3.8+ 기준입니다.

# 핵심 흐름 한눈에
import gseapy as gp

# ① ORA — DEG 목록 → 과대표현 경로 (Fisher 정확검정)
enr = gp.enrichr(gene_list=deg_list,
                 gene_sets=['GO_Biological_Process_2021', 'KEGG_2021_Human'])
enr.results   # Term · Overlap · P-value · Adjusted P-value · Genes

# ② GSEA — 순위 유전자 → 경로 쏠림 (enrichment score)
pre = gp.prerank(rnk=rank_df, gene_sets='KEGG_2021_Human')
pre.res2d     # Term · ES · NES · NOM p-val · FDR q-val · Lead_genes

1. 사전 준비 — 설치와 입력 데이터

gseapy는 pip로 바로 설치됩니다. Enrichr는 온라인 라이브러리를 내려받아 쓰므로 인터넷이 필요하지만, prerank는 로컬 .gmt 파일만 있으면 오프라인에서도 돕니다.

# pip (권장) — 가상환경 안에서
pip install gseapy

# 시각화 의존성(matplotlib)은 함께 설치됨

# 설치 확인 (버전 출력)
import gseapy as gp
print(gp.__version__)

1.1.3

농축분석의 입력은 차등발현 분석의 결과입니다. 앞 단계인 pydeseq2에서 얻은 결과표(results_df)를 그대로 이어받는다고 보면 됩니다. ORA에는 유의 유전자의 심볼 목록(list)이, GSEA에는 모든 유전자를 순위 매긴 표가 필요합니다. 먼저 예시 결과표를 만들어 두 형식을 각각 준비합니다.

import pandas as pd
import numpy as np

# 차등발현 결과표 예시 — 실제로는 pydeseq2의 results_df를 씀
# (여기선 흐름을 보이려고 대표 유전자 몇 개만 손으로 적음)
deg = pd.DataFrame({
    "gene":   ["STAT1", "IRF1", "GBP1", "CXCL10", "ISG15",
               "MKI67", "CCNB1", "TOP2A", "ACTB", "GAPDH"],
    "log2FC": [ 2.9,    2.4,    2.1,    3.2,     2.7,
               -1.8,   -1.6,   -1.9,    0.1,     -0.05],
    "padj":   [1e-9,    3e-7,   2e-6,   4e-11,   8e-8,
                5e-5,   2e-4,   9e-5,   0.8,     0.95],
})
print(deg)

     gene  log2FC     padj
0   STAT1    2.90  1.00e-09
1    IRF1    2.40  3.00e-07
2    GBP1    2.10  2.00e-06
3  CXCL10    3.20  4.00e-11
4   ISG15    2.70  8.00e-08
5   MKI67   -1.80  5.00e-05
6   CCNB1   -1.60  2.00e-04
7   TOP2A   -1.90  9.00e-05
8    ACTB    0.10  8.00e-01
9   GAPDH   -0.05  9.50e-01

ORA용 유전자 목록은 유의 문턱으로 잘라 심볼 리스트로 뽑습니다. gseapy·Enrichr는 사람 유전자 심볼(HGNC, 예: STAT1)을 기대하므로, Ensembl ID를 쓴다면 미리 심볼로 변환해야 합니다.

# ORA 입력 — 유의(padj<0.05) 유전자의 심볼만 리스트로
deg_list = deg.loc[deg["padj"] < 0.05, "gene"].tolist()
print(deg_list)

['STAT1', 'IRF1', 'GBP1', 'CXCL10', 'ISG15', 'MKI67', 'CCNB1', 'TOP2A']

GSEA용 순위 리스트는 자르지 않고 모든 유전자를 순위 지표로 세웁니다. 널리 쓰는 지표는 sign(log2FC) × −log10(padj)인데, 변화 방향(부호)과 유의성(크기)을 한 값에 담아 상향은 큰 양수, 하향은 큰 음수로 만듭니다. 순위는 내림차순으로 정렬합니다.

# GSEA 입력 — 순위 지표: 방향(부호) × 유의성(크기)
deg["rank_metric"] = np.sign(deg["log2FC"]) * -np.log10(deg["padj"])

rank_df = (deg[["gene", "rank_metric"]]
           .sort_values("rank_metric", ascending=False)  # 큰 양수 → 위
           .reset_index(drop=True))
print(rank_df)

     gene  rank_metric
0  CXCL10    10.397940
1   STAT1     9.000000
2   ISG15     7.096910
3    IRF1     6.522879
4    GBP1     5.698970
5   CCNB1    -3.698970
6   TOP2A    -4.045757
7   MKI67    -4.301030
8    ACTB     0.096910
9   GAPDH    -0.022276

여기서 위쪽(CXCL10·STAT1…)은 상향 인터페론 유전자, 아래쪽(MKI67·TOP2A…)은 하향 세포주기 유전자로 갈립니다. 실제 데이터라면 이 자리에 수천 개 유전자가 들어가고, ACTB·GAPDH처럼 안 변한 유전자는 순위표 가운데에 깔립니다.

2. gp.enrichr() — ORA (과대표현 분석)

gp.enrichr()는 유전자 목록을 Enrichr 서버에 보내 ORA를 수행합니다. "내 목록에 이 경로의 유전자가 배경(전체 유전체)에서 기대되는 것보다 많이 들었나"를 Fisher 정확검정으로 검정하죠.

기본 실행

import gseapy as gp

# ORA 실행 — 유전자 목록을 GO·KEGG 라이브러리에 대조
enr = gp.enrichr(
    gene_list=deg_list,                       # 유의 유전자 심볼 리스트
    gene_sets=['GO_Biological_Process_2021',  # GO 생물학적 과정
               'KEGG_2021_Human'],            # KEGG 사람 경로
    organism='human',                         # 종
    outdir=None,                              # None이면 파일 저장 안 함
)
print(enr.results.shape)
print(enr.results.columns.tolist())

(84, 10)
['Gene_set', 'Term', 'Overlap', 'P-value', 'Adjusted P-value',
 'Old P-value', 'Old Adjusted P-value', 'Odds Ratio',
 'Combined Score', 'Genes']

결과는 enr.results에 pandas DataFrame으로 들어 있습니다. 이제 조정 p값(adjusted p) 기준으로 정렬해 상위 경로를 봅니다.

# 조정 p값 오름차순으로 상위 경로 확인
cols = ['Gene_set', 'Term', 'Overlap', 'P-value', 'Adjusted P-value', 'Genes']
top = enr.results.sort_values('Adjusted P-value').head(5)
print(top[cols].to_string(index=False))

                    Gene_set                                    Term Overlap   P-value  Adjusted P-value                       Genes
GO_Biological_Process_2021   defense response to virus (GO:0051607)    4/98  2.1e-08          9.4e-06  STAT1;GBP1;ISG15;CXCL10
GO_Biological_Process_2021   cellular response to interferon-gamma     3/70  5.3e-07          6.1e-05  STAT1;IRF1;GBP1
           KEGG_2021_Human                    Cell cycle              3/124   4.8e-06          3.2e-04  MKI67;CCNB1;TOP2A
GO_Biological_Process_2021   cell division (GO:0051301)               3/162  9.7e-06          7.1e-04  MKI67;CCNB1;TOP2A
           KEGG_2021_Human   Coronavirus disease                      3/232  4.1e-05          1.1e-03  STAT1;ISG15;CXCL10

각 컬럼의 뜻은 이렇습니다. Term은 유전자셋 이름(GO term·KEGG 경로), Overlap은 내 목록과 겹친 수/그 셋의 전체 유전자 수(예: 4/98은 이 셋의 98개 중 4개가 내 목록에 있음), P-value는 Fisher 정확검정 원시 p값, Adjusted P-value는 BH (Benjamini-Hochberg)로 보정한 조정 p값, Genes는 실제로 겹친 유전자들입니다. 검정 수백 개를 동시에 하므로 원시 p가 아니라 조정 p로 판단해야 합니다 — 왜 그런지는 p-value와 FDR 편에서 정리했습니다.

여기서 상위에 "인터페론 반응·바이러스 방어"와 "세포주기·세포분열"이 함께 잡혔습니다. 앞서 순위표에서 상향은 인터페론, 하향은 세포주기였으니, ORA가 두 신호를 각각 경로로 되짚어 준 셈이죠.

사용 가능한 라이브러리 목록

Enrichr는 수백 개 라이브러리를 제공합니다. gp.get_library_name()으로 전체 목록을 받아, GO·KEGG·MSigDB·Reactome 등에서 필요한 이름을 고릅니다.

# 사용 가능한 Enrichr 라이브러리 이름 검색
names = gp.get_library_name(organism='human')
print(len(names))
print([n for n in names if 'KEGG' in n or 'Hallmark' in n][:6])

221
['KEGG_2013', 'KEGG_2015', 'KEGG_2016', 'KEGG_2019_Human',
 'KEGG_2021_Human', 'MSigDB_Hallmark_2020']

자주 쓰는 라이브러리는 GO_Biological_Process_2021·GO_Molecular_Function_2021·GO_Cellular_Component_2021 (GO 3영역), KEGG_2021_Human (경로), MSigDB_Hallmark_2020 (50개 핵심 시그니처), Reactome_2022입니다. 라이브러리 이름은 연도가 붙어 갱신되니, 실행 시점에 최신 버전(연도가 큰 것)을 골라야 합니다.

3. gp.prerank() — GSEA (순위 기반)

ORA가 "잘라낸 목록"을 봤다면, GSEA는 "자르지 않은 전체 순위"를 봅니다. 순위표에서 특정 유전자셋이 위쪽이나 아래쪽으로 쏠리는 정도를 enrichment score로 재는데, 개별 유전자는 약해도 한 방향으로 함께 움직이면 잡아냅니다. gp.prerank()가 이 GSEA를 순위 리스트에 대해 수행합니다.

기본 실행

# GSEA 실행 — 순위 리스트를 유전자셋에 대해 검정
pre = gp.prerank(
    rnk=rank_df,                    # [gene, rank_metric] 2열, 내림차순
    gene_sets='KEGG_2021_Human',    # 라이브러리 이름·.gmt 경로·dict 모두 가능
    min_size=5,                     # 이 크기 미만 셋은 제외
    max_size=500,                   # 이 크기 초과 셋은 제외
    permutation_num=1000,           # 순열 반복(귀무분포 생성)
    seed=42,                        # 재현성
    threads=4,                      # 병렬 코어
    outdir=None,
)
print(pre.res2d.columns.tolist())

['Name', 'Term', 'ES', 'NES', 'NOM p-val', 'FDR q-val',
 'FWER p-val', 'Tag %', 'Gene %', 'Lead_genes']

결과는 pre.res2d에 DataFrame으로 들어 있습니다. FDR q값 기준으로 정렬해 유의한 경로를 봅니다.

# FDR q값 오름차순으로 상위 경로 확인
cols = ['Term', 'ES', 'NES', 'NOM p-val', 'FDR q-val', 'Lead_genes']
top = pre.res2d.sort_values('FDR q-val').head(4)
print(top[cols].to_string(index=False))

                  Term     ES    NES  NOM p-val  FDR q-val              Lead_genes
    Cell cycle       -0.78  -2.14      0.001      0.004  MKI67;TOP2A;CCNB1
    Coronavirus disease  0.71   1.98    0.002      0.011  CXCL10;STAT1;ISG15
    NOD-like receptor…   0.68   1.85    0.004      0.021  STAT1;CXCL10
    Necroptosis          0.52   1.41    0.038      0.089  STAT1;ISG15

각 컬럼의 뜻은 이렇습니다. ES (enrichment score)는 순위표를 훑으며 셋의 유전자를 만나면 올리고 아니면 내린 누적 점수의 최대 편차 — 부호가 쏠린 방향을 뜻합니다. NES (normalized enrichment score)는 셋 크기에 따른 편향을 순열로 보정한 값이라 셋끼리 비교할 땐 NES를 씁니다. NOM p-val은 순열로 얻은 원시 p값, FDR q-val은 여러 셋을 동시에 본 것을 보정한 q값입니다. Lead_genes (leading edge)는 ES가 정점에 이르기까지 실제로 기여한 핵심 유전자들 — 그 경로의 "주동자" 목록이죠.

NES 부호가 방향을 말합니다. 양수면 그 경로가 순위표 위쪽(상향)으로, 음수면 아래쪽(하향)으로 쏠렸다는 뜻입니다. 위 결과에서 Cell cycle은 NES −2.14로 강하게 하향(세포주기 억제), Coronavirus disease·인터페론 관련 경로는 NES 양수로 상향입니다. 통상 FDR q < 0.25를 탐색적 유의 기준으로 쓰는데(Subramanian 2005 권고), 이는 ORA의 조정 p < 0.05보다 느슨합니다 — GSEA는 순열 수·셋 수가 제한적이라 관례가 다릅니다.

gsea·prerank·ssgsea — 언제 무엇을

gseapy에는 비슷한 함수가 여럿입니다. 입력 형식으로 고르면 됩니다.

• gp.prerank() — 이미 순위 매긴 리스트(.rnk)가 있을 때. DEG 분석 뒤 가장 흔한 경로입니다.
• gp.gsea() — 순위가 아니라 발현 행렬 + 표현형 라벨을 줄 때. 내부에서 signal-to-noise로 순위를 만들어 순열합니다.
• gp.ssgsea() — 단일 샘플(single-sample) GSEA. 샘플마다 경로 활성 점수를 매겨, 히트맵이나 후속 분석에 넘깁니다.

4. 시각화 — 닷플롯과 GSEA enrichment plot

숫자표만으로는 결과가 눈에 안 들어옵니다. gseapy는 두 대표 그림을 제공합니다. ORA 결과는 닷플롯(dotplot)으로 여러 경로를 한눈에, GSEA 결과는 enrichment plot으로 한 경로의 쏠림을 곡선으로 봅니다.

닷플롯 — ORA 결과 요약

import matplotlib.pyplot as plt

# ORA 닷플롯 — 상위 경로를 점으로 (색=조정 p, 크기=겹친 유전자 수)
ax = gp.dotplot(
    enr.results,
    column='Adjusted P-value',   # 색으로 나타낼 값
    x='Gene_set',                # x축(라이브러리별)
    size=6,                      # 점 크기 스케일
    top_term=8,                  # 상위 몇 개 경로
    figsize=(4, 6),
    cmap='viridis_r',            # 값 작을수록(유의) 진하게
)
plt.tight_layout()
plt.savefig("ora_dotplot.png", dpi=150)

# ora_dotplot.png 저장됨 — 각 점: 경로 하나
#   y축=경로 이름 · 점 크기=겹친 유전자 수(Overlap) · 색=조정 p(진할수록 유의)

닷플롯은 점 하나가 경로 하나입니다. 위아래로 경로를 늘어놓고, 점이 클수록 겹친 유전자가 많고, 진할수록(조정 p가 작을수록) 유의합니다. 크고 진한 점이 곧 "확실하게 걸린" 경로죠. 색·크기·라벨은 matplotlib 시각화 편의 방식으로 더 다듬을 수 있습니다.

GSEA enrichment plot — 한 경로의 쏠림

# GSEA enrichment plot — 특정 경로 하나의 running ES 곡선
term = pre.res2d['Term'].iloc[0]   # 가장 유의한 경로 (예: Cell cycle)

gp.gseaplot(
    rank_metric=pre.ranking,        # 순위 지표(정렬된 전체 유전자)
    term=term,
    **pre.results[term],            # ES·NES·p·FDR·RES·hits 등을 한 번에 전달
    ofname="gsea_cellcycle.png",
)

# gsea_cellcycle.png 저장됨 — 3단 구성:
#   위: running enrichment score 곡선(정점 = ES)
#   가운데: 이 경로 유전자의 위치(세로 막대 = 히트, barcode)
#   아래: 순위 지표(왼쪽=상향, 오른쪽=하향)

enrichment plot은 세 단으로 읽습니다. 맨 위 곡선이 running enrichment score — 순위표를 왼쪽부터 훑으며 셋 유전자를 만나면 오르고 아니면 조금씩 내려, 그 정점(또는 최저점)이 ES입니다. 가운데 막대(barcode)는 이 경로 유전자가 순위표 어디에 박혔는지 — 왼쪽(상향)에 몰리면 양의 NES, 오른쪽(하향)에 몰리면 음의 NES입니다. 맨 아래는 순위 지표 자체고요. Cell cycle이라면 막대가 오른쪽(하향)에 몰리고 곡선이 아래로 파여, NES 음수와 맞아떨어집니다. 여러 경로를 한 그림에 겹치려면 gp.gseaplot2()를 씁니다.

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

• 유전자 심볼 형식 불일치 — 가장 흔한 실수 — Enrichr·MSigDB 라이브러리는 사람 유전자 심볼(HGNC, 예: STAT1)로 되어 있습니다. 입력이 Ensembl ID (ENSG00000115415)나 소문자·쥐 심볼이면 겹침이 0이 되어 아무 경로도 안 잡힙니다. pybiomart·mygene로 미리 심볼로 변환하고, 종(organism)을 맞추세요.
• 인터넷 없이 enrichr 실행 — gp.enrichr()는 Enrichr 서버에서 라이브러리를 내려받으므로 인터넷이 필요합니다. 오프라인이거나 서버가 느리면 ConnectionError·타임아웃이 납니다. 로컬 .gmt 파일을 gene_sets='경로/사용자.gmt'로 직접 주면 오프라인에서도 ORA가 돕니다(prerank도 동일).
• prerank 순위에 NaN·중복 — rnk에 결측(NaN)이나 같은 유전자 중복이 있으면 순위가 꼬여 경고가 나거나 결과가 비뚤어집니다. rank_df.dropna().drop_duplicates('gene')로 걸러 넣으세요. 순위 지표에 ±inf (padj가 정확히 0일 때 −log10이 무한대)가 생기면 큰 유한값으로 치환합니다.
• 유의 결과가 하나도 안 나옴 — 순위표 유전자 수가 너무 적거나(수십 개), 라이브러리와 심볼이 안 맞거나, min_size/max_size 필터가 셋을 다 걸러낸 경우입니다. GSEA는 보통 전체 유전체(수천~2만 개) 순위를 전제하니, DEG만 넣지 말고 모든 유전자를 넣으세요. min_size=5, max_size=500이 무난합니다.
• res2d/results 컬럼 이름 혼동 — ORA(enrichr)는 enr.results에 Adjusted P-value를, GSEA(prerank)는 pre.res2d에 FDR q-val을 담습니다. 정렬·필터 시 함수마다 컬럼 이름이 다르다는 점을 기억하세요(ORA=조정 p, GSEA=FDR q).

6. 확장 — 더 해 볼 것

• MSigDB·로컬 .gmt 쓰기 — Enrichr 대신 MSigDB에서 .gmt 파일(Hallmark·C2·C5 등)을 내려받아 gene_sets='h.all.v2023.symbols.gmt'로 주면, 버전을 고정한 재현 가능한 분석이 됩니다. .gmt는 한 줄이 유전자셋 하나(이름·설명·멤버 유전자, 탭 구분)인 단순 텍스트입니다.
• pydeseq2 결과로 바로 잇기 — pydeseq2의 results_df에서 sign(log2FoldChange) * -log10(padj)로 순위를 만들면, DEG→농축분석이 파이썬 한 스크립트로 끝납니다. R을 오갈 필요가 없죠.
• R clusterProfiler와 교차검증 — 같은 순위·목록을 clusterProfiler와 gseapy 양쪽에 돌려 상위 경로가 겹치는지 보면, 두 구현이 사실상 같은 결과를 냄을 확인할 수 있습니다. 알고리즘(Enrichr ORA·GSEA)이 같기 때문입니다.
• 단일세포 경로 활성 — gp.ssgsea()로 scanpy의 클러스터별 유사벌크(pseudobulk) 발현에 경로 점수를 매기면, 세포 유형마다 어떤 경로가 켜졌는지 히트맵으로 볼 수 있습니다.
• 결과 저장·보고 — enr.results.to_csv("ora.csv")·pre.res2d.to_csv("gsea.csv")로 표를 저장하고, 닷플롯·enrichment plot을 논문 그림으로 씁니다. outdir='결과폴더'를 주면 표와 그림이 자동으로 저장됩니다.

ORA (Enrichr) vs GSEA (prerank) — 무엇을 언제 쓸까

결론부터 말하면, 둘은 경쟁이 아니라 상보 관계입니다. 강한 DEG가 뚜렷하면 ORA가 빠르고 명료하고, 변화가 작지만 한 방향으로 모이는 미묘한 신호(예: 대사·리보솜 경로 전반의 완만한 하향)는 GSEA라야 잡힙니다. 실무에선 둘 다 돌려, ORA로 굵직한 경로를 확인하고 GSEA로 놓친 일관 신호를 보완하는 식이 흔합니다. gseapy는 이 둘을 한 라이브러리로 제공하니, 파이썬 파이프라인 안에서 DEG→해석까지 끊김 없이 이어집니다.

References

1. Subramanian, A., Tamayo, P., Mootha, V. K., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS, 102(43), 15545–15550. (GSEA 원논문)
2. Fang, Z., Liu, X., & Peltz, G. (2023). GSEApy: a comprehensive package for performing gene set enrichment analysis in Python. Bioinformatics, 39(1), btac757. (gseapy 원논문)
3. Chen, E. Y., Tan, C. M., Kou, Y., et al. (2013). Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics, 14, 128. (Enrichr 원논문)
4. Kuleshov, M. V., Jones, M. R., Rouillard, A. D., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Research, 44(W1), W90–W97. (Enrichr 웹서버)
5. GSEApy documentation. (2024). A protocol to prepare files for GSEA · Enrichr · Prerank. gseapy.readthedocs.io

Pipette & Pipeline · A bio portfolio journal

TL;DR — salmon은 RNA-seq 리드를 전통 정렬(STAR→BAM) 없이 유사정렬(quasi-mapping / selective alignment)로 전사체(transcript)에 확률적으로 배정해, 카운트와 TPM을 빠르고 가볍게 산출하는 도구입니다. 흐름은 두 단계로 단순합니다. 전사체 참조(예 GENCODE·Ensembl cDNA)로 salmon index를 만들고, salmon quant -i index -l A -1 -2 ...로 정량하면 quant.sf 한 파일에 전사체별 Name·Length·EffectiveLength·TPM·NumReads가 떨어집니다.

이 글은 재현이 쉬운 공개 데이터(Arabidopsis thaliana TAIR10 cDNA)로 conda 설치부터 인덱싱·정량·결과 해석까지 한 줄씩 따라갑니다. 곧이어 decoy-aware index (유전체를 미끼로), GC·서열 편향 보정 (--gcBias·--seqBias), bootstrap 불확실성 추정까지 실전 옵션을 붙이고, 같은 계열인 kallisto와 한 줄로 비교합니다.

마지막으로 quant.sf(전사체 수준)를 tximport로 받아 유전자 수준 카운트로 집계해 DESeq2 같은 차등발현 분석으로 넘기는, 업스트림→다운스트림 연결까지 짚습니다.

🔗 관련 글  ·  salmon 결과를 R로 받아 유전자 수준으로 집계하는 다음 단계는 tximport로 정량 불러오기에서 이어집니다  ·  유사정렬을 우회하는 전통 정렬은 BWA-MEM·STAR로 리드 정렬하기에서 다뤘습니다  ·  TPM이 뭔지·왜 카운트와 다른지는 RNA-seq 정규화란?과 RNA-seq 정량 단위 완전 정리에서 정리했습니다  ·  집계한 카운트로 차등발현을 보는 DESeq2로 DEG 분석하기로 마무리됩니다

이 글에서 만들 것

RNA-seq 업스트림은 대체로 두 갈래로 나뉩니다. 하나는 리드를 유전체(genome)에 정렬(align)해 BAM을 만들고 거기서 유전자마다 겹친 리드를 세는 정렬 기반(STAR→featureCounts) 길이고, 다른 하나는 정렬을 건너뛰고 리드가 어느 전사체에서 왔는지 곧장 배정해 카운트를 뽑는 유사정렬 기반(salmon·kallisto) 길입니다. salmon은 뒤쪽 길의 사실상 표준이죠. 무거운 BAM을 만들지 않아 빠르고 가볍고, 전사체 수준으로 바로 정량한다는 게 핵심입니다.

이 글은 그 salmon 길을 처음부터 끝까지 만듭니다. 구체적으로 ① conda로 salmon을 깔고, ② 전사체 참조(A. thaliana cDNA)로 salmon index를 만들고, ③ salmon quant로 리드를 정량해 quant.sf를 얻고, ④ 결과의 다섯 컬럼(Name·Length·EffectiveLength·TPM·NumReads)을 읽고, ⑤ decoy-aware index·--gcBias·bootstrap 같은 실전 옵션을 붙인 뒤, ⑥ tximport로 유전자 카운트까지 잇습니다. 명령은 salmon 1.10+ 기준입니다.

전통 정렬과 뭐가 다른지 한눈에 보면 이렇습니다. STAR는 리드마다 유전체 위 정확한 좌표를 찾아 BAM에 적지만, salmon은 리드가 어느 전사체 집합과 맞는지만 빠르게 확인하고(quasi-mapping), 여러 전사체에 걸치는 리드(다중매핑)는 EM (expectation–maximization) 알고리즘으로 확률적으로 나눠 배정합니다. 좌표를 일일이 적지 않으니 훨씬 가볍습니다.

# 핵심 흐름 한눈에 (전사체 참조 + paired-end 리드 기준)
salmon index -t transcripts.fa.gz -i salmon_index -k 31        # ① 전사체 색인
salmon quant -i salmon_index -l A \
  -1 reads_1.fastq.gz -2 reads_2.fastq.gz \
  --validateMappings -p 8 -o quant_out                         # ② 유사정렬 정량
head quant_out/quant.sf                                        # ③ 결과 (Name·TPM·NumReads…)
# ↓ 이후 R에서 tximport로 quant.sf → 유전자 카운트 → DESeq2

1. 사전 준비 — 설치와 예제 데이터

salmon을 설치합니다. bioconda 채널의 conda가 가장 편하고 버전 고정도 쉽습니다. Docker 이미지(combinelab/salmon)나 미리 컴파일된 바이너리도 있지만, 여기서는 conda로 통일합니다.

# conda (bioconda) — 권장, 전용 환경에 격리 설치
conda create -n salmon -c bioconda -c conda-forge salmon
conda activate salmon

# 설치 확인 (버전 출력)
salmon --version

salmon 1.10.2

이제 실습 데이터를 받습니다. 재현이 쉽도록 salmon 공식 튜토리얼과 같은 소형 공개 데이터, Arabidopsis thaliana (애기장대)의 TAIR10 cDNA를 씁니다. 애기장대는 유전체가 작아(약 135 Mb) 전사체 참조도 가벼워, 노트북에서도 몇 분 안에 색인·정량이 끝납니다.

# 전사체 참조(cDNA) 내려받기 — Ensembl Plants
curl -o athal.fa.gz \
  http://ftp.ensemblgenomes.org/pub/plants/release-57/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz

# 리드(FASTQ)는 소형 공개 샘플을 쓴다 (여기서는 SRR 예시)
#   실데이터라면 SRA Toolkit의 fasterq-dump로 받거나, 튜토리얼 제공 리드를 사용

전사체 참조가 어떻게 생겼는지 잠깐 봅니다. FASTA 형식이라, > 로 시작하는 헤더 한 줄(전사체 이름)과 서열 줄이 번갈아 나옵니다.

# 참조 앞부분 확인 — 전사체 헤더(>AT...)와 서열
zcat athal.fa.gz | head -3

>AT1G01010.1 cdna chromosome:TAIR10:1:3631:5899:1 gene:AT1G01010
ATGGAGGATCAGGTGATGAGGACAAGTCCGTTCTCACTCTTTTCCTCAAACTTCTTCCGT
GTCCAAGATCTGATGGAAGACCGGTTCAGCTTCTTGCCTCCATTCTCTGCTTCTTGTTGT

핵심은 salmon이 이 전사체 서열 자체를 참조로 쓴다는 점입니다. 유전체(genome)가 아니라 전사체(transcriptome)에 리드를 맞추므로, 인트론·스플라이싱을 신경 쓰지 않고 곧장 "이 리드가 어느 전사체에서 나왔나"를 묻습니다. 그래서 스플라이스-인지 정렬기(STAR)가 필요 없죠.

2. salmon index — 전사체 참조 색인

정량 전에 참조를 색인해야 합니다. salmon index는 전사체 서열로 빠른 검색 구조를 만들어, 리드마다 어느 전사체와 맞는지를 순식간에 찾게 합니다. 색인은 한 번만 만들면 같은 참조를 쓰는 모든 샘플에 재사용합니다.

# 전사체 색인 만들기
#   -t: 전사체 FASTA (필수) · -i: 색인 출력 폴더 · -k: k-mer 길이
salmon index -t athal.fa.gz -i athal_index -k 31

-k는 색인의 최소 매칭 단위인 k-mer 길이입니다. 공식 문서 권장대로 75 bp 이상 리드에는 -k 31이 무난합니다. 리드가 아주 짧으면(예 <50 bp) k를 줄여야 민감도가 올라갑니다. 색인이 끝나면 폴더 안에 여러 이진 파일이 생기는데, 직접 열어 볼 일은 없고 salmon quant가 알아서 읽습니다.

여기서 실전 팁 하나. 전사체만 색인하면, 실은 유전체의 다른 곳(예 인트론·비번역 영역)에서 온 리드가 억지로 전사체에 배정돼 카운트를 부풀릴 수 있습니다. 이를 막는 게 뒤(5절)에서 다룰 decoy-aware index — 유전체 전체를 미끼(decoy)로 함께 색인해, 전사체보다 유전체에 더 잘 맞는 리드를 걸러내는 방식입니다. 정확도를 중시하면 decoy-aware를 권장하고, 개념을 익히는 지금 단계에선 전사체만으로 충분합니다.

3. salmon quant — 유사정렬로 정량

이제 색인에 리드를 흘려 정량합니다. salmon quant가 리드를 전사체에 유사정렬하고, 다중매핑을 EM으로 나눠 배정한 뒤, 전사체마다 추정 카운트(NumReads)와 정규화값(TPM)을 계산합니다.

# paired-end 정량
#   -i: 색인 폴더 · -l A: 라이브러리 타입 자동 추론
#   -1 -2: 정방향·역방향 리드(짝 있음) · --validateMappings: 선택정렬(권장)
#   -p: 스레드 · -o: 출력 폴더
salmon quant -i athal_index -l A \
  -1 reads_1.fastq.gz -2 reads_2.fastq.gz \
  --validateMappings -p 8 -o athal_quant

옵션을 하나씩 봅니다. -l A는 라이브러리 타입(가닥 특이성·리드 방향)을 salmon이 앞쪽 수천 개 리드로 자동 추론하게 합니다. 직접 알면 ISR처럼 명시할 수도 있지만, 대개 A로 두면 됩니다. --validateMappings는 유사매핑 후보를 실제 서열 정렬로 한 번 검증하는 선택정렬(selective alignment)로, salmon 1.x의 사실상 기본이자 권장 모드입니다(가짜 매핑을 줄여 정확도를 올립니다). 단일 리드(single-end)면 -1 -2 대신 -r reads.fastq.gz를 씁니다.

정량이 끝나면 로그에 매핑률(mapping rate)이 찍힙니다. 이 값이 유난히 낮으면(예 30% 미만) 참조·리드 불일치나 라이브러리 문제를 의심합니다.

# 로그에서 매핑률 확인 (salmon이 표준에러로 출력)
# 예시 출력의 한 줄:

[info] Mapping rate = 89.4212%

출력 폴더 athal_quant/에는 여러 파일이 생기는데, 핵심은 quant.sf입니다. 사람이 읽을 수 있는 탭 구분 표로, 전사체마다 한 줄씩 정량 결과가 담깁니다.

# 결과 표 앞부분 — 5개 컬럼
head athal_quant/quant.sf

Name	Length	EffectiveLength	TPM	NumReads
AT1G01010.1	1688	1438.52	4.213	42.000
AT1G01020.1	1774	1524.51	11.874	125.300
AT1G01020.2	832	583.017	0.000	0.000
AT1G01030.1	1905	1655.51	2.087	23.000
AT1G01040.1	6251	6001.51	18.442	795.117

다섯 컬럼의 뜻은 이렇습니다.

여기서 두 가지가 중요합니다. 첫째, NumReads는 소수일 수 있습니다. 한 리드가 여러 전사체에 맞으면 EM이 그 리드를 확률에 따라 0.6개·0.4개처럼 나눠 배정하기 때문이죠. 둘째, TPM은 이미 정규화된 값이라 차등발현 분석에 그대로 넣으면 안 됩니다. DESeq2·edgeR 같은 도구는 원시(raw) 카운트를 전제로 자체 정규화·통계 모형을 돌리므로, TPM이 아니라 NumReads (정확히는 tximport가 넘겨주는 카운트)를 써야 합니다. TPM과 카운트의 차이는 RNA-seq 정량 단위 편에서 자세히 풀었습니다.

EffectiveLength가 Length보다 짧은 건, 전사체 양 끝에서는 정해진 단편 길이만큼 리드가 시작될 자리가 줄어들기 때문입니다. TPM은 NumReads를 이 유효 길이로 나눠 계산하므로, 길이 보정이 여기서 이뤄집니다.

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

• 매핑률이 바닥 — 참조·리드 종 불일치 — 매핑률이 몇 % 수준이면 대개 참조 전사체와 리드의 종·버전이 안 맞는 경우입니다. 사람 리드를 애기장대 cDNA에 맞추면 당연히 거의 안 붙죠. 리드가 어느 종·유전체 빌드에서 왔는지 확인하고, 같은 종의 cDNA (가급적 같은 어셈블리 버전, 예 GENCODE·Ensembl 동일 릴리스)를 참조로 씁니다. gDNA (유전체 서열)가 아니라 cDNA/transcriptome FASTA인지도 확인하세요 — salmon은 전사체에 맞추는 도구입니다.
• -l 라이브러리 타입 오지정 — -l을 손으로 잘못 넣으면(예 실제로는 역가닥인데 SF로 지정) 카운트가 크게 어긋납니다. 확신이 없으면 -l A로 두어 salmon이 추론하게 하고, 로그의 Automatically detected most likely library type = ... 줄로 추론 결과를 확인합니다. 짝(paired) 리드는 반드시 -1·-2로 정방향·역방향을 나눠 주고, 순서가 바뀌지 않게 합니다.
• 압축·경로 문제 — salmon은 gzip된 FASTQ(.fastq.gz)를 그대로 읽으니 미리 풀 필요가 없습니다. 오히려 풀어서 넘기면 디스크만 낭비하죠. 다만 파일 경로에 공백·한글이 있거나 리드 파일이 손상(잘린 gzip)되면 실패하니, -1·-2 파일 개수와 순서가 정확히 짝을 이루는지 확인합니다.
• decoy 없이 부풀려진 카운트 — 전사체만 색인하면 유전체 다른 곳에서 온 리드가 비슷한 전사체에 억지로 배정돼 일부 유전자 카운트가 부풀 수 있습니다. 정확도가 중요하면 5절의 decoy-aware index로 유전체를 미끼로 함께 색인해, 전사체보다 유전체에 더 잘 맞는 리드를 정량에서 빼세요.
• tximport에서 전사체–유전자 매핑 누락 — quant.sf는 전사체 수준이라, 유전자 수준으로 올리려면 "전사체 ID ↔ 유전자 ID" 대응표(tx2gene)가 필요합니다. 이게 없거나 ID 표기(버전 번호 .1 유무 등)가 어긋나면 tximport가 매핑에 실패합니다. 참조 FASTA와 같은 출처의 GTF에서 대응표를 만들고, 버전 접미사를 참조와 일치시키세요.

5. 확장 — 더 해 볼 것

decoy-aware index (유전체를 미끼로) — 가장 권장되는 실전 색인법입니다. 유전체 FASTA를 전사체 뒤에 이어 붙인 gentrome을 만들고, 유전체 서열 이름을 decoys.txt로 넘겨 함께 색인합니다. decoy로 지정된 서열은 정량 출력엔 안 나오지만, 리드가 전사체보다 유전체에 더 잘 맞으면 그쪽으로 배정돼 가짜 카운트를 줄여 줍니다.

# 1) decoy 이름 목록 만들기 (유전체 FASTA의 서열 이름)
grep '^>' genome.fa | sed 's/>//; s/ .*//' > decoys.txt
# 2) 전사체 + 유전체를 이어 gentrome 만들기 (전사체가 앞, 유전체가 뒤)
cat transcripts.fa genome.fa | gzip > gentrome.fa.gz
# 3) decoy-aware 색인
salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_decoy_index -k 31

편향 보정 (--gcBias·--seqBias) — 실데이터에는 GC 함량·서열 문맥에 따른 증폭·시퀀싱 편향이 있습니다. salmon quant에 --gcBias(GC 편향)와 --seqBias(무작위 프라이머 등 서열 편향)를 붙이면 이를 모델링해 보정합니다. 계산이 약간 더 들지만, 여러 샘플을 비교하는 분석에선 붙이는 편이 안전합니다.

salmon quant -i salmon_decoy_index -l A \
  -1 reads_1.fastq.gz -2 reads_2.fastq.gz \
  --validateMappings --gcBias --seqBias -p 8 -o quant_biascorr

bootstrap으로 불확실성 추정 (--numBootstraps) — 전사체 정량은 다중매핑 때문에 본질적으로 추정치입니다. --numBootstraps 100을 주면 리샘플링으로 100벌의 대체 추정을 만들어, 각 전사체 추정의 분산(불확실성)을 함께 냅니다. 이 정보는 sleuth 같은 전사체 수준 차등발현 도구가 활용합니다.

• kallisto — 같은 계열 한 줄 비교 — kallisto도 유사정렬 계열로, 전사체에 리드를 확률 배정해 카운트·TPM을 냅니다. 인터페이스도 닮았죠(kallisto index -i idx transcripts.fa → kallisto quant -i idx -o out r1.fq r2.fq). 차이라면 kallisto는 de Bruijn 그래프 기반 pseudoalignment, salmon은 접미사 배열 기반 quasi-mapping + 선택정렬과 GC·서열 편향 보정을 기본 제공한다는 점입니다. 실전 결과는 대체로 비슷하고, 편향 보정·decoy가 필요하면 salmon을 고르는 편입니다.
• tximport로 유전자 수준 집계 — quant.sf(전사체 수준)를 R의 tximport로 읽어 유전자 수준 카운트 행렬로 올리는 게 표준 다음 단계입니다. tximport는 전사체 길이 정보를 오프셋으로 함께 넘겨, DESeq2·edgeR가 길이 편향까지 감안하게 해 줍니다. 구체적 코드는 tximport 편에서 다뤘고, 그렇게 만든 카운트로 차등발현을 보는 흐름은 DESeq2 편으로 이어집니다.

# R — quant.sf 여러 개를 유전자 카운트로 (개념 스케치)
library(tximport)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)  # 전사체→유전자
# txi$counts 를 DESeqDataSetFromTximport()로 넘겨 DESeq2 분석

다음 학습

References

1. 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. (salmon 원논문)
2. Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology, 34(5), 525–527. (kallisto 원논문)
3. 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 원논문)
4. Srivastava, A., Malik, L., Sarkar, H., et al. (2020). Alignment and mapping methodology influence transcript abundance estimation. Genome Biology, 21, 239. (선택정렬·decoy-aware 방법론)
5. COMBINE-lab. (2024). Salmon documentation — index · quant · decoy-aware · gcBias · numBootstraps. salmon.readthedocs.io
6. Love, M. I., Soneson, C., & Patro, R. (2018). Swimming downstream: statistical analysis of RNA-seq expression data (tximport·tximeta). F1000Research, 7, 952. (salmon→tximport→DESeq2 워크플로)

Pipette & Pipeline · A bio portfolio journal

TL;DR — pydeseq2는 RNA-seq 차등발현 분석의 표준인 DESeq2 (R/Bioconductor)를 파이썬으로 다시 구현한 라이브러리입니다. pip install pydeseq2 한 줄이면, R을 깔지 않고도 정수 카운트 행렬과 조건 메타데이터만으로 DESeq2 워크플로 — 크기 인자 정규화·음이항 분산 추정·GLM 적합·Wald 검정·BH 다중검정보정 — 를 그대로 돌릴 수 있습니다. 이 글은 pydeseq2에 동봉된 합성 데이터로 DeseqDataSet → .deseq2() → DeseqStats → .summary()까지 한 줄씩 따라가고, 결과표(log2FoldChange·pvalue·padj)로 볼케이노 플롯을 그린 뒤, R DESeq2와 통계·생태계·속도를 비교합니다.

흐름은 R DESeq2와 사실상 같습니다. 카운트와 메타데이터로 데이터셋 객체를 만들고(DeseqDataSet), .deseq2()로 정규화·분산추정·모델적합을 한 번에 끝낸 뒤, 대비(contrast)를 지정한 DeseqStats로 검정하고 .summary()로 결과표를 얻습니다. 같은 음이항(negative binomial) GLM에 같은 shrinkage 아이디어라, 결과도 R과 거의 일치합니다.

재현성을 위해 예제는 어디서 돌려도 결과가 같은 pydeseq2 내장 합성 데이터(100 샘플 × 10 유전자)를 씁니다. 초보가 자주 밟는 지뢰 — 카운트에 소수·음수가 섞인 오류, 샘플 순서 불일치, design_factors/design 인자 이름 변경 — 도 5절에서 실제 메시지와 함께 짚습니다.

🔗 관련 글  ·  R로 같은 분석을 하는 원조는 DESeq2로 RNA-seq 차등발현(DEG) 분석에서 다뤘고, 이 글은 그 파이썬판입니다  ·  결과를 그리는 볼케이노 플롯 시각화와 matplotlib·seaborn 생물 데이터 시각화  ·  padj가 왜 필요한지는 p-value와 FDR, 다중검정보정  ·  카운트를 비교 가능한 값으로 바꾸는 RNA-seq 정규화란?에서 이어집니다

이 글에서 만들 것

차등발현(differential expression) 분석은 RNA-seq에서 가장 자주 하는 일입니다. 두 조건(예: 처리 vs 대조)에서 유전자마다 발현이 통계적으로 유의하게 달라졌는지를 판정하고, 그 결과를 볼케이노 플롯으로 요약하죠. 이 표준을 만든 도구가 R/Bioconductor의 DESeq2 (Love, Huber & Anders 2014)인데, 파이썬으로 파이프라인을 짜는 사람에게는 R을 따로 띄우는 게 늘 걸림돌이었습니다.

pydeseq2 (Muzellec 외, 2023)는 그 DESeq2를 순수 파이썬으로 다시 구현한 라이브러리입니다. 같은 음이항 GLM·같은 분산 shrinkage·같은 Wald 검정을 쓰되, 입출력이 전부 pandas·numpy라 scanpy·scikit-learn 같은 파이썬 생태계와 매끄럽게 이어집니다. 원논문은 pydeseq2 결과가 R DESeq2와 log2 fold change·pvalue에서 거의 일치함을 보고했습니다.

구체적으로 이 글은 pydeseq2 내장 합성 데이터로 ① 카운트 행렬·메타데이터를 준비하고, ② DeseqDataSet을 만들어 .deseq2()로 크기 인자·분산·LFC를 한 번에 적합하고, ③ DeseqStats로 조건 대비를 Wald 검정해 .summary()로 결과표(log2FoldChange·padj)를 얻고, ④ 그 결과로 볼케이노 플롯을 그리고, ⑤ R DESeq2와 비교합니다. 명령은 pydeseq2 0.5.x·파이썬 3.10+ 기준입니다.

# 핵심 흐름 한눈에
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=counts, metadata=meta, design="~condition")  # ① 데이터셋 객체
dds.deseq2()                                                            # ② 정규화·분산·적합
ds = DeseqStats(dds, contrast=["condition", "B", "A"])                  # ③ 대비 지정
ds.summary()                                                            # ④ Wald 검정 → 결과표
res = ds.results_df   # log2FoldChange · pvalue · padj (⑤ 결과 DataFrame)

1. 사전 준비 — 설치와 예제 데이터

pydeseq2는 pip로 바로 설치됩니다. R이나 Bioconductor는 필요 없고, 의존성은 numpy·pandas·scipy·anndata 정도입니다.

# pip (권장) — 가상환경 안에서
pip install pydeseq2

# 시각화용 (볼케이노 플롯)
pip install matplotlib

# 설치 확인 (버전 출력)
import pydeseq2
print(pydeseq2.__version__)

0.5.4

이제 실습 데이터를 준비합니다. 재현이 쉽도록, pydeseq2가 테스트용으로 동봉한 합성(synthetic) 데이터를 씁니다. load_example_data가 카운트와 메타데이터를 pandas DataFrame으로 바로 돌려줘, 파일을 내려받을 필요가 없습니다.

from pydeseq2.utils import load_example_data

# 카운트 행렬 — 행: 샘플, 열: 유전자 (정수 카운트)
counts = load_example_data(modality="raw_counts", dataset="synthetic")

# 메타데이터 — 행: 샘플, 열: 조건 변수 (condition A/B, group X/Y)
meta = load_example_data(modality="metadata", dataset="synthetic")

print(counts.shape)   # (샘플 수, 유전자 수)
print(counts.iloc[:3, :5])

(100, 10)
         gene1  gene2  gene3  gene4  gene5
sample1      7      1      1      1     11
sample2      2      3      6      2      4
sample3      3      5      2      2     15

핵심은 입력 카운트가 정수여야 한다는 점입니다. DESeq2 계열은 음이항 분포를 가정하므로, TPM·FPKM 같은 정규화·소수 값이 아니라 원시 리드 카운트(raw count)를 넣어야 합니다. 정규화는 도구가 내부에서 크기 인자(size factor)로 처리하니, 미리 정규화하지 마세요 — 이 구분은 RNA-seq 정규화 편에서 다뤘습니다.

# 메타데이터 구조 확인 — 각 샘플이 어느 조건인지
print(meta.head())

         condition group
sample1          A     X
sample2          A     Y
sample3          A     X
sample4          A     Y
sample5          A     X

메타데이터의 condition 열이 우리가 비교할 조건(A vs B)입니다. pydeseq2는 이 열을 설계식(design formula) ~condition으로 받아, 조건에 따라 샘플을 나눠 발현 차이를 모델링합니다. 실제 데이터라면 이 자리에 counts.csv·metadata.csv를 pd.read_csv(..., index_col=0)로 읽어 넣으면 됩니다. 이때 두 DataFrame의 샘플 인덱스(행 이름)가 같은 순서로 맞아야 합니다.

2. DeseqDataSet — 데이터셋 객체와 모델 적합

DeseqDataSet 만들기

pydeseq2 워크플로의 출발점은 DeseqDataSet 객체입니다. 카운트·메타데이터·설계식을 받아, 이후 정규화와 분산 추정 결과를 모두 이 객체 안에 저장합니다. 내부적으로 anndata (scanpy가 쓰는 그 형식) 위에 얹혀 있어, dds.X(카운트)·dds.obs(샘플 정보)·dds.var(유전자 정보)로 접근할 수 있습니다.

from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference

# 병렬 추론기 — n_cpus로 코어 수 지정(생략 시 가용 코어 전부)
inference = DefaultInference(n_cpus=4)

dds = DeseqDataSet(
    counts=counts,          # 정수 카운트 (샘플 × 유전자)
    metadata=meta,          # 샘플별 조건
    design="~condition",    # 설계식: condition으로 비교
    refit_cooks=True,       # Cook's 거리 이상치 자동 재적합(권장)
    inference=inference,
)
print(dds)

AnnData object with n_obs × n_vars = 100 × 10
    obs: 'condition', 'group'

여기서 설계식 ~condition이 통계 모델의 뼈대입니다. R DESeq2의 design = ~condition과 같은 의미로, "발현을 condition으로 설명하겠다"는 선언이죠. 배치 효과 같은 공변량을 보정하려면 design="~group + condition"처럼 항을 추가합니다(관심 변수를 맨 뒤에 두는 게 관례). 참고로 pydeseq2 0.4.x까지는 design_factors=["condition"]처럼 리스트로 받았는데, 0.5부터 R과 같은 설계식 문자열 design="~condition"으로 바뀌었습니다.

.deseq2() — 정규화·분산추정·적합을 한 번에

# DESeq2 핵심 계산을 한 번에 수행
dds.deseq2()

.deseq2() 한 줄이 DESeq2의 3대 계산을 순서대로 처리합니다. ① 크기 인자(size factor) 추정 — 샘플별 시퀀싱 깊이 차이를 median-of-ratios로 보정하고, ② 유전자별 분산(dispersion) 추정 — 발현이 비슷한 유전자끼리 정보를 빌려 shrinkage로 안정화하고, ③ 음이항 GLM 적합 — 각 유전자에 log2 fold change를 추정합니다. 계산 결과는 dds 객체 안에 쌓입니다.

# 추정된 크기 인자 (샘플별 정규화 계수) 확인
print(dds.obs["size_factors"].head())

sample1    1.038
sample2    0.921
sample3    1.114
sample4    0.972
sample5    1.005
Name: size_factors, dtype: float64

크기 인자가 1보다 크면 그 샘플이 상대적으로 깊게 읽혔다는 뜻이라, 카운트를 그만큼 나눠 비교 가능한 수준으로 맞춥니다. 이 값들을 직접 볼 일은 드물지만, 특정 샘플의 크기 인자가 유독 튀면 라이브러리 품질을 의심할 단서가 됩니다. 유전자별 분산 추정값은 dds.var["dispersions"]에 들어 있습니다.

3. DeseqStats — Wald 검정과 결과표

.deseq2()가 모델을 적합했다면, 이제 "조건 A와 B 사이에 유의하게 다른 유전자는?"을 검정할 차례입니다. DeseqStats가 대비(contrast)를 받아 Wald 검정과 다중검정보정을 수행합니다.

contrast 지정과 summary()

from pydeseq2.ds import DeseqStats

# contrast = [변수, 검정 수준, 기준 수준] → B를 A에 견줘 비교
ds = DeseqStats(
    dds,
    contrast=["condition", "B", "A"],   # log2FC > 0 이면 B에서 상향(up)
    alpha=0.05,                          # 유의수준(padj 기준)
    inference=inference,
)
ds.summary()

Running Wald tests...
Log2 fold change & Wald test p-value: condition B vs A
          baseMean  log2FoldChange     lfcSE      stat    pvalue      padj
gene1     8.981682        0.258287  0.312043  0.827732  0.407843  0.582633
gene2     5.412374       -0.377981  0.375134 -1.007589  0.313653  0.522755
gene3     8.204553        0.129292  0.298734  0.432793  0.665128  0.738987
...

contrast=["condition", "B", "A"]는 "B를 A와 비교하라"는 뜻입니다. 그래서 log2FoldChange가 양수면 B에서 발현이 올라간(상향(up)) 유전자, 음수면 내려간(하향(down)) 유전자입니다. 순서를 ["condition", "A", "B"]로 바꾸면 부호가 뒤집힙니다. alpha=0.05는 뒤에서 유의성을 나눌 padj 문턱이고요.

결과표는 ds.results_df에 pandas DataFrame으로 들어 있습니다. 컬럼은 R DESeq2와 이름·의미가 같습니다.

# 결과표를 padj 오름차순으로 정렬해서 보기
res = ds.results_df
print(res.sort_values("padj").head())

       baseMean  log2FoldChange     lfcSE      stat        pvalue          padj
gene8   12.4103        2.8571      0.4126    6.9245    4.37e-12      4.37e-11
gene4    9.8820       -1.9438      0.3852   -5.0462    4.51e-07      2.26e-06
gene2    5.4124       -0.3780      0.3751   -1.0076    0.313653      0.522755
gene1    8.9817        0.2583      0.3120    0.8277    0.407843      0.582633
gene6    7.1055        0.1902      0.3344    0.5688    0.569470      0.632745

각 컬럼의 뜻은 이렇습니다. baseMean은 모든 샘플에서 정규화된 평균 발현량, log2FoldChange는 조건 간 발현비의 log2 값(1이면 2배·2면 4배), lfcSE는 그 표준오차, stat은 Wald 통계량, pvalue는 원시 pvalue, 그리고 padj는 다중검정보정(BH, Benjamini-Hochberg)을 거친 조정 pvalue입니다. 유전자 수천 개를 동시에 검정하면 우연히 낮은 pvalue가 쏟아지므로, 반드시 padj로 판단해야 합니다 — 왜 그런지는 p-value와 FDR 편에서 정리했습니다.

lfc_shrink — fold change 축소(선택)

발현이 낮거나 분산이 큰 유전자는 log2 fold change 추정이 요동칩니다. lfc_shrink는 이런 잡음을 줄여, 볼케이노·MA 플롯에서 신뢰할 만한 크기만 남깁니다(R DESeq2의 lfcShrink와 같은 역할).

# LFC 축소 — coeff는 설계식이 만든 계수 이름
ds.lfc_shrink(coeff="condition[T.B]")

축소는 pvalue를 바꾸지 않고 log2FoldChange 값만 0쪽으로 당깁니다. 유의성 판정은 그대로 두면서 효과 크기 시각화를 안정화하는 셈이죠. 이제 이 결과표로 볼케이노 플롯을 그립니다.

4. 볼케이노 플롯 — 결과 시각화

볼케이노 플롯(volcano plot)은 차등발현 결과를 한 장으로 요약하는 표준 그림입니다. x축에 log2FoldChange(변화 방향·크기), y축에 −log10(padj)(유의성)를 놓아, 오른쪽·왼쪽 위 구석에 몰린 점이 "크게 변하고 통계적으로도 유의한" 관심 유전자입니다. matplotlib만으로 그릴 수 있습니다.

import numpy as np
import matplotlib.pyplot as plt

res = ds.results_df.dropna(subset=["padj"])            # padj 없는 행 제거

# 유의 기준: padj < 0.05 이고 |log2FC| > 1 (2배 이상 변화)
sig = (res["padj"] < 0.05) & (res["log2FoldChange"].abs() > 1)
up   = sig & (res["log2FoldChange"] > 0)               # 상향(up)
down = sig & (res["log2FoldChange"] < 0)               # 하향(down)

x = res["log2FoldChange"]
y = -np.log10(res["padj"])                             # 작을수록 위로

plt.figure(figsize=(6, 5))
plt.scatter(x[~sig], y[~sig], s=12, c="#B8B0A2", label="비유의(ns)")   # 회색: 유의하지 않음
plt.scatter(x[up],   y[up],   s=16, c="#E76F51", label="상향(up)")     # 코랄: 상향
plt.scatter(x[down], y[down], s=16, c="#2E7D6B", label="하향(down)")   # 그린: 하향

plt.axhline(-np.log10(0.05), ls="--", c="#9AA0A6", lw=1)   # padj 0.05 기준선
plt.axvline( 1, ls="--", c="#9AA0A6", lw=1)                # log2FC +1
plt.axvline(-1, ls="--", c="#9AA0A6", lw=1)                # log2FC -1
plt.xlabel("log2 fold change (B vs A)")
plt.ylabel(r"$-\log_{10}$ adjusted $p$")
plt.legend(frameon=False)
plt.tight_layout()
plt.savefig("volcano.png", dpi=150)

# volcano.png 저장됨 — 유의 유전자가 좌우 상단에 강조됨

이 코드는 어느 결과표에나 그대로 적용됩니다. 다만 내장 합성 데이터는 유전자가 10개뿐이라 볼케이노가 성글게 나옵니다. 유전자 수천 개짜리 실제 데이터(예: airway·GEO 카운트)를 넣으면, 아래처럼 전형적인 화산 모양이 나타납니다.

# (실데이터 예시) 유의 유전자 개수 세기
n_up   = int(up.sum())
n_down = int(down.sum())
print(f"상향 {n_up}개 · 하향 {n_down}개 · 전체 검정 {len(res)}개")

상향 214개 · 하향 187개 · 전체 검정 12648개

여기서 상향 214·하향 187은 유전자 수천 개짜리 실데이터를 돌렸을 때의 대표 예시입니다. 볼케이노에서 좌우 상단으로 갈수록 "크게, 그리고 확실하게" 달라진 유전자라, 후속 검증(qPCR)이나 경로 분석(pathway analysis)의 1순위 후보가 됩니다. 점 색·기준선·라벨은 볼케이노 시각화 편과 matplotlib 시각화 편의 방식으로 더 다듬을 수 있습니다.

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

• 카운트가 정수가 아님 — 가장 흔한 실수 — DeseqDataSet에 TPM·FPKM처럼 소수·정규화된 값을 넣으면 ValueError: counts must be non-negative integers 류의 오류가 나거나, 결과가 통계적으로 무의미해집니다. DESeq2 계열은 원시 리드 카운트(raw count)를 전제로 음이항 분포를 적합하므로, featureCounts·htseq-count·salmon의 정수 카운트를 그대로 넣어야 합니다. 정규화는 .deseq2()가 내부에서 크기 인자로 처리합니다.
• 샘플 순서·이름 불일치 — 카운트의 행(샘플)과 메타데이터의 행이 서로 다른 순서거나 이름이 어긋나면, 조건이 엉뚱한 샘플에 붙어 결과가 통째로 틀립니다. counts.index.equals(meta.index)로 두 인덱스가 같은지 먼저 확인하고, 어긋나면 meta = meta.loc[counts.index]로 맞추세요. pydeseq2는 이 정렬을 자동으로 해 주지 않습니다.
• design_factors vs design 인자 이름 — 오래된 튜토리얼을 따라 DeseqDataSet(..., design_factors="condition")으로 쓰면, pydeseq2 0.5+에서는 TypeError: unexpected keyword argument 'design_factors'가 납니다. 0.5부터 R과 같은 설계식 문자열 design="~condition"으로 바뀌었기 때문입니다. 반대로 0.4.x에 design=을 쓰면 인식하지 못하니, pydeseq2.__version__을 먼저 확인하세요.
• 저발현 유전자로 인한 padj NaN — 발현이 너무 낮은 유전자는 독립 필터링(independent filtering)·Cook's 이상치 처리 과정에서 padj가 NaN으로 남습니다. 이는 오류가 아니라 "검정에서 제외됨"이라는 신호입니다. 볼케이노를 그리기 전 res.dropna(subset=["padj"])로 걸러야 −log10 계산에서 에러가 안 납니다. 사전에 counts.loc[:, counts.sum(axis=0) >= 10]처럼 총 카운트가 낮은 유전자를 미리 쳐내면 검정력도 올라갑니다.
• 결과가 비었거나 대비 오류 — contrast=["condition", "B", "A"]에서 B·A가 실제 메타데이터의 condition 값과 정확히 일치해야 합니다. 대소문자·공백만 달라도 KeyError가 나거나 빈 결과가 나옵니다. meta["condition"].unique()로 실제 수준(level) 이름을 확인한 뒤 그대로 넣으세요.

6. 확장 — 더 해 볼 것

• 실제 공개 데이터로 돌리기 — 합성 데이터 대신 GEO·recount3의 실제 카운트(예: airway 데이터셋)를 pd.read_csv로 읽어 넣으면, 수천 유전자짜리 진짜 볼케이노가 나옵니다. 전처리는 저발현 유전자 필터링(sum >= 10)만 하면 충분합니다.
• 다요인 설계·배치 보정 — design="~batch + condition"처럼 공변량을 넣으면, 배치 효과를 보정한 채 조건 차이만 뽑을 수 있습니다. 관심 변수를 설계식 맨 뒤에 두는 게 관례입니다.
• scanpy·단일세포와 잇기 — pydeseq2는 anndata 기반이라, scanpy로 다룬 단일세포(single-cell) 데이터를 유사벌크(pseudobulk)로 합친 뒤 그대로 넘겨 차등발현을 볼 수 있습니다. R을 거치지 않고 파이썬 한 흐름으로 끝나는 게 최대 장점이죠 — 단일세포 개념은 scRNA-seq 편에서 다뤘습니다.
• 결과 저장·경로 분석 — res.to_csv("deg.csv")로 결과표를 저장하고, 유의 유전자 목록을 gseapy(GSEA 파이썬판)나 Enrichr로 넘기면 경로 분석까지 파이썬 안에서 이어집니다.
• R DESeq2와 교차검증 — 같은 카운트를 R DESeq2와 pydeseq2 양쪽에 돌려 log2FoldChange·padj를 산점도로 비교하면, 두 구현이 거의 일직선으로 일치함을 눈으로 확인할 수 있습니다. 방법론이 궁금하면 R DESeq2 편과 나란히 놓고 보세요.

R DESeq2 vs pydeseq2 — 무엇을 언제 쓸까

결론부터 말하면, 결과는 거의 같으니 생태계로 고르면 됩니다. 파이프라인이 이미 파이썬(scanpy·scikit-learn) 위에 있다면 pydeseq2가 R을 오가는 수고를 없애 줍니다. 반대로 apeglm 축소나 복잡한 상호작용 설계처럼 DESeq2 생태계의 성숙한 확장이 필요하면 R이 아직 앞섭니다. pydeseq2는 원논문에서 R DESeq2와의 일치를 정량적으로 검증했으니, 파이썬 사용자에게는 "R을 안 깔아도 되는 DESeq2"로 충분히 신뢰할 만합니다.

References

1. Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. (DESeq2 원논문)
2. Muzellec, B., Teleńczuk, M., Cabeli, V., & Andreux, M. (2023). PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 39(9), btad547. (pydeseq2 원논문)
3. Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B, 57(1), 289–300. (BH 다중검정보정 원논문)
4. PyDESeq2 documentation. (2024). A simple PyDESeq2 workflow · Step-by-step workflow. pydeseq2.readthedocs.io
5. Anders, S., & Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biology, 11(10), R106. (음이항 모델·크기 인자 기초)

Pipette & Pipeline · A bio portfolio journal

TL;DR — NGS 파이프라인에서 정렬 결과(SAM/BAM)와 변이 결과(VCF)를 다루는 두 표준 도구가 samtools (SAM/BAM/CRAM 처리)와 bcftools (VCF/BCF 처리)입니다. 이 글은 samtools에 동봉된 예제(toy.sam·toy.fa)로, SAM→정렬된 BAM→인덱스(.bai)→커버리지 확인→bcftools mpileup+call로 변이 호출→QUAL·DP로 필터링까지 한 줄씩 따라갑니다. 모든 명령·플래그에 한글 주석을 달고, flagstat·coverage 출력 예시까지 그대로 보여 재현 가능하게 구성했습니다.

samtools의 핵심 흐름은 단순합니다. view로 형식을 바꾸거나 FLAG·MAPQ로 거르고, sort로 좌표순 정렬하고, index로 임의 좌표 접근을 위한 색인을 만든 뒤, flagstat·coverage·depth로 정렬 품질과 커버리지를 확인합니다. bcftools는 그 BAM에서 mpileup+call로 변이를 뽑아(GATK의 가벼운 대안), filter 표현식으로 저품질 변이를 걷어냅니다.

재현성을 위해 예제는 어디서 돌려도 결과가 같은 toy.sam을 씁니다. 초보가 가장 많이 밟는 지뢰 — 정렬 안 된 BAM을 인덱싱하다 나는 오류, bcftools mpileup의 -f 참조 누락, @SQ 헤더 불일치 — 도 4절에서 실제 에러 메시지와 함께 짚습니다.

🔗 관련 글  ·  SAM/BAM의 11개 필드·FLAG·CIGAR를 아직 안 봤다면 SAM/BAM 파일이란?을 먼저 읽는 걸 권합니다  ·  변이 결과 포맷은 VCF 파일이란?에서 8개 컬럼·INFO·GT를 정리했습니다  ·  BAM을 만들기 직전 단계인 정렬은 BWA-MEM·STAR로 리드 정렬하기에서, 더 무거운 변이 호출기는 GATK으로 변이 호출하기에서 이어집니다

이 글에서 만들 것

시퀀싱 파이프라인은 대체로 이렇게 흐릅니다. FASTQ (원시 리드)를 품질검사·트리밍하고, 참조 유전체에 정렬해 SAM/BAM을 얻고, 그 정렬을 정리한 뒤 변이를 호출해 VCF를 만듭니다. 이 가운데 정렬 결과를 다루는 도구가 samtools, 변이 결과를 다루는 도구가 bcftools입니다. 정렬 자체(FASTQ→SAM)는 BWA·STAR 편에서 다뤘으니, 이 글은 그다음 — SAM을 정렬된 BAM으로 만들고, 커버리지를 확인하고, 변이를 뽑아 필터링하는 구간 — 을 처음부터 끝까지 만듭니다.

samtools (Heng Li 외, 2009~)는 SAM (Sequence Alignment/Map)·BAM (SAM의 이진 압축본)·CRAM을 읽고 쓰고 거르는 표준 도구입니다. bcftools는 그 짝으로, 변이 정보 포맷인 VCF (Variant Call Format)·BCF (VCF의 이진본)를 다룹니다. 둘 다 htslib이라는 같은 C 라이브러리 위에 얹혀 있어, 2021년 GigaScience 논문 "Twelve years of SAMtools and BCFtools"로 함께 정리됐습니다.

구체적으로 이 글은 samtools 예제 파일 하나로 ① view로 SAM을 BAM으로 바꾸고 FLAG·MAPQ로 거르고, ② sort로 좌표순 정렬한 뒤 index로 .bai 색인을 만들고, ③ flagstat·coverage·depth로 정렬 통계와 커버리지를 보고, ④ faidx로 참조를 색인한 뒤 bcftools mpileup+call로 변이를 호출하고, ⑤ bcftools filter·norm·query로 변이를 걸러 원하는 컬럼만 뽑습니다. 명령은 samtools·bcftools 1.21+ 기준입니다.

# 핵심 흐름 한눈에 (toy 예제 기준)
samtools view -b toy.sam | samtools sort -o toy.sorted.bam   # ① SAM → 정렬된 BAM
samtools index toy.sorted.bam                                # ② .bai 색인 생성
samtools flagstat toy.sorted.bam                             # ③ 정렬 통계
samtools faidx toy.fa                                         # ④ 참조 .fai 색인
bcftools mpileup -f toy.fa toy.sorted.bam \
  | bcftools call -mv -Oz -o calls.vcf.gz                    # ⑤ 변이 호출 (VCF)
bcftools filter -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz      # ⑥ 저품질 변이 필터

1. 사전 준비 — 설치와 예제 데이터

samtools와 bcftools를 설치합니다. 대부분 환경에서 conda (bioconda 채널)가 가장 편하고, 데비안·우분투 계열이면 apt로도 됩니다.

# conda (bioconda) — 권장, 버전 고정이 쉬움
conda install -c bioconda -c conda-forge samtools bcftools

# 또는 데비안/우분투
sudo apt-get update && sudo apt-get install samtools bcftools

# 설치 확인 (버전 출력)
samtools --version | head -1
bcftools --version | head -1

samtools 1.21
bcftools 1.21

이제 실습 데이터를 받습니다. 재현이 쉽도록, 어디서 돌려도 결과가 같은 samtools 공식 예제 toy.sam·toy.fa를 씁니다. 참조 서열 2개(ref 45 bp·ref2 40 bp)와 정렬 레코드 12개(ref에 r001~r004 6개·ref2에 x1~x6 6개)로 이루어진 아주 작은 데이터입니다.

# samtools 저장소에서 예제 파일 내려받기
wget https://raw.githubusercontent.com/samtools/samtools/develop/examples/toy.sam
wget https://raw.githubusercontent.com/samtools/samtools/develop/examples/toy.fa

# 내용 확인 — 헤더(@SQ 2줄) + 정렬 레코드 12개
cat toy.sam

@SQ	SN:ref	LN:45
@SQ	SN:ref2	LN:40
r001	163	ref	7	30	8M4I4M1D3M	=	37	39	TTAGATAAAGAGGATACTG	*	XX:B:S,12561,2,20,112
r002	0	ref	9	30	1S2I6M1P1I1P1I4M2I	*	0	0	AAAAGATAAGGGATAAA	*
r003	0	ref	9	30	5H6M	*	0	0	AGCTAA	*
r004	0	ref	16	30	6M14N1I5M	*	0	0	ATAGCTCTCAGC	*
r003	16	ref	29	30	6H5M	*	0	0	TAGGC	*
r001	83	ref	37	30	9M	=	7	-39	CAGCGCCAT	*
x1	0	ref2	1	30	20M	*	0	0	aggttttataaaacaaataa	????????????????????
x2	0	ref2	2	30	21M	*	0	0	ggttttataaaacaaataatt	?????????????????????
x3	0	ref2	6	30	9M4I13M	*	0	0	ttataaaacAAATaattaagtctaca	??????????????????????????
x4	0	ref2	10	30	25M	*	0	0	CaaaTaattaagtctacagagcaac	?????????????????????????
x5	0	ref2	12	30	24M	*	0	0	aaTaattaagtctacagagcaact	????????????????????????
x6	0	ref2	14	30	23M	*	0	0	Taattaagtctacagagcaacta	???????????????????????

각 줄은 리드 하나의 정렬입니다. 왼쪽부터 리드 이름(QNAME)·FLAG·참조 이름(RNAME)·위치(POS)·MAPQ·CIGAR 순이고, 이 11개 필드의 의미는 SAM/BAM 파일 편에서 자세히 풀었습니다. 여기서 눈여겨볼 건 두 번째 열 FLAG입니다. FLAG는 여러 상태를 비트로 겹쳐 담은 정수인데, 이걸 이해하면 samtools view의 필터가 한결 명확해집니다.

2. samtools view — 형식 변환·FLAG·MAPQ 필터

SAM FLAG — 상태를 비트로 겹쳐 담은 정수

FLAG는 "이 리드가 짝이 있나·정방향인가·중복인가" 같은 여러 예/아니오를 2의 거듭제곱 비트로 합친 값입니다. 예컨대 FLAG가 163이면 1 (짝 있음) + 2 (제대로 짝지음) + 32 (짝이 역방향) + 128 (두 번째 리드)의 합이죠. 자주 쓰는 비트는 이렇습니다.

숫자를 외울 필요는 없습니다. samtools flags 서브명령이 비트↔이름을 바꿔 주거든요.

# FLAG 값 → 어떤 비트가 켜졌는지 이름으로 풀기
samtools flags 0x904

0x904	2308	UNMAP,SECONDARY,SUPPLEMENTARY

0x904(=2308)는 unmapped·secondary·supplementary 세 비트를 합친 값입니다. 이걸 -F(제외)에 넣으면 "정렬 실패·이차·보충 정렬을 빼고 진짜 일차 정렬만" 남길 수 있죠. 실전에서 자주 쓰는 관용구입니다.

view — 형식 변환과 필터

samtools view는 SAM↔BAM 형식 변환과 필터를 동시에 하는 만능 명령입니다. 자주 쓰는 옵션은 -b(BAM 출력), -h(헤더 포함), -q(MAPQ 문턱), -f(비트 요구), -F(비트 제외), -c(개수만)입니다.

# SAM → BAM 변환 (-b), 헤더 포함(-h는 BAM에선 기본 포함)
samtools view -b toy.sam -o toy.bam

# MAPQ ≥ 20 이고, unmapped·secondary·supplementary(0x904) 제외한 리드 개수만 세기
samtools view -c -q 20 -F 0x904 toy.bam

12

toy 데이터의 정렬 12개는 모두 MAPQ 30이고 이차·보충 정렬이 없어, 이 조건에 12개가 다 걸립니다. 문턱을 -q 60으로 올리면 0개가 되죠. -c는 레코드를 출력하지 않고 개수만 돌려줘, "필터 조건에 몇 개가 걸리나"를 빠르게 볼 때 씁니다. -q 20은 MAPQ (정렬 신뢰도, 값이 클수록 그 위치가 유일할 가능성이 높음)가 20 미만인 리드를 버립니다 — MAPQ 개념은 read alignment 편에서 다뤘습니다.

# 특정 좌표 구간만 잘라내기 (색인이 있어야 동작 — 3절에서 index 생성 후)
samtools view -b toy.sorted.bam ref:1-20 -o region.bam

view에 참조이름:시작-끝 형식을 주면 그 구간에 걸친 리드만 뽑습니다. 다만 이 임의 좌표 접근은 BAM이 좌표순으로 정렬되고 색인(.bai)이 있을 때만 동작합니다. 그래서 다음 순서가 sort → index입니다.

3. samtools sort·index·flagstat·coverage — 정렬·색인·통계

sort — 좌표순 정렬

정렬기(BWA 등)가 내놓은 BAM은 리드가 읽힌 순서 그대로라, 참조 좌표순이 아닙니다. 대부분의 하류 도구(변이 호출·색인·구간 추출)는 좌표순 정렬을 전제하므로, samtools sort로 먼저 정렬합니다.

# 좌표순 정렬 (기본값) → toy.sorted.bam
samtools sort toy.bam -o toy.sorted.bam

sort는 기본적으로 참조 위치(leftmost coordinate)순으로 정렬합니다. 리드 이름순이 필요하면(예: fixmate·markdup 전처리) -n을 줍니다. 큰 파일이면 -@로 스레드를, -m으로 스레드당 메모리를 늘려 속도를 올릴 수 있습니다.

index — 임의 좌표 접근을 위한 색인

samtools index는 정렬된 BAM에 .bai 색인을 붙여, 파일 전체를 읽지 않고도 특정 좌표로 바로 점프할 수 있게 합니다. IGV 같은 뷰어나 구간 추출이 이 색인에 기댑니다.

# .bai 색인 생성 → toy.sorted.bam.bai
samtools index toy.sorted.bam

여기서 흔한 함정 하나. 정렬 안 된 BAM을 인덱싱하려 하면 오류가 납니다.

# (오류 예시) 정렬 안 된 toy.bam을 그대로 인덱싱하면
samtools index toy.bam

samtools index: "toy.bam" is not sorted or is malformed. Coordinate-sorted files are required for indexing.

색인은 좌표순 정렬을 반드시 전제하므로, index 전에 sort가 와야 합니다. 이 순서는 4절에서 다시 정리합니다.

flagstat — 정렬 통계 한눈에

samtools flagstat은 FLAG를 집계해 "전체 몇 개 중 몇 개가 정렬됐고, 짝지음 비율은 얼마인가"를 한 화면에 보여 줍니다. 실제 대규모 BAM에서는 이런 모양입니다(paired-end, 약 100만 리드 예시).

samtools flagstat sample.sorted.bam

1005678 + 0 in total (QC-passed reads + QC-failed reads)
1000000 + 0 primary
3210 + 0 secondary
2468 + 0 supplementary
18234 + 0 duplicates
18234 + 0 primary duplicates
998765 + 0 mapped (99.31% : N/A)
993087 + 0 primary mapped (99.31% : N/A)
1000000 + 0 paired in sequencing
500000 + 0 read1
500000 + 0 read2
985432 + 0 properly paired (98.54% : N/A)
991200 + 0 with itself and mate mapped
1887 + 0 singletons (0.19% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

각 줄은 QC통과 + QC실패 형식입니다. 가장 먼저 보는 값은 mapped 비율(정렬률, 여기선 99.31%)과 properly paired 비율(제대로 짝지어 정렬된 비율, 98.54%)입니다. 정렬률이 유난히 낮거나 properly paired가 낮으면 참조 불일치·오염·라이브러리 문제를 의심합니다. duplicates가 많으면 PCR 중복 표시(markdup)를 점검하죠.

coverage·depth — 얼마나 깊게 읽혔나

변이를 믿으려면 그 자리가 충분히 깊게 읽혔는지(커버리지)를 봐야 합니다. samtools coverage는 참조 서열마다 요약을, samtools depth는 위치마다 깊이를 줍니다. toy 데이터는 두 참조 모두에 리드가 걸려 있어(ref r001~r004·ref2 x1~x6) 둘 다 요약이 나옵니다.

# 참조 서열별 커버리지 요약 (한 줄에 서열 하나)
samtools coverage toy.sorted.bam

#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
ref	1	45	6	39	86.67	2.0	40	30
ref2	1	40	6	36	90	3.32	40	30

coverage 컬럼은 그 서열에서 1× 이상 덮인 염기 비율(%), meandepth는 평균 깊이입니다. 실제 데이터라면 이 표에서 커버리지가 낮은 구간을 먼저 살펴, 변이를 믿어도 되는 자리인지 판단합니다. -m을 주면 같은 정보를 ASCII 히스토그램으로 그려 줍니다.

# 위치별 깊이 (CHROM  POS  DEPTH), -a는 깊이 0 위치도 포함
samtools depth toy.sorted.bam | head -5

ref	7	1
ref	8	1
ref	9	2
ref	10	2
ref	11	2

depth는 위치마다 세 컬럼(참조 이름·1-based 위치·깊이)을 내놓습니다. 특정 유전자·엑손의 깊이 프로파일을 그리거나, 최소 깊이 미달 구간을 찾을 때 씁니다.

4. bcftools — mpileup·call로 변이 호출하고 filter로 거르기

이제 정렬된 BAM에서 변이를 뽑습니다. bcftools는 두 단계로 나뉩니다. mpileup이 각 위치의 염기 증거를 모아 유전형 가능도(genotype likelihood)를 계산하고, call이 그걸 받아 실제 변이를 판정합니다.

faidx — 참조 색인 (변이 호출의 전제)

bcftools mpileup은 참조 FASTA를 -f로 반드시 받아야 합니다. 참조가 있어야 어느 자리가 참조와 다른지 알 수 있으니까요. 그리고 그 참조는 samtools faidx로 미리 색인(.fai)해 둬야 임의 좌표를 빠르게 읽습니다.

# 참조 FASTA 색인 → toy.fa.fai 생성
samtools faidx toy.fa
cat toy.fa.fai

ref	45	5	45	46
ref2	40	57	40	41

.fai는 다섯 컬럼(이름·길이·바이트 오프셋·줄당 염기수·줄당 바이트수)입니다. toy.fa는 서열이 한 줄이라 줄당 염기수가 길이와 같죠. 이 색인 덕에 samtools faidx toy.fa ref:1-10처럼 참조의 특정 구간만 즉시 뽑을 수도 있습니다.

mpileup + call — 변이 호출

# mpileup: 위치별 염기 증거 → 유전형 가능도(VCF/BCF)
#   -f: 참조 FASTA (필수)
# call: 가능도 → 변이 판정
#   -m: 다대립 호출기(권장 기본) · -v: 변이 자리만 출력 · -Oz: gzip VCF 출력
bcftools mpileup -f toy.fa toy.sorted.bam \
  | bcftools call -mv -Oz -o calls.vcf.gz

toy 데이터는 리드가 적어 변이가 거의 안 나오지만, 명령의 흐름과 출력 구조는 실제와 같습니다. 아래는 실제 소형 BAM (예: 1000 Genomes 저커버리지 BAM의 한 구간)에서 얻은 것과 같은 형태의 VCF입니다(대표 예시).

# 결과 VCF 훑어보기 (헤더 제외 데이터 줄 앞부분)
bcftools view calls.vcf.gz | grep -v '^##' | head -4

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
20	1001420	.	G	A	221.4	.	DP=34;MQ=59	GT:PL	1/1:255,102,0
20	1004201	.	C	T	18.7	.	DP=6;MQ=42	GT:PL	0/1:52,0,71
20	1008833	.	T	TA	164.2	.	INDEL;DP=27	GT:PL	0/1:198,0,201

VCF의 8개 고정 컬럼(CHROM·POS·ID·REF·ALT·QUAL·FILTER·INFO)에 이어, FORMAT과 샘플별 값이 붙습니다. QUAL은 변이가 진짜일 확률의 Phred 점수(클수록 신뢰), INFO의 DP는 그 자리 총 깊이, FORMAT의 GT는 유전형(0/1 이형접합·1/1 동형접합 변이)입니다. 이 컬럼들의 의미는 VCF 파일 편에서 정리했습니다. 위에서 두 번째 변이(POS 1004201)는 QUAL 18.7·DP 6으로, 다음 단계에서 걸러질 후보입니다.

filter — 표현식으로 저품질 변이 제거

bcftools filter는 -e(제외)·-i(포함) 표현식으로 변이를 거릅니다. 표현식에는 QUAL·INFO 필드를 이름으로 참조할 수 있습니다.

# QUAL이 20 미만 '또는' 깊이(DP)가 10 미만인 변이를 제외
bcftools filter -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz -Oz -o filtered.vcf.gz

# 걸러진 뒤 남은 변이 확인
bcftools view filtered.vcf.gz | grep -v '^##'

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
20	1001420	.	G	A	221.4	.	DP=34;MQ=59	GT:PL	1/1:255,102,0
20	1008833	.	T	TA	164.2	.	INDEL;DP=27	GT:PL	0/1:198,0,201

QUAL<20 || INFO/DP<10 표현식이 QUAL 18.7·DP 6이던 POS 1004201을 정확히 걸러, 신뢰할 만한 변이 둘만 남았습니다. 실전에서는 변이를 삭제하는 대신, -s LowQual로 FILTER 컬럼에 이름표만 붙여(soft filter) 나중에 되살릴 수 있게 하는 방식도 많이 씁니다.

# (대안) 삭제 대신 FILTER 컬럼에 'LowQual' 표시만 (soft filter)
bcftools filter -s LowQual -e 'QUAL<20 || INFO/DP<10' calls.vcf.gz -Oz -o soft.vcf.gz

norm·stats·query — 정규화·요약·추출

변이를 비교·병합하려면 표현을 통일해야 합니다. bcftools norm은 인델을 왼쪽 정렬(left-align)하고, 한 줄에 뭉친 다대립 변이를 -m -로 이대립씩 쪼갭니다.

# 인델 정규화(left-align) + 다대립 → 이대립 분해 (-m -)
bcftools norm -f toy.fa -m - filtered.vcf.gz -Oz -o norm.vcf.gz

# 변이 요약 통계 (SNP·인델 개수 등)
bcftools stats norm.vcf.gz | grep '^SN'

SN	0	number of samples:	1
SN	0	number of records:	2
SN	0	number of SNPs:	1
SN	0	number of indels:	1
SN	0	number of multiallelic sites:	0

bcftools query는 원하는 컬럼만 뽑아 표로 만들 때 씁니다. INFO·고정 컬럼은 %필드로, 샘플별 FORMAT 값은 대괄호 [...]로 감쌉니다.

# 원하는 컬럼만 탭 구분으로 추출 ([%GT]는 샘플별 유전형)
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT]\n' norm.vcf.gz

20	1001420	G	A	221.4	1/1
20	1008833	T	TA	164.2	0/1

\t는 탭, \n은 줄바꿈입니다. 이렇게 뽑은 표는 그대로 R·pandas로 넘겨 하류 분석에 쓰기 좋습니다. query는 VCF를 사람이·스크립트가 다루기 쉬운 평범한 표로 바꿔 주는, 파이프라인 마지막의 단골 도구입니다.

자주 발생하는 에러 & 해결

• 정렬 안 된 BAM 인덱싱 오류 — 가장 흔한 함정 — samtools index가 "... is not sorted or is malformed. Coordinate-sorted files are required for indexing."를 뱉으면, BAM이 좌표순이 아니라는 뜻입니다. 정렬기 출력이나 이름순 정렬(-n) 결과를 그대로 인덱싱하면 이 오류가 납니다. 반드시 samtools sort(좌표순)를 먼저 돌린 뒤 index합니다. 순서는 view → sort → index로 기억하세요.
• bcftools mpileup의 -f 참조 누락 — mpileup에 -f ref.fa를 빼면 변이를 계산할 기준이 없어 실패합니다. 참조 FASTA는 필수이며, 그 참조는 samtools faidx ref.fa로 미리 색인(.fai)돼 있어야 합니다(없으면 mpileup이 자동 생성하려 하지만, 쓰기 권한이 없는 경로면 실패). 정렬에 쓴 참조와 같은 FASTA를 써야 좌표가 맞습니다.
• @SQ 헤더 불일치 — 참조가 다를 때 — BAM 헤더의 @SQ(참조 이름·길이)와 변이 호출에 준 참조 FASTA가 다르면, 좌표가 어긋나거나 the sequence ... not found 류의 오류가 납니다. 정렬·변이 호출·주석 달기까지 같은 참조 유전체 버전(예: GRCh38)을 일관되게 써야 합니다. chr1 vs 1처럼 이름 규칙만 달라도 문제가 되니, 헤더를 samtools view -H로 먼저 확인하세요.
• 인덱스가 낡음(BAM보다 오래됨) — BAM을 다시 만들었는데 .bai가 옛것이면 뷰어·구간 추출이 어긋납니다. BAM을 갱신할 때마다 samtools index를 다시 돌려 색인을 최신으로 유지합니다. 마찬가지로 .fa를 바꾸면 .fai도 다시 만들어야 합니다.
• 큰 염색체 색인 — .csi가 필요할 때 — 참조 서열 하나가 512 Mbp를 넘으면 기본 .bai 색인이 못 담습니다. 이때는 samtools index -c로 .csi 색인을 만듭니다. 사람 유전체는 대개 .bai로 충분하지만, 일부 식물·양서류 유전체에서 마주칩니다.
• 원격 BAM 구간 추출이 느리거나 실패 — samtools는 http·https·ftp·s3 URL을 직접 읽어 samtools view -b URL chr:start-end로 원격 BAM의 한 구간만 잘라올 수 있습니다. 다만 그 URL 옆에 .bai 색인이 함께 있어야 하고, 공개 FTP는 느리거나 경로가 바뀌기도 합니다. 재현이 중요하면 이 글처럼 작은 로컬 예제로 먼저 검증한 뒤 실데이터로 옮기는 편이 안전합니다.

확장 — 더 해 볼 것

markdup으로 PCR 중복 표시 — 라이브러리 증폭 과정에서 생긴 중복 리드는 변이 근거를 부풀립니다. samtools markdup으로 표시(또는 제거)하는데, 전처리 순서가 정해져 있습니다.

samtools sort -n -o namesort.bam input.bam        # 1) 이름순 정렬
samtools fixmate -m namesort.bam fixmate.bam      # 2) 짝 정보 태그 부착(-m 필수)
samtools sort -o positionsort.bam fixmate.bam     # 3) 다시 좌표순 정렬
samtools markdup positionsort.bam markdup.bam     # 4) 중복 표시 (-r이면 제거)

• CRAM으로 저장 공간 줄이기 — BAM보다 더 압축되는 CRAM은 참조를 알아야 복원됩니다. samtools view -T ref.fa -C in.bam -o out.cram으로 변환하면 대용량 저장에서 용량을 크게 아낄 수 있습니다. 참조만 잘 관리하면 무손실입니다.
• 여러 샘플 합동 호출(joint calling) — bcftools mpileup -f ref.fa a.bam b.bam c.bam | bcftools call -mv처럼 BAM 여러 개를 한 번에 넘기면, 샘플들을 함께 보고 변이를 호출합니다. 코호트 분석의 기본이죠. 더 정교한 GVCF 기반 합동 호출은 GATK 편의 방식과 비교됩니다.
• VCF 병합·교집합 — bcftools merge로 여러 샘플 VCF를 하나로 합치고, bcftools isec으로 두 콜셋의 교집합·차집합을 구할 수 있습니다. 서로 다른 호출기(bcftools vs GATK vs DeepVariant) 결과를 대조할 때 유용합니다. 변이 호출기별 차이는 변이 호출 개념 편에서 다뤘습니다.
• 파이프라인으로 묶기 — 위 명령들을 셸 스크립트나 Snakemake·Nextflow 규칙으로 엮으면, FASTQ부터 필터된 VCF까지 재현 가능한 한 흐름으로 자동화할 수 있습니다. 각 단계 출력에 로그·체크섬을 남기면 디버깅이 쉬워집니다.

References

1. Li, H., Handsaker, B., Wysoker, A., et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079. (samtools 원논문)
2. Danecek, P., Bonfield, J. K., Liddle, J., et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008. (samtools·bcftools 통합 논문)
3. Danecek, P., Auton, A., Abecasis, G., et al. (2011). The variant call format and VCFtools. Bioinformatics, 27(15), 2156–2158. (VCF 포맷 원논문)
4. HTSlib project. (2024). samtools manual — view · sort · index · flagstat · coverage · faidx · markdup. htslib.org/doc/samtools
5. HTSlib project. (2024). bcftools manual — mpileup · call · filter · norm · stats · query. htslib.org/doc/bcftools
6. SAMtools/hts-specs. (2024). Sequence Alignment/Map (SAM) · Variant Call Format (VCF) 명세. hts-specs

Pipette & Pipeline · A bio portfolio journal

TL;DR — 발현 행렬 같은 표 형식 수치 데이터를 numpy (파이썬 수치 연산의 표준 라이브러리)로 빠르게 다루고, scipy.stats (과학 계산용 통계 모듈)로 두 군을 비교하는 가설검정(hypothesis test)을 한 번에 끝냅니다. 유전자 8개 × 샘플 12개(대조군 6 · 처리군 6)짜리 합성 발현 데이터를 시드를 고정해 만들고, 배열 연산·브로드캐스팅(broadcasting)·z-점수 정규화부터 독립표본 t검정·만-휘트니 U 검정(Mann-Whitney U test)·상관분석, 그리고 다중검정 보정(multiple testing correction)까지 코드로 따라갑니다.

numpy의 핵심은 반복문 없이 배열 전체를 한 번에 계산하는 벡터화(vectorization)와, 모양이 다른 배열을 자동으로 맞춰 주는 브로드캐스팅입니다. 이 둘만 손에 익으면 발현 행렬의 정규화·요약 통계가 한두 줄로 끝나죠. scipy.stats는 그렇게 정리한 데이터에 검정을 붙여, "이 유전자가 처리군에서 정말 변했는가"를 p값으로 답합니다.

이 글의 코드는 무작위 시드를 고정해 그대로 복사하면 같은 숫자가 나옵니다. 생물정보 초보가 가장 많이 밟는 지뢰 — 축(axis) 혼동·정수 나눗셈·NaN 전파, 그리고 무엇보다 다중검정 미보정도 코드로 짚습니다. 유전자 수천 개를 한꺼번에 검정하면 우연히 p < 0.05가 쏟아지기 때문에, 원시 p값을 그대로 믿으면 안 됩니다.

🔗 관련 글  ·  검정에 넣을 발현 행렬을 먼저 표로 다듬으려면 pandas로 생물 데이터 정제하기에서 이어집니다  ·  같은 t검정·상관·비모수 검정을 R로 하는 법은 R로 기초 통계 검정하기와 1:1로 비교됩니다  ·  왜 RNA-seq 결과에 adjusted p값을 쓰는지 개념은 p-value와 FDR, 다중검정보정이란?에서 풀었습니다  ·  여기서 만든 배열을 모델에 넣어 분류까지 가려면 scikit-learn으로 생물 데이터 머신러닝으로 이어집니다

이 글에서 만들 것

발현 데이터를 받으면 가장 먼저 하는 일이 둘입니다. 하나는 숫자를 정리하는 것 — 샘플마다 스케일을 맞추고, 유전자별 평균·분산을 구하고, 정규화합니다. 다른 하나는 묻는 것 — "이 유전자가 처리군과 대조군에서 다른가". 앞은 numpy의 일, 뒤는 scipy.stats의 일입니다.

numpy (Travis Oliphant, 2006~)는 파이썬에서 다차원 배열(ndarray)을 다루는 표준 라이브러리입니다. pandas·scipy·scikit-learn이 전부 그 위에 얹혀 있어, 사실상 파이썬 과학 계산의 바닥을 깔죠. scipy는 그 배열을 받아 최적화·신호처리·통계 같은 과학 연산을 더하는데, 그중 scipy.stats가 검정과 분포를 담당합니다.

이 글은 합성 발현 행렬 하나로 수치 처리부터 통계 검정까지 처음부터 끝까지 만듭니다. ① numpy로 배열을 만들고 인덱싱·슬라이싱하고, ② 브로드캐스팅과 축(axis) 연산으로 유전자별 요약과 z-점수 정규화를 하고, ③ scipy.stats로 독립표본 t검정과 비모수 만-휘트니 U 검정을 돌리고, ④ 피어슨·스피어만 상관을 보고, ⑤ 마지막으로 다중검정 보정(Benjamini-Hochberg)으로 원시 p값을 adjusted q값으로 바꿉니다. 코드는 numpy 2.5·scipy 1.18 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
import numpy as np
from scipy import stats

rng = np.random.default_rng(42)                 # ① 시드 고정
expr = rng.normal(50, 6, size=(8, 12))          # 유전자 8 × 샘플 12 발현 행렬
z = (expr - expr.mean(1, keepdims=True)) / expr.std(1, keepdims=True)  # ② z-점수(축 연산)
ctrl, treat = expr[:, :6], expr[:, 6:]          # 대조군 6 · 처리군 6
p = np.array([stats.ttest_ind(c, t).pvalue      # ③ 유전자별 t검정
              for c, t in zip(ctrl, treat)])
q = stats.false_discovery_control(p)            # ④ 다중검정 보정(BH → q값)

1. 사전 준비 — 설치와 합성 발현 데이터

numpy와 scipy 두 개만 설치하면 됩니다. scipy는 numpy를 의존성으로 끌어오니 함께 깔립니다.

# pip로 설치 (가상환경 권장)
pip install numpy scipy

# 또는 conda
conda install -c conda-forge numpy scipy

import numpy as np
import scipy
print("numpy", np.__version__)
print("scipy", scipy.__version__)

numpy 2.5.0
scipy 1.18.0

이제 분석할 데이터를 만듭니다. 공개 RNA-seq를 받아 정제하는 과정은 pandas 편에서 다뤘으니, 여기서는 재현이 쉽도록 합성 발현 행렬을 직접 생성합니다. 유전자 8개를 행, 샘플 12개를 열로 두고, 앞 6개는 대조군(control)·뒤 6개는 처리군(treatment)으로 나눕니다. 일부 유전자에만 처리군에서 발현이 오르내리게 효과(effect)를 심어, 검정이 그 차이를 잡아내는지 보겠습니다.

핵심은 np.random.default_rng(42)로 난수 생성기에 시드를 고정하는 것입니다. 이렇게 하면 누가 언제 돌려도 똑같은 숫자가 나와, 결과를 그대로 검증할 수 있습니다.

import numpy as np

rng = np.random.default_rng(42)                  # 난수 생성기에 시드 고정 (재현성)
genes = ["GeneA", "GeneB", "GeneC", "GeneD",
         "GeneE", "GeneF", "GeneG", "GeneH"]      # 유전자 8개
n_per = 6                                          # 군당 샘플 6개

base   = np.full(8, 50.0)                          # 기저 발현 (모든 유전자 50)
effect = np.array([0, 0, 8, 12, 0, -10, 0, 6])    # 처리군에서의 변화량 (일부만)

# 정규분포로 발현값 생성 → 대조군 / 처리군
ctrl  = rng.normal(loc=base[:, None],            scale=6, size=(8, n_per))
treat = rng.normal(loc=(base + effect)[:, None], scale=6, size=(8, n_per))
expr  = np.round(np.hstack([ctrl, treat]), 1)     # 좌우로 붙여 8 × 12 행렬

print("shape:", expr.shape)                        # (유전자 수, 샘플 수)
print(expr[:3, :4])                                # 앞 3 유전자, 앞 4 샘플

shape: (8, 12)
[[51.8 43.8 54.5 55.6]
 [50.8 48.1 49.9 44.9]
 [50.4 56.8 52.8 44.8]]

base[:, None]은 모양이 (8,)인 1차원 배열을 (8, 1)로 바꿔, 유전자별 평균을 열 방향으로 펼치는 트릭입니다(브로드캐스팅은 2절에서 자세히). 이렇게 만든 expr는 유전자 8개 × 샘플 12개짜리 2차원 배열로, 실제 발현 행렬과 같은 모양입니다.

2. numpy — 배열 연산·브로드캐스팅·축 연산·정규화

numpy의 진짜 힘은 반복문 없이 배열 전체를 한 번에 계산하는 데 있습니다. 먼저 인덱싱과 슬라이싱부터 봅니다. 발현 행렬에서 특정 유전자·샘플을 골라내는 일은 거의 매번 합니다.

print(expr[3])           # GeneD(4번째 행) 전체 12 샘플
print(expr[:, 0])        # 첫 샘플의 8 유전자 (열 슬라이싱)
print(expr[3, 6:])       # GeneD의 처리군 6 샘플만 (행+열 슬라이싱)

[55.3 49.7 48.9 45.9 57.3 49.1 59.9 59.2 67.1 60.9 54.3 55.2]
[51.8 50.8 50.4 47.3 51.7 38.5 49.6 54.9 52.6 50.7 47.6 50.6]
[59.9 59.2 67.1 60.9 54.3 55.2]

expr[행, 열] 형태로 두 축을 한꺼번에 자릅니다. :는 "그 축 전부", 6:는 "6번째부터 끝까지"입니다. 파이썬 리스트와 달리 numpy 슬라이싱은 데이터를 복사하지 않고 같은 메모리를 가리키는 뷰(view)라, 큰 행렬에서도 메모리·속도 이점이 큽니다.

브로드캐스팅 — 모양이 다른 배열을 자동으로 맞춤

브로드캐스팅(broadcasting)은 모양이 다른 두 배열을 연산할 때, 작은 쪽을 자동으로 늘려 맞춰 주는 규칙입니다. 예컨대 (3, 1) 배열과 (4,) 배열을 더하면, numpy가 양쪽을 (3, 4)로 펼쳐 모든 조합을 계산합니다.

col = np.array([1, 2, 3])      # 모양 (3,)
row = np.array([10, 20, 30, 40])  # 모양 (4,)

print(col[:, None] + row)       # (3,1) + (4,) → (3,4)로 자동 확장

[[11 21 31 41]
 [12 22 32 42]
 [13 23 33 43]]

생물 데이터에서 이게 왜 중요하냐면, "유전자마다 다른 값을 행렬 전체에 빼거나 나누는" 작업이 브로드캐스팅 한 줄로 끝나기 때문입니다. 유전자별 평균을 빼는 정규화가 대표적인데, 바로 다음에 씁니다.

축(axis) 연산 — 행이냐 열이냐

2차원 배열에서 평균·표준편차를 구할 때, 어느 방향으로 묶을지를 axis로 정합니다. 발현 행렬에서는 이 구분이 의미를 가릅니다. axis=0은 행을 가로질러(즉 샘플마다) 묶고, axis=1은 열을 가로질러(즉 유전자마다) 묶습니다.

print("샘플별 평균 axis=0:", expr.mean(axis=0).round(2)[:4])  # 각 샘플의 8 유전자 평균
print("유전자별 평균 axis=1:", expr.mean(axis=1).round(2))     # 각 유전자의 12 샘플 평균
print("유전자별 표준편차 axis=1:", expr.std(axis=1).round(2))

샘플별 평균 axis=0: [51.72 48.79 50.25 49.21]
유전자별 평균 axis=1: [48.8  50.42 54.22 55.23 50.56 45.75 49.98 51.96]
유전자별 표준편차 axis=1: [5.68 4.52 6.29 5.85 2.89 8.04 3.87 4.07]

여기서 외우면 편한 규칙 하나 — axis는 "없어지는 축"입니다. axis=0을 주면 0번 축(행, 길이 8)이 사라지고 길이 12짜리 샘플별 결과가, axis=1을 주면 1번 축(열, 길이 12)이 사라지고 길이 8짜리 유전자별 결과가 남습니다. 효과를 심은 GeneD (axis=1 평균 55.23)와 발현을 낮춘 GeneF (45.75)가 다른 유전자보다 평균이 벗어나 있는 게 보이죠.

z-점수 정규화 — 유전자별로 평균 0·표준편차 1

서로 다른 유전자를 한 화면에서 비교하려면 스케일을 맞춰야 합니다. z-점수 정규화(z-score normalization)는 각 유전자를 평균 0·표준편차 1로 바꿔, 발현의 절대 크기 대신 "그 유전자 기준으로 얼마나 높/낮은가"를 봅니다. 히트맵을 그릴 때 거의 항상 거치는 단계죠.

mu = expr.mean(axis=1, keepdims=True)   # 유전자별 평균 → 모양 (8, 1) 유지
sd = expr.std(axis=1, keepdims=True)    # 유전자별 표준편차 → (8, 1)
z = (expr - mu) / sd                     # 브로드캐스팅: (8,12) - (8,1) → (8,12)

print(z[:2, :4].round(3))                # 정규화 결과 앞부분
print("유전자별 평균(≈0):", z.mean(axis=1).round(3))
print("유전자별 표준편차(≈1):", z.std(axis=1).round(3))

[[ 0.528 -0.88   1.003  1.197]
 [ 0.085 -0.512 -0.114 -1.22 ]]
유전자별 평균(≈0): [-0.  0.  0.  0.  0. -0. -0.  0.]
유전자별 표준편차(≈1): [1. 1. 1. 1. 1. 1. 1. 1.]

핵심은 keepdims=True입니다. 이걸 빼면 mu의 모양이 (8,)이 되어 (8, 12) - (8,) 연산에서 브로드캐스팅 규칙이 어긋나 에러가 납니다. keepdims=True로 (8, 1)을 유지해야 "유전자마다(행마다) 그 평균을 빼라"가 제대로 동작하죠. 이 한 줄이 반복문으로 짜면 여러 줄이 되는 작업을, 벡터화로 깔끔하게 처리합니다.

3. scipy.stats — 독립표본 t검정·만-휘트니 U 검정

배열이 정리됐으니 이제 질문에 답합니다. "이 유전자가 처리군에서 대조군과 정말 다른가". 두 독립 군의 평균을 비교하는 가장 표준적인 도구가 독립표본 t검정(independent two-sample t-test)입니다. scipy.stats.ttest_ind가 검정통계량(t)과 p값을 한 번에 돌려줍니다.

먼저 효과를 크게 심은 GeneD 하나만 자세히 봅니다.

from scipy import stats

gd_ctrl  = expr[3, :6]    # GeneD 대조군 6 샘플
gd_treat = expr[3, 6:]    # GeneD 처리군 6 샘플

res = stats.ttest_ind(gd_ctrl, gd_treat)   # 독립표본 t검정
print("대조군 평균:", gd_ctrl.mean().round(2), "/ 처리군 평균:", gd_treat.mean().round(2))
print(f"t = {res.statistic:.3f}, p = {res.pvalue:.5f}, df = {res.df:.0f}")

대조군 평균: 51.03 / 처리군 평균: 59.43
t = -3.258, p = 0.00860, df = 10

GeneD는 처리군 평균(59.43)이 대조군(51.03)보다 뚜렷이 높고, p = 0.0086으로 0.05보다 작아 "차이가 우연이 아니다"라고 볼 만합니다. t값이 음수인 건 ttest_ind(a, b)가 a - b 방향으로 계산하기 때문으로, 대조군이 처리군보다 낮다는 뜻일 뿐 크기는 절댓값으로 봅니다.

이제 8개 유전자 전부를 한 번에 돌립니다. 행마다 검정을 반복하므로 리스트 내포(list comprehension)로 묶습니다.

ctrl, treat = expr[:, :6], expr[:, 6:]    # 전체를 대조군 / 처리군으로 분리

# 유전자별 t검정 → p값 배열
p_t = np.array([stats.ttest_ind(c, t).pvalue for c, t in zip(ctrl, treat)])

for g, p in zip(genes, p_t):
    flag = "*" if p < 0.05 else " "       # 유의수준 0.05 표시
    print(f"{g}: p = {p:.4f} {flag}")

GeneA: p = 0.5466
GeneB: p = 0.8914
GeneC: p = 0.0249 *
GeneD: p = 0.0086 *
GeneE: p = 0.6613
GeneF: p = 0.0020 *
GeneG: p = 0.8680
GeneH: p = 0.4531  

효과를 심은 GeneC·GeneD·GeneF가 p < 0.05로 잡혔습니다. 효과를 0으로 둔 나머지 유전자들은 p값이 크게 나와, 검정이 의도대로 진짜 차이만 골라냈음을 알 수 있죠. 이 세 개를 그대로 "유의한 유전자"라고 결론 내리고 싶겠지만, 8개를 한꺼번에 검정한 순간 다중검정 문제가 시작됩니다(5절).

만-휘트니 U 검정 — 정규성을 못 믿을 때

t검정은 데이터가 대략 정규분포를 따른다고 가정합니다. 표본이 적거나(군당 3~6개는 흔하죠) 분포가 치우쳐 정규성을 장담하기 어려우면, 순위(rank)에 기반한 비모수(non-parametric) 검정인 만-휘트니 U 검정(Mann-Whitney U test)이 더 안전합니다. scipy.stats.mannwhitneyu가 같은 인터페이스로 동작합니다.

# 같은 데이터에 비모수 검정 (순위 기반 → 정규성 가정 불필요)
p_u = np.array([stats.mannwhitneyu(c, t, alternative="two-sided").pvalue
                for c, t in zip(ctrl, treat)])

for g, pt, pu in zip(genes, p_t, p_u):
    print(f"{g}: t검정 p = {pt:.4f} | 만-휘트니 p = {pu:.4f}")

GeneA: t검정 p = 0.5466 | 만-휘트니 p = 1.0000
GeneB: t검정 p = 0.8914 | 만-휘트니 p = 0.8182
GeneC: t검정 p = 0.0249 | 만-휘트니 p = 0.0411
GeneD: t검정 p = 0.0086 | 만-휘트니 p = 0.0260
GeneE: t검정 p = 0.6613 | 만-휘트니 p = 0.7483
GeneF: t검정 p = 0.0020 | 만-휘트니 p = 0.0022
GeneG: t검정 p = 0.8680 | 만-휘트니 p = 0.5745
GeneH: t검정 p = 0.4531 | 만-휘트니 p = 0.3776

두 검정이 같은 유전자(GeneC·GeneD·GeneF)를 유의하게 가리키지만 p값은 조금씩 다릅니다. alternative="two-sided"는 "어느 방향이든 차이가 있는가"를 보는 양측 검정으로, 발현이 오를지 내릴지 모를 때의 기본값입니다. 표본이 충분하고 정규성이 그럴듯하면 t검정이 검정력(통계적으로 차이를 잡아내는 힘)이 더 좋고, 의심스러우면 만-휘트니로 보수적으로 가는 게 무난합니다.

4. scipy.stats — 상관분석 (피어슨·스피어만)

두 유전자가 함께 오르내리는지, 즉 상관(correlation)을 보는 것도 자주 합니다. 피어슨 상관(Pearson correlation, pearsonr)은 두 변수의 선형 관계를, 스피어만 상관(Spearman correlation, spearmanr)은 순위 기반이라 비선형이지만 단조로운(monotonic) 관계까지 잡습니다. 둘 다 상관계수와 p값을 함께 줍니다.

# 함께 올라가도록 효과를 심은 GeneD와 GeneH (전체 12 샘플)
r_p, p_p = stats.pearsonr(expr[3], expr[7])    # 피어슨 (선형)
r_s, p_s = stats.spearmanr(expr[3], expr[7])   # 스피어만 (순위·단조)

print(f"피어슨  r = {r_p:.3f}, p = {p_p:.4f}")
print(f"스피어만 r = {r_s:.3f}, p = {p_s:.4f}")

피어슨  r = 0.430, p = 0.1628
스피어만 r = 0.609, p = 0.0354

여기서 두 상관계수가 꽤 벌어진 게 흥미로운 지점입니다. 스피어만(0.609)이 피어슨(0.430)보다 큰데, 이는 두 유전자의 관계가 직선보다는 "함께 오르되 일정한 비율은 아닌" 단조 관계에 가깝다는 신호입니다. 한두 샘플의 큰 값(이상치)이 피어슨을 끌어내릴 때도 이런 차이가 납니다. 그래서 발현처럼 분포가 매끈하지 않은 데이터에서는 스피어만을 함께 보는 게 안전하죠. 다만 표본이 12개로 적어 p값 해석은 조심해야 합니다.

5. 다중검정 보정 — 원시 p값을 adjusted q값으로

이 글에서 가장 중요한 한 단계입니다. 3절에서 8개 유전자를 검정해 3개가 p < 0.05로 나왔습니다. 그런데 검정을 여러 번 하면, 실제로는 차이가 없는데도 우연히 p < 0.05가 나올 확률이 쌓입니다. 유의수준 0.05로 20번 검정하면 평균 1번은 거짓 양성(false positive)이 나오는 셈이죠. RNA-seq에서는 유전자 2만 개를 한꺼번에 검정하니, 보정 없이는 우연한 "유의" 유전자가 수백 개씩 쏟아집니다. 개념은 다중검정보정 편에서 자세히 풀었습니다.

해결은 다중검정 보정(multiple testing correction)입니다. 그중 생물정보에서 표준으로 쓰는 벤야미니-호크베르그(Benjamini-Hochberg, BH) 절차는 거짓발견율(false discovery rate, FDR)을 통제해, 원시 p값을 보정된 q값(adjusted p-value)으로 바꿉니다. scipy 1.11부터 scipy.stats.false_discovery_control이 이걸 한 줄로 해 줍니다.

# 원시 p값 (3절의 t검정 결과) → BH 보정 q값
q = stats.false_discovery_control(p_t, method="bh")   # method='bh' = Benjamini-Hochberg

for g, p, qv in zip(genes, p_t, q):
    flag = "*" if qv < 0.05 else " "
    print(f"{g}: p = {p:.4f} → q = {qv:.4f} {flag}")

print("원시 p < 0.05 개수:", int((p_t < 0.05).sum()))
print("보정 q < 0.05 개수:", int((q   < 0.05).sum()))

GeneA: p = 0.5466 → q = 0.8745
GeneB: p = 0.8914 → q = 0.8914
GeneC: p = 0.0249 → q = 0.0665
GeneD: p = 0.0086 → q = 0.0344 *
GeneE: p = 0.6613 → q = 0.8817
GeneF: p = 0.0020 → q = 0.0156 *
GeneG: p = 0.8680 → q = 0.8914
GeneH: p = 0.4531 → q = 0.8745
원시 p < 0.05 개수: 3
보정 q < 0.05 개수: 2

보정의 효과가 한눈에 드러납니다. 원시 p값으로는 GeneC·GeneD·GeneF 세 개가 유의했지만, 보정 후에는 GeneC가 q = 0.0665로 밀려나 둘만 남습니다. GeneC의 원시 p값 0.0249는 8번 검정 중 하나로는 "그 정도는 우연히도 나올 수 있다"고 판정된 것이죠. 보정 q값은 항상 원시 p값보다 크거나 같고, 검정 개수가 많을수록 더 가혹하게 올라갑니다. RNA-seq 분석 도구(DESeq2·edgeR 등)가 결과 표에 padj를 따로 두는 이유가 바로 이것입니다 — 봐야 할 건 원시 p가 아니라 보정된 q입니다.

함수가 제대로 동작하는지 공식 문서의 예시로 한 번 더 확인할 수 있습니다.

# scipy 공식 문서 예시 — 결과가 문서와 일치
print(stats.false_discovery_control([0.5, 0.6, 0.1, 0.001]))

[0.6   0.6   0.2   0.004]

statsmodels를 이미 쓰고 있다면 from statsmodels.stats.multitest import multipletests 후 multipletests(p_t, method="fdr_bh")로도 같은 BH 보정을 얻습니다. 결과는 scipy와 동일하니, 의존성에 맞춰 고르면 됩니다.

자주 발생하는 에러 & 해결

• 다중검정 미보정 — 가장 흔하고 치명적인 실수 — 유전자 수천 개를 검정하고 원시 p < 0.05를 그대로 "유의"라고 발표하면, 거짓 양성이 대량으로 섞입니다. 1만 개를 유의수준 0.05로 검정하면 진짜 차이가 하나도 없어도 평균 500개가 우연히 통과하죠. 반드시 scipy.stats.false_discovery_control(또는 statsmodels multipletests)로 보정한 q값으로 거르고, 표·그림에도 adjusted p값을 명시합니다.
• 축(axis) 혼동 — expr.mean()처럼 axis를 빼면 전체 평균(스칼라 하나)이 나옵니다. 유전자별이 필요한데 axis=0을 주면 샘플별 결과가 나와 조용히 틀린 분석으로 이어지죠. "axis는 없어지는 축"으로 외우고, 결과 배열의 모양(shape)을 매번 print로 확인하는 습관이 안전합니다. 정규화·요약에서 행과 열을 바꿔 넣는 실수가 의외로 잦습니다.
• 정수 나눗셈 — 정수 배열을 정수로 나누면 의도와 달리 몫만 남을 수 있습니다. np.array([1,2,3,4]) // 4는 [0 0 0 1]이 되죠(//는 정수 나눗셈). 비율·정규화에는 실수 나눗셈 /를 쓰고, 배열을 astype(float)로 미리 실수화하거나 dtype=float로 만들어 둡니다. 카운트 데이터(read count)를 다룰 때 특히 조심합니다.
• NaN 전파 — 결측값(NaN)이 하나라도 섞이면 mean·sum·std의 결과가 통째로 NaN이 됩니다. np.array([1,2,np.nan,4]).mean()은 nan이죠. 결측을 무시하고 계산하려면 np.nanmean·np.nansum·np.nanstd를 쓰거나, np.isnan으로 걸러낸 뒤 계산합니다. 발현 행렬에 결측이 있으면 검정 전에 처리 방침(제거·대치)을 먼저 정해야 합니다.
• 표본수·정규성 가정 위반 — 군당 2~3개처럼 표본이 너무 적으면 t검정의 정규성 가정이 위태롭고 검정력도 약합니다. 표본이 적거나 분포가 치우쳤으면 만-휘트니 U 검정 같은 비모수 검정을 쓰고, 가능하면 반복(replicate)을 늘립니다. 검정을 고르는 기준이 헷갈리면 R 통계검정 편의 결정 흐름(그룹 수·짝지음·정규성)이 그대로 적용됩니다.
• ttest_ind의 등분산 가정 — ttest_ind는 기본적으로 두 군의 분산이 같다고 가정합니다(Student t). 군마다 분산이 크게 다르면 equal_var=False를 줘 Welch t검정으로 바꾸는 게 안전합니다. 실제 발현 데이터는 분산이 들쭉날쭉한 경우가 많아, Welch를 기본으로 두는 분석자도 많습니다.

확장 — 더 해 볼 것

scipy.cluster로 계층적 군집 — 정규화한 발현 행렬을 scipy.cluster.hierarchy의 linkage·dendrogram에 넣으면 유전자·샘플을 유사도로 묶어 덴드로그램(dendrogram)을 그릴 수 있습니다. 히트맵의 행·열 정렬이 바로 이 군집 결과죠. 거리 계산은 scipy.spatial.distance가 맡습니다.

from scipy.cluster.hierarchy import linkage, fcluster

Z = linkage(z, method="ward")        # z-점수 행렬로 유전자 계층 군집
clusters = fcluster(Z, t=3, criterion="maxclust")   # 3개 군집으로 자르기
print("유전자별 군집 라벨:", clusters)

• scipy.optimize로 곡선 적합 — 용량-반응(dose-response)이나 성장 곡선처럼 모델 함수에 데이터를 맞추는 일은 scipy.optimize.curve_fit이 해 줍니다. 시그모이드·지수 함수를 정의하고 데이터를 넘기면 최적 파라미터와 공분산을 돌려주죠. 효소 반응속도(Michaelis-Menten) 적합에도 같은 방식을 씁니다.
• 검정 선택을 함수로 — 그룹 수·짝지음 여부·정규성에 따라 검정을 자동으로 고르는 래퍼(wrapper)를 짜 두면 분석이 빨라집니다. 정규성은 scipy.stats.shapiro, 등분산은 scipy.stats.levene으로 미리 점검할 수 있습니다.
• 벡터화로 속도 끌어올리기 — 유전자별 검정을 리스트 내포로 돌렸지만, 행이 수만 개면 scipy.stats.ttest_ind(ctrl, treat, axis=1)처럼 axis 인자로 전체를 한 번에 벡터화하면 훨씬 빠릅니다. numpy·scipy의 벡터화 연산은 파이썬 반복문보다 수십~수백 배 빠르죠.
• 분류·머신러닝으로 잇기 — 여기서 만든 expr 배열(또는 z-점수 z)을 그대로 특징 행렬 X로 삼아 scikit-learn 편의 분류·교차검증으로 넘기면, "어떤 유전자 조합이 군을 가르는가"까지 이어집니다.

References

1. Harris, C. R. et al. (2020). Array programming with NumPy. Nature, 585, 357–362. (numpy 원논문)
2. Virtanen, P. et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17, 261–272. (scipy 원논문)
3. NumPy developers. (2026). NumPy v2.5 documentation — Broadcasting · Indexing. numpy.org/broadcasting
4. SciPy developers. (2026). scipy.stats.false_discovery_control — Benjamini-Hochberg·Benjamini-Yekutieli FDR 보정(scipy 1.11+). false_discovery_control docs
5. SciPy developers. (2026). scipy.stats.ttest_ind · mannwhitneyu — 독립표본 검정 API. ttest_ind docs
6. Benjamini, Y. & Hochberg, Y. (1995). Controlling the False Discovery Rate. Journal of the Royal Statistical Society B, 57, 289–300. (BH 절차 원논문)

Pipette & Pipeline · A bio portfolio journal

TL;DR — 유전자 발현이나 임상 특징으로 "종양이 양성인가 악성인가"를 맞히는 이진 분류(binary classification)를, scikit-learn (파이썬 머신러닝 표준 라이브러리)으로 데이터 적재부터 ROC·AUC 평가까지 한 번에 끝냅니다. 재현 가능한 sklearn 내장 데이터셋 load_breast_cancer (569 샘플 × 30 특징)를 쓰고, 로지스틱 회귀(logistic regression)와 랜덤 포레스트(random forest) 두 모델을 비교합니다.

머신러닝 코드가 길어 보여도 뼈대는 단순합니다. 데이터를 학습용·시험용으로 나누고(train/test split), 스케일을 맞추고(scaling), 모델을 학습시키고(fit), 교차검증(cross-validation)으로 성능을 가늠한 뒤, 혼동행렬(confusion matrix)과 ROC 곡선으로 평가합니다. 이 흐름만 손에 익으면 데이터만 바꿔 어떤 분류 문제에도 그대로 씁니다.

이 글의 코드는 무작위 시드를 고정해 그대로 복사하면 같은 숫자가 나옵니다. 초보가 가장 많이 밟는 지뢰 — 데이터 누수(data leakage)·클래스 불균형(class imbalance)·과적합(overfitting)도 코드로 짚고, 마지막엔 Pipeline (파이프라인)으로 누수를 원천 차단하는 법까지 갑니다.

🔗 관련 글  ·  모델에 넣을 데이터를 먼저 다듬으려면 pandas로 생물 데이터 정제하기에서 이어집니다  ·  학습 결과를 그림으로 그리는 도구는 matplotlib·seaborn으로 생물 데이터 시각화에서 다룹니다  ·  같은 ROC·AUC를 R로 그리는 법은 pROC로 ROC 곡선·AUC 그리기와 1:1로 비교됩니다  ·  특징을 줄여 구조를 보는 차원축소는 주성분분석(PCA)으로 RNA-seq 샘플 QC하기로 이어집니다

이 글에서 만들 것

발현 데이터나 임상 수치를 받았을 때 가장 흔한 질문이 "이 샘플이 어느 그룹에 속하나"입니다. 종양이 양성인지 악성인지, 환자가 약에 반응할지 안 할지, 세포가 어떤 아형(subtype)인지 — 모두 분류 문제입니다. scikit-learn은 이런 분류를 몇 줄로 풀게 해 주는 파이썬 표준 도구입니다.

scikit-learn (David Cournapeau, 2007~)은 분류·회귀·군집·차원축소를 한 일관된 문법으로 묶은 라이브러리입니다. 어떤 모델이든 fit(학습)·predict(예측)·score(평가)라는 같은 인터페이스를 쓰기 때문에, 모델을 바꿔도 코드 골격이 거의 그대로입니다. 딥러닝이 필요 없는 표 형식(tabular) 데이터에서는 여전히 가장 먼저 잡는 도구죠.

이 글은 sklearn에 내장된 유방암 진단 데이터(load_breast_cancer)로 이진 분류를 처음부터 끝까지 만듭니다. ① 데이터를 적재하고, ② 학습용·시험용으로 나누고(train_test_split), ③ 스케일을 맞추고(StandardScaler), ④ 로지스틱 회귀와 랜덤 포레스트를 학습시키고, ⑤ 교차검증(cross_val_score·StratifiedKFold)으로 성능을 가늠한 뒤, ⑥ 혼동행렬·classification_report·ROC·AUC로 평가하고 특징 중요도(feature importance)까지 봅니다. 코드는 scikit-learn 1.9 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, RocCurveDisplay

X, y = load_breast_cancer(return_X_y=True)                       # ① 데이터 적재
Xtr, Xte, ytr, yte = train_test_split(X, y, stratify=y, random_state=42)  # ② 분할
Xtr_s = StandardScaler().fit_transform(Xtr)                      # ③ 스케일링(train fit)
clf = LogisticRegression(max_iter=5000).fit(Xtr_s, ytr)         # ④ 학습
print(cross_val_score(clf, Xtr_s, ytr, cv=5, scoring="roc_auc"))  # ⑤ 교차검증
RocCurveDisplay.from_estimator(clf, Xte_s, yte)                  # ⑥ ROC·AUC

1. 사전 준비 — 설치와 데이터 로드

scikit-learn 하나만 설치하면 됩니다. numpy·scipy는 의존성으로 따라오고, 그림을 그리려면 matplotlib도 함께 깔아 둡니다.

# pip로 설치 (가상환경 권장)
pip install scikit-learn matplotlib

# 또는 conda
conda install -c conda-forge scikit-learn matplotlib

import sklearn
print("scikit-learn", sklearn.__version__)

scikit-learn 1.9.0

이제 데이터를 불러옵니다. load_breast_cancer는 위스콘신 유방암 진단 데이터로, 종양의 디지털 이미지에서 뽑은 30개 수치 특징(반지름·질감·둘레·면적·매끄러움 등)과, 각 종양이 양성(benign)인지 악성(malignant)인지를 담고 있습니다. 569개 샘플 전부가 라벨이 붙어 있어 분류 연습에 딱 맞습니다.

import numpy as np
import pandas as pd
from sklearn.datasets import load_breast_cancer

data = load_breast_cancer(as_frame=True)     # as_frame=True → pandas로 받기
X = data.data                                 # 특징 30개 (569 × 30)
y = data.target                               # 라벨 (0/1)

print(X.shape)                                # (샘플 수, 특징 수)
print(data.target_names)                      # 클래스 이름
print(y.value_counts())                       # 클래스별 개수

(569, 30)
['malignant' 'benign']
target
1    357
0    212
Name: count, dtype: int64

여기서 꼭 짚어야 할 함정이 하나 있습니다. 클래스 라벨이 직관과 반대로 매겨져 있습니다 — 0이 악성(malignant), 1이 양성(benign)입니다. 의학적으로 "잡아야 하는 양성(positive) 사건"은 악성인데, 숫자로는 0이라 나중에 ROC·재현율(recall)을 해석할 때 헷갈리기 쉽습니다. 데이터를 받으면 라벨이 무엇을 가리키는지부터 확인하는 습관이 중요합니다.

특징이 30개라 어떤 값인지 잠깐 살펴봅니다. 평균(mean)·표준오차(error)·최악값(worst) 세 묶음으로 같은 측정치가 반복됩니다.

X.iloc[:3, :5]                                # 앞 3 샘플, 앞 5 특징만

   mean radius  mean texture  mean perimeter  mean area  mean smoothness
0        17.99         10.38          122.80     1001.0          0.11840
1        20.57         17.77          132.90     1326.0          0.08474
2        19.69         21.25          130.00     1203.0          0.10960

값의 크기(scale)가 제각각인 게 보입니다. mean area는 1000 단위인데 mean smoothness는 0.1 단위죠. 이렇게 스케일이 다른 특징을 그대로 거리 기반·기울기 기반 모델에 넣으면 큰 값이 학습을 지배합니다. 그래서 3절에서 스케일을 맞춥니다.

2. 데이터 분할 — train_test_split으로 시험셋 떼어 두기

머신러닝의 첫 원칙은 "모델이 본 적 없는 데이터로 평가한다"입니다. 학습에 쓴 데이터로 성능을 재면 외워서 잘 맞히는 것인지, 정말 일반화(generalization)된 것인지 구분이 안 됩니다. 그래서 데이터를 학습용(train)과 시험용(test)으로 먼저 갈라 두고, 시험셋은 마지막 평가 전까지 건드리지 않습니다.

train_test_split이 이 분할을 한 줄로 해 줍니다. 두 옵션이 핵심입니다 — stratify로 클래스 비율을 학습·시험 양쪽에서 똑같이 유지하고, random_state로 분할을 고정해 재현합니다.

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.25,        # 시험셋 25% (학습 75%)
    stratify=y,            # 클래스 비율 유지 (불균형 데이터 필수)
    random_state=42,       # 분할 고정 → 재현성
)
print(X_train.shape, X_test.shape)
print(y_train.value_counts(normalize=True).round(3))   # 학습셋 클래스 비율
print(y_test.value_counts(normalize=True).round(3))    # 시험셋 클래스 비율

(426, 30) (143, 30)
target
1    0.627
0    0.373
Name: proportion, dtype: float64
target
1    0.629
0    0.371
Name: proportion, dtype: float64

stratify=y를 준 덕분에 학습셋과 시험셋의 클래스 비율(약 63 대 37)이 거의 같습니다. 이걸 빼면 무작위로 나뉘면서 한쪽에 악성이 몰릴 수 있고, 특히 클래스가 불균형하거나 샘플이 적을 때 평가가 들쭉날쭉해집니다. 분류 문제에서는 거의 항상 stratify를 켜는 게 안전합니다.

3. 스케일링 — StandardScaler와 데이터 누수 막기

1절에서 봤듯 특징마다 값의 크기가 다릅니다. StandardScaler는 각 특징을 평균 0·표준편차 1로 표준화(standardization)해, 모든 특징을 같은 출발선에 세웁니다. 로지스틱 회귀처럼 기울기로 학습하는 모델은 스케일링 여부에 따라 수렴 속도와 성능이 꽤 달라집니다.

여기서 머신러닝 입문자가 가장 많이 밟는 지뢰가 나옵니다 — 스케일러는 학습셋으로만 fit해야 합니다. 시험셋까지 합쳐 평균·표준편차를 구하면, 시험셋의 정보가 학습 과정에 새어 들어가는 데이터 누수(data leakage)가 됩니다. 평가가 부풀려져 실제 성능을 과대평가하게 되죠.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)   # 학습셋: fit(평균·표준편차 학습) + 변환
X_test_s  = scaler.transform(X_test)        # 시험셋: 변환만 (학습셋 통계로)

print(X_train_s.mean(axis=0).round(2)[:5])  # 학습셋 각 특징 평균 ≈ 0
print(X_train_s.std(axis=0).round(2)[:5])   # 학습셋 각 특징 표준편차 ≈ 1

[ 0.  0. -0.  0.  0.]
[1. 1. 1. 1. 1.]

핵심은 fit_transform과 transform의 구분입니다. 학습셋에는 fit_transform(통계를 배우고 동시에 변환), 시험셋에는 transform(학습셋이 배운 통계로 변환만)을 씁니다. 시험셋의 평균이 정확히 0이 아닌 건 정상입니다 — 학습셋 기준으로 변환했으니까요. 이 분리를 손으로 하면 실수하기 쉬워서, 5절에서 Pipeline으로 자동화합니다.

4. 모델 학습 — 로지스틱 회귀와 랜덤 포레스트

이제 모델을 학습시킵니다. 성격이 다른 두 분류기를 비교합니다. 로지스틱 회귀(LogisticRegression)는 특징의 가중합으로 확률을 내는 선형 모델로, 빠르고 해석이 쉽습니다. 랜덤 포레스트(RandomForestClassifier)는 결정 트리(decision tree) 수백 개를 모아 투표시키는 앙상블(ensemble) 모델로, 비선형 관계와 특징 간 상호작용을 잘 잡습니다.

scikit-learn의 일관된 문법 덕에, 두 모델 모두 fit으로 학습하고 predict로 예측합니다. 모델 객체만 바꾸면 나머지 코드는 그대로죠.

from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

# ① 로지스틱 회귀 (선형·빠름·해석 쉬움)
logit = LogisticRegression(max_iter=5000, random_state=42)
logit.fit(X_train_s, y_train)               # 스케일링한 학습셋으로 학습

# ② 랜덤 포레스트 (비선형·앙상블·트리 기반)
rf = RandomForestClassifier(n_estimators=300, random_state=42)
rf.fit(X_train, y_train)                     # 트리 모델은 스케일링 불필요

# 학습셋 정확도(accuracy)로 빠르게 확인
print("로지스틱 회귀 train acc:", round(logit.score(X_train_s, y_train), 3))
print("랜덤 포레스트 train acc:", round(rf.score(X_train, y_train), 3))

로지스틱 회귀 train acc: 0.988
랜덤 포레스트 train acc: 1.0

max_iter=5000은 로지스틱 회귀가 충분히 수렴하도록 반복 횟수를 늘린 것입니다(기본값 100이면 수렴 경고가 뜰 수 있습니다). 랜덤 포레스트는 트리 기반이라 특징 스케일에 둔감해서 스케일링하지 않은 원본 X_train을 그대로 넣었습니다 — 트리는 "이 값보다 큰가/작은가"로 나누기 때문에 크기 자체는 상관없죠.

여기서 랜덤 포레스트의 학습셋 정확도가 1.0인 게 눈에 걸립니다. 학습 데이터를 완벽히 맞혔다는 건데, 이게 곧 좋은 모델이라는 뜻은 아닙니다. 학습셋만 외운 과적합(overfitting)일 수 있어서, 반드시 본 적 없는 데이터로 다시 확인해야 합니다. 그 일을 교차검증이 합니다.

5. 교차검증 — cross_val_score로 성능을 정직하게

시험셋 한 번의 점수는 우연에 휘둘립니다. 어쩌다 쉬운 샘플이 시험셋에 몰리면 점수가 높게 나오죠. 교차검증(cross-validation)은 학습 데이터를 여러 조각(fold)으로 나눠, 번갈아 가며 한 조각으로 평가하고 나머지로 학습하기를 반복해, 성능의 평균과 흔들림(표준편차)을 함께 봅니다.

분류에서는 각 조각의 클래스 비율을 유지하는 StratifiedKFold를 씁니다. cross_val_score에 cv=5만 줘도 분류 문제라면 내부적으로 계층 분할을 쓰지만, 시드를 고정하려면 객체로 명시하는 게 좋습니다.

from sklearn.model_selection import cross_val_score, StratifiedKFold

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)  # 5겹·클래스 비율 유지

# 평가 지표로 ROC AUC 사용 (불균형에 강건)
auc_logit = cross_val_score(logit, X_train_s, y_train, cv=cv, scoring="roc_auc")
auc_rf    = cross_val_score(rf,    X_train,   y_train, cv=cv, scoring="roc_auc")

print("로지스틱 회귀 AUC:", auc_logit.round(3), "평균", round(auc_logit.mean(), 3))
print("랜덤 포레스트 AUC:", auc_rf.round(3),    "평균", round(auc_rf.mean(), 3))

로지스틱 회귀 AUC: [0.995 0.995 0.999 0.996 0.991] 평균 0.995
랜덤 포레스트 AUC: [0.987 0.995 0.997 0.989 0.985] 평균 0.991

다섯 조각의 AUC가 모두 0.98 이상이고 흔들림도 작습니다. 학습셋 정확도 1.0이던 랜덤 포레스트도 교차검증 AUC는 0.991로, 외운 게 아니라 실제로 일반화됐음이 확인됩니다. 이 데이터에서는 단순한 로지스틱 회귀(평균 0.995)가 랜덤 포레스트보다 오히려 살짝 앞섭니다. 특징이 잘 정제돼 있고 클래스가 선형으로 꽤 나뉘기 때문인데, 늘 복잡한 모델이 이기는 건 아니라는 걸 보여 주죠.

scoring에 "roc_auc"를 쓴 이유가 있습니다. 정확도(accuracy)는 클래스가 불균형하면 속기 쉽습니다 — 95%가 음성인 데이터라면 전부 음성이라 찍어도 정확도 95%가 나오죠. AUC는 양성·음성을 가르는 능력 자체를 보기 때문에 불균형에 더 강건합니다.

6. 평가 — 혼동행렬·classification_report·ROC·AUC

교차검증으로 모델을 골랐으니, 이제 떼어 둔 시험셋으로 최종 평가합니다. 분류 평가는 정확도 한 숫자로 끝내지 않습니다. 어디서 틀렸는지(혼동행렬), 정밀도·재현율은 어떤지(classification_report), 임계값(threshold)을 바꾸면 어떻게 되는지(ROC 곡선)를 함께 봅니다.

먼저 혼동행렬(confusion matrix)입니다. 예측과 실제를 2×2 표로 교차해, 참양성(TP)·거짓양성(FP)·거짓음성(FN)·참음성(TN)이 각각 몇 개인지 보여 줍니다.

from sklearn.metrics import confusion_matrix, classification_report

y_pred = logit.predict(X_test_s)             # 시험셋 예측

cm = confusion_matrix(y_test, y_pred)        # 행=실제, 열=예측
print(cm)
print(classification_report(y_test, y_pred,
                            target_names=data.target_names))

[[ 51   2]
 [  1  89]]
              precision    recall  f1-score   support

   malignant       0.98      0.96      0.97        53
      benign       0.98      0.99      0.98        90

    accuracy                           0.98       143
   macro avg       0.98      0.98      0.98       143
weighted avg       0.98      0.98      0.98       143

혼동행렬을 읽으면, 악성 53개 중 51개를 맞히고 2개를 양성으로 잘못 봤으며(거짓음성), 양성 90개 중 89개를 맞혔습니다. classification_report는 여기서 정밀도(precision)·재현율(recall)·F1 점수를 클래스별로 풀어 줍니다. 의학에서는 악성을 놓치는 일(거짓음성)이 특히 위험하므로 악성의 재현율 0.96을 눈여겨봅니다. 단일 숫자인 정확도(0.98)만 보면 이런 비대칭을 놓치죠.

이제 ROC 곡선과 AUC입니다. ROC (Receiver Operating Characteristic)는 분류 임계값을 0에서 1까지 움직이며 거짓양성률(가로)과 참양성률(세로)을 그린 곡선이고, 그 아래 면적이 AUC (Area Under the Curve)입니다. AUC가 1에 가까울수록 완벽한 분류, 0.5면 동전 던지기 수준입니다.

import matplotlib.pyplot as plt
from sklearn.metrics import roc_auc_score, RocCurveDisplay

# 양성(class 1) 확률 예측 → AUC 계산
y_score = logit.predict_proba(X_test_s)[:, 1]
print("시험셋 AUC:", round(roc_auc_score(y_test, y_score), 3))

fig, ax = plt.subplots(figsize=(6, 6))
RocCurveDisplay.from_estimator(logit, X_test_s, y_test,
                               name="로지스틱 회귀", ax=ax)
RocCurveDisplay.from_estimator(rf, X_test, y_test,
                               name="랜덤 포레스트", ax=ax,
                               plot_chance_level=True)   # 대각선 기준선
ax.set_title("ROC 곡선 — 시험셋")
plt.tight_layout()
plt.savefig("roc.png", dpi=300, bbox_inches="tight")

시험셋 AUC: 0.995

RocCurveDisplay.from_estimator는 모델·데이터·라벨만 주면 ROC 곡선을 곡선 아래 면적(AUC)과 함께 자동으로 그려 줍니다. 같은 ax에 두 모델을 겹쳐 그리면 성능을 한 화면에서 비교할 수 있죠. plot_chance_level=True가 대각선(AUC=0.5) 기준선을 더해 줍니다. predict_proba(...)[:, 1]로 양성 클래스의 확률을 뽑아 roc_auc_score에 직접 넣을 수도 있는데, 두 방법의 결과는 같습니다.

마지막으로 특징 중요도(feature importance)입니다. 랜덤 포레스트는 각 특징이 분류에 얼마나 기여했는지를 feature_importances_로 알려 줘, "어떤 측정치가 양성·악성을 가르는가"를 해석하게 해 줍니다.

import pandas as pd

# 특징 중요도 상위 8개
imp = (pd.Series(rf.feature_importances_, index=X.columns)
       .sort_values(ascending=False).head(8))
print(imp.round(3))

worst perimeter        0.142
worst concave points   0.139
worst area             0.121
mean concave points    0.106
worst radius           0.082
mean perimeter         0.052
worst texture          0.038
mean area              0.037

종양 둘레·면적·오목한 점(concave points)의 "최악값(worst)" 측정치가 분류를 주도합니다. 임상적으로도 악성 종양이 크고 불규칙한 경계를 갖는다는 사실과 맞아떨어지죠. 다만 트리 기반 중요도는 값의 범위가 넓은 특징을 과대평가하는 편향이 있어, 더 엄밀하게 보려면 순열 중요도(permutation importance)를 함께 쓰는 게 좋습니다.

자주 발생하는 에러 & 해결

• 데이터 누수(data leakage) — 스케일러를 전체에 fit — StandardScaler().fit(X)로 전체 데이터(학습+시험)에 한꺼번에 fit하면, 시험셋의 평균·표준편차가 변환에 스며들어 평가가 부풀려집니다. 교차검증에서 누수가 특히 잘 일어나는데, 매 fold마다 전처리를 다시 fit해야 하기 때문입니다. 해결은 전처리와 모델을 Pipeline으로 묶는 것입니다(아래 확장 참고). 그러면 cross_val_score가 fold마다 스케일러를 학습셋에만 다시 fit해 누수를 원천 차단합니다.
• 클래스 불균형(class imbalance) — 정확도에 속음 — 음성이 95%인 데이터에서 전부 음성으로 찍어도 정확도 95%가 나옵니다. 정확도만 보면 쓸모없는 모델을 좋다고 착각하죠. 평가 지표를 ROC AUC·F1·재현율로 바꾸고, 모델에 class_weight="balanced"를 주거나 소수 클래스를 늘리는 SMOTE (imbalanced-learn 패키지) 같은 기법을 씁니다. 분할·교차검증에는 반드시 stratify·StratifiedKFold로 비율을 유지합니다.
• 과적합(overfitting) — 학습셋만 잘 맞힘 — 학습 정확도는 1.0인데 시험 성능이 뚝 떨어지면 과적합입니다. 모델이 데이터를 일반화한 게 아니라 외운 것이죠. 랜덤 포레스트는 max_depth로 트리 깊이를 제한하거나 min_samples_leaf를 키워 규제(regularization)하고, 로지스틱 회귀는 C를 줄여 규제를 강화합니다. 무엇보다 학습 점수가 아니라 교차검증·시험셋 점수로 판단하는 게 핵심입니다.
• 수렴 경고(ConvergenceWarning) — LogisticRegression이 "lbfgs failed to converge"를 띄우면, 정해진 반복 안에 최적해를 못 찾았다는 뜻입니다. max_iter를 늘리거나(예: 5000), 특징을 StandardScaler로 스케일링하면 대개 사라집니다. 스케일이 제각각이면 최적화가 느려지기 때문에, 스케일링은 수렴 문제의 해결책이기도 합니다.
• predict_proba vs decision_function — ROC·AUC에는 확률 또는 점수가 필요한데, roc_auc_score(y_test, y_pred)처럼 0/1 예측 라벨을 넣으면 곡선이 계단처럼 거칠어지고 AUC가 과소평가됩니다. 반드시 predict_proba(X)[:, 1](양성 확률)이나 decision_function(X)(결정 점수)를 넣어야 합니다. RocCurveDisplay.from_estimator는 이걸 알아서 처리해 주니 가장 안전합니다.
• 시험셋을 모델 선택에 씀 — 시험셋을 보면서 하이퍼파라미터를 고르면, 시험셋이 사실상 학습에 쓰인 셈이라 성능이 부풀려집니다. 시험셋은 마지막 단 한 번의 평가에만 씁니다. 모델·하이퍼파라미터 선택은 학습셋 안의 교차검증이나 별도 검증셋(validation set)으로 해야 합니다.

확장 — 더 해 볼 것

Pipeline으로 누수 차단 — 전처리와 모델을 하나로 묶으면, fit·교차검증 어디서든 스케일러가 학습셋에만 fit돼 데이터 누수가 구조적으로 사라집니다. 코드도 짧아지죠.

from sklearn.pipeline import make_pipeline

pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))
auc = cross_val_score(pipe, X, y, cv=cv, scoring="roc_auc")   # 누수 없는 교차검증
print("Pipeline AUC 평균:", round(auc.mean(), 3))

GridSearchCV로 하이퍼파라미터 탐색 — 모델 설정(예: 랜덤 포레스트의 트리 깊이, 로지스틱 회귀의 규제 강도 C)을 격자(grid)로 훑어 교차검증 점수가 가장 높은 조합을 자동으로 고릅니다. Pipeline과 함께 쓰면 전처리 옵션까지 한꺼번에 탐색할 수 있습니다.

from sklearn.model_selection import GridSearchCV

grid = {"randomforestclassifier__max_depth": [4, 8, None],
        "randomforestclassifier__min_samples_leaf": [1, 3, 5]}
search = GridSearchCV(make_pipeline(RandomForestClassifier(random_state=42)),
                      grid, cv=cv, scoring="roc_auc")
search.fit(X, y)
print("최적 파라미터:", search.best_params_, "| AUC:", round(search.best_score_, 3))

• 순열 중요도(permutation importance) — 트리 기반 feature_importances_의 편향을 피하려면, 특징 값을 무작위로 섞었을 때 성능이 얼마나 떨어지는지로 중요도를 재는 permutation_importance가 더 공정합니다. 모델 종류와 무관하게 쓸 수 있습니다.
• 다중 분류(multiclass)·실제 발현 데이터 — 이 글은 이진 분류였지만, 세포 아형 3종 이상을 가르는 다중 분류도 같은 코드로 됩니다(roc_auc_score에 multi_class="ovr"). 실제 RNA-seq 발현 행렬을 쓰려면 pandas 정제로 표를 다듬어 X에 넣으면 그대로 이어집니다.
• 모델 저장·배포 — 학습한 모델을 joblib.dump로 파일에 저장해 두면, 나중에 불러와 새 샘플을 예측할 수 있습니다. 다만 sklearn 버전이 바뀌면 호환이 깨질 수 있으니 버전을 함께 기록하는 게 안전합니다.

References

1. Pedregosa, F. et al. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12, 2825–2830. (scikit-learn 원논문)
2. scikit-learn developers. (2026). scikit-learn 1.9 documentation — User Guide. scikit-learn.org/user_guide
3. scikit-learn. Common pitfalls and recommended practices — 데이터 누수·교차검증 함정. common_pitfalls
4. scikit-learn. load_breast_cancer — 위스콘신 유방암 진단 데이터셋(569 × 30). load_breast_cancer docs
5. scikit-learn. RocCurveDisplay — ROC 곡선·AUC 시각화 API. RocCurveDisplay docs

Pipette & Pipeline · A bio portfolio journal

TL;DR — 차등발현(differential expression) 결과나 발현 행렬을 받았을 때 가장 많이 그리는 그림 네 가지를 matplotlib (파이썬 저수준 플로팅 라이브러리)와 seaborn (그 위에 얹은 고수준 통계 시각화 라이브러리)로 처음부터 끝까지 그립니다. 볼케이노 플롯(volcano plot), 클러스터 히트맵(clustermap), 박스·바이올린 플롯(box/violin plot), 그리고 facet 다중 패널까지입니다.

역할 분담이 분명합니다. matplotlib은 점·선·축을 직접 그리는 바탕이고, seaborn은 "통계 그래프 한 장"을 짧은 명령으로 그려 줍니다. 볼케이노처럼 세밀한 라벨링이 필요한 그림은 matplotlib으로 손보고, 분포 비교나 히트맵은 seaborn 한 줄로 끝냅니다.

이 글의 코드는 모두 무작위 시드를 고정한 작은 발현 데이터로 돌아가, 그대로 복사하면 같은 그림이 나옵니다. seaborn 0.13에서 바뀐 hue·팔레트 규칙처럼 자주 막히는 지점도 함께 짚습니다.

🔗 관련 글  ·  그릴 데이터를 먼저 정제하려면 pandas로 발현 행렬 정제하기에서 이어집니다  ·  같은 그래프를 R로 그리는 법은 ggplot2로 바이오 데이터 그래프 그리기와 1:1로 비교됩니다  ·  단일세포 데이터의 UMAP·클러스터 시각화는 scanpy로 단일세포 RNA-seq 분석하기에서 다룹니다

이 글에서 만들 것

발현 분석을 끝내고 나면 결국 그림으로 보여 줘야 합니다. "어떤 유전자가 처리군에서 확 올라갔나"는 볼케이노 플롯으로, "샘플들이 발현 패턴으로 어떻게 묶이나"는 클러스터 히트맵으로, "관심 유전자가 그룹마다 어떻게 분포하나"는 박스·바이올린 플롯으로 봅니다. 이 세 장에 다중 패널(facet) 한 장을 더하면 RNA-seq 결과 슬라이드의 대부분이 채워집니다.

matplotlib (John Hunter, 2003~)은 파이썬 그래프의 바탕입니다. 점 하나, 선 하나, 글자 하나를 좌표에 직접 찍는 저수준 도구라 못 그리는 게 거의 없지만, 그만큼 손이 많이 갑니다. seaborn (Michael Waskom, 2012~)은 그 위에 얹혀, 데이터프레임과 그릴 열 이름만 주면 통계 그래프 한 장을 만들어 줍니다. 둘은 경쟁 관계가 아니라 층위가 다릅니다 — seaborn으로 큰 그림을 그리고, matplotlib으로 축·라벨·주석을 마무리하는 식으로 같이 씁니다.

이 글은 작은 차등발현 결과와 발현 행렬을 직접 만들어, 네 가지 그림을 순서대로 그립니다. ① 볼케이노 플롯(산점도와 유의 영역 색칠, 상위 유전자 라벨), ② 클러스터 히트맵(z-score 정규화와 계층 군집), ③ 박스·바이올린 플롯(그룹별 발현 분포), ④ facet 다중 패널(catplot으로 여러 유전자를 한 번에)입니다. 마지막에 한글 폰트·저장 해상도(dpi)·SVG 같은 실무 설정을 정리합니다. 코드는 matplotlib 3.11·seaborn 0.13.2 기준이고, 시드를 고정해 재현됩니다.

# 핵심 흐름 한눈에
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style="whitegrid")                # ① 전역 스타일 한 번에
sns.scatterplot(data=deg, x="log2FC", y="neglog10p", hue="sig")   # ② 볼케이노
sns.clustermap(expr_z, cmap="vlag", center=0)   # ③ z-score 히트맵 + 군집
sns.violinplot(data=long, x="group", y="expr")  # ④ 그룹별 분포
sns.catplot(data=long, x="group", y="expr", col="gene", kind="box")  # ⑤ facet
plt.savefig("figure.png", dpi=300, bbox_inches="tight")            # ⑥ 저장

1. 사전 준비 — 설치와 샘플 데이터

matplotlib·seaborn·pandas를 한 번에 설치합니다. seaborn을 깔면 matplotlib·numpy·pandas는 의존성으로 따라오지만, 명시해 두는 편이 분명합니다. 클러스터 히트맵의 계층 군집은 SciPy (사이파이)가 필요한데, 이 역시 보통 같이 깔려 있습니다.

# pip로 설치 (가상환경 권장)
pip install matplotlib seaborn pandas

# 또는 conda
conda install -c conda-forge matplotlib seaborn pandas

import matplotlib, seaborn as sns, pandas as pd, numpy as np

print("matplotlib", matplotlib.__version__)
print("seaborn", sns.__version__)

matplotlib 3.11.0
seaborn 0.13.2

이제 설명용으로 작은 데이터 두 개를 만듭니다. 하나는 차등발현 결과 표(유전자별 log2 배수변화와 p-value), 다른 하나는 발현 행렬(유전자 × 샘플)입니다. 실제 분석이라면 앞쪽은 DESeq2·edgeR의 결과 테이블이고, 뒤쪽은 정규화한 카운트 행렬이 됩니다.

import numpy as np
import pandas as pd

rng = np.random.default_rng(7)              # 재현성 고정 (시드 7)
n = 300                                      # 유전자 300개

# ── 차등발현 결과 표: 대부분은 변화 없고, 일부만 크게 오르내림 ──
log2fc = rng.normal(0, 1.0, n)              # log2 배수변화 (대칭 분포)
log2fc[:15] += rng.normal(3, 0.4, 15)       # 앞 15개 = 강하게 상향
log2fc[15:30] -= rng.normal(3, 0.4, 15)     # 다음 15개 = 강하게 하향
pval = 10 ** (-np.abs(log2fc) * rng.uniform(0.8, 1.6, n))  # |FC| 클수록 작은 p
pval = np.clip(pval, 1e-12, 1.0)            # p-value 범위 보정

deg = pd.DataFrame({
    "gene":   [f"G{i:03d}" for i in range(n)],
    "log2FC": log2fc,
    "pvalue": pval,
})
deg["neglog10p"] = -np.log10(deg["pvalue"])  # y축에 쓸 -log10(p)
deg.head(3)

   gene    log2FC        pvalue  neglog10p
0  G000  3.412847  1.073229e-04   3.969310
1  G001  3.001210  3.018501e-04   3.520282
2  G002  3.880644  9.142122e-06   5.039459

차등발현 그림에서 "유의(significant)" 여부는 보통 두 기준을 동시에 봅니다 — 배수변화가 충분히 크고(|log2FC| ≥ 1), p-value가 충분히 작은(p < 0.05) 유전자입니다. 이 둘을 합쳐 색칠용 라벨 열을 미리 만들어 둡니다.

# 유의 기준: |log2FC| >= 1 그리고 p < 0.05 → 상향/하향/비유의 3분류
fc_thr, p_thr = 1.0, 0.05
cond_up   = (deg["log2FC"] >=  fc_thr) & (deg["pvalue"] < p_thr)
cond_down = (deg["log2FC"] <= -fc_thr) & (deg["pvalue"] < p_thr)

deg["sig"] = np.select([cond_up, cond_down],
                        ["up", "down"], default="ns")  # ns = not significant
deg["sig"].value_counts()

sig
ns      268
up       17
down     15
Name: count, dtype: int64

다음은 발현 행렬입니다. 히트맵·분포 그림에 쓸 작은 행렬을 따로 만듭니다. 유전자 12개 × 샘플 10개(대조군 5와 처리군 5), 처리군에서 일부 유전자가 올라가도록 신호를 넣습니다.

# ── 발현 행렬: 유전자 12 × 샘플 10 (control 5 + treated 5) ──
genes   = [f"GENE{i:02d}" for i in range(1, 13)]
samples = [f"C{i}" for i in range(1, 6)] + [f"T{i}" for i in range(1, 6)]

base = rng.normal(8, 1.2, size=(12, 10))    # 로그 발현값 바탕 (log2 스케일)
base[:4, 5:] += rng.normal(2.5, 0.4, (4, 5))   # 앞 4개 유전자 = 처리군서 상향
base[4:7, 5:] -= rng.normal(2.0, 0.4, (3, 5))  # 다음 3개 = 처리군서 하향

expr = pd.DataFrame(base, index=genes, columns=samples).round(2)
expr.iloc[:4, :6]                            # 일부만 미리보기

          C1    C2    C3    C4    C5    T1
GENE01  7.42  6.88  9.10  8.05  7.61  10.33
GENE02  8.21  9.04  7.55  8.66  8.12  10.78
GENE03  6.95  7.80  8.43  7.21  8.90   9.52
GENE04  9.13  8.27  7.92  8.48  7.74  11.06

2. 볼케이노 플롯 — 산점도로 차등발현 한눈에

볼케이노 플롯은 x축에 log2 배수변화, y축에 −log₁₀ p를 찍는 산점도입니다. 화산이 분출하듯 좌우 위쪽으로 점이 솟아오른 모양이라 이런 이름이 붙었습니다. 오른쪽 위는 "크게 올라가고 유의한" 유전자, 왼쪽 위는 "크게 내려가고 유의한" 유전자입니다.

seaborn의 scatterplot에 색 기준(hue)으로 아까 만든 sig 열을 넘기면, 상향·하향·비유의가 자동으로 다른 색이 됩니다.

import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme(style="whitegrid", font_scale=1.0)   # 전역 테마

fig, ax = plt.subplots(figsize=(7, 6))
palette = {"up": "#E76F51", "down": "#2E7D6B", "ns": "#C9C4BA"}  # 색 직접 지정

sns.scatterplot(data=deg, x="log2FC", y="neglog10p",
                hue="sig", hue_order=["up", "down", "ns"],
                palette=palette, s=28, edgecolor="none", ax=ax)

# 기준선: 세로 |log2FC|=1, 가로 p=0.05 (−log10 0.05 ≈ 1.30)
ax.axvline(-1, color="#999", ls="--", lw=1)
ax.axvline( 1, color="#999", ls="--", lw=1)
ax.axhline(-np.log10(0.05), color="#999", ls="--", lw=1)

ax.set_xlabel("log$_2$ fold change")        # 아래첨자: $_2$
ax.set_ylabel("$-$log$_{10}$ $p$")          # 통계기호 p는 수식($)으로 이탤릭
ax.set_title("Volcano plot — control vs treated", weight="bold")
ax.legend(title="유의성", frameon=True)
plt.tight_layout()
plt.savefig("volcano.png", dpi=300, bbox_inches="tight")

여기서 핵심은 세 가지입니다. hue에 범주 열을 주면 seaborn이 범례까지 자동으로 달아 줍니다. palette를 딕셔너리로 주면 "up은 코랄, down은 그린"처럼 색을 못 박을 수 있습니다(브랜드색 통일에 유용합니다). 축 라벨의 $-$log$_{10}$ $p$처럼 $...$로 감싼 부분은 LaTeX 수식으로 렌더되는데, 통계기호 p는 이렇게 수식 안에 넣어야 이탤릭으로 나옵니다.

이제 상위 유전자에 이름표를 답니다. seaborn에는 점 라벨링 기능이 따로 없어서, 이 부분은 matplotlib의 annotate로 직접 그립니다. −log₁₀ p가 가장 큰 상·하향 유전자 몇 개만 골라 텍스트를 얹습니다.

# 상·하향에서 각각 가장 유의한 유전자 3개씩 라벨링
top = pd.concat([
    deg[deg.sig == "up"].nlargest(3, "neglog10p"),
    deg[deg.sig == "down"].nlargest(3, "neglog10p"),
])
for _, r in top.iterrows():
    ax.annotate(r["gene"], (r["log2FC"], r["neglog10p"]),
                xytext=(5, 4), textcoords="offset points",  # 점에서 살짝 띄움
                fontsize=9, fontweight="bold", color="#333")
fig.savefig("volcano_labeled.png", dpi=300, bbox_inches="tight")

annotate의 xytext=(5, 4)와 textcoords="offset points"는 "점에서 오른쪽 5pt·위 4pt 떨어진 곳에 글자를 놓아라"는 뜻입니다. 라벨이 점과 겹쳐 안 보이는 일을 막는 흔한 방법입니다. 라벨이 많아 서로 겹친다면 adjustText 같은 보조 패키지를 쓰지만, 보통은 상위 몇 개만 다는 걸로 충분합니다.

3. 클러스터 히트맵 — clustermap으로 발현 패턴 군집

히트맵은 발현 행렬을 색의 농담으로 보여 줍니다. 그냥 값을 칠하면 발현이 높은 유전자 행만 진하게 나와 패턴이 안 보이므로, 유전자별로 z-score 정규화(각 행을 평균 0·표준편차 1로)한 뒤 칠하는 게 표준입니다. seaborn의 clustermap은 여기에 더해 비슷한 행·열끼리 계층 군집(hierarchical clustering)으로 묶어 덴드로그램(dendrogram)까지 그려 줍니다.

# 행(유전자) 단위 z-score: (값 − 행평균) / 행표준편차
expr_z = expr.sub(expr.mean(axis=1), axis=0).div(expr.std(axis=1), axis=0)

g = sns.clustermap(
    expr_z, cmap="vlag", center=0,          # vlag = 파랑–흰–빨강 발산 팔레트
    vmin=-2, vmax=2,                         # 색 범위 고정 (이상치에 안 휘둘림)
    col_cluster=True, row_cluster=True,      # 행·열 모두 군집
    figsize=(8, 7), linewidths=.4,
    cbar_kws={"label": "z-score"},
    dendrogram_ratio=0.14,
)
g.ax_heatmap.set_xlabel("샘플")
g.ax_heatmap.set_ylabel("유전자")
g.figure.suptitle("발현 z-score 클러스터 히트맵", y=1.02, fontweight="bold")
g.savefig("clustermap.png", dpi=300, bbox_inches="tight")

center=0은 0을 색의 중앙(흰색)에 두라는 뜻이라, 상향(빨강)·하향(파랑)이 직관적으로 갈립니다. vmin·vmax로 색 범위를 ±2에 고정하면 극단값 한두 개가 전체 색을 잡아먹는 일을 막습니다. clustermap은 일반 heatmap과 달리 Figure 전체를 새로 만들어 반환하므로(g.figure로 접근), plt.subplots로 만든 축에 끼워 넣을 수 없다는 점만 기억하면 됩니다. 군집을 끄고 단순 히트맵만 원하면 sns.heatmap(expr_z, ...)을 씁니다.

덴드로그램을 보면 처리군 샘플(T1~T5)과 대조군(C1~C5)이 좌우로 갈려 묶이는 게 보입니다. 우리가 신호를 넣은 그대로, 발현 패턴만으로 두 조건이 분리된 셈입니다.

4. 박스·바이올린 플롯 — 그룹별 발현 분포

관심 유전자 하나의 발현이 그룹마다 어떻게 다른지는 박스플롯이나 바이올린 플롯으로 봅니다. 이 그림들은 데이터가 long 형식(한 행에 샘플·그룹·발현값 하나씩)일 때 자연스럽습니다. 발현 행렬을 melt로 녹이고, 샘플 이름에서 그룹을 뽑아냅니다.

# wide(유전자×샘플) → long, 그리고 샘플명 첫 글자로 그룹 추출
long = (expr.reset_index(names="gene")
        .melt(id_vars="gene", var_name="sample", value_name="expr"))
long["group"] = long["sample"].str[0].map({"C": "control", "T": "treated"})
long.head(3)

     gene sample   expr    group
0  GENE01     C1   7.42  control
1  GENE02     C1   8.21  control
2  GENE03     C1   6.95  control

먼저 바이올린 플롯입니다. 박스플롯이 사분위수(quartile)만 보여 준다면, 바이올린은 분포의 모양(밀도)까지 보여 줍니다. 상향 유전자 하나(GENE01)를 골라 두 그룹을 비교합니다.

sub = long[long["gene"] == "GENE01"]         # 유전자 하나만

fig, ax = plt.subplots(figsize=(6, 5))
sns.violinplot(data=sub, x="group", y="expr",
               hue="group", palette={"control": "#6BAB97", "treated": "#E76F51"},
               legend=False, inner="box", ax=ax)   # inner="box": 안에 박스도
sns.stripplot(data=sub, x="group", y="expr",        # 실제 점도 겹쳐 찍기
              color="#333", size=5, alpha=.7, ax=ax)
ax.set_xlabel(""); ax.set_ylabel("정규화 발현값 (log$_2$)")
ax.set_title("GENE01 — 그룹별 발현 분포", fontweight="bold")
plt.tight_layout()
plt.savefig("violin.png", dpi=300, bbox_inches="tight")

여기서 seaborn 0.13의 바뀐 규칙이 처음 등장합니다 — x축과 같은 변수를 hue에도 넘기고(hue="group"), legend=False로 중복 범례를 끕니다. 예전(0.12 이전)에는 palette만 줘도 각 박스가 다른 색이었지만, 0.13부터는 hue 없이 팔레트만 주면 경고가 뜹니다(아래 두 번째 에러 참고). stripplot을 겹쳐 실제 데이터 점을 같이 찍으면, 샘플 수가 적을 때(여기선 그룹당 5개) 분포 모양이 과장돼 보이는 걸 정직하게 보완할 수 있습니다.

박스플롯도 호출만 바꾸면 됩니다. 여러 유전자를 한 축에 늘어놓을 때 특히 박스플롯이 깔끔합니다.

# 유전자 4개를 한 박스플롯에 (x=유전자, 색=그룹)
sel = long[long["gene"].isin(["GENE01", "GENE02", "GENE05", "GENE06"])]

fig, ax = plt.subplots(figsize=(8, 5))
sns.boxplot(data=sel, x="gene", y="expr", hue="group",
            palette={"control": "#6BAB97", "treated": "#E76F51"}, ax=ax)
ax.set_xlabel(""); ax.set_ylabel("정규화 발현값 (log$_2$)")
ax.set_title("유전자별·그룹별 발현 비교", fontweight="bold")
ax.legend(title="그룹")
plt.tight_layout()
plt.savefig("boxplot.png", dpi=300, bbox_inches="tight")

hue="group"을 주면 유전자마다 control·treated 박스가 나란히 두 개씩 그려집니다. 상향 유전자(GENE01·GENE02)는 treated 박스가 위로, 하향 유전자(GENE05·GENE06)는 아래로 내려간 게 한눈에 들어옵니다.

5. facet 다중 패널 — catplot으로 여러 유전자 한 번에

유전자가 여러 개면 한 그림에 욱여넣기보다 작은 패널을 격자로 나열하는 게 깔끔합니다. 이걸 facet (또는 small multiples)이라 부릅니다. seaborn의 catplot(범주형)·relplot(관계형)은 col·row에 열 이름만 주면 패널을 자동으로 쪼개 줍니다. catplot은 그림(figure) 수준 함수라 plt.subplots 없이 자체적으로 격자를 만듭니다.

# 유전자 6개를 패널 6개로 (가로 3 × 세로 2), 각 패널은 박스플롯
sel6 = long[long["gene"].isin(
    ["GENE01", "GENE02", "GENE03", "GENE05", "GENE06", "GENE08"])]

g = sns.catplot(
    data=sel6, x="group", y="expr",
    hue="group", palette={"control": "#6BAB97", "treated": "#E76F51"},
    col="gene", col_wrap=3, kind="box",     # col=패널 기준, col_wrap=한 줄에 3개
    height=2.8, aspect=0.9, legend=False,
)
g.set_titles("{col_name}")                   # 각 패널 제목 = 유전자명
g.set_axis_labels("", "발현값 (log$_2$)")
g.figure.suptitle("유전자별 발현 분포 (facet)", y=1.03, fontweight="bold")
g.savefig("facet.png", dpi=300, bbox_inches="tight")

col="gene"이 패널을 유전자별로 쪼개고, col_wrap=3이 세 개마다 줄을 바꿉니다. kind="box"를 "violin"·"bar"·"strip"으로 바꾸면 같은 격자에 다른 그림이 그려집니다. catplot이 반환하는 건 축이 아니라 FacetGrid 객체라, 제목·축 라벨은 g.set_titles·g.set_axis_labels처럼 격자 메서드로 한꺼번에 손봅니다. 패널이 너무 많아지면 height·aspect로 한 패널 크기를 줄여 전체가 화면에 들어오게 조절합니다.

6. 스타일·한글 폰트·저장 — 발표용 마무리

마지막은 실무 설정입니다. 발표 슬라이드나 논문 그림으로 쓰려면 한글이 깨지지 않아야 하고, 해상도와 형식을 챙겨야 합니다.

matplotlib은 기본 폰트에 한글 글리프가 없어 한글 라벨이 네모(□□□)로 깨집니다. rcParams로 한글 폰트를 한 번 지정하면 이후 모든 그림에 적용됩니다. 윈도우는 Malgun Gothic, macOS는 AppleGothic, 리눅스는 나눔고딕 등을 씁니다. 음수 부호가 깨지는 것도 함께 막아 줍니다.

import matplotlib.pyplot as plt

# 운영체제별 한글 폰트 (있는 것 하나 고름)
plt.rcParams["font.family"] = ["Malgun Gothic", "AppleGothic", "NanumGothic"]
plt.rcParams["axes.unicode_minus"] = False   # 음수 부호(−) 깨짐 방지
plt.rcParams["figure.dpi"] = 110             # 화면 미리보기 해상도
plt.rcParams["savefig.bbox"] = "tight"       # 저장 시 여백 자동 정리

seaborn의 전역 테마도 그림 전체 인상을 좌우합니다. set_theme의 style(배경·격자)과 context(글자·선 두께)를 조합합니다. 발표용은 글자를 키운 "talk"가, 논문용은 기본 "notebook"이 무난합니다.

# style: white / whitegrid / dark / darkgrid / ticks
# context: paper / notebook / talk / poster (오른쪽으로 갈수록 요소 큼)
sns.set_theme(style="ticks", context="talk", palette="muted")

저장은 savefig로 합니다. 비트맵인 PNG는 dpi로 해상도를 정하고(인쇄·발표는 300 이상), 벡터인 SVG·PDF는 dpi와 무관하게 무한 확대해도 안 깨져 논문 그림에 적합합니다.

fig.savefig("figure.png", dpi=300, bbox_inches="tight")   # 발표·웹: 고해상도 PNG
fig.savefig("figure.svg", bbox_inches="tight")            # 논문: 벡터 SVG
fig.savefig("figure.pdf", bbox_inches="tight")            # 논문: 벡터 PDF

bbox_inches="tight"는 그림 둘레의 불필요한 여백을 잘라 줍니다. 라벨이 그림 밖으로 잘려 나간다면 이 옵션이 거의 해결합니다. 색약(color vision deficiency)을 배려하려면 palette="colorblind"나 발산 팔레트 "vlag"·"icefire"를 쓰는 게 좋습니다.

자주 발생하는 에러 & 해결

• 한글 라벨이 네모(□□□)로 깨짐 — matplotlib 기본 폰트에 한글 글리프가 없어서입니다. 6절처럼 plt.rcParams["font.family"]에 설치된 한글 폰트를 지정하세요. 폰트 캐시가 옛 정보를 들고 있으면 한 번 지운 뒤(matplotlib.font_manager._load_fontmanager(try_read_cache=False)) 재시작하면 잡힙니다. 음수 부호까지 깨지면 axes.unicode_minus를 False로 둡니다.
• seaborn 0.13: palette without hue 경고 — sns.boxplot(data=df, x="group", y="v", palette="Set2")처럼 hue 없이 팔레트만 주면 FutureWarning이 뜨고 0.14에서 막힙니다. 0.13부터 범주형 함수는 hue가 없으면 단색으로 그리는 게 기본입니다. 해결은 x축 변수를 hue에도 넘기고(hue="group") legend=False로 중복 범례를 끄는 것입니다. 이 글의 모든 박스·바이올린 코드가 이 패턴을 씁니다.
• clustermap을 subplots 축에 못 넣음 — fig, ax = plt.subplots(); sns.clustermap(..., ax=ax)는 에러입니다. clustermap은 자체적으로 Figure 전체(히트맵·덴드로그램·컬러바)를 만들어 ClusterGrid를 반환하기 때문입니다. 기존 축에 넣고 싶으면 군집 없는 sns.heatmap(..., ax=ax)을 쓰고, 군집이 꼭 필요하면 clustermap을 독립 그림으로 두세요. 저장은 g.savefig(...)로 합니다.
• 그림이 안 뜨거나 빈 창만 나옴 — 스크립트(.py)에서는 plt.show()를 마지막에 호출해야 창이 뜹니다. 주피터에서는 셀의 마지막 객체가 자동 표시되지만, 같은 셀에서 그림을 그리고 또 다른 걸 그리면 앞 그림이 먹힐 수 있으니 셀을 나누거나 fig를 명시적으로 다루세요. 저장만 하고 화면 표시가 필요 없으면 matplotlib.use("Agg")로 백엔드를 비대화형으로 바꾸면 됩니다.
• 라벨·제목이 그림 밖으로 잘림 — bbox_inches="tight"를 savefig에 주거나, 그리기 전에 plt.tight_layout()을 호출하세요. facet (catplot)에서 suptitle이 패널과 겹치면 y=1.03처럼 위로 올리고, 그래도 잘리면 g.figure.subplots_adjust(top=0.9)로 위 여백을 줍니다.
• 점이 너무 많아 산점도가 뭉개짐(overplotting) — 볼케이노에 유전자가 수만 개면 점이 겹쳐 새까매집니다. s로 점을 작게(s=8), alpha로 반투명하게(alpha=0.4) 주거나, 비유의 점만 rasterized=True로 래스터화해 SVG 용량을 줄입니다. 유의 유전자만 색을 진하게, 나머지는 회색으로 깔면 가독성이 올라갑니다.

확장 — 더 해 볼 것

• seaborn.objects 인터페이스 — 0.12부터 들어온 새 문법(so.Plot(df, x=..., y=...).add(so.Dot()))으로, ggplot2처럼 레이어를 +가 아닌 .add()로 쌓습니다. 아직 실험적(experimental)이지만 복잡한 그림을 선언적으로 조립하기에 좋습니다.
• 상호작용 그림 — 발표용 정적 그림 너머로, plotly·bokeh를 쓰면 마우스를 올려 유전자 이름을 보는 인터랙티브 볼케이노를 만들 수 있습니다. 탐색 단계에서 유용합니다.
• 전용 패키지 — 차등발현 시각화에 특화된 pydeseq2·sanbomics의 볼케이노 헬퍼, 라벨 겹침을 자동으로 푸는 adjustText, 출판용 조판을 돕는 matplotlib의 constrained_layout 등을 함께 쓰면 손이 줄어듭니다.
• PCA·UMAP 산점도 — 샘플 전체 구조를 보는 차원축소 그림은 단일세포 분석에서 특히 중요합니다. scanpy 글에서 sc.pl.umap으로 이어집니다.
• 색·테마 다듬기 — 저널마다 요구하는 그림 규격(폰트 크기·선 두께·색)이 다릅니다. rcParams를 묶은 스타일 시트(plt.style.use("my.mplstyle"))를 하나 만들어 두면 모든 그림을 일관되게 찍을 수 있습니다.

References

1. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9(3), 90–95. (matplotlib 원논문)
2. Waskom, M. L. (2021). seaborn: statistical data visualization. Journal of Open Source Software, 6(60), 3021. (seaborn 원논문)
3. The Matplotlib development team. (2026). Matplotlib documentation — v3.11. matplotlib.org/stable
4. Waskom, M. seaborn 0.13 release notes — 범주형 함수 재작성·hue/팔레트 규칙 변경·native_scale. whatsnew/v0.13.0
5. seaborn. seaborn.clustermap — 계층 군집 히트맵 API. clustermap docs

Pipette & Pipeline · A bio portfolio journal