scRNA基礎(chǔ)分析-1:安裝包焰枢、導(dǎo)入數(shù)據(jù)、過(guò)濾質(zhì)控 - 簡(jiǎn)書(shū)
scRNA基礎(chǔ)分析-2：降維與聚類(lèi) - 簡(jiǎn)書(shū)
scRNA基礎(chǔ)分析-3：鑒定細(xì)胞類(lèi)型 - 簡(jiǎn)書(shū)
scRNA基礎(chǔ)分析-4：細(xì)胞亞類(lèi)再聚類(lèi)卵迂、注釋 - 簡(jiǎn)書(shū)
scRNA基礎(chǔ)分析-5：偽時(shí)間分析 - 簡(jiǎn)書(shū)
scRNA基礎(chǔ)分析-6：富集分析 - 簡(jiǎn)書(shū)
之前已將1222個(gè)細(xì)胞分成10個(gè)cluster，這一步我們鑒定下每個(gè)cluster為什么類(lèi)型細(xì)胞绒净。
library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
scRNA <- readRDS("scRNA.rds")
法1：mark基因鑒定
- 通過(guò)細(xì)胞類(lèi)型特異性表達(dá)的marker基因識(shí)別（常用）见咒，就好像流式細(xì)胞儀用特定的抗體篩選細(xì)胞一樣。
- 最關(guān)鍵的是要找到各個(gè)cluster的顯著高表達(dá)的基因（cluster間基因表達(dá)差異分析）挂疆，主要是利用Seurat包的
FindAllMarkers函數(shù)论颅，通過(guò)設(shè)置其test.use參數(shù)，從而選擇不同的method囱嫩。常用的有以下三種：
wilcox(default)
- 為默認(rèn)參數(shù)恃疯，也是最建議使用的。
diff.wilcox = FindAllMarkers(scRNA)
#大概1-2min
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
#先把基因這一列放在第一列墨闲，然后選取p值小于0.05的行(結(jié)果行數(shù)不變今妄，說(shuō)明都挺好的)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
#0-9共10個(gè)cluster，每個(gè)cluster選取top10高變基因
1-1
MAST
- MASK是專(zhuān)門(mén)針對(duì)單細(xì)胞數(shù)據(jù)差異分析設(shè)計(jì)的
diff.mast = FindAllMarkers(scRNA, test.use = 'MAST')
# 較慢鸳碧，4-5min左右
all.markers = diff.mast %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DESeq2
- 為bulk RNA差異分析的經(jīng)典方法
diff.deseq2 = FindAllMarkers(scRNA, test.use = 'DESeq2', slot = 'counts')
#選用 "DESeq2"盾鳞，還需要將slot參數(shù)設(shè)為"counts"
# 太慢，放棄了
all.markers = diff.deseq2 %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
top10基因熱圖可視化
top10 = CaseMatch(search = as.vector(top10$gene), match = rownames(scRNA))
#這里選取的是wilcox方法挑選的差異基因
plot1 = DoHeatmap(scRNA, features = top10, group.by = "seurat_clusters", group.bar = T, size = 4)
#結(jié)果可以看出cluster間差異還是挺明顯的
1-2
挑選部分感興趣基因可視化
select_genes <- c('LYZ','CD79A','CD8A','CD8B','GZMB','FCGR3A')
#t挑選基因瞻离，最好分布在幾個(gè)cluster里
VlnPlot(scRNA, features = select_genes, pt.size=0, ncol=2)
# 圖1-3
FeaturePlot(scRNA, features = select_genes, reduction = "tsne", label=T, ncol=2)
# 圖1-4
1-3
1-4
人工注釋
在上述工作的基礎(chǔ)上腾仅，就可以進(jìn)行人工注釋cluster細(xì)胞類(lèi)型了。
- 首先要明確自己研究組織可能有哪些細(xì)胞類(lèi)型套利，各種細(xì)胞的marker基因是什么推励；
其中細(xì)胞的marker基因主要通過(guò)相關(guān)領(lǐng)域的文獻(xiàn)收集，也可以通過(guò)專(zhuān)門(mén)的數(shù)據(jù)庫(kù)查找肉迫，作者推薦了兩個(gè)比較常用的數(shù)據(jù)庫(kù)：
CellMarker：http://biocc.hrbmu.edu.cn/CellMarker/index.jsp
PanglaoDB：https://panglaodb.se/index.html -
然后對(duì)比差異分析得到的各個(gè)cluster的顯著高表達(dá)基因验辞，綜合分析就可以判斷細(xì)胞類(lèi)型了。
這里具體就不演示了喊衫，直接展示下作者的鑒定結(jié)果跌造。
1-5
法2：SingleR包注釋
library(SingleR)
refdata <- HumanPrimaryCellAtlasData()
#參考數(shù)據(jù)庫(kù)，等待時(shí)間較長(zhǎng)族购。建議下載成功后壳贪，儲(chǔ)存為Rdata，以后方便使用寝杖。
testdata <- GetAssayData(scRNA, slot="data")
clusters <- scRNA@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main,
# 注意此時(shí)labels參數(shù)為 refdata$label.main违施，與下一節(jié)亞類(lèi)再注釋時(shí)的設(shè)置不同
method = "cluster", clusters = clusters,
assay.type.test = "logcounts", assay.type.ref = "logcounts")
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = F)
#如下為singleR的細(xì)胞cluster鑒定結(jié)果。
2-1
#結(jié)合上述結(jié)果朝墩，給scRNA增添celltype注釋信息
scRNA@meta.data$celltype = "NA"
#先新增列celltype醉拓，值均為NA伟姐，然后利用下一行代碼循環(huán)填充
for(i in 1:nrow(celltype)){
scRNA@meta.data[which(scRNA@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
可視化1
p1 = DimPlot(scRNA, group.by="celltype", label=T, label.size=5, reduction='tsne')
p2 = DimPlot(scRNA, group.by="celltype", label=T, label.size=5, reduction='umap')
p3 = plotc <- p1+p2+ plot_layout(guides = 'collect')
2-2
可視化2
select_genes <- c('LYZ','CD79A','GZMB','FCGR3A')
#vlnplot展示
p1 <- VlnPlot(scRNA, features = select_genes, pt.size=0, group.by="celltype", ncol=2)
#featureplot展示
p2 <- FeaturePlot(scRNA, features = select_genes, reduction = "tsne", label=T, ncol=2)
p3=p1|p2
2-3
#最后保存分析結(jié)果收苏，以供后續(xù)的分析
saveRDS(scRNA, file="scRNA.rds")
