챕터 25 본편: 단일세포 전사체 분석

단일세포 전사체 분석 실습

개요

이 장에서는 단일세포 RNA 시퀀싱(scRNA-seq) 데이터의 분석 파이프라인을 파이썬으로 직접 구현하는 실습을 다룬다. 단일세포 전사체학의 이론적 배경, 기술적 원리, 세포 아틀라스 프로젝트에 대한 내용은 10장 단일 세포 전사체학을 참조한다. 차원 축소와 클러스터링 알고리즘의 수학적 원리에 대한 상세한 내용은 11장 차원 축소와 데이터 분석을 참조한다.

분석 파이프라인 개요

단일세포 전사체 데이터 분석은 다음과 같은 단계로 구성된다:

1. 정량화 및 데이터 생성: Cell x Gene matrix 생성

2. 전처리: 배치 효과 교정, 정규화, 로그 변환

3. 차원 축소: PCA를 통한 초기 차원 축소

4. 그래프 구성: k-최근접 이웃(KNN) 그래프 생성

5. 클러스터링: Leiden 알고리즘을 이용한 세포 유형 분류

6. 시각화: UMAP 임베딩을 통한 2차원 시각화

7. 차등 발현 분석: 클러스터 간 마커 유전자 발견

실습 환경 구성

작업 디렉토리 생성

$ mkdir -p ~/week7
$ cd ~/week7

UV 가상환경 설정

UV를 사용하여 파이썬 가상환경을 생성하고 필요한 패키지를 설치한다.

$ uv venv --python 3.13
$ source .venv/bin/activate
$ uv pip install numpy pandas scikit-learn seaborn matplotlib
$ uv pip install umap-learn leidenalg igraph
$ uv pip install ipykernel

Jupyter 커널 등록

$ python -m ipykernel install --user --name week7 --display-name "week7"

데이터 시뮬레이션

실제 scRNA-seq 데이터를 사용하기 전에, 데이터의 특성을 이해하기 위해 시뮬레이션 데이터를 생성한다.

기본 패키지 및 상수 설정

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import nbinom
 
np.random.seed(42)
 
n_celltypes = 20
n_cells_per_celltype = 100
n_genes = 500
n_total_cells = n_celltypes * n_cells_per_celltype

Negative Binomial 기반 카운트 생성

scRNA-seq 데이터는 음이항 분포(Negative Binomial distribution)를 따르는 것으로 알려져 있다. 이를 시뮬레이션하는 함수를 정의한다.

def simulate_counts(mu0, mu1):
    mu = np.random.uniform(mu0, mu1)
    r = np.random.uniform(2, 4)
    p = r / (r + mu)
    return nbinom.rvs(r, p, size=n_cells_per_celltype)

세포 유형별 발현 데이터 생성

각 세포 유형에 대해 마커 유전자와 배경 유전자를 구분하여 발현량을 생성한다.

total_counts = []
for i in range(n_celltypes):
    n_marker_genes = np.random.randint(5, 50)
    marker_gene_indices = np.random.choice(n_genes, n_marker_genes, replace=False)
    counts = np.zeros((n_cells_per_celltype, n_genes), dtype=int)
    for gene_idx in range(n_genes):
        if gene_idx in marker_gene_indices:
            counts[:, gene_idx] = simulate_counts(10, 1000)
        else:
            counts[:, gene_idx] = simulate_counts(1, 5)
    total_counts.append(counts)
 
mat = np.vstack(total_counts)
np.random.shuffle(mat)
 
df = pd.DataFrame(mat,
    columns=[f"Gene_{i}" for i in range(n_genes)],
    index=[f"Cell_{i}" for i in range(n_total_cells)]
)

데이터 전처리

Mean-Variance 관계 확인

scRNA-seq 데이터의 특징적인 mean-variance 관계를 확인한다.

gene_means = df.mean(axis=0)
gene_vars = df.var(axis=0)
 
sns.scatterplot(x=gene_means, y=gene_vars, s=2)
plt.xlabel('Mean Expression')
plt.ylabel('Variance of Expression')
plt.ylim(0, None)

라이브러리 크기 정규화

세포마다 시퀀싱 깊이가 다르므로, 라이브러리 크기로 정규화한다.

cell_summed = df.sum(axis=1).values
df_normalized = df / cell_summed.reshape(-1, 1) * np.median(cell_summed)

로그 변환

정규화된 데이터에 로그 변환을 적용하여 분산을 안정화한다.

df_log_normalized = np.log1p(df_normalized)

로그 변환 전후의 유전자 발현 분포를 비교한다.

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sns.histplot(df['Gene_100'], bins=50, ax=axes[0])
axes[0].set_title('Raw counts')
sns.histplot(df_log_normalized['Gene_100'], bins=50, ax=axes[1])
axes[1].set_title('Log-normalized')

차원 축소

PCA 수행

scikit-learn의 PCA를 사용하여 차원을 축소한다.

from sklearn.decomposition import PCA
 
pca = PCA(n_components=50)
pca_result = pca.fit_transform(df_log_normalized)
 
sns.scatterplot(x=pca_result[:, 0], y=pca_result[:, 1], s=1)
plt.xlabel('PC1')
plt.ylabel('PC2')

KNN 그래프 생성

UMAP 패키지의 nearest_neighbors 함수를 사용하여 KNN 그래프를 생성한다.

from umap.umap_ import nearest_neighbors
 
knn = nearest_neighbors(pca_result,
                        n_neighbors=15,
                        metric="euclidean",
                        metric_kwds=None,
                        angular=False,
                        random_state=None)

UMAP 임베딩

KNN 그래프를 기반으로 UMAP 임베딩을 수행한다.

from umap import UMAP
 
umap_model = UMAP(n_components=2, min_dist=0.5, spread=1.0, precomputed_knn=knn)
umap_result = umap_model.fit_transform(pca_result)
 
sns.scatterplot(x=umap_result[:, 0], y=umap_result[:, 1], s=1)
plt.xlabel('UMAP1')
plt.ylabel('UMAP2')

클러스터링

Leiden 클러스터링을 위한 그래프 구성

KNN 인덱스를 사용하여 igraph 그래프 객체를 생성한다.

import leidenalg
import igraph as ig
 
indices = knn[0]
 
edges = []
for i in range(indices.shape[0]):
    for j in range(1, indices.shape[1]):
        edges.append((i, indices[i, j]))
 
g = ig.Graph(edges=edges, directed=False)
g.simplify()

Leiden 클러스터링 수행

Leiden 알고리즘을 사용하여 세포를 클러스터링한다.

partition = leidenalg.find_partition(
    g,
    leidenalg.RBConfigurationVertexPartition,
    resolution_parameter=1.0
)
 
labels = np.array(partition.membership)

클러스터링 결과 시각화

UMAP 좌표에 클러스터 라벨을 색상으로 표시한다.

sns.scatterplot(
    x=umap_result[:, 0], y=umap_result[:, 1], s=1,
    hue=labels, palette='tab20', legend=None
)
plt.xlabel('UMAP1')
plt.ylabel('UMAP2')

클러스터별 Heatmap

클러스터링 결과로 세포를 정렬하여 발현 패턴을 확인한다.

sorted_indices = np.argsort(labels)
df_sorted = df_log_normalized.iloc[sorted_indices, :]
 
sns.heatmap(df_sorted, cmap='viridis', cbar=False)

차등 발현 분석

Wilcoxon Rank-Sum 검정

두 클러스터 간의 차등 발현 유전자를 찾기 위해 Wilcoxon rank-sum 검정을 수행한다.

from scipy.stats import ranksums
 
cluster_0_indices = np.where(labels == 0)[0]
cluster_1_indices = np.where(labels == 1)[0]
 
deg_genes = []
p_values = []
adjusted_p_values = []
log2fc = []
 
for gene in df_log_normalized.columns:
    stat, p_value = ranksums(
        df_log_normalized.iloc[cluster_0_indices][gene],
        df_log_normalized.iloc[cluster_1_indices][gene]
    )
    adjusted_p_value = p_value * n_genes  # Bonferroni correction
    if adjusted_p_value < 0.01:
        deg_genes.append(gene)
        p_values.append(p_value)
        adjusted_p_values.append(adjusted_p_value)
        log2fc.append(
            np.log2(df_log_normalized.iloc[cluster_0_indices][gene].mean() + 1) -
            np.log2(df_log_normalized.iloc[cluster_1_indices][gene].mean() + 1)
        )

Volcano Plot

차등 발현 분석 결과를 volcano plot으로 시각화한다.

sns.scatterplot(x=log2fc, y=-np.log10(adjusted_p_values), s=5)
plt.xlabel('Log2 Fold Change')
plt.ylabel('-Log10 Adjusted P-Value')
plt.axhline(y=2, color='red', linestyle='--', alpha=0.5)
plt.axvline(x=-1, color='red', linestyle='--', alpha=0.5)
plt.axvline(x=1, color='red', linestyle='--', alpha=0.5)

상위 DEG Heatmap

상위 차등 발현 유전자의 발현 패턴을 heatmap으로 확인한다.

top_deg_indices = np.argsort(np.abs(log2fc))[-50:]
top_deg_genes = [deg_genes[i] for i in top_deg_indices]
 
df_top_degs = df_sorted[top_deg_genes].iloc[:200]
sns.heatmap(df_top_degs, cmap='viridis', cbar=True)

Scanpy를 이용한 분석

Scanpy는 단일세포 데이터 분석을 위한 파이썬 패키지로, 앞서 다룬 모든 분석 단계를 통합적으로 제공한다.

Scanpy 설치

$ uv pip install scanpy

데이터 로드

h5ad 형식의 데이터를 로드한다.

import scanpy as sc
 
adata = sc.read_h5ad("data.h5ad")

품질 관리

# 미토콘드리아 유전자 비율 계산
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
 
# QC plot
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4)

전처리 및 차원 축소

# 필터링
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
 
# 정규화
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
 
# 고변이 유전자 선택
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]
 
# 스케일링
sc.pp.scale(adata, max_value=10)
 
# PCA
sc.tl.pca(adata, svd_solver='arpack')

클러스터링 및 시각화

# 이웃 그래프 구성
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
 
# UMAP
sc.tl.umap(adata)
 
# Leiden 클러스터링
sc.tl.leiden(adata, resolution=1.5)
 
# 시각화
sc.pl.umap(adata, color=['leiden'])

마커 유전자 발견

# 클러스터별 마커 유전자 찾기
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
 
# Dot plot
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)

실습 과제

실습 25.1: 시뮬레이션 데이터 분석

1. 본문에서 제시된 코드를 사용하여 시뮬레이션 데이터를 생성한다.

2. 정규화, PCA, UMAP, Leiden 클러스터링을 순차적으로 수행한다.

3. resolution 파라미터를 0.5, 1.0, 2.0으로 변경하면서 클러스터 수의 변화를 관찰한다.

실습 25.2: Scanpy 분석 파이프라인

1. UV로 새로운 가상환경을 생성하고 Scanpy를 설치한다.

2. 제공된 h5ad 파일을 로드한다.

• 데이터 경로: /bce/lectures/2025-bioinformatics/data/scrnaseq/brain_small.h5ad

1. Scanpy 튜토리얼을 참조하여 다음을 수행한다:

• 기본 QC plot 및 고변이 유전자(HVG) 선택 plot

• Leiden 클러스터링 결과(resolution 1.5)로 색이 표현된 PCA 및 UMAP

• 각 클러스터별 DEG 5개씩 포함된 dot plot

참고 튜토리얼: https://scverse-tutorials.readthedocs.io/en/latest/notebooks/basic-scrna-tutorial.html

출처

이 본편 글은 원본 docx에 기록된 아래 출처를 기준으로 게시했습니다.

• 원본 문서: https://chaek.org/books/introduction-to-biomedical-informatics
• 원본 저장소: https://github.com/chaek-union/introduction-to-biomedical-informatics
• 원본 파일: practices/25-scrnaseq-analysis.md

문제 풀이

Chapter 25. 단일세포 전사체 분석

0 / 24

Gemini AI 채점

주관식 답안은 Gemini API로 채점합니다. API 키는 이 브라우저에만 저장됩니다.

API KEY 미등록

1. 1. [쉬움] 객관식

   단일세포 전사체 데이터 분석 파이프라인의 단계로 옳지 않은 것을 고르시오.

   선택지

   Cell x Gene matrix 생성 PCA를 통한 초기 차원 축소 Leiden 알고리즘을 이용한 클러스터링 ANNOVAR를 이용한 변이 주석
2. 2. [쉬움] 객관식

   라이브러리 크기 정규화 코드의 목적으로 옳은 것을 고르시오.

   cell_summed = df.sum(axis=1).values
   df_normalized = df / cell_summed.reshape(-1, 1) * np.median(cell_summed)
   

   선택지

   세포마다 다른 시퀀싱 깊이를 보정한다. 모든 유전자를 삭제한다. Leiden 클러스터를 UMAP 좌표로 바꾼다. FASTQ 파일을 BAM 파일로 변환한다.
3. 3. [쉬움] 객관식

   이 장에서 UV 가상환경을 생성하는 명령어로 옳은 것을 고르시오.

   선택지

   `source .venv/bin/activate` `uv venv --python 3.13` `uv pip install numpy pandas matplotlib seaborn` `python -m ipykernel install --user --name week7 --display-name "week7"`
4. 4. [쉬움] 객관식

   Jupyter 커널을 week7 이름으로 등록하는 명령어로 옳은 것을 고르시오.

   선택지

   `uv venv --python 3.13` `source .venv/bin/activate` `python -m ipykernel install --user --name week7 --display-name "week7"` `uv pip install numpy pandas matplotlib seaborn`
5. 5. [보통] 객관식

   본문의 데이터 시뮬레이션에서 설정한 총 세포 수로 옳은 것을 고르시오.

   n_celltypes = 20
   n_cells_per_celltype = 100
   n_total_cells = n_celltypes * n_cells_per_celltype
   

   선택지

   2000 500 1000 100
6. 6. [보통] 객관식

   본문의 시뮬레이션 설정에서 n_genes = 500이고 총 세포 수가 2000이면 df의 기본 형태로 옳은 것은 무엇인가?

   선택지

   500행 × 2000열 100행 × 20열 20행 × 100열 2000행 × 500열
7. 7. [보통] 객관식

   PCA 수행 코드에서 n_components=50의 의미로 옳은 것을 고르시오.

   pca = PCA(n_components=50)
   pca_result = pca.fit_transform(df_log_normalized)
   

   선택지

   50개의 클러스터를 강제로 생성한다. 50개의 주성분으로 차원을 축소한다. 50개의 유전자를 삭제한다. 50개의 FASTQ 파일을 병합한다.
8. 8. [보통] 객관식

   UMAP 임베딩 코드에서 precomputed_knn=knn의 의미로 옳은 것을 고르시오.

   umap_model = UMAP(n_components=2, min_dist=0.5, spread=1.0, precomputed_knn=knn)
   

   선택지

   UMAP 결과를 CSV 파일로 저장한다. Leiden 알고리즘을 비활성화한다. 앞에서 계산한 KNN 정보를 UMAP에 사용한다. 모든 세포 유형을 하나로 합친다.
9. 9. [보통] 객관식

   마커 유전자 탐색에서 본문 코드가 사용하는 기준으로 옳은 것을 고르시오.

   선택지

   각 유전자의 FASTQ 품질 점수를 비교한다. 각 클러스터의 평균 발현을 나머지 세포의 평균 발현과 비교한다. 각 클러스터의 BAM 파일 크기를 비교한다. 각 세포의 염색체 수를 비교한다.
10. 10. [보통] 객관식

    정규화 계산에서 어떤 세포의 총 카운트가 1000, 전체 세포 총 카운트의 중앙값이 500, 특정 유전자의 원래 카운트가 20이면 정규화 값으로 옳은 것은 무엇인가?

    정규화 값 = 원래 카운트 / 세포 총 카운트 × 중앙값
    

    선택지

    40 1000 500 10
11. 11. [어려움] 객관식

    KNN 결과 indices를 이용하여 edges를 만들 때 for j in range(1, indices.shape[1])처럼 1부터 시작하는 이유로 가장 적절한 것은 무엇인가?

    선택지

    첫 번째 이웃은 자기 자신일 수 있어 자기 연결을 피하기 위해서이다. 첫 번째 유전자는 항상 삭제해야 하기 때문이다. 첫 번째 클러스터는 UMAP에서 사용하지 않기 때문이다. 첫 번째 열은 FASTQ 품질 점수이기 때문이다.
12. 12. [어려움] 객관식

    Leiden 클러스터링 수행 코드에서 사용하는 객체와 함수의 연결로 옳은 것을 고르시오.

    선택지

    `DeseqDataSet`과 `dds.fit()` `table_annovar.pl`과 `-vcfinput` `ig.Graph`와 `leidenalg.find_partition` `spatialdata_io.xenium`과 `sdata.tables`
13. 13. [어려움] 객관식

    Wilcoxon Rank-Sum 검정 코드에서 비교하는 대상으로 옳은 것은 무엇인가?

    선택지

    두 이미지 파일의 픽셀 크기 두 FASTQ 파일의 어댑터 서열 두 VCF 파일의 FILTER 상태 두 클러스터의 유전자별 발현값
14. 14. [어려움] 객관식

    Bonferroni correction에서 p_value = 0.00001, n_genes = 500일 때 adjusted_p_value로 옳은 것은 무엇인가?

    선택지

    0.00001 0.005 0.5 0.05
15. 15. [어려움] 객관식

    Scanpy 전처리 및 클러스터링 흐름으로 옳은 것을 고르시오.

    선택지

    `annovar` → `octopus` → `filter_cells` → `pca` `leiden` → `umap` → `filter_cells` → `normalize_total` → `log1p` `filter_cells/filter_genes` → `normalize_total` → `log1p` → `highly_variable_genes` → `scale` → `pca` → `neighbors` → `umap` → `leiden` `rank_genes_groups_dotplot` → `read_h5ad` → `filter_genes` → `normalize_total`
16. 16. [어려움] 객관식

    Scanpy에서 클러스터별 마커 유전자를 찾고 dot plot을 그리는 코드 조합으로 옳은 것을 고르시오.

    선택지

    `sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')` 후 `sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)` `sq.gr.ripley(adata)` 후 `sq.pl.ripley(adata)` `trainer.save_checkpoint()` 후 `load_from_checkpoint()` `annovar/table_annovar.pl` 후 `less myanno.hg19_multianno.vcf`
17. 주관식 1. [쉬움] 주관식 · Gemini 채점

    ~/week7 디렉토리를 만들고 해당 디렉토리로 이동하는 명령어를 작성하시오.
18. 주관식 2. [보통] 주관식 · Gemini 채점

    다음 조건을 만족하는 Python 코드를 작성하시오.

    • df의 각 유전자 평균을 gene_means에 저장한다.
    • df의 각 유전자 분산을 gene_vars에 저장한다.
    • sns.scatterplot으로 평균-분산 관계를 그린다.
19. 주관식 3. [보통] 주관식 · Gemini 채점

    df_normalized에 로그 변환을 적용하여 df_log_normalized에 저장하는 코드를 작성하시오.
20. 주관식 4. [보통] 주관식 · Gemini 채점

    simulate_counts(mu0, mu1) 함수의 핵심 계산 흐름을 작성하시오. mu, r, p = r / (r + mu), nbinom.rvs를 포함할 것.
21. 주관식 5. [보통] 주관식 · Gemini 채점

    KNN 결과 indices = knn[0]를 이용하여 edges 리스트를 만들고 igraph 그래프를 생성하는 코드를 작성하시오.
22. 주관식 6. [어려움] 주관식 · Gemini 채점

    Leiden 클러스터링을 수행하고 labels 배열을 생성하는 코드를 작성하시오.
23. 주관식 7. [어려움] 주관식 · Gemini 채점

    Wilcoxon Rank-Sum 검정에서 cluster 0과 cluster 1의 특정 유전자 발현을 비교하는 코드의 핵심 부분을 작성하시오.
24. 주관식 8. [어려움] 주관식 · Gemini 채점

    Scanpy를 사용하여 data.h5ad를 읽고 필터링, 정규화, 로그 변환, HVG 선택, PCA, 이웃 그래프, UMAP, Leiden 클러스터링까지 수행하는 코드 흐름을 작성하시오.