在完成對樣本的整合、聚類與可視化后，樣本細(xì)胞形成了許多cluster，需要得出這些cluster代表的細(xì)胞種類并將其注釋到圖上。

注釋方法分為手動注釋與軟件注釋兩大類

手動注釋

1. 尋找各cluster特異的高表達(dá)基因（marker）

markers <- FindAllMarkers(object=scRNA_harmony,test.use = "wilcox",only.pos = TRUE,
                          logfc.threshold = 0.25)
#篩選每個cluster特異表達(dá)基因表達(dá)量最高的前十個
all.markers = markers %>% select(gene,everything()) %>% subset(p_val<0.05)
top10 <- all.markers %>% group_by(cluster) %>% top_n(n=10,wt=avg_log2FC)

再查閱文獻(xiàn)資料找出top10 的marker基因?qū)?yīng)的細(xì)胞種類

2.尋找conserved marker

單細(xì)胞測序分析時經(jīng)常會出現(xiàn)不同分組，例如藥物處理與沒有藥物處理的上皮細(xì)胞基因表達(dá)也許存在差異，而這些差異基因不能作為marker，F(xiàn)indConservedMarker可以解決這個問題

#尋找兩個sample中共有的marker基因
#FindConservedMarke()每次只能處理一個cluster，下面以cluster1 為例
marker2 <- FindConservedMarkers(scRNA_harmony,ident.1 = 1,grouping.var = "orig.ident",
                                only.pos = TRUE,
                                min.diff.pct = 0.25,
                                min.pct = 0.25,
                                logfc.threshold = 0.25)

3.利用氣泡圖

#例如T細(xì)胞的marker基因IL7R，CD3D，NKG7
DotPlot(scRNA_harmony,features=c("IL7R","CD3D","NKG7"))

Snipaste_2022-08-11_12-07-01.png

可以根據(jù)氣泡圖看出，cluster2的T細(xì)胞marker基因表達(dá)量較高，所以可以被歸類為T細(xì)胞

4.重命名

用上述方法得出細(xì)胞種類后對cluster進(jìn)行重命名

#以cluster7為例
scRNA_harmony <- RenameIdents(scRNA_harmony,"7"="Osteoblastic")

軟件注釋

1.SingleR package

用已知的數(shù)據(jù)注釋未知數(shù)據(jù)

第一步 加載需要的package并load前文已經(jīng)處理好的scRNA.harmony數(shù)據(jù)

##系統(tǒng)報錯改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())

library(harmony)
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(org.Hs.eg.db)
library(SingleR)
library(scRNAseq)
load("../scRNA_harmony")

第二步 準(zhǔn)備參考數(shù)據(jù)與需要被注釋的數(shù)據(jù)

利用下載的已知的人類細(xì)胞圖譜的數(shù)據(jù)注釋未知的scRNA_harmony的cluster

#下載參考數(shù)據(jù)
refdata = HumanPrimaryCellAtlasData()

#提取出scRNA的轉(zhuǎn)錄表達(dá)數(shù)據(jù)
testdata <- GetAssayData(scRNA_harmony,slot="data")
#提取每個細(xì)胞的cluster信息
clusters <- scRNA_harmony@meta.data$seurat_clusters

#開始使用SingleR進(jìn)行分析
cellpred <- SingleR(test=testdata,ref=refdata,labels=refdata$label.main,
                    method="cluster",cluster=clusters,
                   assay.type.test = "logcounts",assay.type.ref = "logcounts")
#制作細(xì)胞類型的注釋文件
celltype = data.frame(ClusterID = rownames(cellpred),celltype=cellpred$labels,stringsAsFactors = FALSE)

第三步 將注釋結(jié)果加入metadata中

方法一

scRNA_harmony@meta.data$celltype="NA"
for(i in 1:nrow(celltype)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$seurat_clusters==celltype$ClusterID[i]),'celltype']
  <- celltype$celltype[i]
}

• 先在metadata中加入一列名為celltype的新列，本列所有數(shù)據(jù)都設(shè)置為NA
• 利用for循環(huán)，將metadata中每個cluster所含的細(xì)胞所在的行找出，并將該細(xì)胞對應(yīng)的celltype列賦上SingleR的標(biāo)記結(jié)果

方法二

scRNA_harmony@meta.data$SingleR = celltype[match(clusters,celltype$ClusterID),'celltype']

• 在metadata中加入名為SingleR的新列
• match(clusters,celltype$ClusterID)代表celltype中各個clusters所在的行數(shù)
• 再將對應(yīng)行數(shù)的celltype賦給SingleR列

2 Garnett方法

SingleR package注釋法中參考數(shù)據(jù)來源與package自帶的數(shù)據(jù)，因此有較大的局限性，如果要自定義參考數(shù)據(jù)則需要采用Garnett法。
Garnett法是利用基于機(jī)器學(xué)習(xí)訓(xùn)練分類器，利用訓(xùn)練好的分類器對細(xì)胞類型進(jìn)行分類。其本質(zhì)類似于氣泡圖注釋法，利用分類器對氣泡圖進(jìn)行更準(zhǔn)確的判斷。

第一步 制作marker file（txt文件）

格式如下
> B cells
expressed:CD19,MS4A1,CD79A,ACTN,ACTB
not expressed:
subtype of:

• 第一行是">"后接細(xì)胞種類
• 第二行是選定的marker基因
• 第三行是選定的陰性marker
• 第四行是細(xì)胞亞型
• 前兩行是必填，但是后兩行可以不填
  本次示例采用下載的hsPBMC marker file

第二步 創(chuàng)建CDS對象 優(yōu)化marker file

優(yōu)化的目的是使marker file更適合所分析的單細(xì)胞數(shù)據(jù)

#安裝garnett之前要先安裝monocle3
devtools::install_github("cole-trapnell-lab/monocle3")
library(monocle3)
devtools::install_github('cole-trapnell-lab/garnett',ref='monocle3')
library(garnett)
#加載scRNA_harmony數(shù)據(jù)
load("scRNA_harmony.rdata")
pbmc <- scRNA_harmony
##創(chuàng)建CDS對象
data <- GetAssayData(pbmc,assay='RNA',slot="counts")
cell_metadata <- pbmc@meta.data
#gene_annotation實(shí)際上沒有任何實(shí)際意義，僅僅是為了滿足輸入要求
gene_annotation <- data.frame(gene_short_name=rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
                         cell_metadata = cell_metadata,
                         gene_metadata=gene_annotation)
#對CDS對象進(jìn)行歸一化與降維處理，preprocess_cds函數(shù)相當(dāng)于seurat中NormalizeData+ScaleData+RunPCA
cds <- preprocess_cds(cds,num_dim=30)
#優(yōu)化marker file
marker_check <- check_markers(cds,"hsPBMC_markers.txt",
                              db=org.Hs.eg.db,
                              cds_gene_id_type = "SYMBOL",
                              marker_file_gene_id_type = "SYMBOL")
plot_markers(marker_check)

得到的結(jié)果如圖所示

示例.png

根據(jù)圖片，不在database中的marker和high overlap（即特異性不高）的marker都需要刪除替換

第三步 得到分類器

分類器一般有兩種獲得途徑，可以上Garne官網(wǎng)下載已經(jīng)訓(xùn)練好的分類器（資源較少，一般找不到所需分類器），另一種方法就是自行訓(xùn)練

#使用marker file和cds對象訓(xùn)練分類器
pbmc_classifier <-train_cell_classifier(cds=cds,
                                        marker_file = "hsPBMS_markers.txt",
                                        db=org.Hs.eg.db,
                                        cds_gene_id_type = "SYMBOL",
                                        num_unknown = 50,
                                        marker_file_gene_id_type = "SYMBOL")
saveRDS(pbmc_classifier,"my_classifier.rds")

第四步 使用訓(xùn)練器進(jìn)行注釋

#讀取marker file
hsPBMC <- readRDS("hsPBMC.rds")
pData(cds)$garnett_cluster <- pData(cds)$seurat_clusters
cds <- classify_cells(cds,
                      hsPBMC,
                      db = org.Hs.eg.db,
                      cluster_extend = TRUE,
                      cds_gene_id_type = "SYMBOL")
#提取分類結(jié)果
cds.meta <- subset(pData(cds),select=c('cell_type','cluster_ext_type')) %>% as.data.frame()
#將結(jié)果返回給seurat對象
pbmc <- AddMetaData(pbmc,metadata = cds.meta)

3 nnls（非負(fù)最小二乘回歸）法

計(jì)算需要被注釋的數(shù)據(jù)與已知的參考數(shù)據(jù)中哪些類型相關(guān)性更大

公式.png

其中Ta表示需要被注釋的數(shù)據(jù)集A中的細(xì)胞的基因表達(dá)，Mb表示參考數(shù)據(jù)集中的細(xì)胞基因表達(dá)

##系統(tǒng)報錯改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

讀取數(shù)據(jù)并進(jìn)行歸一化降維處理

dir=c("BC2/","BC21/")
names(dir)=c("sample2","sample21")
#讀取數(shù)據(jù)并進(jìn)行歸一化降維處理
scRNA.list <- list()
for(i in 1:length(dir)){
  counts=Read10X(data.dir = dir[i])
  scRNA.list[[i]]<- CreateSeuratObject(counts,min.cells = 3,min.features = 2000)
  
}

seu.obj <- scRNA.list[[1]]
seu.obj <- seu.obj %>% NormalizeData(verbose=FALSE) %>% FindVariableFeatures(selection.method='vst') %>% ScaleData(verbose=FALSE)%>% RunPCA(pc.genes=seu.obj@var.genes,npcs=100,verbose=FALSE)%>% FindNeighbors(dims=1:10) %>% FindClusters(resolution=0.5) %>% RunUMAP(dims=1:10)

seu.obj2 = scRNA.list[[2]]
seu.obj2 <- seu.obj2 %>% NormalizeData(verbose=FALSE) %>% FindVariableFeatures(selection.method='vst')%>%ScaleData(verbose=FALSE) %>% RunPCA(pc.genes=seu.obj2@var.genes,npcs=100,verbose=FALSE) %>%FindNeighbors(dims=1:10) %>% FindClusters(resolution=0.5) %>% RunUMAP(dims=1:10)

#找到兩個數(shù)據(jù)集共有的基因
shared_gene <- intersect(rownames(seu.obj),row.names(seu.obj2))

#計(jì)算所有基因在每個cluster中的表達(dá)量之和
seu.obj.mat <- AggregateExpression(seu.obj,assays="RNA",features=shared_gene)$RNA
#除以測序深度進(jìn)行normalization
seu.obj.mat <- seu.obj.mat / rep(colSums(seu.obj.mat),each=nrow(seu.obj.mat))
seu.obj.mat <- log10(seu.obj.mat * 100000+1)

seu.obj.mat2 <- AggregateExpression(seu.obj2,assays="RNA",features=shared_gene)$RNA
seu.obj.mat2 <- seu.obj.mat2 / rep(colSums(seu.obj.mat2),each=nrow(seu.obj.mat2))
seu.obj.mat2 <- log10(seu.obj.mat2 * 100000+1)

篩選基因

根據(jù)目標(biāo)數(shù)據(jù)集（A數(shù)據(jù)集中選定的cluster）中給定的細(xì)胞類型與整體細(xì)胞類型（數(shù)據(jù)集A）中位數(shù)的倍數(shù)變化，選擇前200個，得到list1；然后根據(jù)目標(biāo)數(shù)據(jù)集中給定的細(xì)胞類型與其他細(xì)胞類型（去除選定的cluster的數(shù)據(jù)集A）最大值的倍數(shù)變化，選擇前200個，得到list2，然后將兩個list進(jìn)行合并

#以cluster2為為例
cluster <- 2
seu.obj.gene <- seu.obj.mat[,cluster+1]
#得出list1
gene_fc <- seu.obj.gene / apply(seu.obj.mat,1,median)
gene_list1 <- names(sort(gene_fc,decreasing=TRUE)[1:200])
#得出list2
gene_fc <- seu.obj.gene/apply(seu.obj.mat[,-(cluster+1)],1,max)
gene_list2 <- names(sort(gene_fc,decreasing=TRUE)[1:200])
#合并兩個list
gene_list <- unique(c(gene_list,gene_list2))

進(jìn)行相關(guān)性計(jì)算

#A數(shù)據(jù)集中的某個細(xì)胞類型
Ta <- seu.obj.mat[gene_list,cluster+1]
#B數(shù)據(jù)集中的所有細(xì)胞類型
Mb <- seu.obj.mat2[gene_list,]
library(lsei)
solv <- nnls(Mb,Ta)
corr <- solv$x
corr

運(yùn)行結(jié)果中B數(shù)據(jù)集里與cluster2相關(guān)性最強(qiáng)的cluster（根據(jù)運(yùn)行結(jié)果可知是cluster6）可以視作與cluster同一種類的細(xì)胞

進(jìn)行驗(yàn)證

由結(jié)果可知，A數(shù)據(jù)集中的cluster2與B數(shù)據(jù)集中的cluster6相關(guān)性最高，接下來嘗試找出B數(shù)據(jù)集中的cluster6與A數(shù)據(jù)集中的哪個cluster相關(guān)性最高（即將上述過程反向操作）

cluster=6
seu.obj.gene <- seu.obj.mat2[,cluster+1]

gene_fc <- seu.obj.gene / apply(seu.obj.mat2,1,median)
gene_list1 <- names(sort(gene_fc,decreasing=TRUE)[1:200])

gene_fc <- seu.obj.gene/apply(seu.obj.mat[,-(cluster+1)],1,max)
gene_list2 <- names(sort(gene_fc,decreasing=TRUE)[1:200])

gene_list <- unique(c(gene_list,gene_list2))   

Ta <- seu.obj.mat2[gene_list,cluster+1]
Mb <- seu.obj.mat[gene_list,]
solv <- nnls(Mb,Ta)
corr <- solv$x
corr

批量計(jì)算相關(guān)系數(shù)

上述是計(jì)算單個cluster相關(guān)系數(shù)的pipeline
接下來對數(shù)據(jù)集中多個cluster進(jìn)行批量計(jì)算

library(lsei)
list1 <- list()
for (cluster in seq(1,ncol(seu.obj.mat))) {
  seu.obj.gene <- seu.obj.mat[,cluster]
  
  gene_fc <- seu.obj.gene / apply(seu.obj.mat,1,median)
  gene_list1 <- names(sort(gene_fc,decreasing=TRUE)[1:200])
  
  gene_fc <- seu.obj.gene/apply(seu.obj.mat[,-cluster],1,max)
  gene_list2 <- names(sort(gene_fc,decreasing=TRUE)[1:200])

  gene_list <- unique(c(gene_list,gene_list2))   

  Ta <- seu.obj.mat[gene_list,cluster]
  Mb <- seu.obj.mat2[gene_list,]
  solv <- nnls(Mb,Ta)
  corr <- solv$x
  corr
  list1[[cluster]] <- corr
}
#用a預(yù)測b
list2 <- list()
for (cluster in seq(1,ncol(seu.obj.mat2))) {
  seu.obj.gene <- seu.obj.mat2[,cluster]
  
  gene_fc <- seu.obj.gene / apply(seu.obj.mat2,1,median)
  gene_list1 <- names(sort(gene_fc,decreasing=TRUE)[1:200])
  gene_fc <- seu.obj.gene/apply(seu.obj.mat2[,-cluster],1,max)
  gene_list2 <- names(sort(gene_fc,decreasing=TRUE)[1:200])

  gene_list <- unique(c(gene_list,gene_list2))   
  Ta <- seu.obj.mat2[gene_list,cluster]
  Mb <- seu.obj.mat[gene_list,]
  solv <- nnls(Mb,Ta)
  corr <- solv$x
  corr
  list2[[cluster]] <- corr
} 
#合并結(jié)果
mat1 <- do.call(rbind,list1)
mat2 <- do.call(rbind,list2)

#計(jì)算系數(shù)
beta <- 2 * (mat1+0.01) * t(mat2+0.01)
row.names(beta) <- paste0("C",1:nrow(beta)-1)
colnames(beta)<- paste0("C",1:ncol(beta)-1)

進(jìn)行可視化

#可視化
pheatmap::pheatmap(beta,cluster_rows = FALSE,cluster_cols = FALSE)

熱圖色塊顏色越深對應(yīng)的兩個cluster的相關(guān)性越高

根據(jù)熱圖來優(yōu)化聚類結(jié)果

library(psych)
colnames(seu.obj@meta.data)
Idents(seu.obj)="seurat_clusters"
#計(jì)算表達(dá)矩陣每個cluster的平均值
exp = AverageExpression(seu.obj)
#畫出相關(guān)性熱圖
corrda <- corr.test(exp$RNA,exp$RNA,method="spearman")
pheatmap::pheatmap(corrda$r)
#畫出聚類圖
DimPlot(seu.obj,label=T)

根據(jù)相關(guān)性熱圖可以判斷聚類圖上距離較近的cluster是否可以歸于同一個cluster
示例：

聚類圖.png

相關(guān)性熱圖.png

根據(jù)熱圖聚類情況可知，cluster0,6,9的相關(guān)性很強(qiáng)，在聚類圖中距離也很相近，可以被歸為一個cluster