隨著單細(xì)胞技術(shù)的發(fā)展和普及伶唯，越來越多的課題組開始使用單細(xì)胞技術(shù)進(jìn)行自己的課題研究恤左。然而蔓肯，在單細(xì)胞的注釋上遭铺，商業(yè)公司提供的流程化注釋往往難以滿足實(shí)際的需求丽柿，甚至得到牛頭不對馬嘴的結(jié)果。這期講帶大家實(shí)現(xiàn)全方位的單細(xì)胞注釋指南魂挂。

1. 引言

單細(xì)胞RNA測序（scRNA-seq）是一種先進(jìn)的分子生物學(xué)技術(shù)甫题，可以在單個細(xì)胞水平上評估基因表達(dá)。在過去的幾年里涂召，scRNA-seq在細(xì)胞注釋中扮演了重要角色坠非，它已經(jīng)成為揭示細(xì)胞異質(zhì)性和發(fā)展過程中轉(zhuǎn)錄組動態(tài)的有力工具。目前果正，常見的scRNA-seq注釋方法包括基于表達(dá)模式的細(xì)胞標(biāo)記麻顶、參考基因表達(dá)和聚類特征等。
然而舱卡，盡管scRNA-seq可以在單細(xì)胞的水平上揭示基因的表達(dá)辅肾，但是目前在細(xì)胞注釋上依然存在不小的挑戰(zhàn)。首先轮锥，技術(shù)本身的特點(diǎn)如擴(kuò)增偏差矫钓、噪音和低覆蓋率等，可能導(dǎo)致結(jié)果的偏差和不確定性舍杜。其次新娜，數(shù)據(jù)的處理和分析也是一個復(fù)雜而關(guān)鍵的步驟，如何準(zhǔn)確地鑒定和分類不同的細(xì)胞類型仍然是一個挑戰(zhàn)既绩。此外概龄，由于細(xì)胞的異質(zhì)性和動態(tài)性，以及在特定生理狀態(tài)下轉(zhuǎn)錄組的變化饲握，細(xì)胞狀態(tài)鑒定和動態(tài)分析也是進(jìn)一步研究的難點(diǎn)私杜。

2.常見注釋策略

細(xì)胞注釋按照是否有人工參與可以分為自動注釋和人工注釋。

2.1 自動注釋

最為常見的就是使用SingleR進(jìn)行救欧，其根據(jù)參考數(shù)據(jù)集的表達(dá)矩陣衰粹，和輸入數(shù)據(jù)表達(dá)矩陣的相似性來識別細(xì)胞類型。其使用非常的方便笆怠，是絕大多數(shù)商業(yè)公司提供的全自動流程方案铝耻。然而，缺點(diǎn)是顯而易見的蹬刷。首先瓢捉，SingleR包含的參考數(shù)據(jù)集有限频丘，難以滿足實(shí)際樣本的需求，參考的有限導(dǎo)致注釋出現(xiàn)牛頭不對馬腳的結(jié)果泡态。
另一類策略則是根據(jù)Marker基因進(jìn)行椎镣，這里常用的軟件有cellassign和scCATCH等，他們的思路是類似的兽赁，首先需要我們輸入?yún)⒖嫉膍arker基因，marker基因可以從cell marker2.0冷守、sctype或已經(jīng)發(fā)表文獻(xiàn)等獲得刀崖。這里推薦大家使用scCATCH，其自帶了cell marker的數(shù)據(jù)庫對新手較為友好拍摇。通過marker基因進(jìn)行細(xì)胞分群的鑒定可以特異性的根據(jù)物種亮钦、組織的特點(diǎn)，注釋到更為合理的結(jié)果充活。然而蜂莉，這也并非“萬全之策”，由于樣本之間混卵、組織之間甚至是批次之間映穗，分析策略之間等一系列因素的差異，適用于別人的marker不見得就一定適用于自己的樣本幕随。

2.2 人工注釋

人工注釋蚁滋，顧名思義，就是自己去注釋赘淮。如何自己去注釋呢辕录？手動注釋在很大程度上比較類似基于marker基因注釋的方法∩倚叮基于參考文獻(xiàn)和以往經(jīng)驗(yàn)得到的marker基因的表達(dá)豐度進(jìn)行走诞。

3.實(shí)戰(zhàn)

3.1 SingleR自動注釋

使用的數(shù)據(jù)是經(jīng)典的pbmc3k （https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）
安裝好 R并安裝R包
Seurat、SingleR蛤高、tidyverse等蚣旱。
首先我們導(dǎo)入數(shù)據(jù)，pbmc解壓出來的數(shù)據(jù)是包括三個文件barcodes.tsv genes.tsv matrix.mtx戴陡。首先是進(jìn)行導(dǎo)入數(shù)據(jù) -> 質(zhì)控處理 -> 標(biāo)準(zhǔn)化 -> 降緯 -> 分群

library(Seurat)
library(tidyverse)
library(SingleR)

pbmc <- Read10X('~/pbmc/filtered_gene_bc_matrices/hg19') #這里填寫解壓后的目錄
pbmc <- CreateSeuratObject(pbmc) #轉(zhuǎn)化為Seurat對象
Pattern = '^MT-' #線粒體基因的名字姻锁，根據(jù)實(shí)際去匹配，人是^MT-開頭猜欺，小鼠^mt-位隶，大鼠^Mt-
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = Pattern)
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt <  15) #進(jìn)行細(xì)胞過濾，這里是要求每個細(xì)胞里面至少有200個的基因开皿，但不大于8000個基因涧黄，線粒體基因含量不大于15%篮昧。這里的過濾就根據(jù)實(shí)際來，沒有絕對的標(biāo)準(zhǔn)笋妥。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000) #進(jìn)行l(wèi)og歸一化
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000) #尋找高變基因
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes) #根據(jù)基因數(shù)進(jìn)行數(shù)據(jù)的縮放
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) #進(jìn)行PCA降緯懊昨，特征是我們前面找到的高變基因
pbmc <- FindNeighbors(pbmc, dims = 1:30) #尋找臨近的特征點(diǎn)
pbmc <- FindClusters(object = pbmc, verbose = T, resolution = 0.7) #分群，關(guān)鍵在于resolution調(diào)整分辨率春宣，在0～1之間酵颁，約大就分越多分群

這樣我們就得到就將細(xì)胞進(jìn)行了初步的分群，隨后開始導(dǎo)入SingleR的參數(shù)數(shù)據(jù)集月帝。
也有幾種方法躏惋，一種我們手動下載，一種是安裝``celldex```包導(dǎo)入嚷辅， 如 celldex::HumanPrimaryCellAtlasData簿姨，或者導(dǎo)入別人已經(jīng)有注釋的rds,方法也非常簡單。只需要

ref <- celldex::HumanPrimaryCellAtlasData() #用celldex包導(dǎo)入簸搞，第一次需要一點(diǎn)時間下載扁位，該數(shù)包還有其他更多的參考數(shù)據(jù)集，可以詳細(xì)看手冊
ref <- readRDS('pbmc.rds') #用別人注釋好的rds趁俊，需要先準(zhǔn)備域仇，沒有準(zhǔn)備的話用方法一
pbmc_anno <- SingleR(test = pbmc@assays$RNA@data, ref = ref, labels = ref$label.main, cluster = pbmc$seurat_clusters)  #這里我們選擇按照分群去注釋

這個時候我們可以看到，已經(jīng)注釋好了寺擂。主要有T細(xì)胞和B細(xì)胞殉簸，單核細(xì)胞和血小板等我們進(jìn)行一下可視化展示

image.png

pbmc <- RunUMAP(object = pbmc, dims = 1:30) #首先進(jìn)行UMAP降緯
pbmc$CellType <- pbmc_anno$labels[match(pbmc$seurat_clusters, row.names(pbmc_anno))] #加上SingleR的注釋結(jié)果
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T) 
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T) 
p1+p2
ggsave('CellAnnot_umap.png', height = 4.5, width = 12, dpi = 300)

CellAnnot_umap.png

我們對比一下官方給出的手動注釋結(jié)果，可以看到已經(jīng)是基本接近了沽讹。

image.png

3.2 scCATCH注釋

這是我個人更推薦的一種方法般卑，首先主要安裝 scCATCH包，我們接著上面的代碼爽雄。

library(scCATCH)
obj <- createscCATCH(data = pbmc@assays$RNA@data, cluster = as.character(pbmc$seurat_clusters) ) #創(chuàng)建CATCH對象
obj <- findmarkergene(obj, marker = cellmatch, if_use_custom_marker = F ,use_method = "2",logfc = 0.4, cell_min_pct = 0.25, pvalue = 0.05, cancer = 'Normal', species = 'Human', tissue = 'Blood') #這里需要注意的關(guān)鍵參數(shù)是marker我們使用了自帶的cellmatch蝠检， 指定好癌癥的類型，物種和組織挚瘟，我們可以通過unique(cellmatch$tissue)去查看自帶的組織等叹谁。
obj <- findcelltype(obj)
obj@celltype$cell_type

我們看看效果

marker <- tibble( cluster = parse_number(obj@celltype$cluster), cell_tpye =  obj@celltype$cell_type) 
pbmc$CellType <- marker$cell_tpye[match(pbmc$seurat_clusters,marker$cluster)]
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T) 
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T,label.size=2) 
p1+p2
ggsave('CellAnnot_umap2.png', height = 4.5, width = 12, dpi = 300)

CellAnnot_umap2.png

效果也是很不錯。那么如果我們需要參考別人的結(jié)果如何進(jìn)行呢乘盖？這里還是只需要借助scCATCH包即可焰檩。我們需要準(zhǔn)備以下樣式的表格，包含以下的行頭订框，比如我們使用pbmc3k提供的marker
celltype gene subtype1 subtype2 subtype3 pmid

對應(yīng)代碼如下：

cell_marker <- tribble(
~celltype, ~gene,  ~subtype1,~subtype2,~subtype3, ~pmid,
  "T cell", "IL7R",  "Naive CD4+", NA, NA, NA,
  "T cell", "CCR7",  "Naive CD4+", NA, NA, NA,
  "Mono", "CD14",  "CD14+", NA, NA, NA,
  "Mono", "LYZ",  "CD14+", NA, NA, NA,
  "T cell", "IL7R",  "Memory CD4+", NA, NA, NA,
  "T cell", "S100A4",  "Memory CD4+", NA, NA, NA,
  "B cell", "MS4A1",  NA, NA, NA, NA,
  "T cell", "CD8A",  "CD8+", NA, NA, NA,
  "Mono", "FCGR3A",  "FCGR3A+", NA, NA, NA,
  "Mono", "MS4A7",  "FCGR3A+", NA, NA, NA,
  "NK", "GNLY",  NA, NA, NA, NA,
  "NK", "NKG7",  NA, NA, NA, NA,
  "DC", "FCER1A",  NA, NA, NA, NA,
  "DC", "CST3",  NA, NA, NA, NA,
  "Platelet", "PPBP",  NA, NA, NA, NA,
)
obj <- findmarkergene(obj, marker = cell_marker, if_use_custom_marker = T , use_method = "1",logfc = 0.4, cell_min_pct = 0.25, pvalue = 0.05)
obj <- findcelltype(obj)
marker <- tibble( cluster = parse_number(obj@celltype$cluster), cell_tpye =  obj@celltype$cell_type) 
pbmc$CellType <- marker$cell_tpye[match(pbmc$seurat_clusters,marker$cluster)]
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T) 
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T,label.size=2) 
p1+p2
ggsave('CellAnnot_umap3.png', height = 4.5, width = 12, dpi = 300)

CellAnnot_umap3.png