注釋是很關(guān)鍵的一步，這幾天先把注釋完全搞懂∏捎椋現(xiàn)在有幾個(gè)問題需要解決：

1碉怔、搞清楚不同注釋方法（手動(dòng)、網(wǎng)站禁添、用已發(fā)表的文章里面的去注釋同一測序樣品）撮胧；
2、如何提取細(xì)胞亞群老翘，進(jìn)行亞群注釋

單細(xì)胞數(shù)據(jù)經(jīng)過降維處理后芹啥，得到的tesn和umap圖就是根據(jù)表達(dá)情況，把表達(dá)相似的細(xì)胞聚在一起（也就是不同的cluster）铺峭。注釋這一步其實(shí)就是給這些表達(dá)相似的細(xì)胞群起名字墓怀。

參考：

【生信純干貨】 生信分析手把手教學(xué) 第六課
第五講-基于Seurat的細(xì)胞注釋分析
單細(xì)胞數(shù)據(jù)分析系列講座（三）：細(xì)胞類型注釋

注釋原則：

image.png

image.png

image.png

1、手動(dòng)注釋

注釋網(wǎng)站 CellMarker
就是先搜索組織 再定位組織里面的細(xì)胞 打開之后卫键，給到你某個(gè)組織中特定的細(xì)胞的 Marker 基因傀履。然后用 Vlnplot 看marker在不同的群里面的特異性表達(dá)情況。 然后慢慢篩選莉炉、調(diào)整钓账。整個(gè)過程需要對測序的細(xì)胞比較熟悉，尤其是腫瘤細(xì)胞絮宁。

marker 展示方法：

image.png

image.png

image.png

image.png

image.png

image.png

2梆暮、自動(dòng)注釋（內(nèi)置參考集，也就是celldex包）-人工檢查

這一種注釋方法比較常用绍昂，最常用的是SingleR啦粹。使用時(shí)，需要下載celldex參考數(shù)據(jù)包窘游。

使用內(nèi)置參考進(jìn)行注釋（最簡便的）
使用SingleR的最簡單方法是使用內(nèi)置參考對細(xì)胞進(jìn)行注釋卖陵。celldex包通過專用的檢索功能提供了7個(gè)參考數(shù)據(jù)集（主要來自大量RNA-seq或微陣列數(shù)據(jù)）。
singleR自帶的7個(gè)參考數(shù)據(jù)集张峰，需要聯(lián)網(wǎng)才能下載泪蔫，其中5個(gè)是人類數(shù)據(jù)，2個(gè)是小鼠的數(shù)據(jù)：
BlueprintEncodeData Blueprint (Martens and Stunnenberg 2013) and Encode (The ENCODE Project Consortium 2012) （人）
DatabaseImmuneCellExpressionData The Database for Immune Cell Expression(/eQTLs/Epigenomics)(Schmiedel et al. 2018)（人）
HumanPrimaryCellAtlasData the Human Primary Cell Atlas (Mabbott et al. 2013)（人）
MonacoImmuneData, Monaco Immune Cell Data - GSE107011 (Monaco et al. 2019)（人）
NovershternHematopoieticData Novershtern Hematopoietic Cell Data - GSE24759（人）
ImmGenData the murine ImmGen (Heng et al. 2008) （鼠）
MouseRNAseqData a collection of mouse data sets downloaded from GEO (Benayoun et al. 2019).鼠）

pred<- SingleR(test = norm_count,   #輸入表達(dá)矩陣
               ref = ref,           #參考數(shù)據(jù)
               clusters = scRNA$originalexp_snn_res.0.2,
               labels = ref$label.main)
head(pred)
table(pred$labels)
scRNA$singleR_cluster = pred$labels[match(scRNA$originalexp_snn_res.0.2,
                                          rownames(pred))]
table(scRNA$singleR_cluster)

3喘批、自動(dòng)注釋（已發(fā)表文章作為參考集）-人工檢查

前兩種注釋方法介紹的帖子非常多撩荣，這里主要關(guān)注一下自己用已發(fā)表的文章制作參考集铣揉，然后用singleR做注釋的方法。 參考復(fù)現(xiàn)的是全網(wǎng)介紹最詳細(xì)的 TOP生物信息 的帖子餐曹，參考帖子鏈接如下：

單細(xì)胞輔助注釋工具-SingleR
SingleR如何使用自定義的參考集
使用SingleR包進(jìn)行單細(xì)胞類型注釋分析
【單細(xì)胞】SingleR注釋細(xì)胞類型
學(xué)習(xí)對scRNA-seq數(shù)據(jù)注釋的R包SingleR
Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer

先復(fù)現(xiàn)TOP菌的代碼逛拱，然后找一篇文章自己做一次。

image.png

演示用到的數(shù)據(jù)集來自2020年發(fā)表在Nature Genetics的一篇結(jié)直腸癌研究台猴。文中用到了韓國患者(SMC)和比利時(shí)患者(KUL3)兩套數(shù)據(jù)集朽合，兩套數(shù)據(jù)平行分析，相互印證饱狂。

TOP菌只取了髓系細(xì)胞曹步，用KUL3數(shù)據(jù)集中的髓系細(xì)胞為參考，來注釋SMC里面的髓系細(xì)胞休讳。我這里用的是全部的細(xì)胞做注釋讲婚。

23例韓國大腸癌患者的單細(xì)胞（GSE132465 ）
6例比利時(shí)大腸癌患者的單細(xì)胞（GSE144735）

代碼如下：

對應(yīng)TOP菌的 0.R

library(Seurat)
library(tidyverse)
library(data.table)

setwd(dir="D:/SingleR自定義-TOP-0801")

### 表達(dá)數(shù)據(jù)和注釋可以到NCBI下載 ##########
#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132465
#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144735


#1、這一步是全部的細(xì)胞俊柔，TOP菌只用了髓系細(xì)胞筹麸，比利時(shí)人KUL3的注釋文件
KUL3.anno1=read.table("D:/SingleR自定義-TOP-0801/下載的原始數(shù)據(jù)/GSE144735_processed_KUL3_CRC_10X_annotation.txt/GSE144735_processed_KUL3_CRC_10X_annotation.txt",header = T,sep = "\t",stringsAsFactors = F)
KUL3.all.anno=KUL3.anno1%>%filter(Cell_subtype!="Unknown" & Cell_subtype!="Proliferating")

#2、讀取原始數(shù)據(jù)  fread和read.table是一樣的
mat=fread("D:/SingleR自定義-TOP-0801/下載的原始數(shù)據(jù)/GSE144735_processed_KUL3_CRC_10X_raw_UMI_count_matrix .txt/GSE144735_processed_KUL3_CRC_10X_raw_UMI_count_matrix.txt")
mat=as.data.frame(mat)
rownames(mat)=mat$Index
mat$Index=NULL

#3雏婶、制作KUL3文件
KUL3.all.mat=mat[,KUL3.all.anno$Index]
rm(list = c("KUL3.anno1","mat"))
###################

#4物赶、KUL3的Seurat標(biāo)準(zhǔn)流程
BLS <- CreateSeuratObject(counts = KUL3.all.mat, project = "BLS") %>%
  Seurat::NormalizeData() %>%
  FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>% 
  ScaleData()
BLS <- RunPCA(BLS, npcs = 50, verbose = FALSE)

BLS@meta.data$Index=rownames(BLS@meta.data)
BLS@meta.data=inner_join(BLS@meta.data,KUL3.all.anno,by="Index")
rownames(BLS@meta.data)=BLS@meta.data$Index

BLS <- FindNeighbors(BLS, dims = 1:20)
BLS <- FindClusters(BLS, resolution = 0.5)

BLS <- RunUMAP(BLS, dims = 1:20)
BLS <- RunTSNE(BLS, dims = 1:20)

DimPlot(BLS, reduction = "tsne", group.by = "Cell_subtype", pt.size=2)
DimPlot(BLS, reduction = "tsne", group.by = "Class", pt.size=2)

saveRDS(BLS,file = "KUL3_BLS.rds")


#######################

#5、開始做韓國人的留晚，
SMC.anno1=read.table("D:/SingleR自定義-TOP-0801/下載的原始數(shù)據(jù)/GSE132465_GEO_processed_CRC_10X_cell_annotation.txt/GSE132465_GEO_processed_CRC_10X_cell_annotation.txt",header = T,sep = "\t",stringsAsFactors = F)
SMC.all.anno=SMC.anno1%>%filter(Cell_subtype!="Unknown" & Cell_subtype!="Proliferating")

mat=fread("D:/SingleR自定義-TOP-0801/下載的原始數(shù)據(jù)/GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix .txt/GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix.txt",select = c("Index",SMC.all.anno$Index))
mat=as.data.frame(mat)
rownames(mat)=mat$Index
mat$Index=NULL

SMC.all.mat=mat[,SMC.all.anno$Index]
rm(list = c("SMC.anno1","mat")
###################
Han <- CreateSeuratObject(counts = SMC.all.mat, project = "Han") %>%
  Seurat::NormalizeData() %>%
  FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>% 
  ScaleData()
Han <- RunPCA(Han, npcs = 50, verbose = FALSE)

對應(yīng)TOP菌的 1.R

library(Seurat)
library(tidyverse)
library(SummarizedExperiment)
library(scuttle)

KUL3_BLS=readRDS("KUL3_BLS.rds")
KUL3_BLS_count=KUL3_BLS[["RNA"]]@counts  #提取counts矩陣

#準(zhǔn)備注釋信息
pdata=KUL3_BLS@meta.data[,c("Index","Cell_subtype")]
rownames(pdata)=pdata$Index
pdata$Index=NULL
colnames(pdata)="ref_label"

KUL3_BLS_SE <- SummarizedExperiment(assays=list(counts=KUL3_BLS_count),colData = pdata)
KUL3_BLS_SE <- logNormCounts(KUL3_BLS_SE) #Compute log-transformed normalized expression values #參考集一定是logcounts的數(shù)據(jù)

saveRDS(KUL3_BLS_SE,"KUL3_BLS_SE.ref.rds")

對應(yīng)TOP菌的 2.R

library(tidyverse)
library(Seurat)
library(SingleR)
library(scuttle)
library(SummarizedExperiment)

KUL3_BLS_SE=readRDS("KUL3_BLS_SE.ref.rds")  #導(dǎo)入?yún)⒖技?SMC_Han=readRDS("SMC_Han.rds")  #待查詢的數(shù)據(jù)集酵紫，韓國人的

#下面開始是構(gòu)建待查詢數(shù)據(jù)集的SE對象
SMC_Han_count=SMC_Han[["RNA"]]@counts  # 導(dǎo)出表達(dá)矩陣counts

# 下面這三行是把帶查詢和參考數(shù)據(jù)集 取交集，注釋的方法依賴的就是交集基因
common_gene <- intersect(rownames(SMC_Han_count), rownames(KUL3_BLS_SE))
KUL3_BLS_SE <- KUL3_BLS_SE[common_gene,]
SMC_Han_count <- SMC_Han_count[common_gene,]

SMC_Han_SE <- SummarizedExperiment(assays=list(counts=SMC_Han_count)) #創(chuàng)建SE對象
SMC_Han_SE <- logNormCounts(SMC_Han_SE)   # 對創(chuàng)建SE對象倔丈，做logNormCounts處理

#開始SingleR注釋： labels參數(shù)的意思是:用參考集里面的說明labels去做注釋,
#由于制作的參考集只有ref_label（用$符號就可以看）一個(gè)憨闰，如果用SingleR自帶的包的話状蜗，可能有2個(gè)需五、3個(gè)labels。例如SingleR的 main大類轧坎、fine小類宏邮。  
singleR_res <- SingleR(test = SMC_Han_SE, ref = KUL3_BLS_SE, labels = KUL3_BLS_SE$ref_label) 

anno_df <- as.data.frame(singleR_res$labels) #將labels轉(zhuǎn)換成數(shù)據(jù)框
anno_df$Index <- rownames(singleR_res)  #給數(shù)據(jù)框加行名
colnames(anno_df)[1] <- 'ref_label_from_KUL3'  # 換列名


#把得到的注釋信息整合到seurat對象里面
SMC_Han@meta.data=SMC_Han@meta.data%>%inner_join(anno_df,by="Index")   # 屬性數(shù)據(jù)庫和anno_df有公共列Index
rownames(SMC_Han@meta.data)=SMC_Han@meta.data$Index                    #屬性數(shù)據(jù)框是沒有行名的，再把Index賦給行名

DimPlot(SMC_Han, reduction = "tsne", group.by = "Cell_subtype", pt.size=2)+
DimPlot(SMC_Han, reduction = "tsne", group.by = "ref_label_from_KUL3", pt.size=2)

最后結(jié)果：

可以看到用比利時(shí)人KUL3注釋的缸血，和作者手動(dòng)注釋的結(jié)果差不多蜜氨。 TOP菌提到的一個(gè)很關(guān)鍵的因素就是：比利時(shí)人也就是參考集的細(xì)胞比較少，然后要注釋的韓國人組的細(xì)胞多捎泻，這樣可能不太好飒炎。最好是參考集的細(xì)胞數(shù)比待注釋的細(xì)胞多情況下比較準(zhǔn)確。

還有一個(gè)要解決的問題就是：看了很多帖子說可以同時(shí)參考幾個(gè)參考集笆豁，然后去注釋郎汪，這樣更準(zhǔn)確赤赊。但是沒找到詳細(xì)的參考資料。

作者手動(dòng)注釋的圖

用KUL3比利時(shí)人注釋的