學習目標：
- 了解如何確定單個亞群的標記
- 理解聚類和標記識別的迭代過程
單細胞RNA-seq標記鑒定
現(xiàn)在我們已經確定了我們想要的集群，我們可以繼續(xù)進行標記識別，這將使我們能夠驗證某些集群的身份，并幫助推測任何未知集群的身份。
目標：
- 為了確定每個集群的基因標記
- 使用標記來識別每個集群的細胞類型
- 為了確定是否需要基于細胞類型標記重新聚類，可能需要合并或拆分聚類
挑戰(zhàn)：
- 對結果過度解讀
- 結合不同類型的標記進行識別
建議：
- 把結果看作是需要驗證的假設?？浯蟮膒值可能導致結果的過度解釋(基本上每個細胞都被用作一個復制)。得分高的值可信度高。識別每個聚類的所有保存在條件之間的標記
- 識別特定簇之間表達有差異的標記
通過聚類分析產生了以下聚類：
請記住，我們在聚類分析中遇到了以下問題：
- 集群 7 和 20 的細胞類型標識是什么？
- 對應于相同細胞類型的簇是否具有生物學意義的差異？是否存在這些細胞類型的亞群中？
- 通過識別這些簇的其他標記基因，我們能否對這些細胞類型身份獲得更高的可信度？
我們可以使用 Seurat 探索幾種不同類型的識別標記，以獲得這些問題的答案。每個都有自己的優(yōu)點和缺點：
-
對每個集群的所有標記進行識別：該分析將每個集群與其他所有集群進行比較，并輸出差異表達/存在的基因。
- 用于識別未知簇并提高對假設細胞類型的可信度。
-
識別每個簇的保守標記：該分析首先尋找在每個條件下差異表達/存在的基因，然后找出所有條件下在簇中保守的基因。這些基因可以幫助確定集群的身份。
- 用于多個條件以識別不同條件下保守的細胞類型標記。
-
特定簇之間的標記識別：該分析探討了特定簇之間差異表達的基因。
- 用于確定從上述分析中似乎代表相同細胞類型（即具有相似標記）的簇之間的基因表達差異。
識別每個簇的所有標記
這種類型的分析通常推薦用于評估單個樣品組/條件。使用 FindAllMarkers() 函數，我們將每個集群與所有其他集群進行比較，以識別潛在的標記基因。每個集群中的細胞被視為重復，本質上通過一些統(tǒng)計檢驗進行差異表達分析。
注意：默認是 Wilcoxon Rank Sum檢驗，但也可以選擇其他檢驗方法。
該FindAllMarkers()函數有三個重要參數，它們?yōu)榇_定基因是否為標記提供閾值：
-
logfc.threshold：相對于所有其他聚類組合中基因平均表達的最小log2倍變化。默認值為 0.25。-
缺點：
- 如果平均logfc不滿足閾值，可能會錯過那些在感興趣的集群內的一小部分細胞中表達的細胞標記，而在其他集群中不會
- 由于不同細胞類型的代謝輸出略有差異，可能會返回大量代謝/核糖體基因，這對于區(qū)分細胞類型身份沒有什么用
-
缺點：
-
min.diff.pct：表達該基因的細胞百分數與所有其他基因組合中表達基因的細胞百分數之間的最小差異百分數。- 缺點：可能會錯過那些在所有細胞中表達的細胞標記物，但在這個特定的細胞類型中高度上調
-
min.pct：只檢測在兩個種群中最小部分細胞中檢測到的基因。這意味著通過不測試那些很少表達的基因來加快功能。默認是0.1。- 缺點：由于不是所有的基因都能在所有的細胞中被檢測到(即使是表達出來的)，如果設置一個非常高的值可能會導致許多假陰性
您可以使用這些參數的任意組合，這取決于您想要的嚴格/寬松程度。此外，默認情況下，這個函數將返回陽性和陰性表達變化的基因，添加參數only.pos來選擇只保留陽性表達的基因。為每個集群查找標記的代碼如下所示。
## DO NOT RUN THIS CODE ##
# Find markers for every cluster compared to all remaining cells, report only the positive ones
markers <- FindAllMarkers(object = seurat_integrated,
only.pos = TRUE,
logfc.threshold = 0.25)
注意：此命令可能需要很長時間才能運行，因為它正在針對其他所有細胞處理每個單獨的集群。
在所有條件下鑒定保守標記
由于我們的數據集中有不同條件下的樣本，因此我們最好的選擇是找到保守標記。此函數在內部按樣本組/條件分離出細胞，然后針對單個指定簇與所有其他簇（或第二個簇，如果指定）執(zhí)行差異基因表達測試。計算每個條件下基因的 p 值，然后使用R 包中的 MetaDE 薈萃分析方法跨組組合。
在我們開始我們的標記識別之前，我們將明確設置我們的默認檢測，我們想要使用標準化數據，而不是整合數據。
DefaultAssay(seurat_integrated) <- "RNA"
注意：默認檢測應該已經是
RNA，因為我們在之前的聚類質量控制課程中設置了它。但我們鼓勵您運行上面的這行代碼，以確保萬一運行環(huán)境在您的分析上游某處發(fā)生了更改。請注意，RNA原始計數和歸一化計數存儲在檢測的counts和data插槽。默認情況下，查找標記的函數將使用標準化數據。
函數FindConservedMarkers(), 具有以下結構：
FindConservedMarkers() 句法：
## DO NOT RUN ##
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：元數據中的變量（列標題），用于指定將細胞分成組
對于我們的分析，只使用大于 0.25 的對數倍數變化閾值。返回每個集群的陽性標記。
讓我們在一個集群上測試它，看看它是如何工作的：
cluster0_conserved_markers <- FindConservedMarkers(seurat_integrated,
ident.1 = 0,
grouping.var = "sample",
only.pos = TRUE,
logfc.threshold = 0.25)
FindConservedMarkers()函數的輸出是一個矩陣，其中包含按我們指定的聚類的基因 ID 列出的推定標記的排名列表，以及相關的統(tǒng)計數據。請注意，為每個組（在我們的例子中，Ctrl 和 Stim）計算相同的統(tǒng)計數據集，最后兩列對應于兩個組的組合 p 值。我們在下面描述了其中的一些列：
- 基因：基因符號
- 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 值
注意：由于每個細胞都被視為一個重復，這將導致每個組內的 p 值膨脹！一個基因可能具有較低 的p 值 < 1e-50，但這并不能轉化為高度可靠的標記基因。
在查看輸出時，我們建議在pct.1和pct.2之間尋找表達差異很大的標記，以及更大的折疊變化。例如，如果pct.1= 0.90 和pct.2= 0.80，它可能不會像標記一樣令人興奮。但是，如果pct.2改為 = 0.1，則更大的差異會更令人信服。此外，感興趣的是表達標記的大多數細胞是否在我感興趣的簇中。如果pct.1低，例如 0.3，則可能不那么有趣。如上所述，這兩個也是在運行函數時可能包含的參數。
添加基因注釋
添加基因注釋列信息會很有幫助。為此，我們將data使用下面提供的代碼加載位于您文件夾中的注釋文件：
annotations <- read.csv("data/annotation.csv")
注意：如果您有興趣了解我們如何獲得此注釋文件，請查看鏈接的材料。
首先，我們將帶有基因標識符的行名變成它自己的列。然后我們將這個注釋文件與我們的結果合并FindConservedMarkers()：
# Combine markers with gene descriptions
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)
練習
在上一課中，我們通過檢查已知細胞標記物 FCGR3A 和 MS4A7 的表達，將簇 9 確定為 FCGR3A+ 單核細胞。使用FindConservedMarkers()函數找到簇 9 的保守標記。你觀察到什么？您是否將 FCGR3A 和 MS4A7 視為簇 9 中的高表達基因？
在多樣本上運行
該FindConservedMarkers()函數一次 接受一個集群，我們可以運行該函數次數與集群數量一樣多。然而，這不是很有效。相反，我們將首先創(chuàng)建一個函數來查找包含我們想要包含的所有參數的保守標記。我們還將添加幾行代碼來修改輸出。我們的函數將：
- 運行
FindConservedMarkers()函數 - 使用
rownames_to_column()函數將行名稱傳輸到列 - 合并注釋
- 使用
cbind()函數創(chuàng)建集群 ID 列
# Create function to get conserved markers for any given cluster
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, .)
}
現(xiàn)在我們已經創(chuàng)建了這個函數，我們可以將它用作適當map函數的參數。我們希望map函數的輸出是一個數據幀，每個集群輸出通過行綁定在一起，我們將使用該map_dfr()函數。
map 家族語法：
## DO NOT RUN ##
map_dfr(inputs_to_function, name_of_function)
現(xiàn)在，讓我們嘗試使用此函數來查找未識別細胞類型的簇的保守標記：簇 7 和簇 20。
# Iterate function across desired clusters
conserved_markers <- map_dfr(c(7,20), get_conserved)
為所有集群尋找標記
對于您的數據，您可能希望在所有集群上運行此函數，在這種情況下，您可以輸入
0:20而不是c(7,20); 但是，運行需要相當長的時間。此外，當您在所有集群上運行此函數時，在某些情況下，您的集群可能沒有足夠的細胞用于特定組- 并且您的函數將失敗。對所有的集群，您需要使用FindAllMarkers().
評估標記基因
我們想使用這些基因列表來查看我們可以識別這些簇與哪些細胞類型相關。讓我們看看每個集群的頂部的基因，看看這是否給我們任何提示。我們可以通過兩個組的平均倍數變化查看前 10 個標記，對于每個集群進行快速閱讀：
# Extract top 10 markers per cluster
top10 <- conserved_markers %>%
mutate(avg_fc = (ctrl_avg_log2FC + stim_avg_log2FC) /2) %>%
group_by(cluster_id) %>%
top_n(n = 10,
wt = avg_fc)
# Visualize top 10 markers per cluster
View(top10)
我們看到許多熱休克和 DNA 損傷基因出現(xiàn)在第7 簇中。根據這些標記，這些很可能是受壓或垂死的細胞。然而，如果我們更詳細地探索這些細胞的質量指標（即覆蓋在集群上的 mitoRatio 和 nUMI），我們并沒有真正看到支持該論點的數據。如果我們仔細查看標記基因列表，我們還會看到一些 T 細胞相關基因和激活標記. 這些有可能被激活（細胞毒性）T 細胞。有廣泛的研究支持熱休克蛋白與反應性 T 細胞在慢性炎癥中誘導抗炎細胞因子的關聯(lián)。在這個集群中，我們需要對免疫細胞有更深入的了解，才能真正梳理結果并得出最終結論。
對于第20 組，富集的基因似乎沒有一個讓我們感興趣。我們經常查看pct.1與pct.2良好標記基因相比差異較大的基因。例如，我們可能對基因 TPSB2 感興趣，該基因顯示簇中表達該基因的細胞比例很大，但其他簇中表達該基因的細胞很少。如果我們“谷歌”TPSB2，我們會找到GeneCards 網站。
“β類胰蛋白酶似乎是肥大細胞中表達的主要同工酶，而在嗜堿性粒細胞中，α-類胰蛋白酶占主導地位。類胰蛋白酶被認為是哮喘和其他過敏性和炎癥性疾病發(fā)病機制的介質?！?/p>
因此，簇20可能代表肥大細胞。肥大細胞是免疫系統(tǒng)的重要細胞，屬于造血譜系。研究已經確定肥大細胞特征顯著富含絲氨酸蛋白酶，例如TPSAB1 和 TPSB2，這兩種蛋白酶都出現(xiàn)在我們的保守標記列表中。另一個不是絲氨酸蛋白酶但是已知肥大細胞特異性基因并出現(xiàn)在我們列表中的基因是 FCER1A（編碼 IgE 受體的一個亞基）。此外，我們看到 GATA1 和 GATA2 出現(xiàn)在我們的列表中，它們不是肥大細胞標記基因，但在肥大細胞中大量表達，是已知的轉錄因子調節(jié)各種肥大細胞特異性基因。
可視化標記基因
為了更好地了解集群 20的細胞類型身份，我們可以使用該FeaturePlot()函數按集群探索不同已識別標記的表達。
# Plot interesting marker gene expression for cluster 20
FeaturePlot(object = seurat_integrated,
features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
我們還可以使用小提琴圖探索特定標記的表達范圍：
小提琴圖類似于箱線圖，不同之處在于它們還顯示了數據在不同值下的概率密度，通常由核密度評估器平滑。小提琴圖比普通箱線圖提供更多信息。箱線圖僅顯示匯總統(tǒng)計數據，例如均值/中位數和四分位距，而小提琴圖則顯示了數據的完整分布。當數據分布是多峰（多于一個峰）時，這種差異特別有用。在這種情況下，小提琴圖顯示了不同峰值的存在、它們的位置和相對幅度。
# Vln plot - cluster 20
VlnPlot(object = seurat_integrated,
features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"))
這些結果和圖可以幫助我們確定這些簇的身份，或驗證我們之前探索過預期細胞類型的規(guī)范標記后假設的身份。
識別每個簇的基因標記
我們分析的最后一組問題：對應于相同細胞類型的簇是否具有生物學意義的差異？有時返回的標記列表不能充分分離某些集群。例如，我們之前已經將集群 0、2、4、10 和 18 確定為 CD4+ T細胞，但是這些細胞集群之間是否存在生物學相關的差異？我們可以使用該FindMarkers()函數來確定兩個特定簇之間差異表達的基因。
我們可以嘗試所有比較的組合，但我們將從簇 2 與所有其他 CD4+ T 細胞簇開始：
# Determine differentiating markers for CD4+ T cell
cd4_tcells <- FindMarkers(seurat_integrated,
ident.1 = 2,
ident.2 = c(0,4,10,18))
# Add gene symbols to the DE table
cd4_tcells <- cd4_tcells %>%
rownames_to_column(var = "gene") %>%
left_join(y = unique(annotations[, c("gene_name", "description")]),
by = c("gene" = "gene_name"))
# Reorder columns and sort by padj
cd4_tcells <- cd4_tcells[, c(1, 3:5,2,6:7)]
cd4_tcells <- cd4_tcells %>%
dplyr::arrange(p_val_adj)
# View data
View(cd4_tcells)
在這些頂級基因中，CREM 基因作為激活標記脫穎而出。我們知道另一個激活標志物是 CD69，而幼稚或記憶細胞的標志物包括 SELL 和 CCR7 基因。有趣的是，SELL 基因也位居榜首。讓我們使用這些新的細胞狀態(tài)標記直觀地探索激活狀態(tài)：
| 細胞狀態(tài) | 標記 |
|---|---|
| 幼稚T細胞 | CCR7，賣出 |
| 活化的 T 細胞 | CREM，CD69 |
# Plot gene markers of activated and naive/memory T cells
FeaturePlot(seurat_integrated,
reduction = "umap",
features = c("CREM", "CD69", "CCR7", "SELL"),
label = TRUE,
order = TRUE,
min.cutoff = 'q10',
repel = TRUE
)
由于初始狀態(tài)和激活狀態(tài)的標記都顯示在標記列表中，因此可視化表達是有幫助的?；谶@些圖，集群 0 和 2 似乎是可靠的初始 T 細胞。然而，對于激活的 T 細胞，很難說。我們可能會說簇 4 和 18 是激活的 T 細胞，但 CD69 表達不如 CREM 明顯。我們將標記初始細胞并將剩余的簇標記為 CD4+ T 細胞。
現(xiàn)在利用所有這些信息，我們可以推測不同簇的細胞類型，并用細胞類型標簽繪制細胞。
| 集群 ID | 細胞類型 |
|---|---|
| 0 | 幼稚或記憶 CD4+ T 細胞 |
| 1 | CD14+ 單核細胞 |
| 2 | 幼稚或記憶 CD4+ T 細胞 |
| 3 | CD14+ 單核細胞 |
| 4 | CD4+ T 細胞 |
| 5 | CD8+ T 細胞 |
| 6 | B細胞 |
| 7 | 應激細胞/活化的T細胞 |
| 8 | 自然殺傷細胞 |
| 9 | FCGR3A+ 單核細胞 |
| 10 | CD4+ T 細胞 |
| 11 | B細胞 |
| 12 | 自然殺傷細胞 |
| 13 | CD8+ T 細胞 |
| 14 | CD14+ 單核細胞 |
| 15 | 常規(guī)樹突狀細胞 |
| 16 | 巨核細胞 |
| 17 | B細胞 |
| 18 | CD4+ T 細胞 |
| 19 | 漿細胞樣樹突狀細胞 |
| 20 | 肥大細胞 |
然后我們可以將集群的身份重新分配給這些細胞類型：
# Rename all identities
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")
# Plot the UMAP
DimPlot(object = seurat_integrated,
reduction = "umap",
label = TRUE,
label.size = 3,
repel = TRUE)
如果我們想刪除潛在的壓力細胞，我們可以使用以下subset()函數：
# Remove the stressed or dying cells
seurat_subset_labeled <- subset(seurat_integrated,
idents = "Stressed cells / Activated T cells", invert = TRUE)
# Re-visualize the clusters
DimPlot(object = seurat_subset_labeled,
reduction = "umap",
label = TRUE,
label.size = 3,
repel = TRUE)
現(xiàn)在我們想要保存我們最終標記的 Seurat 對象和輸出sessionInfo()：
# Save final R object
write_rds(seurat_integrated,
file = "results/seurat_labelled.rds")
# Create and save a text file with sessionInfo
sink("sessionInfo_scrnaseq_Oct2020.txt")
sessionInfo()
sink()
您可以在此鏈接中找到有關該
sink()功能的更多信息。
現(xiàn)在我們已經定義了我們的集群和每個集群的標記，我們有幾個不同的選擇：
