Scanpy tutorials - preprocessing and clustering ( Re-assess quality control and cell filtering ,Differentially-expressed Genes as Markers) tutorial 코드 설명

Re-assess quality control and cell filtering

sc.pl.umap(
    adata,
    color=["leiden", "predicted_doublet", "doublet_score"],
    # increase horizontal space between panels
    wspace=0.5,
    size=3,
)

 

 

sc.pl.umap() : UMAP 을 통해 차원 축소된 데이터를 시각화 한다.

color=["leiden", "predicted_doublet", "doublet_score"] : 각각의 UMAP 플롯에서 색상을 지정하는 기준이다. 첫번째 그래프는 클러스터별로 색깔을 입힌 것이다. doublet score 를 보며 threshold 를 조절해볼 수 있고 클러스터가 더 명확하게 구분되도록 조절할 수 있다.

# 'predicted_doublet'가 False인 세포들만 남기기
adata_filtered = adata[adata.obs["predicted_doublet"] == False]

 

doublet ture 세포는 제거하고, adata_filtered 객체를 사용하여 UMAP 을 다시 시각화하면 이중세포가 제거된 상태로 결과를 확인할 수 있다.

sc.pl.umap(
    adata,
    color=["leiden", "log1p_total_counts", "pct_counts_mt", "log1p_n_genes_by_counts"],
    wspace=0.5,
    ncols=2,
)

 

UMAP (Uniform Manifold Approximation and Projection) - GPT 설명

UMAP은 비선형 차원 축소 기법으로, 고차원 데이터의 구조를 저차원(주로 2D)에서 효과적으로 표합니다. (PCA는 선형 차원 축소)

PCA 후 UMAP: 고차원 데이터를 UMAP에 바로 적용할 수도 있지만, 데이터가 매우 고차원(수천 개 이상의 변수)일 경우, PCA로 먼저 차원을 줄인 후 UMAP을 적용하는 것이 일반적입니다. 이렇게 하면 계산 비용이 줄어들고, 노이즈가 줄어들어 더 안정적인 결과를 얻을 수 있습니다.

 

 

for res in [0.02, 0.5, 2.0]:
    sc.tl.leiden(
        adata, key_added=f"leiden_res_{res:4.2f}", resolution=res, flavor="igraph"
    )

 

데이터의 클러스터 구조를 다양한 수준에서 확인하고 생물학적인 의미를 고려해 최적의 해상도를 찾는다.

sc.tl.leiden() : Scanpy의 Leiden 클러스터링을 여러 해상도 설정으로 반복 수행한다.

flavor="igraph": igraph 라이브러리를 사용하여 Leiden 알고리즘을 실행한다.

 

sc.pl.umap(
    adata,
    color=["leiden_res_0.02", "leiden_res_0.50", "leiden_res_2.00"],
    legend_loc="on data",
)

 세포를 클러스터로 그룹화한 후에 각 클러스터의 세포 유형을 특정하기 위해 마커 유전자(marker genes) 찾는다. 자동으로 찾을 수 있는데 튜토리얼에 따로 안나와 있으니 그 과정은 다음 포스트에서 정리할 예정이다. 각 resolution 별로 분석을 진행해서 내가 보고자하는 cell type 이 제일 두드러지고, 알려진 reference와 비교해서 납득이되게끔 만들면 된다.

 

*UMAP은 고차원 데이터를 저차원으로 시각화하는 방식으로, 동일한 샘플이라도 시각화 모양 매번 다를 수 있다. 내가한 분석과 reference 와 비교할 때는 주요 세포 유형이 비슷하게 구별되는지를 보면 된다.

marker_genes = {
    "CD14+ Mono": ["FCN1", "CD14"],
    "CD16+ Mono": ["TCF7L2", "FCGR3A", "LYN"],
    # Note: DMXL2 should be negative
    "cDC2": ["CST3", "COTL1", "LYZ", "DMXL2", "CLEC10A", "FCER1A"],
    "Erythroblast": ["MKI67", "HBA1", "HBB"],
    # Note HBM and GYPA are negative markers
    "Proerythroblast": ["CDK6", "SYNGR1", "HBM", "GYPA"],
    "NK": ["GNLY", "NKG7", "CD247", "FCER1G", "TYROBP", "KLRG1", "FCGR3A"],
    "ILC": ["ID2", "PLCG2", "GNLY", "SYNE1"],
    "Naive CD20+ B": ["MS4A1", "IL4R", "IGHD", "FCRL1", "IGHM"],
    # Note IGHD and IGHM are negative markers
    "B cells": [
        "MS4A1",
        "ITGB1",
        "COL4A4",
        "PRDM1",
        "IRF4",
        "PAX5",
        "BCL11A",
        "BLK",
        "IGHD",
        "IGHM",
    ],
    "Plasma cells": ["MZB1", "HSP90B1", "FNDC3B", "PRDM1", "IGKC", "JCHAIN"],
    # Note PAX5 is a negative marker
    "Plasmablast": ["XBP1", "PRDM1", "PAX5"],
    "CD4+ T": ["CD4", "IL7R", "TRBC2"],
    "CD8+ T": ["CD8A", "CD8B", "GZMK", "GZMA", "CCL5", "GZMB", "GZMH", "GZMA"],
    "T naive": ["LEF1", "CCR7", "TCF7"],
    "pDC": ["GZMB", "IL3RA", "COBLL1", "TCF4"],
}

 

Marker gene을 찾아서 딕셔너리 형태로 저장한다. 후에 어떻게 자동화해서 찾는지 알아보자.

sc.pl.dotplot(adata, marker_genes, groupby="leiden_res_0.02", standard_scale="var")

 

마커 유전자 발현 패턴을 시각화해서 각 클러스터의 세포 유형을 확인하는 과정이다. sc.pl.dotplot() 함수는 각 클러스터에서 마커 유전자의 발현 수준을 직관적으로 보여준다.

 

해상도 설정을 0.02로 해서 5개 클러스터(0~4) 로 출력됨. 해상도를 높이면 클러스터 숫자가 늘어난다. 클러스터 조정 방식은 Leiden 클러스터링 알고리즘이 데이터 내에서 밀접한 관계를 가지는 세포들을 더 큰 그룹으로 묶도록 재구성한다. 

지금 나눠진 모양을 봤을때, 클러스터가 깔끔하게 나눠져 보이지 않는다. 1번 클러스터의 경우 CD14+ Mono, Naive CD20, B cells 등 여러 세포의 마커가 겹쳐보이기 때문에 해상도를 높여 클러스터 수를 늘릴 필요가 있다. 아니면, "넓은 범주의 세포 계통(lineage)"으로 클러스터링 한다면 의미있는 데이터로 해석할 수 있다.

adata.obs["cell_type_lvl1"] = adata.obs["leiden_res_0.02"].map(
    {
        "0": "Lymphocytes",
        "1": "Monocytes",
        "2": "Erythroid",
        "3": "B Cells",
    }
)

 

adata.obs 에 세포계통으로 labeling 한 컬럼을 추가한다. 

sc.pl.dotplot(adata, marker_genes, groupby="leiden_res_0.50", standard_scale="var")

 

해상도를 0.5로 올려 더 세밀하게 클러스팅한다. 아까보다 클러스터수가 늘었고 세포를 좀 더 고수준으로 분리할 수 있다.

standard_scale = 'var' 은 각 유전자의 발현값을 0에서 1사이의 범위로 표준화 한다. 유전자 간 발현 수준을 비교할 때 특정 유전자가 지나치게 큰 값으로 표현되지 않도록 조정한다. 

 

Differentially-expressed Genes as Markers

# Obtain cluster-specific differentially expressed genes
sc.tl.rank_genes_groups(adata, groupby="leiden_res_0.50", method="wilcoxon")

 

• sc.tl.rank_genes_groups 함수는 각 클러스터에서 차별적으로 발현되는 유전자들(differentially-expressed genes)을 찾기 위한 함수이다.
• Wilcoxon 또는 t-test와 같은 통계 검정을 사용하여, 특정 클러스터에서 다른 클러스터에 비해 유의미하게 발현이 높은 유전자들을 식별한다.
• 이 과정을 통해 **클러스터별로 특이적으로 발현되는 유전자(마커 유전자)**를 찾을 수 있으며, 이 마커 유전자를 바탕으로 세포 유형을 유추한다.

sc.pl.rank_genes_groups_dotplot(
    adata, groupby="leiden_res_0.50", standard_scale="var", n_genes=5
)

 

• n_genes=5: 각 클러스터에서 상위 5개의 차별 발현 유전자를 시각화한다.
• dot plot을 통해 각 클러스터에서 발현양 상위 5개 유전자를 다른 클러스터와 시각적으로 비교할 수 있다.

sc.get.rank_genes_groups_df(adata, group="7").head(5)

 

7번 클러스터의 발현양 상위5개 유전자를 df으로 반환한다. names 는 유전자 이름, scores는 통계적 검정(t-test, Wilcoxon test)에서 계산된 값으로 높은 점수일수록 해당 유전자가 해당 그룹에서 더 강하게 발현한다. logfoldchanges 는 7번 클러스터와 나머지 그룹 간 발현량의 비율 차이를 로그로 변환한 것

* 위의 경우에는 MT 유전자가 DEG로 나왔는데, 원래는 MT 유전자들은 제거하고 해야한다.

dc_cluster_genes = sc.get.rank_genes_groups_df(adata, group="7").head(5)["names"]
sc.pl.umap(
    adata,
    color=[*dc_cluster_genes, "leiden_res_0.50"],
    legend_loc="on data",
    frameon=False,
    ncols=3,
)

 

sc.get.rank_genes_groups_df(adata, group="7").head(5)["names"]: adata 객체에서 그룹 "7"에 대한 상위 5개의 차등 발현 유전자 이름을 추출하여 dc_cluster_genes라는 변수에 저장한다. 

 

 

sc.pl.umap 함수는 UMAP을 사용하여 각 유전자의 발현을 시각화한다. 클러스터 7번에서 DEG를 찾은 것이기 때문에 마지막 그래프의 7번 위치를 다른 그래프들에 겹쳐 보면, 발현양이 높은 노란색으로 표시되는걸 확인할 수 있다. 이 5개 유전자가 7번 클러스터의 특성을 나타낸다. 어떤 세포 유형(cell type)에 해당하는지 특정하는데 사용된다.