단일세포 분석에서 우리는 원하는 클러스터를 식별한 후, 마커 유전자를 확인하여 특정 클러스터의 정체성을 검증하고 미지의 클러스터의 정체성을 추측하는 데 도움을 얻을 수 있습니다.
학습 목표
- 개별 클러스터의 마커 유전자를 식별하는 방법 학습
- 클러스터링과 마커 식별 간의 반복적인 과정 수행
목표
- 각 클러스터의 유전자 마커 결정
- 마커를 사용하여 각 클러스터의 세포 유형 식별
- 세포 유형 마커를 기반으로 재클러스터링 필요성 판단 (클러스터 병합 또는 분할)
도전 과제
- 결과 과대 해석 가능성
- 다양한 마커 유형 결합을 통한 식별 필요
권장 사항
- 결과를 검증이 필요한 가설로 간주. 과도한 p값은 결과 과대 해석으로 이어질 수 있음. 가장 신뢰도 높은 상위 마커를 우선적으로 분석. 각 조건에서 각 클러스터의 보수적 마커 식별.
- 특정 클러스터 간 차등 발현 마커 식별.
이전 클러스터링 분석 결과는 다음과 같습니다:
클러스터링 분석에서 다음과 같은 문제에 직면했습니다:
- 클러스터 7과 20의 세포 유형은 무엇인가?
- 동일한 세포 유형에 해당하는 클러스터는 생물학적으로 의미 있는 차이를 가지는가? 이러한 세포 유형에는 아류가 존재하는가?
- 이러한 클러스터의 다른 마커 유전자를 식별함으로써 세포 유형에 대한 식별 결과를 검증할 수 있는가?
우리는 이러한 질문에 답하기 위해 Seurat를 사용하여 여러 유형의 마커를 탐색할 수 있습니다. 각 방법은 고유한 장점과 단점을 가집니다:
-
각 클러스터의 모든 마커 식별: 이 분석은 각 클러스터를 다른 모든 클러스터와 비교하고 차등 발현 유전자를 출력합니다. 미지의 클러스터를 식별하고 가설적 세포 유형에 대한 신뢰도를 높이는 데 사용됩니다.
-
각 클러스터의 보수적 마커 식별: 이 분석은 먼저 각 조건에서 차등 발현하는 유전자를 찾은 다음 모든 조건에서 클러스터에서 보존되는 유전자를 보고합니다. 이러한 유전자는 클러스터의 정체성을 결정하는 데 도움이 됩니다. 여러 조건에서 보존된 세포 유형 마커를 식별하는 데 적합합니다.
-
특정 클러스터 간의 마커 식별: 이 분석은 특정 클러스터 간의 차등 발현 유전자를 탐색합니다. 동일한 세포 유형을 나타내는 것처럼 보이는 클러스터(즉, 유사한 마커를 가짐) 간의 유전자 발현 차이를 결정하는 데 사용됩니다.
5. 각 클러스터의 모든 마커 식별
단일 샘플을 평가할 때는 일반적으로 이 유형의 분석을 사용하는 것이 좋습니다. FindAllMarkers() 함수를 사용하여 각 클러스터를 다른 모든 클러스터와 비교하여 잠재적인 마커 유전자를 식별합니다. 각 클러스터의 세포는 반복으로 간주되며, 기본적으로 일부 통계 테스트를 통해 차등 발현 분석을 수행합니다.
참고: 기본값은 Wilcoxon 순위 합 검정이지만 다른 옵션도 사용 가능합니다.
FindAllMarkers() 함수에는 유전자가 마커인지 결정하는 임계값을 제공하는 세 가지 중요한 매개변수가 있습니다:
-
logfc.threshold: 클러스터 내 유전자 평균 발현량 대비 다른 모든 클러스터에서의 평균 발현량의 최소 log2 배 변화량. 기본값은 0.25입니다.
-
단점:
-
관심 클러스터 내에서 소수의 세포에서 발현되지만 다른 클러스터에서는 발현되지 않는 세포 마커를 놓칠 수 있음
-
다양한 세포 유형의 대사 출력이 약간 다르기 때문에 세포 유형 정체성을 구별하는 데 유용하지 않은 대량의 대사/리보솜 유전자를 반환할 수 있음
-
min.diff.pct: 클러스터에서 발현되는 유전자의 세포 백분율과 다른 모든 클러스터에서 발현되는 유전자의 세포 백분율 간의 최소 백분율 차이.
-
단점: 모든 세포에서 발현되지만 이 특정 세포 유형에서 고도로 상향 조절된 세포 마커를 놓칠 수 있음
-
min.pct: 두 집단 중 하나에서만 검출된 유전자만 테스트. 거의 발현되지 않는 유전자를 테스트하지 않음으로써 속도를 높이려는 목적. 기본값은 0.1입니다.
-
단점: 매우 높은 값으로 설정하면 많은 거짓음성이 발생할 수 있음, 모든 유전자가 모든 세포에서 검출되지는 않기 때문(발현되더라도)
원하는 엄격성에 따라 이러한 매개변수를 임의로 조합하여 사용할 수 있습니다. 또한, 기본적으로 이 함수는 양성 및 음성 발현 변화를 모두 보여줍니다. 일반적으로 only.pos 매개변수를 추가하여 긍정적인 변화만 유지하도록 선택합니다. 각 클러스터에 대한 마커를 찾는 코드는 다음과 같습니다.
# 각 클러스터와 모든 나머지 세포를 비교하여 마커 찾기, 양성 결과만 보고
markers <- FindAllMarkers(object = seurat_integrated,
only.pos = TRUE,
logfc.threshold = 0.25)
6. 모든 조건에서 보수적 마커 식별
데이터셋에 서로 다른 조건을 나타내는 샘플이 있으므로 가장 좋은 선택은 보수적인 마커를 찾는 것입니다. 이 방법은 내부적으로 샘플 그룹/조건별로 세포를 분리한 다음, 단일 지정 클러스터에 대해 차등 유전자 발현 테스트를 수행합니다. 각 조건의 유전자 수준 p값을 계산한 후, MetaDE R 패키지의 메타 분석 방법을 사용하여 그룹 간에 결합합니다.
마커 식별을 시작하기 전에 기본 분석을 명시적으로 설정하고, 표준화 데이터를 사용하도록 지정합니다.
DefaultAssay(seurat_integrated) <- "RNA"
FindConservedMarkers() 함수는 다음과 같은 구조를 가집니다:
FindConservedMarkers(seurat_integrated,
ident.1 = cluster,
grouping.var = "sample",
only.pos = TRUE,
min.diff.pct = 0.25,
min.pct = 0.25,
logfc.threshold = 0.25)
이전에 FindAllMarkers() 함수에 대해 설명한 일부 매개변수를 인식할 수 있을 것입니다. 이는 내부적으로 해당 함수를 사용하여 각 그룹에서 먼저 마커를 찾기 때문입니다. 여기서는 FindConservedMarkers() 사용 시 제공되는 추가 매개변수를 나열합니다:
- ident.1: 이 함수는 한 번에 하나의 클러스터만 평가합니다. 여기서 관심 클러스터를 지정합니다.
- grouping.var: 메타데이터의 변수(열 제목), 세포를 그룹으로 지정
분석을 위해 상대적으로 느슨한 기준을 사용하고, logfc.threshold를 0.25보다 큰 값만 사용합니다. 또한 각 클러스터의 긍정적 마커만 반환하도록 지정합니다.
하나의 클러스터에서 작동 방식을 살펴보겠습니다:
cluster0_conserved_markers <- FindConservedMarkers(seurat_integrated,
ident.1 = 0,
grouping.var = "sample",
only.pos = TRUE,
logfc.threshold = 0.25)
FindConservedMarkers() 함수의 출력은 지정한 클러스터별로 정렬된 추정 마커 목록과 관련 통계 데이터를 포함하는 행렬입니다. 각 그룹(예: Ctrl 및 Stim)에 대해 동일한 통계 데이터 세트가 계산되며, 마지막 두 열은 두 그룹의 결합 p값에 해당합니다. 아래에서 이러한 열 중 일부를 설명합니다:
- gene: 유전자 기호
- condition_p_val: 조건에 대한 다중 테스트 보정되지 않은 p값
- condition_avg_logFC: 조건의 평균 로그 배 변화량. 양수 값은 해당 유전자가 클러스터에서 더 높게 발현됨을 나타냅니다.
- condition_pct.1: 클러스터에서 유전자가 검출된 세포 백분율
- condition_pct.2: 다른 클러스터에서 평균적으로 검출된 해당 유전자의 백분율
- condition_p_val_adj: 조건의 보정 p값, 데이터셋의 모든 유전자를 사용한 bonferroni 보정 기반
- max_pval: 각 그룹/조건에 대해 계산된 p값의 최대 p값
- minimump_p_val: 결합 p값
출력을 확인할 때, pct.1과 pct.2 간의 발현 차이가 크고 배 변화량이 큰 마커를 찾는 것이 좋습니다. 예를 들어, pct.1 = 0.90이고 pct.2 = 0.80이면 올바른 마커가 아닐 수 있습니다. 하지만 pct.2가 0.1이면 더 큰 차이가 더 설득력이 있을 것입니다. 또한, 관심 클러스터에 있는 대부분의 발현 마커 세포가 있을 경우(예: pct.1이 낮게, 0.3처럼) 올바른 마커가 아닐 수 있습니다. 위에서 언급했듯이, 이 두 가지는 함수 실행 시 포함할 수 있는 매개변수이기도 합니다.
6.1. 유전자 주석 추가
유전자 주석 정보가 포함된 열을 추가하면 도움이 됩니다. 이를 위해, 데이터 폴더에 있는 주석 파일을 로드하는 코드를 사용합니다:
annotations <- read.csv("data/annotation.csv")
먼저 유전자 식별자가 포함된 행 이름을 별도의 열로 변환합니다. 그런 다음 이 주석 파일을 FindConservedMarkers()의 결과와 병합합니다:
# 마커에 유전자 설명 추가
cluster0_ann_markers <- cluster0_conserved_markers %>%
rownames_to_column(var="gene") %>%
left_join(y = unique(annotations[, c("gene_name", "description")]),
by = c("gene" = "gene_name"))
View(cluster0_ann_markers)
6.2. 여러 샘플에서 실행
FindConservedMarkers() 함수는 한 번에 하나의 클러스터만 허용하므로, 클러스터 수만큼 이 함수를 여러 번 실행할 수 있습니다. 하지만 이는 효율적이지 않습니다. 대신, 먼저 모든 매개변수를 포함하는 보수적 마커를 찾는 함수를 만듭니다. 또한 출력을 수정하기 위한 몇 줄의 코드를 추가합니다. 단계는 다음과 같습니다:
- FindConservedMarkers() 함수 실행
- rownames_to_column() 함수를 사용하여 행 이름을 열로 전송
- 주석 병합
- cbind() 함수를 사용하여 클러스터 ID 열 생성
# 주어진 클러스터의 보수적 마커를 가져오는 함수 생성
get_conserved <- function(cluster){
FindConservedMarkers(seurat_integrated,
ident.1 = cluster,
grouping.var = "sample",
only.pos = TRUE) %>%
rownames_to_column(var = "gene") %>%
left_join(y = unique(annotations[, c("gene_name", "description")]),
by = c("gene" = "gene_name")) %>%
cbind(cluster_id = cluster, .)
}
이제 이 함수를 만들었으므로, 적절한 매핑 함수의 매개변수로 사용할 수 있습니다. map 시리즈 함수의 출력이 데이터 프레임이고 각 클러스터 출력이 행으로 결합되도록 하려면 map_dfr() 함수를 사용합니다.
# 예시
map_dfr(inputs_to_function, name_of_function)
이제 이 함수를 사용하여 미지의 세포 유형 클러스터의 보수적 마커를 찾아보겠습니다: 클러스터 7과 클러스터 20.
conserved_markers <- map_dfr(c(7,20), get_conserved)
6.3. 마커 유전자 평가
이러한 유전자 목록을 사용하여 이러한 클러스터가 식별하는 세포 유형을 확인하고 싶습니다. 각 클러스터의 상위 유전자를 살펴보겠습니다. 두 그룹의 평균 배 변화량을 통해 상위 10개 마커를 빠르게 확인할 수 있습니다:
# 각 클러스터에서 상위 10개 마커 추출
top10 <- conserved_markers %>%
mutate(avg_fc = (ctrl_avg_log2FC + stim_avg_log2FC) /2) %>%
group_by(cluster_id) %>%
top_n(n = 10,
wt = avg_fc)
# 각 클러스터의 상위 10개 마커 시각화
View(top10)
클러스터 7에서 많은 열 shock 및 DNA 손상 유전자가 나타납니다. 이러한 마커를 기반으로, 이것들은 스트레스 또는 죽어가는 세포일 가능성이 높습니다. 그러나 세포의 품질 지표(즉, mitoRatio 및 nUMI에 표시된 클러스터)를 자세히 살펴보면 해당 주장을 뒷받침하는 데이터를 실제로 볼 수 없습니다. 마커 유전자 목록을 자세히 살펴보면, T 세포 관련 유전자와 활성화 마커도 발견할 수 있습니다. 이들은 활성화된(세포독성) T 세포일 수 있습니다. 열 shock 단백질이 만성 염증에서 반응성 T 세포에 의해 항염증성 사이토카인이 유도된다는 것을 뒷받침하는 많은 연구가 있습니다. 이 클러스터는 면역 세포에 대한 더 깊은 이해가 필요하여 결과를 정리하고 최종 결론을 도출할 수 있습니다.
클러스터 20의 경우, 풍부한 유전자는 공통 관계를 가지지 않는 것처럼 보입니다. pct.1과 pct.2 간의 차이가 큰 유전자를 확인하여 좋은 마커 유전자를 얻을 수 있습니다. 예를 들어, 대부분의 세포에서 발현되지만 다른 클러스터에서 발현되는 세포가 거의 없는 TPSB2 유전자에 관심이 있을 수 있습니다.
따라서 클러스터 20은 비만 세포를 나타낼 수 있습니다. 비만 세포는 면역 시스템의 중요한 세포이며 혈액 생성 계통에 속합니다. 연구는 비만 세포 특징이 세린 프로테아제(예: TPSAB1 및 TPSB2)가 풍부하다는 것을 확인했습니다. 이러한 프로테아제는 모두 우리의 보수적 마커 목록에 나타납니다. 또 다른 세린 프로테아제가 아니지만 알려진 비만 세포 특이적 유전자이며 목록에 나타나는 유전자는 FCER1A(면역글로불린 E 수용체의 하나의 아단위를 코딩)입니다. 또한, GATA1 및 GATA2가 목록에 나타나는데, 이들은 비만 세포 마커 유전자는 아니지만 비만 세포에서 대량으로 발현되며 다양한 비만 세포 특이적 유전자를 조절하는 것으로 알려진 전사 인자입니다.
6.4. 마커 유전자 시각화
클러스터 20의 세포 유형 정체성을 더 잘 이해하기 위해, FeaturePlot() 함수를 사용하여 클러스터별로 다양한 식별된 마커의 발현을 탐색할 수 있습니다.
# 클러스터 20의 마커 유전자 발현 플롯
FeaturePlot(object = seurat_integrated,
features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
특정 마커의 발현 범위를 탐색하기 위해 바이올린 플롯도 사용할 수 있습니다:
# 바이올린 플롯 - 클러스터 20
VlnPlot(object = seurat_integrated,
features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"))
이러한 결과와 차트는 우리가 이러한 클러스터의 정체성을 결정하거나, 이전에 탐색한 예상 세포 유형의 표준 마커를 기반으로 가설적 정체성을 검증하는 데 도움이 될 수 있습니다.
7. 각 클러스터의 마커 식별
분석의 마지막 그룹 문제는 동일한 세포 유형에 해당하는 클러스터가 생물학적으로 의미 있는 차이를 가지는지 여부에 관한 것입니다. 때로는 반환된 마커 목록이 특정 클러스터를 충분히 분리하지 못할 수 있습니다. 예를 들어, 이전에 0, 2, 4, 10 및 18번 클러스터를 CD4+ T 세포로 식별했지만, 이 세포 클러스터 간에는 생물학적으로 관련된 차이가 있는가? FindMarkers() 함수를 사용하여 두 특정 클러스터 간의 차등 발현 유전자를 결정할 수 있습니다.
모든 비교 조합을 시도할 수 있지만, 클러스터 2와 다른 모든 CD4+ T 세포 클러스터부터 시작하겠습니다:
# CD4+ T 세포의 분화 마커 결정
cd4_tcells <- FindMarkers(seurat_integrated,
ident.1 = 2,
ident.2 = c(0,4,10,18))
# DE 테이블에 유전자 기호 추가
cd4_tcells <- cd4_tcells %>%
rownames_to_column(var = "gene") %>%
left_join(y = unique(annotations[, c("gene_name", "description")]),
by = c("gene" = "gene_name"))
# 열 재정렬 및 padj 기준 정렬
cd4_tcells <- cd4_tcells[, c(1, 3:5,2,6:7)]
cd4_tcells <- cd4_tcells %>%
dplyr::arrange(p_val_adj)
# 데이터 확인
View(cd4_tcells)
이러한 상위 유전자에서, CREM 유전자가 활성화 마커로 두드러집니다. 우리는 또 다른 활성화 마커가 CD69이고, 미성숙 또는 기억 세포의 마커에는 SELL 및 CCR7 유전자가 포함된다는 것을 알고 있습니다. 흥미롭게도, SELL 유전자도 상위에 있습니다. 이러한 새로운 세포 상태 마커를 사용하여 활성화 상태를 시각적으로 탐색해 보겠습니다:
# 활성화 및 기억 T 세포의 유전자 마커 플롯
FeaturePlot(seurat_integrated,
reduction = "umap",
features = c("CREM", "CD69", "CCR7", "SELL"),
label = TRUE,
order = TRUE,
min.cutoff = 'q10',
repel = TRUE
)
초기 상태와 활성화 상태의 마커가 모두 마커 목록에 나타나므로 발현을 시각화하는 것이 도움이 됩니다. 이러한 플롯을 기반으로, 클러스터 0과 2는 신선한 T 세포로 신뢰할 수 있습니다. 그러나 활성화된 T 세포에 대해서는 말하기 어렵습니다. 우리는 클러스터 4와 18이 활성화된 T 세포라고 말할 수 있지만, CD69의 발현은 CREM만큼 명확하지 않습니다. 우리는 미성숙 세포를 표시하고 나머지 클러스터를 CD4+ T 세포로 표시하겠습니다.
이제 이 모든 정보를 바탕으로, 다양한 클러스터의 세포 유형을 추측하고 세포 유형 레이블이 지정된 세포를 플롯할 수 있습니다.
그런 다음 이러한 세포 유형에 클러스터의 정체성을 다시 할당할 수 있습니다:
# 모든 정체성 이름 변경
seurat_integrated <- RenameIdents(object = seurat_integrated,
"0" = "Naive or memory CD4+ T cells",
"1" = "CD14+ monocytes",
"2" = "Naive or memory CD4+ T cells",
"3" = "CD14+ monocytes",
"4" = "CD4+ T cells",
"5" = "CD8+ T cells",
"6" = "B cells",
"7" = "Stressed cells / Activated T cells",
"8" = "NK cells",
"9" = "FCGR3A+ monocytes",
"10" = "CD4+ T cells",
"11" = "B cells",
"12" = "NK cells",
"13" = "CD8+ T cells",
"14" = "CD14+ monocytes",
"15" = "Conventional dendritic cells",
"16" = "Megakaryocytes",
"17" = "B cells",
"18" = "CD4+ T cells",
"19" = "Plasmacytoid dendritic cells",
"20" = "Mast cells")
# UMAP 시각화
DimPlot(object = seurat_integrated,
reduction = "umap",
label = TRUE,
label.size = 3,
repel = TRUE)
잠재적인 스트레스 세포를 제거하려면 subset() 함수를 사용할 수 있습니다:
# 죽어가는 세포 제거
seurat_subset_labeled <- subset(seurat_integrated,
idents = "Stressed cells / Activated T cells", invert = TRUE)
# 클러스터 재시각화
DimPlot(object = seurat_subset_labeled,
reduction = "umap",
label = TRUE,
label.size = 3,
repel = TRUE)
이제 클러스터와 각 클러스터의 마커를 정의했으므로, 다음과 같은 다양한 선택지가 있습니다:
-
식별된 세포 유형의 마커를 실험적으로 검증.
-
세포 아류 발견을 위한 세포 유형 하위 집합 탐색
-
조건 ctrl과 stim 간 차등 발현 분석 수행
-
세포 유형 또는 세포 상태 간의 상황을 결정하려는 경우, 궤적 분석 또는 계보 추적 수행:
-
분화 과정
-
시간에 따른 발현 변화
-
발현 과정에서의 세포 상태 변화