注釋是很關(guān)鍵的一步，這幾天先把注釋完全搞懂。現(xiàn)在有幾個問題需要解決：
1、搞清楚不同注釋方法（手動、網(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ì)胞類型注釋
注釋原則：
1、手動注釋
注釋網(wǎng)站 CellMarker
就是先搜索組織 再定位組織里面的細(xì)胞 打開之后，給到你某個組織中特定的細(xì)胞的 Marker 基因。然后用 Vlnplot 看marker在不同的群里面的特異性表達(dá)情況。 然后慢慢篩選、調(diào)整。整個過程需要對測序的細(xì)胞比較熟悉，尤其是腫瘤細(xì)胞。
marker 展示方法：
2、自動注釋（內(nèi)置參考集，也就是celldex包）-人工檢查
這一種注釋方法比較常用，最常用的是SingleR。使用時，需要下載celldex參考數(shù)據(jù)包。
使用內(nèi)置參考進(jìn)行注釋（最簡便的）
使用SingleR的最簡單方法是使用內(nèi)置參考對細(xì)胞進(jìn)行注釋。celldex包通過專用的檢索功能提供了7個參考數(shù)據(jù)集（主要來自大量RNA-seq或微陣列數(shù)據(jù)）。
singleR自帶的7個參考數(shù)據(jù)集，需要聯(lián)網(wǎng)才能下載，其中5個是人類數(shù)據(jù)，2個是小鼠的數(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、自動注釋（已發(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菌的代碼，然后找一篇文章自己做一次。
演示用到的數(shù)據(jù)集來自2020年發(fā)表在Nature Genetics的一篇結(jié)直腸癌研究。文中用到了韓國患者(SMC)和比利時患者(KUL3)兩套數(shù)據(jù)集，兩套數(shù)據(jù)平行分析，相互印證。
TOP菌只取了髓系細(xì)胞，用KUL3數(shù)據(jù)集中的髓系細(xì)胞為參考，來注釋SMC里面的髓系細(xì)胞。我這里用的是全部的細(xì)胞做注釋。
23例韓國大腸癌患者的單細(xì)胞（GSE132465 ）
6例比利時大腸癌患者的單細(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ì)胞，比利時人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（用$符號就可以看）一個，如果用SingleR自帶的包的話，可能有2個、3個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é)果：
可以看到用比利時人KUL3注釋的，和作者手動注釋的結(jié)果差不多。 TOP菌提到的一個很關(guān)鍵的因素就是：比利時人也就是參考集的細(xì)胞比較少，然后要注釋的韓國人組的細(xì)胞多，這樣可能不太好。最好是參考集的細(xì)胞數(shù)比待注釋的細(xì)胞多情況下比較準(zhǔn)確。
還有一個要解決的問題就是：看了很多帖子說可以同時參考幾個參考集，然后去注釋，這樣更準(zhǔn)確。但是沒找到詳細(xì)的參考資料。
