Seurat 라이브러리 치트시트(Cheatsheet)

Tip

scRNA-seq

Bioinformatics

Seurat

Cheatsheet

Author

Taeyoon Kim

Published

November 12, 2024

Modified

November 12, 2024

Seurat 은 scRNA-seq 분석에서 가장 널리 사용되는 R 패키지로, 다양한 기능을 통해 데이터를 처리하고 시각화합니다. 이 문서는 주요 명령어와 기능을 요약하여 사용자가 데이터를 효율적으로 다룰 수 있도록 지원합니다. 특히, 데이터 로드, 전처리, 클러스터링 및 시각화와 같은 일반적인 작업에 대한 예제를 제공합니다. Seurat v5.0.1 버전을 사용한 환경에서 테스트되었습니다. 다양한 명령어와 튜토리얼은 공식 페이지 에서 확인하세요.

1 파일 불러오기

scRNA-seq 정보를 저장할 수 있는 파일에는 여러 형식이 있으며, 모든 형식을 여기에 모두 나열할 수는 없습니다. 다만 보편적인 형식은 아래와 같습니다.

SMARTseq2 에서는 행렬과 샘플 메타데이터를 csv, tsv 또는 txt 형식으로 저장합니다.
10X Chromium 데이터의 기본 형식 중 하나는 .tsv 로 저장된 유전자 및 세포에 대한 주석이 포함된 압축된 희소 행렬 파일 .mtx 를 사용합니다.
현재는 10X Chromium 에서 많이 사용되는 형식은 HDF5 파일입니다. HDF5 는 속도가 빠르고 확장 가능하며 한 번에 사용할 데이터의 일부만을 불러오는 기능도 제공합니다. 또한 메타데이터를 동일한 파일에 저장하며 데이터를 이진 압축된 희소 행렬 형식으로 저장합니다.

# csv, tsv, txt 등의 형식 읽기
raw_matrix <- read.delim(
    file = "data/folder_sample1.csv",
    row.names = 1
)

# 희소 행렬로 변환
sparse_matrix <- Matrix::Matrix(
    data   = raw_matrix,
    sparse = TRUE
)

# mtx 형식에서 읽고 변환하기 (폴더에 있는 파일들을 읽음)
sparse_matrix <- Seurat::Read10X(
    data.dir = "data/folder_sample1"
)

# h5 형식에서 읽고 변환하기
sparse_matrix <- Seurat::Read10X_h5(
    filename = "data/matrix_file.h5",
    use.names = TRUE
)

메타 데이터는 데이터 프레임 형식으로 가져와야 합니다.

# csv 또는 tsv 형식에서 메타데이터 읽기
metadata <- read.delim(
    file = "data/metadata.csv",
    sep = ",",
    row.names = 1
)

2 Seurat 객체 만들기

개발자들은 데이터 분석 프로세스를 조금 더 간단하게 만들기 위해 데이터를 간결한 형식으로 저장하는 객체를 만들었습니다. Seurat 에서는 CreateSeuratObject 함수를 사용해 희소 행렬 또는 카운트가 포함된 데이터를 입력으로 객체를 생성합니다.

SeuratObject <- CreateSeuratObject(
    counts = sparse_matrix,
    assay = "RNA",
    project = "SAMPLE1",
    meta.data = metadata
)

Seurat 객체는 객체 이름 뒤에 @ 또는 $ 문자를 사용하여 내용에 쉽게 접근할 수 있습니다.

@ 기호는 assays, meta.data, graphs 및 reduction 슬롯을 포함한 모든 분석 결과에 접근할 수 있습니다.
$ 기호를 사용하면 Seurat 객체의 메타데이터 열에 직접 액세스할 수 있습니다. 즉, SeuratObject$column1 은 SeuratObject@meta.data$column1 과 동일한 결과를 출력합니다.

기본적으로 분석 데이터는 RNA 이라는 이름의 assay 슬롯에 저장됩니다. 분석을 수행할 때 추가 assays 가 생성되며 예를 들어 유전자의 원시 데이터 (counts), 정규화한 데이터 (data), 스케일링/회귀된 데이터 (scale.data) 및 유전자의 분산에 관한 정보 (var) 들이 있습니다. Seurat 객체의 모든 슬롯에 대한 자세한 설명은 Seurat 위키 를 참고하세요.

2.1 메타데이터에 새로운 열 추가

새로운 열을 추가하려면 $ 기호를 사용하여 벡터를 메타데이터 열에 할당할 수 있습니다.

SeuratObject$NEW_COLUMN_NAME <- setNames(
    colnames(SeuratObject),
    vector_NEW_DATA
)

또는 AddMetaData 함수를 사용하여 열을 추가할 수 있습니다.

SeuratObject <- AddMetaData(
    object   = SeuratObject,
    metadata = vector_NEW_DATA,
    col.name = "NEW_COLUMN_NAME"
)

2.2 Seurat 객체 하위 집합 만들기

일부 데이터를 따로 때어내는 방법은 R 의 데이터프레임 문법을 사용하는 방법과 Seurat 함수를 사용하는 방법이 있습니다. 먼저 데이터 프레임 문법은 아래와 같습니다.

# Feature 하위 집합
SeuratObject[vector_FEATURES_TO_USE, ]
# 세포 하위 집합
SeuratObject[, vector_CELLS_TO_USE]

Seurat 에서 제공하는 함수는 아래와 같습니다.

# Identity 클래스를 기반으로 Seurat 객체를 하위 집합으로 추출합니다. 자세한 내용은 ?SubsetData를 참조하세요.
subset(x = SeuratObject, idents = "Bcell")
subset(x = SeuratObject, idents = c("CD4T", "CD8T"), invert = TRUE)

# 유전자/특성의 발현 수준을 기준으로 하위 집합을 만듭니다.
subset(x = SeuratObject, subset = MS4A1 > 3)

# 기준의 조합에 따라 하위 집합을 만듭니다.
subset(x = SeuratObject, subset = MS4A1 > 3 & PC1 > 5)
subset(x = SeuratObject, subset = MS4A1 > 3, idents = "B 세포")

# 객체 메타 데이터의 값에 따라 하위 집합을 만듭니다.
subset(x = SeuratObject, subset = orig.ident == "Replicate1")

# 각 identity 클래스당 세포 수를 다운샘플링합니다.
subset(x = SeuratObject, downsample = 100)

2.3 시각화하기

가장 일반적으로 사용되는 시각화 함수는 바이올린 플롯과 UMAP(차원 축소 후의 산점도) 입니다.

Warning

UMAP 을 그리기 전에 반드시 먼저 차원 축소를 계산해야 합니다.

연속 변수에 대하여 바이올린 플롯 시각화하는 코드

VlnPlot(
    object = SeuratObject,
    group.by = "orig.ident",
    features = c("percent_mito"),
    pt.size = 0.1,
    ncol = 4,
    y.max = 100
) + NoLegend()

연속 변수에 대한 UMAP 을 그리는 방법 (예시: 유전자 발현량):

FeaturePlot(
    object = SeuratObject,
    features = c("FEATURE_1", "FEATURE_2", "FEATURE_3"),
    reduction = "umap",
    dims = c(1, 2),
    order = TRUE,
    pt.size = 0.1,
    ncol = 3
)

범주형 변수에 대한 UMAP 을 그리는 방법 (예시: 세포 유형):

DimPlot(
    object = SeuratObject,
    group.by = c("DATASET"),
    reduction = "umap",
    dims = c(1, 2),
    pt.size = 0.1,
    label = TRUE,
    ncol = 3
)

이것 말고도 다양한 시각화 함수가 있으니 자세한 것은 공식 문서 를 참고하세요.

2.4 데이터 합치기

두 개이상의 데이터셋을 결합하려면:

merged_obj <- merge(
    x = SeuratObject1,
    y = c(SeuratObject2, SeuratObject3, SeuratObject4),
    add.cell.ids = c("Dataset1", "Dataset2", "Dataset3", "Dataset4")
)

3 전체 Feature 의 수 확인하기

일반적인 방식은 적어도 200 개의 유전자가 있는 세포와 유전자는 적어도 3 개 이상의 세포에서 발현된 것만 남기는 것입니다. 다만 이 수치는 데이터를 준비할 때 사용된 라이브러리 기술에 매우 의존적입니다. 따라서 주의해서 사용해야 합니다. Seurat 에서는 nCount_RNA: 세포 당 UMIs(Unique Molecular Identifiers) 의 수; nFeature_RNA: 각 세포에서 감지된 유전자의 수로 저장하고 있습니다. 시각화를 통해 유전자 수와 UMI 의 갯수를 살펴보는 코드는 아래와 같습니다.

VlnPlot(
    SeuratObject,
    group.by = "orig.ident",
    features = c("nFeature_RNA", "nCount_RNA"),
    pt.size = 0.1
) + NoLegend()

4 유전자 관련 품질 관리

단일 세포에서 가장 많이 감지되는 유전자는 주로 미토콘드리아 유전자 (MT-), 리보솜 유전자 (RPL 및 RPS), 그리고 기타 항상 발현되는 구조 단백질들 (예: ACTB, TMSB4X, B2M, EEF1A1) 입니다. 이러한 유전자 데이터는 세포 품질 지표 계산, 데이터 정규화, 배치 효과 보정 등 단일 세포 데이터 분석의 여러 목적에 활용될 수 있습니다.¹

4.1 미토콘드리아 유전자 비율

높은 비율의 미토콘드리아 유전자는 일반적으로 손상된 세포를 의미합니다. 이는 세포질 속의 RNA 가 쉽게 소실되는 반면, 미토콘드리아 속의 RNA 는 비교적 덜 손실되기 때문입니다. 다음은 미토콘드리아 유전자 비율을 계산하고 메타데이터 테이블에 추가하는 예시입니다.

# 미토콘드리아 유전자의 비율 계산
SeuratObject <- PercentageFeatureSet(
    object = SeuratObject,
    pattern = "^MT",
    assay = "RNA",
    col.name = "percent_mito"
)

Note

유전자 관련 품질 관리에 대한 설명은 사람 (human) 에 해당하는 내용으로 다른 종의 데이터를 처리하는 경우에는 코드가 수정되어야 합니다. 예시로 생쥐 (mouse) 의 경우 미토콘드리아 (mt-), 리보소말 (Rpl) 입니다.

4.2 리보솜 유전자 비율

리보솜 유전자는 모든 세포에서 가장 많이 발현되는 유전자 중 하나이며, 미토콘드리아 유전자와는 반대로 미토콘드리아 수에 반비례하는 경향이 있습니다. 즉, 미토콘드리아의 수가 많을수록 리보솜 유전자의 감지 수준은 낮아집니다 (다만, 이 관계는 비선형적입니다). 리보솜 단백질에서 유래하는 유전자 발현 비율을 계산하는 방법은 미토콘드리아 유전자 비율을 계산하는 방법과 유사합니다.

# 리보솜 유전자의 비율 계산
SeuratObject <- PercentageFeatureSet(
    SeuratObject,
    pattern  = "^RP[SL]",
    col.name = "percent_ribo"
)

4.3 유전자 유형 및 염색체 위치 정보

RNA 시퀀싱에서 얻은 유전자는 여러 유형으로 구분할 수 있습니다. 이에는 단백질 발현 정보를 포함하는 Coding 유전자, 단백질 정보를 포함하지 않는 Non-coding 유전자, BCR 및 TCR 과 관련된 VDJ 영역 유전자 (세포 분화에 대한 정보를 포함), 그리고 siRNA 등이 포함됩니다. 이러한 정보를 활용하면 원하는 분석 유형에 따라 관심 없는 특정 유전자 범주를 필터링할 수 있습니다.

특히, scRNA-seq 는 일반적으로 poly-A 를 사용하여 RNA 를 농축하기 때문에 데이터에 존재하는 유전자의 약 80-90% 가 단백질 코딩 유전자입니다. 또한, 유전자의 염색체 위치 정보를 수집하면 성 염색체에 의한 효과를 식별하는 데 유용할 수 있습니다.

# 사람 유전자로 변경하기 위해 먼저 ENSEMBL에서 마우스 유전자 주석 검색
library(biomaRt)

# ENSEMBL 마트 설정
mart <- biomaRt::useMart(
    biomart = "ensembl",
    dataset = "hsapiens_gene_ensembl",
    host= "www.ensembl.org"
)

# 선택된 속성으로 마우스 유전자 주석 검색
annot <- biomaRt::getBM(
    mart = mart,
    attributes = c("external_gene_name", "gene_biotype", "chromosome_name")
)

유전자 이름을 해당 생물 유형과 일치시킵니다.

# 유전자 이름과 해당 염색체 위치 일치
item <- annot[match(rownames(SeuratObject), annot[, 1]), "chromosome_name"]

# NA 값을 "unknown"으로 대체
item[is.na(item)] <- "unknown"

# 유효하지 않은 염색체 값을 "other"로 대체
item[!item %in% as.character(c(1:23, "X", "Y", "MT"))] <- "other"

예를 들어 단백질 코딩 유전자에만 초점을 맞추려면 다음과 같이 수행 할 수 있습니다.

# SeuratObject의 차원 확인
dim(SeuratObject)

# 유전자 주석에서 단백질 코딩 유전자 선택
sel <- annot[match(rownames(SeuratObject), annot[, 1]), 2] == "protein_coding"

# 선택된 유전자 이름 가져오기
genes_use <- rownames(SeuratObject)[sel]
genes_use <- as.character(na.omit(genes_use))

# SeuratObject에서 단백질 코딩 유전자만 필터링
SeuratObject <- SeuratObject[genes_use, ]

# 필터링 후 SeuratObject의 차원 확인
dim(SeuratObject)

5 세포 주기에 대한 정보

각각의 세포가 어떤 주기에 해당하는지는 특정 유전자의 평균 발현을 계산하여 예측합니다. Seurat 은 CellCycleScoring 함수를 통해 이것은 쉽게 계산할 수 있습니다 (사람 유전자에만 해당). 실행 방법은 다음과 같습니다.

# CellCycleScoring를 실행하기 전에 데이터를 정규화하고 로그 변환해야 합니다.
SeuratObject <- NormalizeData(SeuratObject)

# Cell Cycle 점수 계산
SeuratObject <- CellCycleScoring(
    object = SeuratObject,
    g2m.features = cc.genes$g2m.genes,
    s.features = cc.genes$s.genes
)

6 PCA 플롯에 메타데이터 시각화

메타데이터를 PCA 플롯에서 시각화하는 것은 각각의 메타데이터 변수를 개별적으로 분석하는 것보다 쉽습니다. 모든 정보를 하나의 플롯으로 결합하기 위해 주성분 분석 (PCA) 을 사용할 것입니다. 상위 주성분 (PC1 과 PC2) 에 대한 시각화는 데이터 세트 간의 차이가 어떻게 반영되는지 직관적으로 보여줍니다.

# 선택된 메타데이터 매개 변수를 사용하여 PCA 계산
metadata_use <- grep(
    "perc", colnames(SeuratObject@meta.data), value = TRUE
)

metadata_use <- c(
    "nCount_RNA", "nFeature_RNA", "S.Score", "G2M.Score",
    metadata_use
)

# PCA 실행
PC <- prcomp(
    SeuratObject@meta.data[, metadata_use],
    center = TRUE,
    scale. = TRUE
)

# 객체에 PCA(메타데이터에서 실행됨) 추가
SeuratObject@reductions[["pca_metadata"]] <- CreateDimReducObject(
    embeddings = PC$x,
    key = "metadataPC_",
    assay = "RNA"
)

# 메타데이터에서 실행된 PCA 플롯
DimPlot(
    SeuratObject,
    reduction = "pca_metadata",
    dims = c(1, 2),
    group.by = "orig.ident"
)

Note

제로 카운트 유전자와 모든 값이 0 인 데이터 열은 PCA 와 같은 분석에 영향을 줄 수 있어서 이런 데이터를 제외하기 위해 아래 코드를 사용합니다.

# 제로 카운트 유전자 제거
SeuratObject <- SeuratObject[Matrix::rowSums(SeuratObject) > 0, ]

# 모든 0인 메타데이터 열 제외
metadata_use <- metadata_use[
    colSums(SeuratObject@meta.data[, metadata_use] != 0) > 0
]

7 데이터 정규화 및 회귀

추가 분석을 수행하기 전에 서로 다른 시퀀싱 깊이를 고려하여 데이터를 정규화해야 합니다. 일반적으로 로그 변환을 통해 데이터를 정규화하며, 세포 주기, 미토콘드리아 비율, 감지된 유전자 수와 같은 데이터 품질 관리를 위한 메타데이터는 제거하는 것이 좋습니다.

7.1 데이터 정규화

scRNA-seq 에서 가장 일반적인 정규화 방법은 로그 정규화입니다. 이 방법은 각 유전자의 카운트를 모든 유전자 카운트의 합 (라이브러리 크기라고도 함) 으로 나누어 라이브러리 크기의 차이를 보상합니다. 그 후, 결과에 상수를 곱하여 모든 세포가 동일한 시퀀싱 깊이를 가지도록 합니다. 대부분의 bulk RNA-seq 에서는 상수가 보통 1e6 이며, 이로 인해 CPM(백만 단위 카운트) 이 생성됩니다. 그러나 단일 세포 라이브러리 크기는 그보다 훨씬 작기 때문에 1e3 에서 1e4(10,000 단위 카운트) 까지 사용됩니다.

\[ NormCounts = \frac{(GeneCounts \times 10000)}{LibrarySize} \]

라이브러리 크기에 대해 보정된 값은 로그 변환하여 로그 정규 분포로 만듭니다.

\[ logNormCounts = \ln(NormCounts + 1) \]

Seurat에서는 NormalizeData 함수를 사용하여 이 과정을 수행할 수 있습니다.

# 분산이 0인 유전자 제거
SeuratObject <- SeuratObject[Matrix::rowSums(SeuratObject) > 0, ]

# 데이터 정규화 (LogNormalize 사용)
SeuratObject <- NormalizeData(
    object = SeuratObject,
    scale.factor = 10000,
    normalization.method = "LogNormalize"
)

7.2 데이터 회귀 (Regression)

데이터에서 혼란 요인을 제거하기 위해 회귀 분석을 수행할 수 있습니다. 예를 들어, 세포 주기나 품질 지표 (예: 미토콘드리아 비율 또는 감지된 유전자 수) 와 같은 혼란 요인을 제거할 수 있습니다. Seurat에서는 ScaleData 함수를 사용하여 이러한 회귀를 수행할 수 있습니다.

SeuratObject <- ScaleData(
    SeuratObject,
    vars.to.regress = c("nFeature_RNA", "percent_mito"),
    verbose = FALSE
)

이렇게 하면 데이터에 회귀가 적용되며, 이제 데이터가 준비되었으므로 추가 분석을 진행할 수 있습니다.

8 highly variable features 선택하기

대규모 데이터 분석에서 중요한 단계는 샘플 간에 현저한 차이를 보이는 features(유전자, 전사체, 단백질, 대사물질 등) 를 식별하는 것입니다. 예를 들어, 다양한 종류의 상피 세포가 포함된 데이터셋을 가정해 보겠습니다. 상피 세포의 하위 유형을 구별하기 위해 다음과 같은 전략을 사용할 수 있습니다:

모든 상피 세포에서 공통적으로 발현되며 수준이 거의 동일한 유전자만 사용합니다.
상피 세포에서 감지되지 않는 유전자만 사용합니다.
상피 세포 간에 발현 수준이 크게 다른 유전자만 사용합니다.
데이터셋의 모든 유전자를 사용합니다.

상피 세포 간 발현이 크게 다른 유전자를 사용하는 것이 가장 효과적인 방법이며, 그 다음으로 모든 유전자를 사용하는 것이 좋습니다. 반면, 상피 세포 간에 균일하게 발현되는 유전자나 상피 세포에서 감지되지 않는 유전자를 사용하는 것은 하위 유형을 구별하는 데 충분한 정보를 제공하지 않을 수 있습니다.그러나 단일 세포 분석에서는 보통 세포의 상피 세포 유형을 사전에 알지 못합니다 (그것이 발견하려는 대상이기 때문입니다). 따라서 이 작업을 수행하기 위해 다른 방법이 필요합니다. 일반적인 접근법은 각 샘플 전체에서 유전자의 분산을 기준으로 유전자를 순위 지정하는 것입니다. 분산이 높은 유전자는 세포를 구분하는 데 효과적일 가능성이 높습니다.일반적으로 발현 수준이 높은 유전자는 자연스럽게 더 높은 변동성을 나타내므로, 각 유전자의 평균 로그 발현량으로 변동성을 정규화하는 것이 중요합니다.

# Variable feature 찾기
SeuratObject <- FindVariableFeatures(
    object = SeuratObject,
    nfeatures = 3000,
    selection.method = "vst",
    verbose = FALSE,
    assay = "RNA"
)

# 상위 20개의 Variable feature 시각화
top20 <- head(VariableFeatures(SeuratObject), 20)

LabelPoints(
    plot = VariableFeaturePlot(SeuratObject),
    points = top20,
    repel = TRUE
)

9 스케일링 및 중심화

9.1 (선형) 스케일링 및 중심화

각 유전자는 서로 다른 발현 수준을 가지므로, 높은 발현 값을 가진 유전자는 자연스럽게 더 높은 변동성을 보입니다. 이는 후속 분석에서 잘 포착될 수 있습니다. 따라서 각 유전자에 대해 유사한 가중치를 부여하는 것이 중요합니다. 일반적으로 PCA 를 수행하기 전에 각 유전자를 중심화하고 스케일링하는 것이 권장됩니다. 이러한 정밀한 스케일링 방법은 Z- 점수 정규화라고 하며, PCA, 클러스터링 및 히트맵 시각화에 매우 유용합니다.

또한, 세포 주기, 시퀀싱 깊이, 미토콘드리아 비율 등과 같은 데이터셋에서 원하지 않는 변동 소스를 제거하기 위해 회귀 분석을 사용할 수 있습니다. 이는 이러한 매개변수를 모형의 공변량으로 사용하여 일반화된 선형 회귀 (GLM) 를 수행함으로써 달성됩니다. 이후 모델의 잔차는 “회귀된 데이터”로 간주됩니다.

SeuratObject <- ScaleData(
    object = SeuratObject,
    vars.to.regress = c("nCount_RNA", "percent_mito", "nFeatures_RNA"),
    model.use = "linear",
    assay = "RNA",
    do.scale = TRUE,
    do.center = TRUE
)

Seurat 은 기본적으로 Seurat 객체에 있는 모든 변수 기능에 대해 스케일링을 실행합니다. 변수 기능이 있을 경우 이러한 스케일링 단계의 성능을 크게 향상시킬 수 있습니다. 그러나 결과는 변수 유전자 자체에 기반하여 PCA, 차원 축소 또는 클러스터링과 같은 후속 분석의 결과가 변경되지 않습니다. 전체 유전자 수에 대해 데이터를 스케일링하려는 경우 FindVariableFeatures 를 사용한 후 함수 호출에서 이를 지정할 수 있습니다.

SeuratObject <- ScaleData(
    object = SeuratObject,
    vars.to.regress = c("nCount_RNA", "percent_mito", "nFeatures_RNA"),
    model.use  = "linear",
    assay = "RNA",
    do.scale = TRUE,
    do.center = TRUE,
    features = rownames(SeuratObject)
)

9.2 (포아송) 스케일링 및 중심화

위의 선형 스케일링과 중심화는 로그 - 선형 데이터 분포를 가정하기 때문에, RNA-seq 데이터 (단일 세포 포함) 는 보다 음이항 분포에 가까울 수 있으며, 이로 인해 변이를 올바르게 회귀하지 못할 수 있습니다. 위의 절차의 대안으로 “포아송” 또는 “음이항” 분포를 사용하여 원시 UMI 카운트 데이터에 대해 실행할 수 있습니다. 이는 포아송 모델을 사용하여 유전자별 GLM 회귀를 수행하는 방법입니다.

SeuratObject <- ScaleData(
    object = SeuratObject,
    vars.to.regress = c("nCount_RNA", "mito.percent", "nFeatures_RNA"),
    model.use = "poisson",
    assay = "RNA",
    do.scale = TRUE,
    do.center = TRUE
)

10 SCTransform

포아송 분포를 기반으로 한 스케일링과 중심화는 경우에 따라 데이터를 과적합시킬 수 있습니다. 이러한 문제를 극복하기 위해 유사한 풍부도를 가진 유전자 간의 정보를 통합하여 회귀 모델에서 안정된 매개변수 추정을 얻을 수 있습니다. 이 방법을 scTransform 이라고 하며, 간단히 말해 제한된 음이항 모델을 사용하여 유전자별 일반 선형 모델 회귀를 수행하는 방식입니다.

SeuratObject <- SCTransform(
    object = SeuratObject,
    assay = "RNA",
    vars.to.regress = "percent_mito",
    new.assay.name = "sctransform",
    do.center = TRUE,
    verbose = FALSE
)

11 차원 축소하기

scRNA-seq 데이터는 수천 개의 유전자로 구성된 전체 유전자 발현 공간을 포함하고 있으며, 이 공간은 상당한 노이즈를 포함하고 있어 시각화가 어렵습니다. 따라서 대부분의 scRNA-seq 분석에서는 데이터의 일부 변동을 제거하기 위해 PCA(주성분 분석) 또는 ICA(독립 성분 분석) 와 같은 방법으로 시작합니다.단순한 scRNA-seq 데이터셋에서는 몇 가지 세포 유형만 포함되어 있는 경우 PCA 가 데이터를 2 차원 또는 3 차원으로 시각화하는 데 충분할 수 있습니다. 그러나 데이터의 복잡성이 증가함에 따라 비선형 차원 축소 방법을 사용하여 데이터를 2 차원으로 투영해야 할 필요가 있습니다. 이러한 비선형 차원 축소 방법에는 t-SNE, UMAP, 확산 맵 등이 있습니다.어떤 차원 축소 방법을 선택할지는 간단한 문제가 아니며, 여러 가지 적절한 답변이 있을 수 있는 복잡한 문제입니다. 더 자세한 내용을 알고 싶다면 다음 논문 을 참고하시기 바랍니다.

11.1 PCA

주성분 분석 (PCA) 은 데이터를 새로운 좌표계로 변환하는 직교 선형 변환입니다. 이 과정에서 데이터의 가장 큰 분산이 첫 번째 주성분에 위치하고, 두 번째로 큰 분산은 두 번째 주성분에 위치하도록 합니다. 이러한 과정을 통해 데이터의 분산을 가장 잘 설명하는 방식으로 내부 구조를 드러냅니다. 일반적으로 이 동작은 데이터의 차원을 줄여 변환된 데이터의 차원을 축소하는 방식으로 이해될 수 있습니다.

# 주성분 분석 실행
SeuratObject <- RunPCA(
    object = SeuratObject,
    assay = "RNA",
    npcs = 100,
    verbose = FALSE
)

11.2 t-SNE

t- 분포 확률적 이웃 임베딩 (t-distributed stochastic neighbor embedding, t-SNE) 은 고차원 데이터를 저차원 공간으로 임베딩하여 시각화하는 데 적합한 비선형 차원 축소 기술입니다. 이 방법은 각 고차원 객체를 유사한 객체와 가까운 점으로 모델링하고, 서로 다른 객체는 먼 확률적 점으로 모델링하여 두 개 또는 세 개의 차원으로 표현합니다.²

# tSNE 실행
SeuratObject <- RunTSNE(
    object = SeuratObject,
    reduction = "pca",
    perplexity = 30,
    max_iter = 1000,
    theta = 0.5,
    eta = 200,
    exaggeration_factor = 12,
    dims.use = 1:50,
    verbose = FALSE,
    num_threads = 0
)

11.3 UMAP

UMAP(Uniform Manifold Approximation and Projection) 은 t-SNE 와 유사하게 시각화를 위한 차원 축소에 사용할 수 있는 기술이며, 일반적인 비선형 차원 축소에도 활용될 수 있습니다.³

UMAP 알고리즘은 실제 데이터에 적용할 수 있는 실용적이고 확장 가능한 방법입니다. 이 알고리즘은 시각화 품질에서 t-SNE 와 경쟁할 만큼 우수하며, 더 나은 실행 시간 성능을 바탕으로 전역 구조를 더 잘 보존한다고 알려져 있습니다. 또한 UMAP 은 임베딩 차원에 대한 계산 제한이 없어 머신러닝을 위한 일반적인 목적의 차원 축소 기술로 널리 사용될 수 있습니다.⁴

SeuratObject <- RunUMAP(
    object = SeuratObject,
    reduction = "pca",
    dims = 1:50,
    n.components = 2,
    n.neighbors = 20,
    repulsion.strength = 1,
    verbose = FALSE,
    n.epochs = 200,
    metric = "euclidean",
    seed.use = 42,
    reduction.name = "umap"
)

11.4 확산 맵 (Diffusion Maps)

확산 맵 (Diffusion Maps, DM) 은 데이터 집합을 유클리드 공간 (일반적으로 저차원) 으로 임베딩하는 차원 축소 기법 중 하나입니다. 이 기법은 데이터에 대한 확산 연산자의 고유 벡터와 고유값을 기반으로 임베딩을 계산합니다. 임베딩된 공간에서의 점들 사이의 유클리드 거리는 해당 점들을 중심으로 하는 확률 분포 간의 “확산 거리”와 동일합니다. 이 방법은 주어진 데이터에서 랜덤 워크를 수행할 때, 인접한 데이터 포인트로 이동하는 것이 멀리 떨어진 데이터 포인트로 이동하는 것보다 더 자주 발생한다는 기본적인 관찰에 기반하여 작동합니다.확산 맵은 주로 데이터가 추출된 기저 매니폴드를 발견하는 데 중점을 두고 있으며, 이는 주성분 분석 (PCA) 이나 다차원 스케일링 (MDS) 과 같은 선형 차원 축소 방법과는 다른 접근 방식입니다.⁵

# 추가 라이브러리 로드
library(destiny)

# destiny 패키지를 사용하여 확산 맵 실행
dm <- DiffusionMap(
    data = SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:50],
    k = 20,
    n_eigs = 20)

# DM 임베딩에서 셀 이름 수정
rownames(dm@eigenvectors) <- colnames(SeuratObject)

# SeuratObject에 DM 임베딩 추가
SeuratObject@reductions[["dm"]] <- CreateDimReducObject(
    embeddings = dm@eigenvectors,                                
    key = "DC_",
    assay = "RNA"
)

11.5 독립 성분 분석 (ICA)

독립 성분 분석 (Independent Component Analysis, ICA) 은 다변량 신호를 가산성 하위 구성 요소로 분리하는 계산 방법입니다. 이 기법은 하위 구성 요소가 가우시안이 아닌 신호이며 서로 통계적으로 독립적이라는 가정을 기반으로 합니다. ICA 는 블라인드 소스 분리 (blind source separation) 의 특별한 경우로 간주됩니다.⁶

SeuratObject <- RunICA(
    object = SeuratObject,
    assay = "pca",
    nics = 20,
    reduction.name = "ica"
)

12 그래프 만들기

scRNA-seq 데이터의 클러스터링을 전체 발현 데이터 혹은 PCA 공간 (선형 거리를 제공) 에서 하는 대신 그래프 (셀 간 유사성에 대한 비선형 표현) 를 먼저 만들고 클러스터링을 하는 것이 매우 효과적입니다. 그래프는 모든 셀 (노드 또는 정점) 을 나타내며, 일정한 유사성 기준에 따라 이들 사이에 에지를 연결합니다. 예를 들어, 50 개의 주요 구성 요소를 사용하여 PCA 공간에서 서로 X 거리 이하인 모든 셀 사이에 에지를 갖는 그래프를 구성할 수 있습니다.

12.1 KNN

KNN 은 “K Nearest Neighbors”의 약자로, 데이터 마이닝 및 기계 학습 분야에서 기본적이고 인기 있는 주제입니다. KNN 그래프는 두 정점 p 와 q 가 연결되는 그래프로, 그 거리가 K 번째로 작은 거리 중 하나인 경우에만 연결됩니다. 이때 쌍별 거리 측정치는 해밍 거리, 코사인 거리, 유클리드 거리 등 다양한 방법을 사용할 수 있습니다. 본 논문에서는 벡터 간 유사성을 측정하기 위해 유클리드 거리를 사용합니다. KNN 그래프 데이터 구조는 데이터 마이닝에서 많은 이점을 제공합니다. 예를 들어, 수십억 개의 데이터 세트에서는 KNN 그래프를 오프라인으로 사전 구축하여 인덱스로 사용하는 것이 여러 번의 온라인 KNN 검색보다 훨씬 효율적입니다.

SeuratObject <- FindNeighbors(
    SeuratObject,
    assay = "RNA",
    compute.SNN = FALSE, # FALSE로 설정하면 KNN 그래프만 계산됩니다.
    reduction = "pca",
    dims = 1:50,
    graph.name = "nn",
    prune.SNN = 1/15,
    k.param = 20,
    force.recalc = TRUE
)

KNN 그래프는 모든 셀 간의 연결이 1 로 표시된 행렬로 나타납니다. 이를 가중치가 없는 그래프라고 하며, Seurat 의 기본 설정입니다. 그러나 일부 셀 간의 연결은 다른 연결보다 더 중요할 수 있으며, 이 경우 그래프의 가중치는 0 에서 최대 거리까지 다양해집니다. 일반적으로 거리가 작을수록 두 점이 가까워지고, 그들의 연결이 강해집니다. 이러한 경우를 가중치가 있는 그래프라고 합니다.

가중치가 있는 그래프와 가중치가 없는 그래프 모두 클러스터링에 적합하지만, 대규모 데이터 세트 (100,000 개 이상의 세포) 에서는 가중치가 없는 그래프에서 클러스터링하는 것이 더 빠릅니다.

library(pheatmap)

pheatmap(
    SeuratObject@graphs$nn[1:100, 1:100],
    col = c("white", "black"),
    border_color = "grey90",
    legend = FALSE,
    cluster_rows = FALSE,
    cluster_cols = FALSE,
    fontsize = 2
)

12.2 SNN

추가로 KNN 그래프 외에도 모든 두 점 간에 공유된 최근접 이웃의 수를 계산합니다. 이를 통해 “공유된 최근접 이웃” 그래프를 형성합니다. 이 과정에서는 두 점 간의 각 링크의 가중치를 해당 점들이 공유하는 이웃의 수로 대체합니다. 다시 말해, 이는 최근접 이웃 그래프 내에서 모든 두 점 사이의 길이가 2 인 경로의 수를 나타냅니다.공유된 최근접 이웃 그래프가 생성되면 모든 점 쌍이 비교됩니다. 두 점이 T 이상의 이웃을 공유하는 경우, 즉 공유된 최근접 이웃 그래프에서 가중치가 임계값 T 보다 많은 링크를 가지면, 두 점과 그들이 속한 모든 클러스터가 병합됩니다. 다시 말해, 클러스터는 임계값을 사용하여 희소화한 후 공유된 최근접 이웃 그래프에서의 연결 요소입니다.

SeuratObject <- FindNeighbors(
    SeuratObject,
    assay = "RNA",
    compute.SNN = TRUE, # KNN 그래프와 SNN 그래프 둘 다 계산합니다.
    reduction = "pca",
    dims = 1:50,
    graph.name = "SNN",
    prune.SNN = 1/15,
    k.param = 20,
    force.recalc = TRUE
)

13 여러 데이터셋 합치기 및 배치 보정하기

기존의 일괄 보정 방법은 주로 대량 RNA-seq 을 위해 설계되었습니다. 따라서 이러한 방법을 scRNA-seq 데이터에 적용할 때는 각 일괄 처리 내의 세포 모집단 구성이 동일하다고 가정합니다. 일괄 간의 평균 유전자 발현 사이에 시스템적인 차이가 있을 경우, 이는 기술적인 차이로 귀속됩니다. 그러나 실제로는 scRNA-seq 연구에서 일괄 처리별로 모집단 구성이 동일하지 않을 수 있습니다. 각 일괄 처리에 동일한 세포 유형이 존재한다고 하더라도, 세포 유형의 풍부도는 세포 배양, 조직 추출, 분리 및 분류 등의 미세한 차이에 따라 달라질 수 있습니다. 따라서 일괄 차단 요인에 대한 추정 계수는 순수하게 기술적인 것이 아니라 구성의 차이로 인한 비생물학적 요소를 포함하게 됩니다. 이러한 계수를 기반으로 한 일괄 보정은 세포의 발현 프로파일을 부정확하게 표현할 수 있으며, 보정을 수행하지 않은 경우보다 결과가 더 나쁠 수 있습니다.⁷

많은 다른 scRNA-seq 방법론과 마찬가지로, 데이터셋에 대해 어떤 일괄 보정 기술을 사용해야 할지 선택하기 어려울 수 있습니다.⁸

13.1 Scanorama 사용하기

Note

최근 여러 방법들이 데이터셋을 성공적으로 통합할 수 있다는 것을 보여주고 있습니다. 그러나 이러한 접근법들은 데이터셋에 적어도 하나의 공통 세포 유형이 존재한다고 가정하기 때문에, 세포 구성이 상당히 다른 데이터셋을 결합할 경우 과도한 보정이 발생할 수 있습니다.

Scanorama 의 접근 방식은 두 데이터 세트 간의 유사한 요소를 찾는 상호 최근접 이웃 일치 기술을 일반화하여 여러 데이터 세트 간에 유사한 요소를 탐색합니다. 이 방식은 휴대폰으로 파노라마 사진을 촬영하는 기술과 유사합니다. Scanorama 는 자동으로 유사한 전사 프로필을 가진 세포를 식별하고 이를 일괄 보정 및 통합에 활용합니다. 이 방법은 강력한 특징을 가지고 있으며, 데이터셋의 크기나 출처에 크게 영향을 받지 않으며, 모든 데이터셋이 적어도 하나의 공통 세포 모집단을 공유할 필요가 없습니다.⁹

library(reticulate)
scanorama <- import("scanorama")

# 오브젝트를 배치 별로 분할
SeuratObject.list <- SplitObject(
    object   = SeuratObject,
    split.by = "BATCH"
)

# 어레이 및 제네릭 데이터 만들기
assaylist <- list()
genelist  <- list()

for (i in seq_along(SeuratObject.list)) {
    assaylist[[i]] <- t(
        as.matrix(GetAssayData(SeuratObject.list[[i]], "data"))
    )
    genelist[[i]] <- rownames(SeuratObject.list[[i]])
}

# 데이터 통합 및 보정
integrated.data <- scanorama$integrate(assaylist, genelist)
corrected.data <- scanorama$correct(
    assaylist,
    genelist,
    return_dense = TRUE
)
integrated.corrected.data <- scanorama$correct(
    assaylist,
    genelist,
    return_dimred = TRUE,
    return_dense = TRUE
)

# 데이터 결합 및 이름 지정
intdata <- lapply(integrated.corrected.data[[2]], t)
panorama <- do.call(cbind, intdata)
rownames(panorama) <- as.character(integrated.corrected.data[[3]])
colnames(panorama) <- unlist(sapply(assaylist, rownames))

# 차원 축소 결과 결합
intdimred <- do.call(rbind, integrated.corrected.data[[1]])
colnames(intdimred) <- paste0("PC_", 1:100)

# Seurat에서 Elbow Plot을 그리기 위해 표준 편차 추가
stdevs <- apply(intdimred, MARGIN = 2, FUN = sd)

# 새로운 어레이 오브젝트 만들기
SeuratObject@assays[["pano"]] <- CreateAssayObject(
    data = panorama,
    min.cells = 0,
    min.features = 0
)

# PCA 결과 추가
SeuratObject[["pca_scanorama"]] <- CreateDimReducObject(
    embeddings = intdimred,
    stdev = stdevs,
    key = "PC_",
    assay = "pano"
)

13.2 Harmony 사용하기

Harmony 는 주어진 저차원 임베딩 (예: 주성분 분석) 을 기반으로 시작합니다. Harmony 는 먼저 이 임베딩을 사용하여 셀을 여러 데이터 세트 클러스터로 그룹화합니다. 부드러운 전이를 고려하기 위해 소프트 클러스터링을 사용하여 셀을 잠재적으로 여러 클러스터에 할당합니다. 이러한 클러스터는 이산적인 세포 유형을 식별하기 위한 것이 아니라, 세포 상태 간의 부드러운 전환을 반영하는 대리 변수로 기능합니다. 클러스터링이 완료되면 각 데이터셋에는 클러스터별 중심점이 생성됩니다. 이 중심점은 클러스터별 선형 보정 계수를 계산하는 데 사용됩니다. 클러스터가 세포 유형과 상태에 해당하므로, 클러스터별 보정 계수는 개별 세포 유형 및 세포 상태에 특정한 보정 계수로 작용합니다. 이를 통해 Harmony 는 미묘한 세포 표현형에 민감한 간단한 선형 조정 함수를 학습합니다. 마지막으로, 각 세포는 이러한 클러스터의 가중 평균을 할당받고, 세포별 선형 요소에 의해 보정됩니다. 각 세포가 여러 클러스터에 속할 수 있기 때문에, 각 세포에는 잠재적으로 고유한 보정 계수가 부여됩니다. Harmony 는 셀 클러스터 할당이 안정될 때까지 이 네 단계를 반복합니다. Harmony 는 대규모 데이터셋으로 확장 가능하며, 다양한 인구 집단과 세부적인 하위 집단을 동시에 식별할 수 있는 유연성을 갖추고 있습니다. 또한 복잡한 실험 설계를 수용할 수 있으며, 다양한 모드 간 통합이 가능한 강력한 알고리즘입니다.¹⁰

# 추가 라이브러리 로드
library(harmony)
library(SeuratWrappers)

# Harmony 실행
SeuratObject <- RunHarmony(
  SeuratObject,
  group.by.vars = "batch"
)

14 클러스터링 하기

scRNA-seq 연구의 주요 목표 중 하나는 유사한 세포들을 클러스터로 그룹화하는 것입니다. 클러스터를 정의하는 방법은 여러 가지가 있으며, K- 평균, 계층적 클러스터링, 친화 전파와 같은 전통적인 접근 방식이 있지만, 최근 대부분의 파이프라인에서는 그래프를 사용하여 클러스터를 상호 연결된 세포 그룹으로 정의합니다 (위의 섹션 참조).그래프에서 클러스터를 정의하기 위해 Louvain, Leiden, Infomap 과 같은 다양한 커뮤니티 탐지 알고리즘이 사용됩니다. 이들 알고리즘의 주요 아이디어는 데이터셋 내에서 서로 더 많은 엣지 (높은 유사도의 연결) 를 가진 세포 그룹을 찾는 것입니다.모든 클러스터링 방법은 세포 간의 페어와이즈 거리를 정의하는 것으로 시작합니다. 일반적으로 이는 데이터를 스케일링하고 변수를 선택한 후 PCA 공간에서 세포 간의 유클리드 거리를 계산함으로써 수행됩니다. 따라서 클러스터링 결과는 데이터 전처리 단계에서의 선택에 크게 의존하게 됩니다. 특히 변수 유전자 선택과 PCA 에서 포함할 주성분의 수가 중요한 요소입니다.그래프 기반 방법에서는 PCA 공간에서의 거리를 사용하여 이웃하는 세포 사이에 엣지가 있는 그래프를 생성합니다. 이때 KNN 그래프의 이웃 수와 SNN 구성을 위한 가지치기 설정과 같은 여러 매개변수가 클러스터링 결과에 영향을 미칩니다.¹¹

14.1 루브앵 (Louvain)

루브앵 (Louvain) 방법은 Blondel 이 개발한 대규모 네트워크에서 커뮤니티를 추출하는 기법입니다.¹² 이 방법은 네트워크 노드 수에 대해 $O(n.log2n)$ 의 시간 복잡도를 가지는 탐욕적 최적화 방법입니다. 최적화하는 값은 모듈러리티 (modularity) 로 이는 커뮤니티 내 링크 밀도를 커뮤니티 간 링크와 비교하여 측정하는 범위 내 값입니다. 이 값을 최적화하는 것은 이론적으로 주어진 네트워크의 노드를 가장 잘 그룹화하는 결과를 도출하지만 모든 노드를 그룹으로 나누는 모든 가능한 반복 과정을 거치는 것은 비실용적이므로 휴리스틱 알고리즘이 사용됩니다.

SeuratObject <- FindClusters(
  object = SeuratObject,
  graph.name = "SNN",
  resolution = 0.8,
  algorithm = 1 # 1로 설정하면 Louvain을 의미합니다.
)

클러스터 수는 resolution 매개변수를 사용하여 제어할 수 있으며 값이 높을수록 더 많은 (작은) 클러스터가 생성됩니다.

14.2 라이든 (Leiden)

라이든 알고리즘은 반복적으로 적용되어 모든 커뮤니티의 하위 집합이 지역적으로 최적화된 분할로 수렴합니다.¹³ 이 알고리즘은 빠른 로컬 이동 접근 방식을 기반으로 하여 루브앵 알고리즘보다 더 빠르게 실행됩니다. 라이든 알고리즘은 다음과 같은 세 단계로 구성됩니다:

노드의 로컬 이동: 네트워크 내의 각 노드는 주변 노드와의 연결을 고려하여 위치를 조정합니다. 이 과정에서 각 노드는 자신의 이웃 중에서 가장 높은 모듈러리티를 제공하는 위치로 이동하게 됩니다. 이를 통해 각 노드가 속한 커뮤니티가 더욱 명확해집니다.
분할의 정제: 초기 클러스터링 결과를 바탕으로, 각 커뮤니티 내에서 노드의 재배치를 통해 더 나은 분할을 찾습니다. 이 단계에서는 커뮤니티 간의 경계를 명확히 하고, 중복된 연결을 줄이는 과정이 포함됩니다.
정제된 분할을 바탕으로 네트워크를 집계: 정제된 분할을 사용하여 새로운 집계 네트워크를 생성합니다. 이때, 정제되지 않은 분할을 기반으로 초기 집계 네트워크의 분할을 생성하여, 전체 네트워크 구조를 반영합니다.

이러한 과정을 통해 라이든 알고리즘은 데이터의 구조를 효과적으로 파악하고, 커뮤니티 간의 관계를 명확히 하여 최적화된 결과를 도출합니다. 라이든 알고리즘은 대규모 네트워크에서도 효율적으로 작동하며, 다양한 분야에서 널리 사용되고 있습니다.

SeuratObject <- FindClusters(
  object = SeuratObject,
  graph.name = "SNN",
  resolution = 0.8,
  algorithm = 4 # 4로 설정하면 Leiden을 의미합니다.
)

클러스터 수는 resolution 매개변수를 사용하여 제어할 수 있으며 값이 높을수록 더 많은 (작은) 클러스터가 생성됩니다.

14.3 계층적 클러스터링

계층적 클러스터링 (Hierarchical Clustering, HC) 은 클러스터의 계층 구조를 구축하는 클러스터 분석 방법입니다. 이 방법은 일반적으로 두 가지 전략으로 나뉘며, 각각 합병적 (agglomerative) 또는 분할적 (divisive) 입니다. 병합과 분할은 보통 탐욕적인 방식으로 결정됩니다. 계층적 클러스터링의 결과는 일반적으로 덴드로그램 (dendrogram) 으로 시각화됩니다.

계층적 병합 클러스터링 (Hierarchical Agglomerative Clustering, HAC) 의 표준 알고리즘은 $O(n^3)$ 의 시간 복잡도를 가지며 $O(n^2)$ 의 메모리를 필요로 하기 때문에, 중간 규모의 데이터 세트에서도 느릴 수 있습니다. 병합적인 경우에는 어떤 클러스터를 결합할 것인지, 분할적인 경우에는 어떤 클러스터를 분할할 것인지를 결정하기 위해 관측값 집합 간의 불일치 측정이 필요합니다. 대부분의 계층적 클러스터링 방법에서는 적절한 메트릭 (관측값 쌍 간의 거리 측정) 과 연결 기준을 사용하여 관측값 집합의 불일치를 관측값 쌍 간 거리의 함수로 정의합니다.

R 의 기본 stats 패키지에는 모든 샘플 쌍 사이의 거리를 계산하는 dist() 함수가 포함되어 있습니다. dist() 에서 사용할 수 있는 거리 방법은 다음과 같습니다:

유클리드 거리 (euclidean)
최대 거리 (maximum)
맨하탄 거리 (manhattan)
캔버라 거리 (canberra)
이진 거리 (binary)
민코프스키 거리 (minkowski)

또한, 이미 셀 간 거리 정보를 포함하는 그래프 (KNN 또는 SNN) 를 사용하여 직접 계층적 클러스터링을 수행할 수도 있습니다. 그러나 그래프의 거리는 반전되어 있으므로 (0 은 멀리 떨어져 있음을 나타내고 1 은 가깝다는 것을 나타냄), 그래프의 최대값 (인접 SNN 의 경우 1) 을 뺀 후 0 이 가까운 거리를 나타내고 1 이 먼 거리를 나타내도록 조정해야 합니다.

샘플 간의 거리를 계산한 후에는 본격적인 계층적 클러스터링을 진행할 수 있습니다. 이를 위해 hclust() 함수를 사용하며, 위에서 생성한 거리 객체를 활용하여 간단히 실행할 수 있습니다. 사용 가능한 방법에 대한 자세한 내용은 HC for networks 에서 확인할 수 있습니다.

# PCA에서 HC 실행
h_pca <- hclust(
    d = dist(
        SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:30],
        method = "euclidean"
    ),
    method = "ward.D2"
)

# 그래프에서 HC 실행
h_graph <- hclust(
    d = 1 - as.dist(SeuratObject@graphs$SNN),
    method = "ward.D2"
)

클러스터 계층이 정의되면, 다음 단계는 특정 클러스터에 속하는 샘플을 식별하는 것입니다.

plot(h, labels = FALSE)

세포에 대한 덴드로그램을 확인한 후, 고정된 높이 (threshold) 에서 나무 (tree) 를 자를 수 있습니다. 이 작업은 cutree 함수를 사용하여 수행할 수 있습니다. 또한, 해상도에 따라 다른 수준에서 나무를 자르는 것도 가능합니다.

# 덴드로그램 높이를 기반으로 나무 자르기
SeuratObject$HC_res <- cutree(
    tree = h,
    k = 18
)

# 클러스터 수를 기반으로 나무 자르기
SeuratObject$HC_res <- cutree(
    tree = h,
    h = 3
)

# 각 클러스터에 속하는 셀 수 확인
table(SeuratObject$HC_res)

클러스터의 수는 높이 (h) 를 기준으로 하거나 직접 k 매개변수를 설정하여 조절할 수 있습니다.

14.4 K- 평균 클러스터링

K- 평균 클러스터링 (k-means clustering) 은 신호 처리에서 시작된 벡터 양자화 방법으로, n 개의 관측값을 k 개의 클러스터로 분할하여 각 관측값이 평균 (클러스터 중심) 과 가장 가까운 클러스터에 속하도록 하는 것을 목표로 합니다. 이 평균은 클러스터의 프로토타입 역할을 합니다. 그러나 알고리즘은 전역 최적값으로 수렴하는 것을 보장하지 않으며, 결과는 초기 클러스터 설정에 따라 달라질 수 있습니다. 일반적으로 알고리즘은 빠르게 실행되므로 서로 다른 초기 조건으로 여러 번 실행하는 것이 일반적입니다. K- 평균 클러스터링은 주로 비슷한 공간 범위의 클러스터 (모두 같은 크기) 를 찾지만, 기대값 최대화 메커니즘 덕분에 클러스터가 다양한 형태를 가질 수 있습니다.

K- 평균은 많은 응용 분야에서 널리 사용되는 클러스터링 알고리즘입니다. R 에서는 kmeans 함수를 사용하여 쉽게 적용할 수 있습니다. 일반적으로 이 알고리즘은 표현 데이터의 축소된 차원 표현 (대부분 PCA 를 사용함) 에 적용되며, 이는 저차원 거리의 해석 가능성 덕분입니다. K- 평균 클러스터링을 수행할 때는 클러스터의 수를 미리 정의해야 하며, 결과는 클러스터 중심의 초기화에 따라 달라지므로 일반적으로 다양한 시작 구성으로 k-means 를 실행하는 것이 권장됩니다 (nstart 인수를 통해 설정 가능).

set.seed(42)  # 재현성을 위해 설정합니다.

# K-means 클러스터링 실행
SeuratObject$kmeans_12 <- kmeans(
    x = SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:50],
    centers = 12,
    iter.max = 50,
    nstart = 10
)$cluster

클러스터의 수는 centers 매개변수를 사용하여 설정할 수 있습니다. 그렇다면 어떤 클러스터링 해상도를 선택해야 할까요?클러스터링 해상도를 결정하는 것은 간단한 문제가 아니며, 여러 방법을 통해 접근할 수 있습니다. 합리적인 클러스터 수를 정하는 것은 복잡한 과정으로, 샘플에 대한 생물학적 지식과 조사가 필요합니다.차별적 유전자 발현 (Differential gene expression) 분석은 이 과정에 큰 도움이 될 수 있습니다. 만약 두 클러스터가 동일한 차별적 유전자 발현을 보이고 구별할 만한 유전자가 없다면, 이들을 개별 클러스터로 분리하는 것은 바람직하지 않을 수 있습니다. 또한, 품질 관리 지표 (QC metrics) 를 통해 생성된 클러스터를 검토하여, 일부 클러스터가 낮은 품질의 세포나 이중체로 구성되지 않았는지 확인하는 것이 중요합니다.

14.5 클러스터링 결과 비교하기

클러스터링 결과를 다른 클러스터링 기법과 비교하는 것은 매우 유용합니다. 이를 다음과 같은 방식으로 수행할 수 있습니다: 위 코드는 비교 대상 클래스에 속하는 세포 수에 대한 정보를 포함하는 테이블을 생성합니다. 이 데이터를 막대 그래프로 시각화하여 클러스터링 결과를 쉽게 비교할 수 있습니다.

# 비교 테이블 생성
comparison_table <- table(
    list(
        SeuratObject@meta.data[,"METADATA_FACTOR_1"],
        SeuratObject@meta.data[,"METADATA_FACTOR_2"]
    )
)

# 데이터를 백분율로 변환
comparison_table <- t(t(comparison_table) / colSums(comparison_table)) * 100

# 막대 그래프 그리기
barplot(
    comparison_table,
    col = 1:nrow(comparison_table),  # 색상 설정
    border = NA,
    las = 2                         # 축 라벨 방향 설정
)

15 차등발현 (Differentially Expressed) 분석

클러스터가 정의된 후, 각 클러스터를 특징짓는 유전자를 찾는 것이 종종 유익합니다. scRNA-seq 에서의 차별적 발현 (DE) 분석 방법은 일반적인 RNA-seq 방법과 다소 다릅니다. 한편으로는 대량 RNA-seq 에 비해 샘플 (개별 세포) 의 수가 훨씬 많지만, 반면에 scRNA-seq 데이터는 노이즈가 많고 드롭아웃에 취약하여 분석이 복잡해집니다.¹⁴

15.1 DEG(Differential expression gene) 찾기

차별적 발현 유전자 (DEG) 는 종종 “마커 유전자”로 언급되지만, 대부분의 DE 테스트는 한 그룹의 세포와 다른 그룹의 세포 간에 발현이 더 높은 유전자를 감지하도록 설계되었다는 점을 인식해야 합니다. DEG 가 자동으로 세포 유형을 위한 고유한 마커가 되지는 않습니다.Seurat 패키지에는 DE 분석을 위한 여러 가지 테스트가 구현되어 있으며, 일부는 scRNA-seq 에 특화되어 있고, 다른 일부는 벌크 RNA-seq 에서 사용됩니다.

wilcox: Wilcoxon 순위 합 검정을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별
bimod: 단일 세포 유전자 발현에 대한 우도 비율 검정 ¹⁵
roc: ROC 분석을 사용하여 유전자 발현의 ’마커’를 식별합니다. 각 유전자에 대해 해당 유전자를 기반으로 구축된 분류기를 평가하여 두 그룹의 세포를 분류합니다. AUC 를 사용하여 각 유전자에 대한 분류기의 성능을 평가하며, AUC 값이 1 인 경우 해당 유전자의 발현 값만으로 두 그룹을 완벽하게 구분할 수 있음을 의미합니다. AUC 값이 0 인 경우 완벽한 분류가 이루어지지만 반대 방향입니다. AUC 값이 0.5 인 경우 해당 유전자가 두 그룹을 구분하는 데 예측력이 없음을 나타냅니다. 이 방법은 예상되는 차별적 발현 유전자의 순위를 지정된 예측력$(abs(AUC-0.5) * 2)$ 행렬로 반환합니다.
t: 두 그룹의 세포 간 차이가 있는 유전자 식별을 위한 Student 의 t-test 사용
negbinom: 음이항 일반화 선형 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별 (UMI 기반 데이터 세트에만 사용)
poisson: 포아송 일반화 선형 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별 (UMI 기반 데이터 세트에만 사용)
LR: 로지스틱 회귀 모델을 사용하여 차이가 있는 유전자 식별. 각 기능을 개별적으로 사용하여 그룹 소속을 예측하는 로지스틱 회귀 모델을 구성하고 이를 가능도 비율 검정과 비교합니다.
MAST: scRNA-seq 데이터에 특화된 허들 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별. DE 테스트를 실행하기 위해 MAST 패키지를 사용합니다.
DESeq2: DESeq2 를 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별. DESeq2 는 음이항 분포를 사용하는 모델에 기초하며, 이 테스트는 세포 그룹 간 평균 차이 (또는 감지율의 백분율) 에 따른 유전자의 사전 필터링을 지원하지 않습니다. 그러나 유전자는 세포 그룹 모두에서 최소 감지율 (min.pct) 기준으로 사전 필터링될 수 있습니다. 이 방법을 사용하려면 DESeq2 를 설치하십시오.

Seurat 객체의 모든 클러스터에 대해 각 클러스터와 모든 다른 세포 간의 DE 예측을 실행하려면 FindAllMarkers 함수를 사용하세요.

markers <- FindAllMarkers(
    object = SeuratObject,
    assay = "RNA",
    logfc.threshold = 0.25,
    test.use = "wilcox",
    slot = "data",
    min.pct = 0.1,
    min.diff.pct = -Inf,
    only.pos = FALSE,
    max.cells.per.ident = Inf,
    latent.vars = NULL,
    min.cells.feature = 3,
    min.cells.group = 3,
    pseudocount.use = 1,
    return.thresh = 0.01
)

아래에 나타낸 여러 가지 기준을 통해 특정 유전자를 포함하거나 출력을 필터링할 수 있습니다. 예를 들어, 상향 조절된 유전자만 테스트하면 분석 속도를 높일 수 있습니다:

only.pos = TRUE 를 사용하여 상향 조절된 유전자만 테스트
유전자가 발현된 최소 세포 수를 설정하는 min.cells.feature
유전자가 발현된 클러스터의 최소 세포 수를 나타내는 min.cells.group
클러스터에서 min.pct 이상 발현된 유전자만 테스트
두 그룹 간의 발현 비율 차이가 min.pct.diff 인 유전자만 테스트
p.value < return.thresh 인 유전자만 반환
logFoldchange > logfc.threshold 인 유전자만 반환

또한 고려해야 할 몇 가지 중요한 사항이 있습니다:

객체에 여러 분석이 포함된 경우, 올바른 분석 결과에서 차등발현 분석을 실행해야 합니다. 예를 들어, 데이터가 통합된 경우에는 Batch 를 고려해야 합니다.

cell_selection <- SeuratObject[, SeuratObject$seurat_clusters == 4]
cell_selection <- SetIdent(cell_selection, value = "Batch")

# 차별적 발현 계산
DGE_cell_selection <- FindAllMarkers(
    object = cell_selection,
    logfc.threshold = 0.2,
    test.use = "wilcox",
    min.pct = 0.1,
    min.diff.pct = -Inf,
    only.pos = FALSE,
    max.cells.per.ident = Inf,
    latent.vars = NULL,
    min.cells.feature = 3,
    min.cells.group = 3,
    pseudocount.use = 1,
    return.thresh = 0.01
)

또한 FindMarkers 함수를 사용하여 두 세트의 세포 간의 차이를 테스트할 수 있으며, 각 그룹에 대한 셀 이름을 cells 및 cells 로 지정할 수 있습니다.

15.2 DEG 시각화하기

DE 테스트를 실행한 후, 다양한 방법으로 유전자를 시각화하고 싶을 수 있습니다. 그러나 먼저 각 클러스터에서 상위 DEG 를 추출해야 합니다. 각 클러스터별로 상위 5 개의 유전자를 선택하는 방법은 다음과 같습니다:

top5 <- markers %>% group_by(cluster) %>% top_n(-5, p_val_adj)

이를 히트맵으로 시각화할 수 있습니다:

DoHeatmap(
    object = SeuratObject,
    features = as.character(unique(top5$gene)),
    group.by = "seurat_clusters",
    assay = "RNA"
)

또는 각 클러스터가 평균 발현으로 색상별로 나타내고, 해당 유전자를 발현하는 세포의 비율로 크기별로 나타내는 점 도플롯을 사용할 수 있습니다:

# coord_flip() 함수를 위해 필요
library(ggplot2) 

DotPlot(
    object = SeuratObject,
    features = as.character(unique(top5$gene)),
    group.by = "seurat_clusters",
    assay = "RNA"
) + coord_flip()

또한 각 유전자에 대해 바이올린 플롯을 그릴 수 있습니다:

VlnPlot(
    object = SeuratObject,
    features = as.character(unique(top5$gene)),
    ncol = 3,
    group.by = "seurat_clusters",
    assay = "RNA"
)

바이올린 플롯을 배치별로 분할할 수도 있습니다. 이것은 감지된 DEG 가 단일 배치로 인해 주도되는 것이 아닌지 확인하는 데 매우 유용할 수 있습니다.

VlnPlot(
    object = SeuratObject,
    features = as.character(unique(top5$gene)),
    ncol = 3,
    group.by = "seurat_clusters",
    assay = "RNA",
    split.by = "Batch"
)

16 기능적 풍부도 분석 (Functional enrichment anaysis)

enrichR 을 사용해 DEG 에 대하여 기능적 풍부도 분석을 할 수 있습니다.

# 추가 패키지 불러오기
library(enrichR)

# 사용 가능한 데이터베이스를 확인
enrichR::listEnrichrDbs()

enrich_results <- enrichr(
    genes = top5$gene[top5$cluster == "1"], # 관심있는 클러스터 이름
    databases = "KEGG_2019_Human" )[[1]] # 사용할 데이터베이스

enrichR 에 포함된 유용한 데이터베이스 목록은 다음과 같습니다:

GO_Biological_Process_2017b
KEGG_2019_Human
KEGG_2019_Mouse
WikiPathways_2019_Human
WikiPathways_2019_Mouse

결과를 시각화하기 위해 간단한 막대 그래프를 사용할 수 있습니다. 예를 들어:

library(rafalib)  # 추가 패키지 로드

# 그래프 설정
rafalib::mypar(1, 1, mar = c(3, 20, 2, 1))

# 막대 그래프 그리기
barplot(
    height = -log10(enrich_results$P.value)[20:1],
    names.arg = enrich_results$Term[20:1],
    horiz = TRUE,
    las = 2,
    border = FALSE,
    cex.names = 0.8  # 이름 크기 설정 (예시 값)
)

# 기준선 추가
abline(v = c(-log10(0.05)), lty = 2)  # 유의성 기준선
abline(v = 0, lty = 1)                 # 0 기준선

17 경로 분석 (trajectory analysis)

탐색적 데이터 분석 및 연구 데이터에는 다양한 방법으로 수행되는 데이터 처리 단계가 포함됩니다. 이러한 데이터 처리 단계는 경로 분석에 필요한 정보를 준비하는 데 도움을 줍니다. 주요 단계는 다음과 같습니다:

경로 분석을 수행할 차원 축소 기법 선택 (PCA > MDS > DM > UMAP > t-SNE)
세포 클러스터링 정보
경로 분석을 위한 “최적” 방법 선택은 일반적으로 다른 scRNA-seq 분석과 마찬가지로 간단하지 않은 과정입니다.¹⁶

17.1 Slingshot 사용하기

슬링샷 (Slingshot) 은 차원 축소, 클러스터링 및 주요 곡선 (Principal curves) 과 같은 다른 경로 분석 방법에서 발견되는 여러 인기 있는 구성 요소를 포함하고 있어 가장 전형적인 방법 중 하나입니다. 슬링샷은 특정 차원 축소 방법에 제한되지 않으며, 어떤 방법을 사용할지 결정할 때 고려해야 할 중요한 세 가지 사항이 있습니다:

데이터의 전체 복잡성을 포착하기에 충분한 차원을 사용해야 합니다. 이는 특히 PCA 및 MDS 와 같은 선형 차원 축소에서 중요하며, 경로의 위상이 이분화보다 복잡할 경우 더욱 그렇습니다. 두 차원에서 경로가 명확하게 보이지 않을 때도 여전히 다차원에서 경로를 시각화할 수 있습니다.
일부 차원 축소 방법은 고밀도 영역에서 거리를 확대할 수 있습니다. t-SNE 와 UMAP 은 이러한 문제를 겪을 수 있습니다. 우리가 사용하는 데이터셋에서는 이러한 문제가 발생하지 않지만, 생물학적 모집단을 무작위로 샘플링할 경우 문제가 발생할 수 있습니다.
일부 차원 축소 방법은 셀의 그룹화를 강제하고 연속성을 제거하려고 합니다. t-SNE 와 UMAP 도 이러한 문제를 일부 가지고 있으며, 이로 인해 실제로는 그렇지 않은 셀들이 함께 묶일 수 있습니다.

이러한 이유로 MDS 와 확산 맵이 종종 선호되는 선택지입니다. 그러나 샘플이 균형 잡힌 경우 t-SNE 와 UMAP 도 효과적으로 작동할 수 있습니다.

dimred <- SeuratObject@reductions[['pca']]@cell.embeddings  # 사용될 embedding 지정

그 답은 간단합니다. 클러스터링은 문제를 단순화하며, 경로 상에서 거의 동일한 위치에 있는 세포 그룹을 찾아냅니다. 이러한 그룹은 이후 세포가 어떤 경로를 통해 이동하고 분기되는지를 파악하는 데 연결될 수 있습니다. 이상적으로 TI 를 위한 클러스터링 방법은 밀도를 기반으로 셀을 그룹화해서는 안 됩니다. 예를 들어 DBSCAN 과 같은 밀도 기반 클러스터링은 정의상 모든 세포가 고밀도 영역을 통해 연결되지만 유전자 발현에서는 차이가 발생할 수 있습니다. 여기서는 Seurat의 표준 클러스터링 방법인 라이든 (Leiden) 클러스터링을 사용할 것입니다.

clustering <- factor(
    SeuratObject@meta.data[, "seurat_clusters"]
)

이제부터 이 클러스터링 결과를 경로 분석에 적용할 것입니다. 전체 프로세스는 단일 함수인 slingshot 을 사용하여 수행할 수 있으며, 이는 경로 분석을 위한 두 가지 주요 단계를 간단히 묶어 놓은 래퍼 (wrapper) 입니다. 프로세스의 첫 번째 단계는 세포 계통 (cell lineages) 을 정의하고, 그 다음 데이터를 통해 경로를 정의하는 곡선을 맞추는 것입니다. 이러한 단계들은 아래에서 설명합니다.

17.1.1 Slingshot 으로 세포 계통 정의하기

# 기본 설정으로 Slingshot 세포 계통 식별을 실행
library(slingshot)

lineages <- getLineages(
     data = dimred,
     clusterLabels = clustering
)

# Lineages를 시각화.
par(mfrow = c(1,2))
plot(
     dimred[, 1:2],
     col = clustering,
     cex = ,
     pch = 16
)

여기서 우리는 경로 분석의 한 가지 중요한 문제를 살펴볼 수 있습니다: 경로 분석이 어디에서 시작되는지를 알 수 없다는 점입니다. 추가 정보가 없이는 이를 파악하기가 거의 불가능합니다. 경로 분석의 시작점과 끝점을 정의하기 위해서는 사전 생물학적 정보가 필요합니다.

# Cluster 5에서 Slingshot을 시작하기.
lineages <- getLineages(
    data = dimred,
    clusterLabels = clustering,
    start.clus = "5"
)

lineages

17.1.2 주요 경로 정의하기

시작 클러스터를 정한 뒤에는 Slingshot 을 사용하여 경로 분석을 수행할 수 있습니다. Slingshot 의 알고리즘은 초기 경로를 반복적으로 수정하여 데이터 포인트와 더 잘 일치되도록 합니다. 구체적으로 Slingshot 은 두 가지 기능을 사용합니다.

각 “계통”마다 주요 경로를 찾습니다: 여기서 “계통”이란 특정 시작 클러스터에서 특정 종단 클러스터로 이어지는 클러스터 집합을 의미합니다. 즉, 각 계통은 세포들이 어떤 경로를 따라 이동하는지를 나타내며, 이를 통해 세포의 발달 과정을 이해할 수 있습니다.
동일한 클러스터 집합을 가진 계통은 주요 곡선이 겹치는 클러스터 주위에 묶이도록 제약됩니다: 이는 서로 연결된 세포 그룹이 같은 경로를 공유하도록 하여, 분석의 일관성을 높이는 역할을 합니다. 이렇게 하면 각 계통이 서로 다른 경로를 따라 이동하더라도, 동일한 클러스터 내에서의 관계를 유지할 수 있습니다.

Slingshot 의 getCurves() 함수는 실행하는 데 시간이 오래 걸리므로 각 계통에서 사용할 셀의 수를 줄여 곡선 적합 프로세스의 속도를 높일 수 있습니다. 이상적으로는 모든 셀을 사용하는 것이 좋지만 여기서는 approx_points 를 300 으로 설정하여 속도를 높였습니다.

# 곡선 데이터 가져오기
curves <- getCurves(
    lineages,
    approx_points = 300,
    thresh = 0.01,
    stretch = NULL,  # stretch 인수에 NULL 설정 (예시)
    allow.breaks = FALSE,
    shrink = NULL    # shrink 인수에 NULL 설정 (예시)
)

# 곡선 데이터 확인
curves

# 차원 축소 결과 플롯
plot(
    dimred,
    col = clustering,
    asp = 1,
    pch = 16
)

# 곡선 추가
lines(
    curves,
    lwd = 3,
    col = 'black'
)

17.2 경로 분석를 따라 변하는 유전자를 찾기 (Differential expression along trajectories)

경로 분석을 해석하는 주요 방법 중 하나는 경로를 따라 변하는 유전자를 찾는 것입니다. 경로 분석에서 차별적 발현 (differential expression) 을 정의하는 방법은 여러 가지가 있습니다:

특정 경로를 따라 발현 변화: 예를 들어, 세포가 유사한 시간에 따라 발현이 어떻게 변화하는지를 분석합니다.
가지 사이의 발현 차이: 서로 다른 경로 (또는 가지) 간의 유전자 발현 차이를 비교합니다.
분기점에서의 발현 변화: 세포가 분기점에 도달했을 때 발현이 어떻게 변화하는지를 살펴봅니다.
경로 분석의 어느 지점에서든 발현 변화: 경로의 모든 지점에서 유전자 발현의 변화를 평가합니다.

tradeSeq 는 경로 분석 결과에서 차별적 발현 유전자를 찾는 도구 중 하나입니다. 이 도구는 일반화된 가법 모형 (Generalized Additive Models, GAMs) 을 사용하여 경로를 따라 유전자 발현을 조정하고 경로의 여러 지점 사이에서 특정 계수가 통계적으로 다른지를 테스트합니다. GAM 모델을 적합시키는 과정은 시간이 많이 걸립니다. 따라서 여기에서는 매우 엄격한 필터링을 통해 많은 유전자를 제거하였습니다. 실제로는 모든 유전자 데이터를 사용하는 것이 좋습니다.

# 추가 패키지 로드
BiocParallel::register(BiocParallel::SerialParam())
library(tradeSeq)

# 연산 속도를 높이기 위해 일부 유전자 제거
counts <- as.matrix(SeuratObject[["RNA"]]$counts)
dim(counts)  # 원본 카운트 행렬 차원 확인

# 필터링된 카운트 행렬 생성
filt_counts <- counts[rowSums(counts > 5) > ncol(counts) / 100, ]
dim(filt_counts)  # 필터링 후 카운트 행렬 차원 확인

# 감마 분포 fitting
sce <- fitGAM(
    counts = as.matrix(filt_counts),
    sds = curves
)

# 결과 플롯
plotGeneCount(
    curve = curves,
    counts = filt_counts,
    clusters = clustering,
    models = sce
)

17.2.1 두 시간점 사이에서 변하는 유전자

원하는 경우 사용자 정의 유사 시간 (Pseudotime) 값을 설정하여 두 시간점 사이에서 변하는 유전자를 찾을 수 있습니다. 예를 들어, 전구 세포 표시자와 같은 특정 유전자를 분석할 수 있습니다. 기본적으로는 시작점과 끝점 사이의 차이를 살펴보는 방식으로 진행됩니다.

# Pseudotime 시작과 끝 간의 연관성 테스트
pseudotime_start_end_association <- startVsEndTest(
    sce,
    pseudotimeValues = c(0, 1)
)

# feature_id 추가
pseudotime_start_end_association$feature_id <- rownames(pseudotime_start_end_association)

# p-value가 0.05 미만인 feature_id 필터링 및 선택
feature_id <- pseudotime_start_end_association %>%
    filter(pvalue < 0.05) %>%
    top_n(1, waldStat) %>%
    pull(feature_id)

# 차별 발현 플롯 그리기
plot_differential_expression(feature_id)

17.2.2 경로가 갈라지는 사이에서 변하는 유전자

두 가지 사이에서 차별적으로 발현되는 유전자는 특히 흥미로운 연구 대상입니다. 이러한 유전자 중 일부는 이전의 두 시간점 사이에서 변하는 유전자 분석에서 이미 발견되었을 수 있습니다. “경로가 갈라지는 사이에서 변하는 유전자”를 정의하는 여러 가지 방법이 있으며 방법에 따라 함수가 다릅니다.¹⁷

diffEndTest: 각 가지의 끝점에서 발현 차이를 분석합니다.
earlyDETest: 가지의 분기점에서 발현 차이를 평가합니다.
patternTest: 경로 분석의 어느 지점에서든 발현 차이를 탐색합니다. 이 마지막 함수는 두 가지 간의 유사 시간이 정렬되어 있어야 함을 요구합니다.

# 차별적 끝 연관성 테스트 수행
different_end_association <- diffEndTest(sce)
different_end_association$feature_id <- rownames(different_end_association)

# p-value가 0.05 미만인 feature_id 필터링 및 정렬
feature_id <- different_end_association %>%
    filter(pvalue < 0.05) %>%
    arrange(desc(waldStat)) %>%
    dplyr::slice(1) %>%
    pull(feature_id)

# 차별 발현 플롯 그리기
plot_differential_expression(feature_id)

# 가지점 연관성 테스트 수행
branch_point_association <- earlyDETest(sce)
branch_point_association$feature_id <- rownames(branch_point_association)

# p-value가 0.05 미만인 feature_id 필터링 및 정렬
feature_id <- branch_point_association %>%
    filter(pvalue < 0.05) %>%
    arrange(desc(waldStat)) %>%
    dplyr::slice(1) %>%
    pull(feature_id)

# 차별 발현 플롯 그리기
plot_differential_expression(feature_id)

더 자세한 것은 공식 문서 를 확인하십시오.

18 더 읽을 거리

scRNA-seq 분야는 방대하고 끊임없이 발전하고 있습니다. 이로 인해 다양한 방법, 도구, 소프트웨어가 존재하여 유사한 작업을 수행하는 데 활용될 수 있습니다. 아래는 이 분야에 대한 더 깊은 이해를 돕기 위해 탐색할 수 있는 몇 가지 추가 자료입니다.

18.1 온라인 강의

18.2 프레젠테이션 및 강의

18.3 도구 및 파이프라인

18.4 웹사이트

19 Reference

Footnotes

https://www.pnas.org/content/116/20/9775)↩︎
https://distill.pub/2016/misread-tsne, http://jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf↩︎
https://umap-learn.readthedocs.io/en/latest↩︎
https://arxiv.org/abs/1802.03426↩︎
https://projecteuclid.org/euclid.ajm/1253539864↩︎
https://en.wikipedia.org/wiki/Independent_component_analysis↩︎
https://www.nature.com/articles/nbt.4091↩︎
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9↩︎
https://www.nature.com/articles/s41587-019-0113-3↩︎
https://www.nature.com/articles/s41592-019-0619-00↩︎
https://f1000research.com/articles/7-1141↩︎
https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/pdf↩︎
https://www.nature.com/articles/s41598-019-41695-z.pdf↩︎
https://www.nature.com/articles/nmeth.4612↩︎
McDavid et al., Bioinformatics, 2013↩︎
https://www.nature.com/articles/s41587-019-0071-9↩︎
Cannoodt, Robrecht, Wouter Saelens, and Yvan Saeys. 2016.”Computational Methods for Trajectory Inference from Single-Cell Transcriptomics.”European Journal of Immunology 46 (11): 2496–2506.↩︎

--- title: "Seurat 라이브러리 치트시트(Cheatsheet)" author: "Taeyoon Kim" categories: [R, Tip, scRNA-seq, Bioinformatics, Seurat, Cheatsheet] draft: false date: "2024-11-12" date-modified: last-modified --- ![](Cheatsheet_Seurat.png){width=100% fig-align="center"} [Seurat](https://satijalab.org/seurat/) 은 scRNA-seq 분석에서 가장 널리 사용되는 R 패키지로, 다양한 기능을 통해 데이터를 처리하고 시각화합니다. 이 문서는 주요 명령어와 기능을 요약하여 사용자가 데이터를 효율적으로 다룰 수 있도록 지원합니다. 특히, 데이터 로드, 전처리, 클러스터링 및 시각화와 같은 일반적인 작업에 대한 예제를 제공합니다. Seurat v5.0.1 버전을 사용한 환경에서 테스트되었습니다. 다양한 명령어와 튜토리얼은 [공식 페이지](https://satijalab.org/seurat/essential_commands.html) 에서 확인하세요. # 파일 불러오기 scRNA-seq 정보를 저장할 수 있는 파일에는 여러 형식이 있으며, 모든 형식을 여기에 모두 나열할 수는 없습니다. 다만 보편적인 형식은 아래와 같습니다. 1. `SMARTseq2` 에서는 행렬과 샘플 메타데이터를 `csv`, `tsv` 또는 `txt` 형식으로 저장합니다. 2. `10X Chromium` 데이터의 기본 형식 중 하나는 `.tsv` 로 저장된 유전자 및 세포에 대한 주석이 포함된 압축된 희소 행렬 파일 `.mtx` 를 사용합니다. 3. 현재는 `10X Chromium` 에서 많이 사용되는 형식은 `HDF5` 파일입니다. `HDF5` 는 속도가 빠르고 확장 가능하며 한 번에 사용할 데이터의 일부만을 불러오는 기능도 제공합니다. 또한 메타데이터를 동일한 파일에 저장하며 데이터를 이진 압축된 희소 행렬 형식으로 저장합니다. ```r # csv, tsv, txt 등의 형식 읽기 raw_matrix <- read.delim( file = "data/folder_sample1.csv", row.names = 1 ) # 희소 행렬로 변환 sparse_matrix <- Matrix::Matrix( data = raw_matrix, sparse = TRUE ) # mtx 형식에서 읽고 변환하기 (폴더에 있는 파일들을 읽음) sparse_matrix <- Seurat::Read10X( data.dir = "data/folder_sample1" ) # h5 형식에서 읽고 변환하기 sparse_matrix <- Seurat::Read10X_h5( filename = "data/matrix_file.h5", use.names = TRUE ) ``` 메타 데이터는 데이터 프레임 형식으로 가져와야 합니다. ```r # csv 또는 tsv 형식에서 메타데이터 읽기 metadata <- read.delim( file = "data/metadata.csv", sep = ",", row.names = 1 ) ``` # Seurat 객체 만들기 개발자들은 데이터 분석 프로세스를 조금 더 간단하게 만들기 위해 데이터를 간결한 형식으로 저장하는 객체를 만들었습니다. **Seurat** 에서는 `CreateSeuratObject` 함수를 사용해 희소 행렬 또는 카운트가 포함된 데이터를 입력으로 객체를 생성합니다. ```R SeuratObject <- CreateSeuratObject( counts = sparse_matrix, assay = "RNA", project = "SAMPLE1", meta.data = metadata ) ``` Seurat 객체는 객체 이름 뒤에 `@` 또는 `$` 문자를 사용하여 내용에 쉽게 접근할 수 있습니다. - `@` 기호는 `assays`, `meta.data`, `graphs` 및 `reduction` 슬롯을 포함한 모든 분석 결과에 접근할 수 있습니다. - `$` 기호를 사용하면 **Seurat** 객체의 메타데이터 열에 직접 액세스할 수 있습니다. 즉, `SeuratObject$column1` 은 `SeuratObject@meta.data$column1` 과 동일한 결과를 출력합니다. 기본적으로 분석 데이터는 RNA 이라는 이름의 `assay` 슬롯에 저장됩니다. 분석을 수행할 때 추가 `assays` 가 생성되며 예를 들어 유전자의 원시 데이터 (counts), 정규화한 데이터 (data), 스케일링/회귀된 데이터 (scale.data) 및 유전자의 분산에 관한 정보 (var) 들이 있습니다. **Seurat** 객체의 모든 슬롯에 대한 자세한 설명은 [Seurat 위키](https://github.com/satijalab/seurat/wiki/Assay) 를 참고하세요. ## 메타데이터에 새로운 열 추가 새로운 열을 추가하려면 `$` 기호를 사용하여 벡터를 메타데이터 열에 할당할 수 있습니다. ```R SeuratObject$NEW_COLUMN_NAME <- setNames( colnames(SeuratObject), vector_NEW_DATA ) ``` 또는 `AddMetaData` 함수를 사용하여 열을 추가할 수 있습니다. ```R SeuratObject <- AddMetaData( object = SeuratObject, metadata = vector_NEW_DATA, col.name = "NEW_COLUMN_NAME" ) ``` ## Seurat 객체 하위 집합 만들기 일부 데이터를 따로 때어내는 방법은 `R` 의 데이터프레임 문법을 사용하는 방법과 **Seurat** 함수를 사용하는 방법이 있습니다. 먼저 데이터 프레임 문법은 아래와 같습니다. ```R # Feature 하위 집합 SeuratObject[vector_FEATURES_TO_USE, ] # 세포 하위 집합 SeuratObject[, vector_CELLS_TO_USE] ``` **Seurat** 에서 제공하는 함수는 아래와 같습니다. ```r # Identity 클래스를 기반으로 Seurat 객체를 하위 집합으로 추출합니다. 자세한 내용은 ?SubsetData를 참조하세요. subset(x = SeuratObject, idents = "Bcell") subset(x = SeuratObject, idents = c("CD4T", "CD8T"), invert = TRUE) # 유전자/특성의 발현 수준을 기준으로 하위 집합을 만듭니다. subset(x = SeuratObject, subset = MS4A1 > 3) # 기준의 조합에 따라 하위 집합을 만듭니다. subset(x = SeuratObject, subset = MS4A1 > 3 & PC1 > 5) subset(x = SeuratObject, subset = MS4A1 > 3, idents = "B 세포") # 객체 메타 데이터의 값에 따라 하위 집합을 만듭니다. subset(x = SeuratObject, subset = orig.ident == "Replicate1") # 각 identity 클래스당 세포 수를 다운샘플링합니다. subset(x = SeuratObject, downsample = 100) ``` ## 시각화하기 가장 일반적으로 사용되는 시각화 함수는 바이올린 플롯과 UMAP(차원 축소 후의 산점도) 입니다. :::{.callout-warning} UMAP 을 그리기 전에 반드시 먼저 차원 축소를 계산해야 합니다. ::: 연속 변수에 대하여 바이올린 플롯 시각화하는 코드 ```R VlnPlot( object = SeuratObject, group.by = "orig.ident", features = c("percent_mito"), pt.size = 0.1, ncol = 4, y.max = 100 ) + NoLegend() ``` 연속 변수에 대한 UMAP 을 그리는 방법 (예시: 유전자 발현량): ```R FeaturePlot( object = SeuratObject, features = c("FEATURE_1", "FEATURE_2", "FEATURE_3"), reduction = "umap", dims = c(1, 2), order = TRUE, pt.size = 0.1, ncol = 3 ) ``` 범주형 변수에 대한 UMAP 을 그리는 방법 (예시: 세포 유형): ```R DimPlot( object = SeuratObject, group.by = c("DATASET"), reduction = "umap", dims = c(1, 2), pt.size = 0.1, label = TRUE, ncol = 3 ) ``` 이것 말고도 다양한 시각화 함수가 있으니 자세한 것은 [공식 문서](https://satijalab.org/seurat/articles/visualization_vignette) 를 참고하세요. ## 데이터 합치기 두 개이상의 데이터셋을 결합하려면: ```R merged_obj <- merge( x = SeuratObject1, y = c(SeuratObject2, SeuratObject3, SeuratObject4), add.cell.ids = c("Dataset1", "Dataset2", "Dataset3", "Dataset4") ) ``` # 전체 Feature 의 수 확인하기 일반적인 방식은 적어도 200 개의 유전자가 있는 세포와 유전자는 적어도 3 개 이상의 세포에서 발현된 것만 남기는 것입니다. 다만 이 수치는 데이터를 준비할 때 사용된 라이브러리 기술에 매우 의존적입니다. 따라서 주의해서 사용해야 합니다. **Seurat** 에서는 `nCount_RNA`: 세포 당 UMIs(Unique Molecular Identifiers) 의 수; `nFeature_RNA`: 각 세포에서 감지된 유전자의 수로 저장하고 있습니다. 시각화를 통해 유전자 수와 UMI 의 갯수를 살펴보는 코드는 아래와 같습니다. ```r VlnPlot( SeuratObject, group.by = "orig.ident", features = c("nFeature_RNA", "nCount_RNA"), pt.size = 0.1 ) + NoLegend() ``` # 유전자 관련 품질 관리 단일 세포에서 가장 많이 감지되는 유전자는 주로 미토콘드리아 유전자 (`MT-`), 리보솜 유전자 (`RPL` 및 `RPS`), 그리고 기타 항상 발현되는 구조 단백질들 (예: ACTB, TMSB4X, B2M, EEF1A1) 입니다. 이러한 유전자 데이터는 세포 품질 지표 계산, 데이터 정규화, 배치 효과 보정 등 단일 세포 데이터 분석의 여러 목적에 활용될 수 있습니다.[^1] ## 미토콘드리아 유전자 비율 높은 비율의 미토콘드리아 유전자는 일반적으로 손상된 세포를 의미합니다. 이는 세포질 속의 RNA 가 쉽게 소실되는 반면, 미토콘드리아 속의 RNA 는 비교적 덜 손실되기 때문입니다. 다음은 미토콘드리아 유전자 비율을 계산하고 메타데이터 테이블에 추가하는 예시입니다. ```r # 미토콘드리아 유전자의 비율 계산 SeuratObject <- PercentageFeatureSet( object = SeuratObject, pattern = "^MT", assay = "RNA", col.name = "percent_mito" ) ``` :::{.callout-note} 유전자 관련 품질 관리에 대한 설명은 사람 (human) 에 해당하는 내용으로 다른 종의 데이터를 처리하는 경우에는 코드가 수정되어야 합니다. 예시로 생쥐 (mouse) 의 경우 미토콘드리아 (`mt-`), 리보소말 (`Rpl`) 입니다. ::: ## 리보솜 유전자 비율 리보솜 유전자는 모든 세포에서 가장 많이 발현되는 유전자 중 하나이며, 미토콘드리아 유전자와는 반대로 미토콘드리아 수에 반비례하는 경향이 있습니다. 즉, 미토콘드리아의 수가 많을수록 리보솜 유전자의 감지 수준은 낮아집니다 (다만, 이 관계는 비선형적입니다). 리보솜 단백질에서 유래하는 유전자 발현 비율을 계산하는 방법은 미토콘드리아 유전자 비율을 계산하는 방법과 유사합니다. ```r # 리보솜 유전자의 비율 계산 SeuratObject <- PercentageFeatureSet( SeuratObject, pattern = "^RP[SL]", col.name = "percent_ribo" ) ``` ## 유전자 유형 및 염색체 위치 정보 RNA 시퀀싱에서 얻은 유전자는 여러 유형으로 구분할 수 있습니다. 이에는 단백질 발현 정보를 포함하는 Coding 유전자, 단백질 정보를 포함하지 않는 Non-coding 유전자, BCR 및 TCR 과 관련된 VDJ 영역 유전자 (세포 분화에 대한 정보를 포함), 그리고 siRNA 등이 포함됩니다. 이러한 정보를 활용하면 원하는 분석 유형에 따라 관심 없는 특정 유전자 범주를 필터링할 수 있습니다. 특히, scRNA-seq 는 일반적으로 `poly-A` 를 사용하여 RNA 를 농축하기 때문에 데이터에 존재하는 유전자의 약 80-90% 가 단백질 코딩 유전자입니다. 또한, 유전자의 염색체 위치 정보를 수집하면 성 염색체에 의한 효과를 식별하는 데 유용할 수 있습니다. ```r # 사람 유전자로 변경하기 위해 먼저 ENSEMBL에서 마우스 유전자 주석 검색 library(biomaRt) # ENSEMBL 마트 설정 mart <- biomaRt::useMart( biomart = "ensembl", dataset = "hsapiens_gene_ensembl", host= "www.ensembl.org" ) # 선택된 속성으로 마우스 유전자 주석 검색 annot <- biomaRt::getBM( mart = mart, attributes = c("external_gene_name", "gene_biotype", "chromosome_name") ) ``` 유전자 이름을 해당 생물 유형과 일치시킵니다. ```r # 유전자 이름과 해당 염색체 위치 일치 item <- annot[match(rownames(SeuratObject), annot[, 1]), "chromosome_name"] # NA 값을 "unknown"으로 대체 item[is.na(item)] <- "unknown" # 유효하지 않은 염색체 값을 "other"로 대체 item[!item %in% as.character(c(1:23, "X", "Y", "MT"))] <- "other" ``` 예를 들어 단백질 코딩 유전자에만 초점을 맞추려면 다음과 같이 수행 할 수 있습니다. ```r # SeuratObject의 차원 확인 dim(SeuratObject) # 유전자 주석에서 단백질 코딩 유전자 선택 sel <- annot[match(rownames(SeuratObject), annot[, 1]), 2] == "protein_coding" # 선택된 유전자 이름 가져오기 genes_use <- rownames(SeuratObject)[sel] genes_use <- as.character(na.omit(genes_use)) # SeuratObject에서 단백질 코딩 유전자만 필터링 SeuratObject <- SeuratObject[genes_use, ] # 필터링 후 SeuratObject의 차원 확인 dim(SeuratObject) ``` # 세포 주기에 대한 정보 각각의 세포가 어떤 주기에 해당하는지는 특정 유전자의 평균 발현을 계산하여 예측합니다. **Seurat** 은 `CellCycleScoring` 함수를 통해 이것은 쉽게 계산할 수 있습니다 (사람 유전자에만 해당). 실행 방법은 다음과 같습니다. ```r # CellCycleScoring를 실행하기 전에 데이터를 정규화하고 로그 변환해야 합니다. SeuratObject <- NormalizeData(SeuratObject) # Cell Cycle 점수 계산 SeuratObject <- CellCycleScoring( object = SeuratObject, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes ) ``` # PCA 플롯에 메타데이터 시각화 메타데이터를 PCA 플롯에서 시각화하는 것은 각각의 메타데이터 변수를 개별적으로 분석하는 것보다 쉽습니다. 모든 정보를 하나의 플롯으로 결합하기 위해 주성분 분석 (PCA) 을 사용할 것입니다. 상위 주성분 (`PC1` 과 `PC2`) 에 대한 시각화는 데이터 세트 간의 차이가 어떻게 반영되는지 직관적으로 보여줍니다. ```r # 선택된 메타데이터 매개 변수를 사용하여 PCA 계산 metadata_use <- grep( "perc", colnames(SeuratObject@meta.data), value = TRUE ) metadata_use <- c( "nCount_RNA", "nFeature_RNA", "S.Score", "G2M.Score", metadata_use ) # PCA 실행 PC <- prcomp( SeuratObject@meta.data[, metadata_use], center = TRUE, scale. = TRUE ) # 객체에 PCA(메타데이터에서 실행됨) 추가 SeuratObject@reductions[["pca_metadata"]] <- CreateDimReducObject( embeddings = PC$x, key = "metadataPC_", assay = "RNA" ) # 메타데이터에서 실행된 PCA 플롯 DimPlot( SeuratObject, reduction = "pca_metadata", dims = c(1, 2), group.by = "orig.ident" ) ``` :::{.callout-note} 제로 카운트 유전자와 모든 값이 0 인 데이터 열은 PCA 와 같은 분석에 영향을 줄 수 있어서 이런 데이터를 제외하기 위해 아래 코드를 사용합니다. ::: ```r # 제로 카운트 유전자 제거 SeuratObject <- SeuratObject[Matrix::rowSums(SeuratObject) > 0, ] # 모든 0인 메타데이터 열 제외 metadata_use <- metadata_use[ colSums(SeuratObject@meta.data[, metadata_use] != 0) > 0 ] ``` # 데이터 정규화 및 회귀 추가 분석을 수행하기 전에 서로 다른 시퀀싱 깊이를 고려하여 데이터를 정규화해야 합니다. 일반적으로 로그 변환을 통해 데이터를 정규화하며, 세포 주기, 미토콘드리아 비율, 감지된 유전자 수와 같은 데이터 품질 관리를 위한 메타데이터는 제거하는 것이 좋습니다. ## 데이터 정규화 scRNA-seq 에서 가장 일반적인 정규화 방법은 로그 정규화입니다. 이 방법은 각 유전자의 카운트를 모든 유전자 카운트의 합 (라이브러리 크기라고도 함) 으로 나누어 라이브러리 크기의 차이를 보상합니다. 그 후, 결과에 상수를 곱하여 모든 세포가 동일한 시퀀싱 깊이를 가지도록 합니다. 대부분의 bulk RNA-seq 에서는 상수가 보통 `1e6` 이며, 이로 인해 CPM(백만 단위 카운트) 이 생성됩니다. 그러나 단일 세포 라이브러리 크기는 그보다 훨씬 작기 때문에 `1e3` 에서 `1e4`(10,000 단위 카운트) 까지 사용됩니다. $$ NormCounts = \frac{(GeneCounts \times 10000)}{LibrarySize} $$ 라이브러리 크기에 대해 보정된 값은 로그 변환하여 로그 정규 분포로 만듭니다. $$ logNormCounts = \ln(NormCounts + 1) $$ **Seurat**에서는 `NormalizeData` 함수를 사용하여 이 과정을 수행할 수 있습니다. ```r # 분산이 0인 유전자 제거 SeuratObject <- SeuratObject[Matrix::rowSums(SeuratObject) > 0, ] # 데이터 정규화 (LogNormalize 사용) SeuratObject <- NormalizeData( object = SeuratObject, scale.factor = 10000, normalization.method = "LogNormalize" ) ``` ## 데이터 회귀 (Regression) 데이터에서 혼란 요인을 제거하기 위해 회귀 분석을 수행할 수 있습니다. 예를 들어, 세포 주기나 품질 지표 (예: 미토콘드리아 비율 또는 감지된 유전자 수) 와 같은 혼란 요인을 제거할 수 있습니다. **Seurat**에서는 `ScaleData` 함수를 사용하여 이러한 회귀를 수행할 수 있습니다. ```r SeuratObject <- ScaleData( SeuratObject, vars.to.regress = c("nFeature_RNA", "percent_mito"), verbose = FALSE ) ``` 이렇게 하면 데이터에 회귀가 적용되며, 이제 데이터가 준비되었으므로 추가 분석을 진행할 수 있습니다. # highly variable features 선택하기 대규모 데이터 분석에서 중요한 단계는 샘플 간에 현저한 차이를 보이는 features(유전자, 전사체, 단백질, 대사물질 등) 를 식별하는 것입니다. 예를 들어, 다양한 종류의 상피 세포가 포함된 데이터셋을 가정해 보겠습니다. 상피 세포의 하위 유형을 구별하기 위해 다음과 같은 전략을 사용할 수 있습니다: 1. 모든 상피 세포에서 공통적으로 발현되며 수준이 거의 동일한 유전자만 사용합니다. 2. 상피 세포에서 감지되지 않는 유전자만 사용합니다. 3. 상피 세포 간에 발현 수준이 크게 다른 유전자만 사용합니다. 4. 데이터셋의 모든 유전자를 사용합니다. 상피 세포 간 발현이 크게 다른 유전자를 사용하는 것이 가장 효과적인 방법이며, 그 다음으로 모든 유전자를 사용하는 것이 좋습니다. 반면, 상피 세포 간에 균일하게 발현되는 유전자나 상피 세포에서 감지되지 않는 유전자를 사용하는 것은 하위 유형을 구별하는 데 충분한 정보를 제공하지 않을 수 있습니다.그러나 단일 세포 분석에서는 보통 세포의 상피 세포 유형을 사전에 알지 못합니다 (그것이 발견하려는 대상이기 때문입니다). 따라서 이 작업을 수행하기 위해 다른 방법이 필요합니다. 일반적인 접근법은 각 샘플 전체에서 유전자의 분산을 기준으로 유전자를 순위 지정하는 것입니다. 분산이 높은 유전자는 세포를 구분하는 데 효과적일 가능성이 높습니다.일반적으로 발현 수준이 높은 유전자는 자연스럽게 더 높은 변동성을 나타내므로, 각 유전자의 평균 로그 발현량으로 변동성을 정규화하는 것이 중요합니다. ```r # Variable feature 찾기 SeuratObject <- FindVariableFeatures( object = SeuratObject, nfeatures = 3000, selection.method = "vst", verbose = FALSE, assay = "RNA" ) # 상위 20개의 Variable feature 시각화 top20 <- head(VariableFeatures(SeuratObject), 20) LabelPoints( plot = VariableFeaturePlot(SeuratObject), points = top20, repel = TRUE ) ``` # 스케일링 및 중심화 ## (선형) 스케일링 및 중심화 각 유전자는 서로 다른 발현 수준을 가지므로, 높은 발현 값을 가진 유전자는 자연스럽게 더 높은 변동성을 보입니다. 이는 후속 분석에서 잘 포착될 수 있습니다. 따라서 각 유전자에 대해 유사한 가중치를 부여하는 것이 중요합니다. 일반적으로 PCA 를 수행하기 전에 각 유전자를 중심화하고 스케일링하는 것이 권장됩니다. 이러한 정밀한 스케일링 방법은 Z- 점수 정규화라고 하며, PCA, 클러스터링 및 히트맵 시각화에 매우 유용합니다. 또한, 세포 주기, 시퀀싱 깊이, 미토콘드리아 비율 등과 같은 데이터셋에서 원하지 않는 변동 소스를 제거하기 위해 회귀 분석을 사용할 수 있습니다. 이는 이러한 매개변수를 모형의 공변량으로 사용하여 일반화된 선형 회귀 (GLM) 를 수행함으로써 달성됩니다. 이후 모델의 잔차는 "회귀된 데이터"로 간주됩니다. ```r SeuratObject <- ScaleData( object = SeuratObject, vars.to.regress = c("nCount_RNA", "percent_mito", "nFeatures_RNA"), model.use = "linear", assay = "RNA", do.scale = TRUE, do.center = TRUE ) ``` Seurat 은 기본적으로 Seurat 객체에 있는 모든 변수 기능에 대해 스케일링을 실행합니다. 변수 기능이 있을 경우 이러한 스케일링 단계의 성능을 크게 향상시킬 수 있습니다. 그러나 결과는 변수 유전자 자체에 기반하여 PCA, 차원 축소 또는 클러스터링과 같은 후속 분석의 결과가 변경되지 않습니다. 전체 유전자 수에 대해 데이터를 스케일링하려는 경우 `FindVariableFeatures` 를 사용한 후 함수 호출에서 이를 지정할 수 있습니다. ```r SeuratObject <- ScaleData( object = SeuratObject, vars.to.regress = c("nCount_RNA", "percent_mito", "nFeatures_RNA"), model.use = "linear", assay = "RNA", do.scale = TRUE, do.center = TRUE, features = rownames(SeuratObject) ) ``` ## (포아송) 스케일링 및 중심화 위의 선형 스케일링과 중심화는 로그 - 선형 데이터 분포를 가정하기 때문에, RNA-seq 데이터 (단일 세포 포함) 는 보다 음이항 분포에 가까울 수 있으며, 이로 인해 변이를 올바르게 회귀하지 못할 수 있습니다. 위의 절차의 대안으로 "포아송" 또는 "음이항" 분포를 사용하여 원시 UMI 카운트 데이터에 대해 실행할 수 있습니다. 이는 포아송 모델을 사용하여 유전자별 GLM 회귀를 수행하는 방법입니다. ```R SeuratObject <- ScaleData( object = SeuratObject, vars.to.regress = c("nCount_RNA", "mito.percent", "nFeatures_RNA"), model.use = "poisson", assay = "RNA", do.scale = TRUE, do.center = TRUE ) ``` # SCTransform 포아송 분포를 기반으로 한 스케일링과 중심화는 경우에 따라 데이터를 과적합시킬 수 있습니다. 이러한 문제를 극복하기 위해 유사한 풍부도를 가진 유전자 간의 정보를 통합하여 회귀 모델에서 안정된 매개변수 추정을 얻을 수 있습니다. 이 방법을 `scTransform` 이라고 하며, 간단히 말해 제한된 음이항 모델을 사용하여 유전자별 일반 선형 모델 회귀를 수행하는 방식입니다. ```r SeuratObject <- SCTransform( object = SeuratObject, assay = "RNA", vars.to.regress = "percent_mito", new.assay.name = "sctransform", do.center = TRUE, verbose = FALSE ) ``` # 차원 축소하기 scRNA-seq 데이터는 수천 개의 유전자로 구성된 전체 유전자 발현 공간을 포함하고 있으며, 이 공간은 상당한 노이즈를 포함하고 있어 시각화가 어렵습니다. 따라서 대부분의 scRNA-seq 분석에서는 데이터의 일부 변동을 제거하기 위해 PCA(주성분 분석) 또는 ICA(독립 성분 분석) 와 같은 방법으로 시작합니다.단순한 scRNA-seq 데이터셋에서는 몇 가지 세포 유형만 포함되어 있는 경우 PCA 가 데이터를 2 차원 또는 3 차원으로 시각화하는 데 충분할 수 있습니다. 그러나 데이터의 복잡성이 증가함에 따라 비선형 차원 축소 방법을 사용하여 데이터를 2 차원으로 투영해야 할 필요가 있습니다. 이러한 비선형 차원 축소 방법에는 t-SNE, UMAP, 확산 맵 등이 있습니다.어떤 차원 축소 방법을 선택할지는 간단한 문제가 아니며, 여러 가지 적절한 답변이 있을 수 있는 복잡한 문제입니다. 더 자세한 내용을 알고 싶다면 다음 [논문](https://www.biorxiv.org/content/10.1101/120378v4) 을 참고하시기 바랍니다. ## PCA 주성분 분석 (PCA) 은 데이터를 새로운 좌표계로 변환하는 직교 선형 변환입니다. 이 과정에서 데이터의 가장 큰 분산이 첫 번째 주성분에 위치하고, 두 번째로 큰 분산은 두 번째 주성분에 위치하도록 합니다. 이러한 과정을 통해 데이터의 분산을 가장 잘 설명하는 방식으로 내부 구조를 드러냅니다. 일반적으로 이 동작은 데이터의 차원을 줄여 변환된 데이터의 차원을 축소하는 방식으로 이해될 수 있습니다. ```R # 주성분 분석 실행 SeuratObject <- RunPCA( object = SeuratObject, assay = "RNA", npcs = 100, verbose = FALSE ) ``` ## t-SNE t- 분포 확률적 이웃 임베딩 (t-distributed stochastic neighbor embedding, t-SNE) 은 고차원 데이터를 저차원 공간으로 임베딩하여 시각화하는 데 적합한 비선형 차원 축소 기술입니다. 이 방법은 각 고차원 객체를 유사한 객체와 가까운 점으로 모델링하고, 서로 다른 객체는 먼 확률적 점으로 모델링하여 두 개 또는 세 개의 차원으로 표현합니다.[^2] ```R # tSNE 실행 SeuratObject <- RunTSNE( object = SeuratObject, reduction = "pca", perplexity = 30, max_iter = 1000, theta = 0.5, eta = 200, exaggeration_factor = 12, dims.use = 1:50, verbose = FALSE, num_threads = 0 ) ``` ## UMAP UMAP(Uniform Manifold Approximation and Projection) 은 t-SNE 와 유사하게 시각화를 위한 차원 축소에 사용할 수 있는 기술이며, 일반적인 비선형 차원 축소에도 활용될 수 있습니다.[^3] UMAP 알고리즘은 실제 데이터에 적용할 수 있는 실용적이고 확장 가능한 방법입니다. 이 알고리즘은 시각화 품질에서 t-SNE 와 경쟁할 만큼 우수하며, 더 나은 실행 시간 성능을 바탕으로 전역 구조를 더 잘 보존한다고 알려져 있습니다. 또한 UMAP 은 임베딩 차원에 대한 계산 제한이 없어 머신러닝을 위한 일반적인 목적의 차원 축소 기술로 널리 사용될 수 있습니다.[^4] ```r SeuratObject <- RunUMAP( object = SeuratObject, reduction = "pca", dims = 1:50, n.components = 2, n.neighbors = 20, repulsion.strength = 1, verbose = FALSE, n.epochs = 200, metric = "euclidean", seed.use = 42, reduction.name = "umap" ) ``` ## 확산 맵 (Diffusion Maps) 확산 맵 (Diffusion Maps, DM) 은 데이터 집합을 유클리드 공간 (일반적으로 저차원) 으로 임베딩하는 차원 축소 기법 중 하나입니다. 이 기법은 데이터에 대한 확산 연산자의 고유 벡터와 고유값을 기반으로 임베딩을 계산합니다. 임베딩된 공간에서의 점들 사이의 유클리드 거리는 해당 점들을 중심으로 하는 확률 분포 간의 "확산 거리"와 동일합니다. 이 방법은 주어진 데이터에서 랜덤 워크를 수행할 때, 인접한 데이터 포인트로 이동하는 것이 멀리 떨어진 데이터 포인트로 이동하는 것보다 더 자주 발생한다는 기본적인 관찰에 기반하여 작동합니다.확산 맵은 주로 데이터가 추출된 기저 매니폴드를 발견하는 데 중점을 두고 있으며, 이는 주성분 분석 (PCA) 이나 다차원 스케일링 (MDS) 과 같은 선형 차원 축소 방법과는 다른 접근 방식입니다.[^5] ```r # 추가 라이브러리 로드 library(destiny) # destiny 패키지를 사용하여 확산 맵 실행 dm <- DiffusionMap( data = SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:50], k = 20, n_eigs = 20) # DM 임베딩에서 셀 이름 수정 rownames(dm@eigenvectors) <- colnames(SeuratObject) # SeuratObject에 DM 임베딩 추가 SeuratObject@reductions[["dm"]] <- CreateDimReducObject( embeddings = dm@eigenvectors, key = "DC_", assay = "RNA" ) ``` ## 독립 성분 분석 (ICA) 독립 성분 분석 (Independent Component Analysis, ICA) 은 다변량 신호를 가산성 하위 구성 요소로 분리하는 계산 방법입니다. 이 기법은 하위 구성 요소가 가우시안이 아닌 신호이며 서로 통계적으로 독립적이라는 가정을 기반으로 합니다. ICA 는 블라인드 소스 분리 (blind source separation) 의 특별한 경우로 간주됩니다.[^6] ```r SeuratObject <- RunICA( object = SeuratObject, assay = "pca", nics = 20, reduction.name = "ica" ) ``` # 그래프 만들기 scRNA-seq 데이터의 클러스터링을 전체 발현 데이터 혹은 PCA 공간 (선형 거리를 제공) 에서 하는 대신 그래프 (셀 간 유사성에 대한 비선형 표현) 를 먼저 만들고 클러스터링을 하는 것이 매우 효과적입니다. 그래프는 모든 셀 (노드 또는 정점) 을 나타내며, 일정한 유사성 기준에 따라 이들 사이에 에지를 연결합니다. 예를 들어, 50 개의 주요 구성 요소를 사용하여 PCA 공간에서 서로 X 거리 이하인 모든 셀 사이에 에지를 갖는 그래프를 구성할 수 있습니다. ## KNN KNN 은 "K Nearest Neighbors"의 약자로, 데이터 마이닝 및 기계 학습 분야에서 기본적이고 인기 있는 주제입니다. KNN 그래프는 두 정점 `p` 와 `q` 가 연결되는 그래프로, 그 거리가 K 번째로 작은 거리 중 하나인 경우에만 연결됩니다. 이때 쌍별 거리 측정치는 해밍 거리, 코사인 거리, 유클리드 거리 등 다양한 방법을 사용할 수 있습니다. 본 논문에서는 벡터 간 유사성을 측정하기 위해 유클리드 거리를 사용합니다. KNN 그래프 데이터 구조는 데이터 마이닝에서 많은 이점을 제공합니다. 예를 들어, 수십억 개의 데이터 세트에서는 KNN 그래프를 오프라인으로 사전 구축하여 인덱스로 사용하는 것이 여러 번의 온라인 KNN 검색보다 훨씬 효율적입니다. ```r SeuratObject <- FindNeighbors( SeuratObject, assay = "RNA", compute.SNN = FALSE, # FALSE로 설정하면 KNN 그래프만 계산됩니다. reduction = "pca", dims = 1:50, graph.name = "nn", prune.SNN = 1/15, k.param = 20, force.recalc = TRUE ) ``` KNN 그래프는 모든 셀 간의 연결이 1 로 표시된 행렬로 나타납니다. 이를 가중치가 없는 그래프라고 하며, Seurat 의 기본 설정입니다. 그러나 일부 셀 간의 연결은 다른 연결보다 더 중요할 수 있으며, 이 경우 그래프의 가중치는 0 에서 최대 거리까지 다양해집니다. 일반적으로 거리가 작을수록 두 점이 가까워지고, 그들의 연결이 강해집니다. 이러한 경우를 가중치가 있는 그래프라고 합니다. 가중치가 있는 그래프와 가중치가 없는 그래프 모두 클러스터링에 적합하지만, 대규모 데이터 세트 (100,000 개 이상의 세포) 에서는 가중치가 없는 그래프에서 클러스터링하는 것이 더 빠릅니다. ```r library(pheatmap) pheatmap( SeuratObject@graphs$nn[1:100, 1:100], col = c("white", "black"), border_color = "grey90", legend = FALSE, cluster_rows = FALSE, cluster_cols = FALSE, fontsize = 2 ) ``` ## SNN 추가로 KNN 그래프 외에도 모든 두 점 간에 공유된 최근접 이웃의 수를 계산합니다. 이를 통해 "공유된 최근접 이웃" 그래프를 형성합니다. 이 과정에서는 두 점 간의 각 링크의 가중치를 해당 점들이 공유하는 이웃의 수로 대체합니다. 다시 말해, 이는 최근접 이웃 그래프 내에서 모든 두 점 사이의 길이가 2 인 경로의 수를 나타냅니다.공유된 최근접 이웃 그래프가 생성되면 모든 점 쌍이 비교됩니다. 두 점이 T 이상의 이웃을 공유하는 경우, 즉 공유된 최근접 이웃 그래프에서 가중치가 임계값 T 보다 많은 링크를 가지면, 두 점과 그들이 속한 모든 클러스터가 병합됩니다. 다시 말해, 클러스터는 임계값을 사용하여 희소화한 후 공유된 최근접 이웃 그래프에서의 연결 요소입니다. ```r SeuratObject <- FindNeighbors( SeuratObject, assay = "RNA", compute.SNN = TRUE, # KNN 그래프와 SNN 그래프 둘 다 계산합니다. reduction = "pca", dims = 1:50, graph.name = "SNN", prune.SNN = 1/15, k.param = 20, force.recalc = TRUE ) ``` # 여러 데이터셋 합치기 및 배치 보정하기 기존의 일괄 보정 방법은 주로 대량 RNA-seq 을 위해 설계되었습니다. 따라서 이러한 방법을 scRNA-seq 데이터에 적용할 때는 각 일괄 처리 내의 세포 모집단 구성이 동일하다고 가정합니다. 일괄 간의 평균 유전자 발현 사이에 시스템적인 차이가 있을 경우, 이는 기술적인 차이로 귀속됩니다. 그러나 실제로는 scRNA-seq 연구에서 일괄 처리별로 모집단 구성이 동일하지 않을 수 있습니다. 각 일괄 처리에 동일한 세포 유형이 존재한다고 하더라도, 세포 유형의 풍부도는 세포 배양, 조직 추출, 분리 및 분류 등의 미세한 차이에 따라 달라질 수 있습니다. 따라서 일괄 차단 요인에 대한 추정 계수는 순수하게 기술적인 것이 아니라 구성의 차이로 인한 비생물학적 요소를 포함하게 됩니다. 이러한 계수를 기반으로 한 일괄 보정은 세포의 발현 프로파일을 부정확하게 표현할 수 있으며, 보정을 수행하지 않은 경우보다 결과가 더 나쁠 수 있습니다.[^7] 많은 다른 scRNA-seq 방법론과 마찬가지로, 데이터셋에 대해 어떤 일괄 보정 기술을 사용해야 할지 선택하기 어려울 수 있습니다.[^8] ## Scanorama 사용하기 :::{.callout-note} 최근 여러 방법들이 데이터셋을 성공적으로 통합할 수 있다는 것을 보여주고 있습니다. 그러나 이러한 접근법들은 데이터셋에 적어도 하나의 공통 세포 유형이 존재한다고 가정하기 때문에, 세포 구성이 상당히 다른 데이터셋을 결합할 경우 과도한 보정이 발생할 수 있습니다. ::: `Scanorama` 의 접근 방식은 두 데이터 세트 간의 유사한 요소를 찾는 상호 최근접 이웃 일치 기술을 일반화하여 여러 데이터 세트 간에 유사한 요소를 탐색합니다. 이 방식은 휴대폰으로 파노라마 사진을 촬영하는 기술과 유사합니다. `Scanorama` 는 자동으로 유사한 전사 프로필을 가진 세포를 식별하고 이를 일괄 보정 및 통합에 활용합니다. 이 방법은 강력한 특징을 가지고 있으며, 데이터셋의 크기나 출처에 크게 영향을 받지 않으며, 모든 데이터셋이 적어도 하나의 공통 세포 모집단을 공유할 필요가 없습니다.[^9] ```r library(reticulate) scanorama <- import("scanorama") # 오브젝트를 배치 별로 분할 SeuratObject.list <- SplitObject( object = SeuratObject, split.by = "BATCH" ) # 어레이 및 제네릭 데이터 만들기 assaylist <- list() genelist <- list() for (i in seq_along(SeuratObject.list)) { assaylist[[i]] <- t( as.matrix(GetAssayData(SeuratObject.list[[i]], "data")) ) genelist[[i]] <- rownames(SeuratObject.list[[i]]) } # 데이터 통합 및 보정 integrated.data <- scanorama$integrate(assaylist, genelist) corrected.data <- scanorama$correct( assaylist, genelist, return_dense = TRUE ) integrated.corrected.data <- scanorama$correct( assaylist, genelist, return_dimred = TRUE, return_dense = TRUE ) # 데이터 결합 및 이름 지정 intdata <- lapply(integrated.corrected.data[[2]], t) panorama <- do.call(cbind, intdata) rownames(panorama) <- as.character(integrated.corrected.data[[3]]) colnames(panorama) <- unlist(sapply(assaylist, rownames)) # 차원 축소 결과 결합 intdimred <- do.call(rbind, integrated.corrected.data[[1]]) colnames(intdimred) <- paste0("PC_", 1:100) # Seurat에서 Elbow Plot을 그리기 위해 표준 편차 추가 stdevs <- apply(intdimred, MARGIN = 2, FUN = sd) # 새로운 어레이 오브젝트 만들기 SeuratObject@assays[["pano"]] <- CreateAssayObject( data = panorama, min.cells = 0, min.features = 0 ) # PCA 결과 추가 SeuratObject[["pca_scanorama"]] <- CreateDimReducObject( embeddings = intdimred, stdev = stdevs, key = "PC_", assay = "pano" ) ``` ## Harmony 사용하기 Harmony 는 주어진 저차원 임베딩 (예: 주성분 분석) 을 기반으로 시작합니다. Harmony 는 먼저 이 임베딩을 사용하여 셀을 여러 데이터 세트 클러스터로 그룹화합니다. 부드러운 전이를 고려하기 위해 소프트 클러스터링을 사용하여 셀을 잠재적으로 여러 클러스터에 할당합니다. 이러한 클러스터는 이산적인 세포 유형을 식별하기 위한 것이 아니라, 세포 상태 간의 부드러운 전환을 반영하는 대리 변수로 기능합니다. 클러스터링이 완료되면 각 데이터셋에는 클러스터별 중심점이 생성됩니다. 이 중심점은 클러스터별 선형 보정 계수를 계산하는 데 사용됩니다. 클러스터가 세포 유형과 상태에 해당하므로, 클러스터별 보정 계수는 개별 세포 유형 및 세포 상태에 특정한 보정 계수로 작용합니다. 이를 통해 Harmony 는 미묘한 세포 표현형에 민감한 간단한 선형 조정 함수를 학습합니다. 마지막으로, 각 세포는 이러한 클러스터의 가중 평균을 할당받고, 세포별 선형 요소에 의해 보정됩니다. 각 세포가 여러 클러스터에 속할 수 있기 때문에, 각 세포에는 잠재적으로 고유한 보정 계수가 부여됩니다. Harmony 는 셀 클러스터 할당이 안정될 때까지 이 네 단계를 반복합니다. Harmony 는 대규모 데이터셋으로 확장 가능하며, 다양한 인구 집단과 세부적인 하위 집단을 동시에 식별할 수 있는 유연성을 갖추고 있습니다. 또한 복잡한 실험 설계를 수용할 수 있으며, 다양한 모드 간 통합이 가능한 강력한 알고리즘입니다.[^10] ```r # 추가 라이브러리 로드 library(harmony) library(SeuratWrappers) # Harmony 실행 SeuratObject <- RunHarmony( SeuratObject, group.by.vars = "batch" ) ``` # 클러스터링 하기 scRNA-seq 연구의 주요 목표 중 하나는 유사한 세포들을 클러스터로 그룹화하는 것입니다. 클러스터를 정의하는 방법은 여러 가지가 있으며, K- 평균, 계층적 클러스터링, 친화 전파와 같은 전통적인 접근 방식이 있지만, 최근 대부분의 파이프라인에서는 그래프를 사용하여 클러스터를 상호 연결된 세포 그룹으로 정의합니다 (위의 섹션 참조).그래프에서 클러스터를 정의하기 위해 Louvain, Leiden, Infomap 과 같은 다양한 커뮤니티 탐지 알고리즘이 사용됩니다. 이들 알고리즘의 주요 아이디어는 데이터셋 내에서 서로 더 많은 엣지 (높은 유사도의 연결) 를 가진 세포 그룹을 찾는 것입니다.모든 클러스터링 방법은 세포 간의 페어와이즈 거리를 정의하는 것으로 시작합니다. 일반적으로 이는 데이터를 스케일링하고 변수를 선택한 후 PCA 공간에서 세포 간의 유클리드 거리를 계산함으로써 수행됩니다. 따라서 클러스터링 결과는 데이터 전처리 단계에서의 선택에 크게 의존하게 됩니다. 특히 변수 유전자 선택과 PCA 에서 포함할 주성분의 수가 중요한 요소입니다.그래프 기반 방법에서는 PCA 공간에서의 거리를 사용하여 이웃하는 세포 사이에 엣지가 있는 그래프를 생성합니다. 이때 KNN 그래프의 이웃 수와 SNN 구성을 위한 가지치기 설정과 같은 여러 매개변수가 클러스터링 결과에 영향을 미칩니다.[^11] ## 루브앵 (Louvain) 루브앵 (Louvain) 방법은 Blondel 이 개발한 대규모 네트워크에서 커뮤니티를 추출하는 기법입니다.[^12] 이 방법은 네트워크 노드 수에 대해 $O(n.log2n)$ 의 시간 복잡도를 가지는 탐욕적 최적화 방법입니다. 최적화하는 값은 모듈러리티 (modularity) 로 이는 커뮤니티 내 링크 밀도를 커뮤니티 간 링크와 비교하여 측정하는 범위 내 값입니다. 이 값을 최적화하는 것은 이론적으로 주어진 네트워크의 노드를 가장 잘 그룹화하는 결과를 도출하지만 모든 노드를 그룹으로 나누는 모든 가능한 반복 과정을 거치는 것은 비실용적이므로 휴리스틱 알고리즘이 사용됩니다. ```r SeuratObject <- FindClusters( object = SeuratObject, graph.name = "SNN", resolution = 0.8, algorithm = 1 # 1로 설정하면 Louvain을 의미합니다. ) ``` 클러스터 수는 `resolution` 매개변수를 사용하여 제어할 수 있으며 값이 높을수록 더 많은 (작은) 클러스터가 생성됩니다. ## 라이든 (Leiden) 라이든 알고리즘은 반복적으로 적용되어 모든 커뮤니티의 하위 집합이 지역적으로 최적화된 분할로 수렴합니다.[^13] 이 알고리즘은 빠른 로컬 이동 접근 방식을 기반으로 하여 루브앵 알고리즘보다 더 빠르게 실행됩니다. 라이든 알고리즘은 다음과 같은 세 단계로 구성됩니다: 1. **노드의 로컬 이동**: 네트워크 내의 각 노드는 주변 노드와의 연결을 고려하여 위치를 조정합니다. 이 과정에서 각 노드는 자신의 이웃 중에서 가장 높은 모듈러리티를 제공하는 위치로 이동하게 됩니다. 이를 통해 각 노드가 속한 커뮤니티가 더욱 명확해집니다. 2. **분할의 정제**: 초기 클러스터링 결과를 바탕으로, 각 커뮤니티 내에서 노드의 재배치를 통해 더 나은 분할을 찾습니다. 이 단계에서는 커뮤니티 간의 경계를 명확히 하고, 중복된 연결을 줄이는 과정이 포함됩니다. 3. **정제된 분할을 바탕으로 네트워크를 집계**: 정제된 분할을 사용하여 새로운 집계 네트워크를 생성합니다. 이때, 정제되지 않은 분할을 기반으로 초기 집계 네트워크의 분할을 생성하여, 전체 네트워크 구조를 반영합니다. 이러한 과정을 통해 라이든 알고리즘은 데이터의 구조를 효과적으로 파악하고, 커뮤니티 간의 관계를 명확히 하여 최적화된 결과를 도출합니다. 라이든 알고리즘은 대규모 네트워크에서도 효율적으로 작동하며, 다양한 분야에서 널리 사용되고 있습니다. ```R SeuratObject <- FindClusters( object = SeuratObject, graph.name = "SNN", resolution = 0.8, algorithm = 4 # 4로 설정하면 Leiden을 의미합니다. ) ``` 클러스터 수는 `resolution` 매개변수를 사용하여 제어할 수 있으며 값이 높을수록 더 많은 (작은) 클러스터가 생성됩니다. ## 계층적 클러스터링 계층적 클러스터링 (Hierarchical Clustering, HC) 은 클러스터의 계층 구조를 구축하는 클러스터 분석 방법입니다. 이 방법은 일반적으로 두 가지 전략으로 나뉘며, 각각 합병적 (agglomerative) 또는 분할적 (divisive) 입니다. 병합과 분할은 보통 탐욕적인 방식으로 결정됩니다. 계층적 클러스터링의 결과는 일반적으로 덴드로그램 (dendrogram) 으로 시각화됩니다. 계층적 병합 클러스터링 (Hierarchical Agglomerative Clustering, HAC) 의 표준 알고리즘은 $O(n^3)$ 의 시간 복잡도를 가지며 $O(n^2)$ 의 메모리를 필요로 하기 때문에, 중간 규모의 데이터 세트에서도 느릴 수 있습니다. 병합적인 경우에는 어떤 클러스터를 결합할 것인지, 분할적인 경우에는 어떤 클러스터를 분할할 것인지를 결정하기 위해 관측값 집합 간의 불일치 측정이 필요합니다. 대부분의 계층적 클러스터링 방법에서는 적절한 메트릭 (관측값 쌍 간의 거리 측정) 과 연결 기준을 사용하여 관측값 집합의 불일치를 관측값 쌍 간 거리의 함수로 정의합니다. `R` 의 기본 `stats` 패키지에는 모든 샘플 쌍 사이의 거리를 계산하는 `dist()` 함수가 포함되어 있습니다. `dist()` 에서 사용할 수 있는 거리 방법은 다음과 같습니다: - 유클리드 거리 (euclidean) - 최대 거리 (maximum) - 맨하탄 거리 (manhattan) - 캔버라 거리 (canberra) - 이진 거리 (binary) - 민코프스키 거리 (minkowski) 또한, 이미 셀 간 거리 정보를 포함하는 그래프 (KNN 또는 SNN) 를 사용하여 직접 계층적 클러스터링을 수행할 수도 있습니다. 그러나 그래프의 거리는 반전되어 있으므로 (0 은 멀리 떨어져 있음을 나타내고 1 은 가깝다는 것을 나타냄), 그래프의 최대값 (인접 SNN 의 경우 1) 을 뺀 후 0 이 가까운 거리를 나타내고 1 이 먼 거리를 나타내도록 조정해야 합니다. 샘플 간의 거리를 계산한 후에는 본격적인 계층적 클러스터링을 진행할 수 있습니다. 이를 위해 `hclust()` 함수를 사용하며, 위에서 생성한 거리 객체를 활용하여 간단히 실행할 수 있습니다. 사용 가능한 방법에 대한 자세한 내용은 [HC for networks](https://en.wikipedia.org/wiki/Hierarchical_clustering_of_networks) 에서 확인할 수 있습니다. ```r # PCA에서 HC 실행 h_pca <- hclust( d = dist( SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:30], method = "euclidean" ), method = "ward.D2" ) # 그래프에서 HC 실행 h_graph <- hclust( d = 1 - as.dist(SeuratObject@graphs$SNN), method = "ward.D2" ) ``` 클러스터 계층이 정의되면, 다음 단계는 특정 클러스터에 속하는 샘플을 식별하는 것입니다. ```R plot(h, labels = FALSE) ``` 세포에 대한 덴드로그램을 확인한 후, 고정된 높이 (threshold) 에서 나무 (tree) 를 자를 수 있습니다. 이 작업은 `cutree` 함수를 사용하여 수행할 수 있습니다. 또한, 해상도에 따라 다른 수준에서 나무를 자르는 것도 가능합니다. ```R # 덴드로그램 높이를 기반으로 나무 자르기 SeuratObject$HC_res <- cutree( tree = h, k = 18 ) # 클러스터 수를 기반으로 나무 자르기 SeuratObject$HC_res <- cutree( tree = h, h = 3 ) # 각 클러스터에 속하는 셀 수 확인 table(SeuratObject$HC_res) ``` 클러스터의 수는 높이 (h) 를 기준으로 하거나 직접 k 매개변수를 설정하여 조절할 수 있습니다. ## K- 평균 클러스터링 K- 평균 클러스터링 (k-means clustering) 은 신호 처리에서 시작된 벡터 양자화 방법으로, n 개의 관측값을 k 개의 클러스터로 분할하여 각 관측값이 평균 (클러스터 중심) 과 가장 가까운 클러스터에 속하도록 하는 것을 목표로 합니다. 이 평균은 클러스터의 프로토타입 역할을 합니다. 그러나 알고리즘은 전역 최적값으로 수렴하는 것을 보장하지 않으며, 결과는 초기 클러스터 설정에 따라 달라질 수 있습니다. 일반적으로 알고리즘은 빠르게 실행되므로 서로 다른 초기 조건으로 여러 번 실행하는 것이 일반적입니다. K- 평균 클러스터링은 주로 비슷한 공간 범위의 클러스터 (모두 같은 크기) 를 찾지만, 기대값 최대화 메커니즘 덕분에 클러스터가 다양한 형태를 가질 수 있습니다. K- 평균은 많은 응용 분야에서 널리 사용되는 클러스터링 알고리즘입니다. `R` 에서는 `kmeans` 함수를 사용하여 쉽게 적용할 수 있습니다. 일반적으로 이 알고리즘은 표현 데이터의 축소된 차원 표현 (대부분 PCA 를 사용함) 에 적용되며, 이는 저차원 거리의 해석 가능성 덕분입니다. K- 평균 클러스터링을 수행할 때는 클러스터의 수를 미리 정의해야 하며, 결과는 클러스터 중심의 초기화에 따라 달라지므로 일반적으로 다양한 시작 구성으로 k-means 를 실행하는 것이 권장됩니다 (`nstart` 인수를 통해 설정 가능). ```R set.seed(42) # 재현성을 위해 설정합니다. # K-means 클러스터링 실행 SeuratObject$kmeans_12 <- kmeans( x = SeuratObject@reductions[["pca"]]@cell.embeddings[, 1:50], centers = 12, iter.max = 50, nstart = 10 )$cluster ``` 클러스터의 수는 `centers` 매개변수를 사용하여 설정할 수 있습니다. 그렇다면 어떤 클러스터링 해상도를 선택해야 할까요?클러스터링 해상도를 결정하는 것은 간단한 문제가 아니며, 여러 방법을 통해 접근할 수 있습니다. 합리적인 클러스터 수를 정하는 것은 복잡한 과정으로, 샘플에 대한 생물학적 지식과 조사가 필요합니다.차별적 유전자 발현 (Differential gene expression) 분석은 이 과정에 큰 도움이 될 수 있습니다. 만약 두 클러스터가 동일한 차별적 유전자 발현을 보이고 구별할 만한 유전자가 없다면, 이들을 개별 클러스터로 분리하는 것은 바람직하지 않을 수 있습니다. 또한, 품질 관리 지표 (QC metrics) 를 통해 생성된 클러스터를 검토하여, 일부 클러스터가 낮은 품질의 세포나 이중체로 구성되지 않았는지 확인하는 것이 중요합니다. ## 클러스터링 결과 비교하기 클러스터링 결과를 다른 클러스터링 기법과 비교하는 것은 매우 유용합니다. 이를 다음과 같은 방식으로 수행할 수 있습니다: 위 코드는 비교 대상 클래스에 속하는 세포 수에 대한 정보를 포함하는 테이블을 생성합니다. 이 데이터를 막대 그래프로 시각화하여 클러스터링 결과를 쉽게 비교할 수 있습니다. ```r # 비교 테이블 생성 comparison_table <- table( list( SeuratObject@meta.data[,"METADATA_FACTOR_1"], SeuratObject@meta.data[,"METADATA_FACTOR_2"] ) ) # 데이터를 백분율로 변환 comparison_table <- t(t(comparison_table) / colSums(comparison_table)) * 100 # 막대 그래프 그리기 barplot( comparison_table, col = 1:nrow(comparison_table), # 색상 설정 border = NA, las = 2 # 축 라벨 방향 설정 ) ``` # 차등발현 (Differentially Expressed) 분석 클러스터가 정의된 후, 각 클러스터를 특징짓는 유전자를 찾는 것이 종종 유익합니다. scRNA-seq 에서의 차별적 발현 (DE) 분석 방법은 일반적인 RNA-seq 방법과 다소 다릅니다. 한편으로는 대량 RNA-seq 에 비해 샘플 (개별 세포) 의 수가 훨씬 많지만, 반면에 scRNA-seq 데이터는 노이즈가 많고 드롭아웃에 취약하여 분석이 복잡해집니다.[^14] ## DEG(Differential expression gene) 찾기 차별적 발현 유전자 (DEG) 는 종종 "마커 유전자"로 언급되지만, 대부분의 DE 테스트는 한 그룹의 세포와 다른 그룹의 세포 간에 발현이 더 높은 유전자를 감지하도록 설계되었다는 점을 인식해야 합니다. DEG 가 자동으로 세포 유형을 위한 고유한 마커가 되지는 않습니다.Seurat 패키지에는 DE 분석을 위한 여러 가지 테스트가 구현되어 있으며, 일부는 scRNA-seq 에 특화되어 있고, 다른 일부는 벌크 RNA-seq 에서 사용됩니다. - **wilcox**: Wilcoxon 순위 합 검정을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별 - **bimod**: 단일 세포 유전자 발현에 대한 우도 비율 검정 [^15] - **roc**: ROC 분석을 사용하여 유전자 발현의 '마커'를 식별합니다. 각 유전자에 대해 해당 유전자를 기반으로 구축된 분류기를 평가하여 두 그룹의 세포를 분류합니다. AUC 를 사용하여 각 유전자에 대한 분류기의 성능을 평가하며, AUC 값이 1 인 경우 해당 유전자의 발현 값만으로 두 그룹을 완벽하게 구분할 수 있음을 의미합니다. AUC 값이 0 인 경우 완벽한 분류가 이루어지지만 반대 방향입니다. AUC 값이 0.5 인 경우 해당 유전자가 두 그룹을 구분하는 데 예측력이 없음을 나타냅니다. 이 방법은 예상되는 차별적 발현 유전자의 순위를 지정된 **예측력**$(abs(AUC-0.5) * 2)$ 행렬로 반환합니다. - **t**: 두 그룹의 세포 간 차이가 있는 유전자 식별을 위한 Student 의 t-test 사용 - **negbinom**: 음이항 일반화 선형 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별 (UMI 기반 데이터 세트에만 사용) - **poisson**: 포아송 일반화 선형 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별 (UMI 기반 데이터 세트에만 사용) - **LR**: 로지스틱 회귀 모델을 사용하여 차이가 있는 유전자 식별. 각 기능을 개별적으로 사용하여 그룹 소속을 예측하는 로지스틱 회귀 모델을 구성하고 이를 가능도 비율 검정과 비교합니다. - **MAST**: scRNA-seq 데이터에 특화된 허들 모델을 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별. DE 테스트를 실행하기 위해 `MAST` 패키지를 사용합니다. - **DESeq2**: DESeq2 를 사용하여 두 그룹의 세포 간 차이가 있는 유전자 식별. DESeq2 는 음이항 분포를 사용하는 모델에 기초하며, 이 테스트는 세포 그룹 간 평균 차이 (또는 감지율의 백분율) 에 따른 유전자의 사전 필터링을 지원하지 않습니다. 그러나 유전자는 세포 그룹 모두에서 최소 감지율 (`min.pct`) 기준으로 사전 필터링될 수 있습니다. 이 방법을 사용하려면 [DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html) 를 설치하십시오. **Seurat** 객체의 모든 클러스터에 대해 각 클러스터와 모든 다른 세포 간의 DE 예측을 실행하려면 `FindAllMarkers` 함수를 사용하세요. ```r markers <- FindAllMarkers( object = SeuratObject, assay = "RNA", logfc.threshold = 0.25, test.use = "wilcox", slot = "data", min.pct = 0.1, min.diff.pct = -Inf, only.pos = FALSE, max.cells.per.ident = Inf, latent.vars = NULL, min.cells.feature = 3, min.cells.group = 3, pseudocount.use = 1, return.thresh = 0.01 ) ``` 아래에 나타낸 여러 가지 기준을 통해 특정 유전자를 포함하거나 출력을 필터링할 수 있습니다. 예를 들어, 상향 조절된 유전자만 테스트하면 분석 속도를 높일 수 있습니다: - `only.pos = TRUE` 를 사용하여 상향 조절된 유전자만 테스트 - 유전자가 발현된 최소 세포 수를 설정하는 `min.cells.feature` - 유전자가 발현된 클러스터의 최소 세포 수를 나타내는 `min.cells.group` - 클러스터에서 `min.pct` 이상 발현된 유전자만 테스트 - 두 그룹 간의 발현 비율 차이가 `min.pct.diff` 인 유전자만 테스트 - `p.value < return.thresh` 인 유전자만 반환 - `logFoldchange > logfc.threshold` 인 유전자만 반환 또한 고려해야 할 몇 가지 중요한 사항이 있습니다: >객체에 여러 분석이 포함된 경우, 올바른 분석 결과에서 차등발현 분석을 실행해야 합니다. 예를 들어, 데이터가 통합된 경우에는 `Batch` 를 고려해야 합니다. ```r cell_selection <- SeuratObject[, SeuratObject$seurat_clusters == 4] cell_selection <- SetIdent(cell_selection, value = "Batch") # 차별적 발현 계산 DGE_cell_selection <- FindAllMarkers( object = cell_selection, logfc.threshold = 0.2, test.use = "wilcox", min.pct = 0.1, min.diff.pct = -Inf, only.pos = FALSE, max.cells.per.ident = Inf, latent.vars = NULL, min.cells.feature = 3, min.cells.group = 3, pseudocount.use = 1, return.thresh = 0.01 ) ``` 또한 `FindMarkers` 함수를 사용하여 두 세트의 세포 간의 차이를 테스트할 수 있으며, 각 그룹에 대한 셀 이름을 `cells` 및 `cells` 로 지정할 수 있습니다. ## DEG 시각화하기 DE 테스트를 실행한 후, 다양한 방법으로 유전자를 시각화하고 싶을 수 있습니다. 그러나 먼저 각 클러스터에서 상위 DEG 를 추출해야 합니다. 각 클러스터별로 상위 5 개의 유전자를 선택하는 방법은 다음과 같습니다: ```r top5 <- markers %>% group_by(cluster) %>% top_n(-5, p_val_adj) ``` 이를 히트맵으로 시각화할 수 있습니다: ```r DoHeatmap( object = SeuratObject, features = as.character(unique(top5$gene)), group.by = "seurat_clusters", assay = "RNA" ) ``` 또는 각 클러스터가 평균 발현으로 색상별로 나타내고, 해당 유전자를 발현하는 세포의 비율로 크기별로 나타내는 점 도플롯을 사용할 수 있습니다: ```r # coord_flip() 함수를 위해 필요 library(ggplot2) DotPlot( object = SeuratObject, features = as.character(unique(top5$gene)), group.by = "seurat_clusters", assay = "RNA" ) + coord_flip() ``` 또한 각 유전자에 대해 바이올린 플롯을 그릴 수 있습니다: ```r VlnPlot( object = SeuratObject, features = as.character(unique(top5$gene)), ncol = 3, group.by = "seurat_clusters", assay = "RNA" ) ``` 바이올린 플롯을 배치별로 분할할 수도 있습니다. 이것은 감지된 DEG 가 단일 배치로 인해 주도되는 것이 아닌지 확인하는 데 매우 유용할 수 있습니다. ```r VlnPlot( object = SeuratObject, features = as.character(unique(top5$gene)), ncol = 3, group.by = "seurat_clusters", assay = "RNA", split.by = "Batch" ) ``` # 기능적 풍부도 분석 (Functional enrichment anaysis) `enrichR` 을 사용해 DEG 에 대하여 기능적 풍부도 분석을 할 수 있습니다. ```r # 추가 패키지 불러오기 library(enrichR) # 사용 가능한 데이터베이스를 확인 enrichR::listEnrichrDbs() enrich_results <- enrichr( genes = top5$gene[top5$cluster == "1"], # 관심있는 클러스터 이름 databases = "KEGG_2019_Human" )[[1]] # 사용할 데이터베이스 ``` `enrichR` 에 포함된 유용한 데이터베이스 목록은 다음과 같습니다: - GO_Biological_Process_2017b - KEGG_2019_Human - KEGG_2019_Mouse - WikiPathways_2019_Human - WikiPathways_2019_Mouse 결과를 시각화하기 위해 간단한 막대 그래프를 사용할 수 있습니다. 예를 들어: ```r library(rafalib) # 추가 패키지 로드 # 그래프 설정 rafalib::mypar(1, 1, mar = c(3, 20, 2, 1)) # 막대 그래프 그리기 barplot( height = -log10(enrich_results$P.value)[20:1], names.arg = enrich_results$Term[20:1], horiz = TRUE, las = 2, border = FALSE, cex.names = 0.8 # 이름 크기 설정 (예시 값) ) # 기준선 추가 abline(v = c(-log10(0.05)), lty = 2) # 유의성 기준선 abline(v = 0, lty = 1) # 0 기준선 ``` # 경로 분석 (trajectory analysis) 탐색적 데이터 분석 및 연구 데이터에는 다양한 방법으로 수행되는 데이터 처리 단계가 포함됩니다. 이러한 데이터 처리 단계는 경로 분석에 필요한 정보를 준비하는 데 도움을 줍니다. 주요 단계는 다음과 같습니다: 1. 경로 분석을 수행할 차원 축소 기법 선택 (PCA > MDS > DM > UMAP > t-SNE) 2. 세포 클러스터링 정보 3. 경로 분석을 위한 "최적" 방법 선택은 일반적으로 다른 scRNA-seq 분석과 마찬가지로 간단하지 않은 과정입니다.[^16] ## Slingshot 사용하기 슬링샷 (Slingshot) 은 차원 축소, 클러스터링 및 주요 곡선 (Principal curves) 과 같은 다른 경로 분석 방법에서 발견되는 여러 인기 있는 구성 요소를 포함하고 있어 가장 전형적인 방법 중 하나입니다. 슬링샷은 특정 차원 축소 방법에 제한되지 않으며, 어떤 방법을 사용할지 결정할 때 고려해야 할 중요한 세 가지 사항이 있습니다: 1. 데이터의 전체 복잡성을 포착하기에 충분한 차원을 사용해야 합니다. 이는 특히 PCA 및 MDS 와 같은 선형 차원 축소에서 중요하며, 경로의 위상이 이분화보다 복잡할 경우 더욱 그렇습니다. 두 차원에서 경로가 명확하게 보이지 않을 때도 여전히 다차원에서 경로를 시각화할 수 있습니다. 2. 일부 차원 축소 방법은 고밀도 영역에서 거리를 확대할 수 있습니다. t-SNE 와 UMAP 은 이러한 문제를 겪을 수 있습니다. 우리가 사용하는 데이터셋에서는 이러한 문제가 발생하지 않지만, 생물학적 모집단을 무작위로 샘플링할 경우 문제가 발생할 수 있습니다. 3. 일부 차원 축소 방법은 셀의 그룹화를 강제하고 연속성을 제거하려고 합니다. t-SNE 와 UMAP 도 이러한 문제를 일부 가지고 있으며, 이로 인해 실제로는 그렇지 않은 셀들이 함께 묶일 수 있습니다. 이러한 이유로 MDS 와 확산 맵이 종종 선호되는 선택지입니다. 그러나 샘플이 균형 잡힌 경우 t-SNE 와 UMAP 도 효과적으로 작동할 수 있습니다. ```r dimred <- SeuratObject@reductions[['pca']]@cell.embeddings # 사용될 embedding 지정 ``` 그 답은 간단합니다. 클러스터링은 문제를 단순화하며, 경로 상에서 거의 동일한 위치에 있는 세포 그룹을 찾아냅니다. 이러한 그룹은 이후 세포가 어떤 경로를 통해 이동하고 분기되는지를 파악하는 데 연결될 수 있습니다. 이상적으로 TI 를 위한 클러스터링 방법은 밀도를 기반으로 셀을 그룹화해서는 안 됩니다. 예를 들어 DBSCAN 과 같은 밀도 기반 클러스터링은 정의상 모든 세포가 고밀도 영역을 통해 연결되지만 유전자 발현에서는 차이가 발생할 수 있습니다. 여기서는 **Seurat**의 표준 클러스터링 방법인 라이든 (Leiden) 클러스터링을 사용할 것입니다. ```r clustering <- factor( SeuratObject@meta.data[, "seurat_clusters"] ) ``` 이제부터 이 클러스터링 결과를 경로 분석에 적용할 것입니다. 전체 프로세스는 단일 함수인 `slingshot` 을 사용하여 수행할 수 있으며, 이는 경로 분석을 위한 두 가지 주요 단계를 간단히 묶어 놓은 래퍼 (wrapper) 입니다. 프로세스의 첫 번째 단계는 세포 계통 (cell lineages) 을 정의하고, 그 다음 데이터를 통해 경로를 정의하는 곡선을 맞추는 것입니다. 이러한 단계들은 아래에서 설명합니다. ### Slingshot 으로 세포 계통 정의하기 ```r # 기본 설정으로 Slingshot 세포 계통 식별을 실행 library(slingshot) lineages <- getLineages( data = dimred, clusterLabels = clustering ) # Lineages를 시각화. par(mfrow = c(1,2)) plot( dimred[, 1:2], col = clustering, cex = , pch = 16 ) ``` 여기서 우리는 경로 분석의 한 가지 중요한 문제를 살펴볼 수 있습니다: 경로 분석이 어디에서 시작되는지를 알 수 없다는 점입니다. 추가 정보가 없이는 이를 파악하기가 거의 불가능합니다. 경로 분석의 시작점과 끝점을 정의하기 위해서는 사전 생물학적 정보가 필요합니다. ```r # Cluster 5에서 Slingshot을 시작하기. lineages <- getLineages( data = dimred, clusterLabels = clustering, start.clus = "5" ) lineages ``` ### 주요 경로 정의하기 시작 클러스터를 정한 뒤에는 Slingshot 을 사용하여 경로 분석을 수행할 수 있습니다. Slingshot 의 알고리즘은 초기 경로를 반복적으로 수정하여 데이터 포인트와 더 잘 일치되도록 합니다. 구체적으로 Slingshot 은 두 가지 기능을 사용합니다. - **각 "계통"마다 주요 경로를 찾습니다**: 여기서 "계통"이란 특정 시작 클러스터에서 특정 종단 클러스터로 이어지는 클러스터 집합을 의미합니다. 즉, 각 계통은 세포들이 어떤 경로를 따라 이동하는지를 나타내며, 이를 통해 세포의 발달 과정을 이해할 수 있습니다. - **동일한 클러스터 집합을 가진 계통은 주요 곡선이 겹치는 클러스터 주위에 묶이도록 제약됩니다**: 이는 서로 연결된 세포 그룹이 같은 경로를 공유하도록 하여, 분석의 일관성을 높이는 역할을 합니다. 이렇게 하면 각 계통이 서로 다른 경로를 따라 이동하더라도, 동일한 클러스터 내에서의 관계를 유지할 수 있습니다. Slingshot 의 `getCurves()` 함수는 실행하는 데 시간이 오래 걸리므로 각 계통에서 사용할 셀의 수를 줄여 곡선 적합 프로세스의 속도를 높일 수 있습니다. 이상적으로는 모든 셀을 사용하는 것이 좋지만 여기서는 `approx_points` 를 300 으로 설정하여 속도를 높였습니다. ```r # 곡선 데이터 가져오기 curves <- getCurves( lineages, approx_points = 300, thresh = 0.01, stretch = NULL, # stretch 인수에 NULL 설정 (예시) allow.breaks = FALSE, shrink = NULL # shrink 인수에 NULL 설정 (예시) ) # 곡선 데이터 확인 curves # 차원 축소 결과 플롯 plot( dimred, col = clustering, asp = 1, pch = 16 ) # 곡선 추가 lines( curves, lwd = 3, col = 'black' ) ``` ## 경로 분석를 따라 변하는 유전자를 찾기 (Differential expression along trajectories) 경로 분석을 해석하는 주요 방법 중 하나는 경로를 따라 변하는 유전자를 찾는 것입니다. 경로 분석에서 차별적 발현 (differential expression) 을 정의하는 방법은 여러 가지가 있습니다: 1. **특정 경로를 따라 발현 변화**: 예를 들어, 세포가 유사한 시간에 따라 발현이 어떻게 변화하는지를 분석합니다. 2. **가지 사이의 발현 차이**: 서로 다른 경로 (또는 가지) 간의 유전자 발현 차이를 비교합니다. 3. **분기점에서의 발현 변화**: 세포가 분기점에 도달했을 때 발현이 어떻게 변화하는지를 살펴봅니다. 4. **경로 분석의 어느 지점에서든 발현 변화**: 경로의 모든 지점에서 유전자 발현의 변화를 평가합니다. [tradeSeq](https://github.com/statOmics/tradeSeq) 는 경로 분석 결과에서 차별적 발현 유전자를 찾는 도구 중 하나입니다. 이 도구는 일반화된 가법 모형 (Generalized Additive Models, GAMs) 을 사용하여 경로를 따라 유전자 발현을 조정하고 경로의 여러 지점 사이에서 특정 계수가 통계적으로 다른지를 테스트합니다. GAM 모델을 적합시키는 과정은 시간이 많이 걸립니다. 따라서 여기에서는 매우 엄격한 필터링을 통해 많은 유전자를 제거하였습니다. 실제로는 모든 유전자 데이터를 사용하는 것이 좋습니다. ```r # 추가 패키지 로드 BiocParallel::register(BiocParallel::SerialParam()) library(tradeSeq) # 연산 속도를 높이기 위해 일부 유전자 제거 counts <- as.matrix(SeuratObject[["RNA"]]$counts) dim(counts) # 원본 카운트 행렬 차원 확인 # 필터링된 카운트 행렬 생성 filt_counts <- counts[rowSums(counts > 5) > ncol(counts) / 100, ] dim(filt_counts) # 필터링 후 카운트 행렬 차원 확인 # 감마 분포 fitting sce <- fitGAM( counts = as.matrix(filt_counts), sds = curves ) # 결과 플롯 plotGeneCount( curve = curves, counts = filt_counts, clusters = clustering, models = sce ) ``` ### 두 시간점 사이에서 변하는 유전자 원하는 경우 사용자 정의 유사 시간 (Pseudotime) 값을 설정하여 두 시간점 사이에서 변하는 유전자를 찾을 수 있습니다. 예를 들어, 전구 세포 표시자와 같은 특정 유전자를 분석할 수 있습니다. 기본적으로는 시작점과 끝점 사이의 차이를 살펴보는 방식으로 진행됩니다. ```R # Pseudotime 시작과 끝 간의 연관성 테스트 pseudotime_start_end_association <- startVsEndTest( sce, pseudotimeValues = c(0, 1) ) # feature_id 추가 pseudotime_start_end_association$feature_id <- rownames(pseudotime_start_end_association) # p-value가 0.05 미만인 feature_id 필터링 및 선택 feature_id <- pseudotime_start_end_association %>% filter(pvalue < 0.05) %>% top_n(1, waldStat) %>% pull(feature_id) # 차별 발현 플롯 그리기 plot_differential_expression(feature_id) ``` ### 경로가 갈라지는 사이에서 변하는 유전자 두 가지 사이에서 차별적으로 발현되는 유전자는 특히 흥미로운 연구 대상입니다. 이러한 유전자 중 일부는 이전의 두 시간점 사이에서 변하는 유전자 분석에서 이미 발견되었을 수 있습니다. "경로가 갈라지는 사이에서 변하는 유전자"를 정의하는 여러 가지 방법이 있으며 방법에 따라 함수가 다릅니다.[^17] - **diffEndTest**: 각 가지의 끝점에서 발현 차이를 분석합니다. - **earlyDETest**: 가지의 분기점에서 발현 차이를 평가합니다. - **patternTest**: 경로 분석의 어느 지점에서든 발현 차이를 탐색합니다. 이 마지막 함수는 두 가지 간의 유사 시간이 정렬되어 있어야 함을 요구합니다. ```r # 차별적 끝 연관성 테스트 수행 different_end_association <- diffEndTest(sce) different_end_association$feature_id <- rownames(different_end_association) # p-value가 0.05 미만인 feature_id 필터링 및 정렬 feature_id <- different_end_association %>% filter(pvalue < 0.05) %>% arrange(desc(waldStat)) %>% dplyr::slice(1) %>% pull(feature_id) # 차별 발현 플롯 그리기 plot_differential_expression(feature_id) # 가지점 연관성 테스트 수행 branch_point_association <- earlyDETest(sce) branch_point_association$feature_id <- rownames(branch_point_association) # p-value가 0.05 미만인 feature_id 필터링 및 정렬 feature_id <- branch_point_association %>% filter(pvalue < 0.05) %>% arrange(desc(waldStat)) %>% dplyr::slice(1) %>% pull(feature_id) # 차별 발현 플롯 그리기 plot_differential_expression(feature_id) ``` 더 자세한 것은 [공식 문서](https://statomics.github.io/tradeSeq/articles/tradeSeq.html) 를 확인하십시오. # 더 읽을 거리 scRNA-seq 분야는 방대하고 끊임없이 발전하고 있습니다. 이로 인해 다양한 방법, 도구, 소프트웨어가 존재하여 유사한 작업을 수행하는 데 활용될 수 있습니다. 아래는 이 분야에 대한 더 깊은 이해를 돕기 위해 탐색할 수 있는 몇 가지 추가 자료입니다. ## 온라인 강의 - [Exercises for the NBIS scRNA-seq course](https://nbisweden.github.io/workshop-scRNAseq/exercises.html) - [The GitHub repository for the NBIS scRNA-seq course](https://github.com/NBISweden/workshop-scRNAseq) - [Single cell RNA-seq course at from Hemberg lab](https://scrnaseq-course.cog.sanger.ac.uk/website/index.html) - [Single cell RNA-seq course in Python](https://chanzuckerberg.github.io/scRNA-python-workshop/intro/about) - [Single cell RNA-seq course at Broad](https://broadinstitute.github.io/2019_scWorkshop/) ## 프레젠테이션 및 강의 - [Single-cell methodologies](https://nbisweden.github.io/workshop-scRNAseq/lectures/single_cell_methodologies_Karolina_Wallenborg_2020.pdf) - [Quaility controls](https://nbisweden.github.io/workshop-scRNAseq/lectures/scRNAseq_QC_Asa_Bjorklund_2020.pdf) - [Normalisation](https://nbisweden.github.io/workshop-scRNAseq/lectures/scRNAseq_normalization_Asa_Bjorklund_2020.pdf) - [Dimensionality reduction](https://nbisweden.github.io/workshop-scRNAseq/lectures/dimensionality_reduction_paulo_czarnewski.pdf) - [Batch correction and data integration](https://nbisweden.github.io/workshop-scRNAseq/lectures/data_integration_paulo_czarnewski_2020.pdf) - [Clustering techniques](https://nbisweden.github.io/workshop-scRNAseq/lectures/scRNAseq_clustering_Asa_Bjorklund_2020.pdf) - [Differential expression analysis](https://nbisweden.github.io/workshop-scRNAseq/lectures/differenatial_expression_olga_dethlefsen_2020.pdf) - [Trajectory inference analyses](https://nbisweden.github.io/workshop-scRNAseq/lectures/trajectory_inference_analysis_paulo_czarnewski.pdf) - [All of the above in a YouTube playlist](https://www.youtube.com/playlist?list=PLjiXAZO27elC_xnk7gVNM85I2IQl5BEJN) ## 도구 및 파이프라인 - [Repository listing many scRNA-seq tools](https://github.com/seandavi/awesome-single-cell) - [Bitbucket repository for QC report scripts](https://bitbucket.org/asbj/qc-summary_scrnaseq) - [Bitbucket repository for an NBIS scRNA-seq pipeline](https://bitbucket.org/scilifelab-lts/lts-workflows-sm-scrnaseq) - [A catalogue of many scRNA-seq tools](https://www.scrna-tools.org/) ## 웹사이트 - [SingleCellExperiment objects for many datasets](https://hemberg-lab.github.io/scRNA.seq.datasets/) - [Conquer datasets - many different datasets based on a salmon pipeline](http://imlspenticton.uzh.ch:3838/conquer/) - [The Human Cell Atlas project](https://www.humancellatlas.org/) - [The EBI Single-cell expression atlas](https://www.ebi.ac.uk/gxa/sc) # Reference [^1]: https://www.pnas.org/content/116/20/9775) [^2]: https://distill.pub/2016/misread-tsne, http://jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf [^3]: https://umap-learn.readthedocs.io/en/latest [^4]: https://arxiv.org/abs/1802.03426 [^5]: https://projecteuclid.org/euclid.ajm/1253539864 [^6]: https://en.wikipedia.org/wiki/Independent_component_analysis [^7]: https://www.nature.com/articles/nbt.4091 [^8]: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9 [^9]: https://www.nature.com/articles/s41587-019-0113-3 [^10]: https://www.nature.com/articles/s41592-019-0619-00 [^11]: https://f1000research.com/articles/7-1141 [^12]: https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/pdf [^13]: https://www.nature.com/articles/s41598-019-41695-z.pdf [^14]: https://www.nature.com/articles/nmeth.4612 [^15]: McDavid et al., Bioinformatics, 2013 [^16]: https://www.nature.com/articles/s41587-019-0071-9 [^17]: Cannoodt, Robrecht, Wouter Saelens, and Yvan Saeys. 2016."Computational Methods for Trajectory Inference from Single-Cell Transcriptomics."European Journal of Immunology 46 (11): 2496–2506.