參考：
細(xì)胞類型注釋
經(jīng)典marker手動注釋
#在文獻(xiàn)里面經(jīng)常用到的經(jīng)典marker
在文獻(xiàn)里面經(jīng)常出現(xiàn)
celltype_marker=c(
"EPCAM",#上皮細(xì)胞 epithelial
"PECAM1",#內(nèi)皮細(xì)胞 endothelial
"COL3A1",#成纖維細(xì)胞 fibroblasts
"CD163","AIF1",#髓系細(xì)胞 myeloid
"CD79A",#B細(xì)胞
"JCHAIN",#漿細(xì)胞 plasma cell
"CD3D","CD8A","CD4",#T細(xì)胞
"GNLY","NKG7",#NK細(xì)胞
"PTPRC"#免疫細(xì)胞
)
結(jié)合對應(yīng)的關(guān)系赴涵，初步確認(rèn)細(xì)胞類型如下:
VlnPlot(test.seu,features = celltype_marker,pt.size = 0,ncol = 2)
ggsave(filename = "marker.png",device = "png",width = 44,height = 33,units = "cm")
image.png
NK細(xì)胞和T細(xì)胞不是很能區(qū)分開拖叙，一些文獻(xiàn)也是直接把這兩個當(dāng)做一個大群來做的横媚，此處我根據(jù)GNLY基因確定第5個cluster為NK細(xì)胞继效，是否正確我們后面再看早像。
免疫細(xì)胞
0: B_cell
4: Plasma_cell
1 2: T_cell
5: NK_cell
13: Unknown
非免疫細(xì)胞
3 6 7 8 10 11 12: Epithelial
14: Endothelial
9: Fibroblasts
doublet
15: Doublet (Myeloid+CD4)
以上根據(jù)經(jīng)典的marker初步確定了細(xì)胞大類桨踪，接著我們找差異基因，看看找出來的每一個cluster的差異基因是不是和前面鑒定的類型一致掀亩。
- 找差異基因
markers <- FindAllMarkers(test.seu, logfc.threshold = 0.25, min.pct = 0.1,
only.pos = TRUE, test.use = "wilcox")
markers_df = markers %>% group_by(cluster) %>% top_n(n = 500, wt = avg_logFC)
#根據(jù)avg_logFC這一列把前500個差異基因取出來舔哪，小于500個時，有多少保留多少槽棍。
#不必在意我們保留多少個差異基因尸红，實(shí)際用到的，也就前面十幾（幾十）個基因刹泄。
write.table(markers_df,file="test.seu_0.5_logfc0.25_markers500.txt",quote=F,sep="\t",row.names=F,col.names=T)
- FindAllMarkers()這個函數(shù)會將cluster中的某一類和剩下的那些類比較，來找差異基因怎爵。與其對應(yīng)的有個FindMarkers()函數(shù)特石，可以指定哪兩個cluster對比。
- logfc.threshold表示logfc的閾值鳖链，這里有兩個地方需要注意：一是Seurat里面的logfc計算公式很特別姆蘸，并不是我們平常在bulk里面那樣算均值，相除芙委，求log逞敷；二是如果想畫火山圖，這個閾值可以設(shè)為0灌侣，不然最后畫出來的火山圖中間會缺一段推捐，我們完全可以把全部基因先拿出來，畫火山圖侧啼，再根據(jù)p value, logFC這些閾值自己過濾牛柒。
- min.pct表示基因在多少細(xì)胞中表達(dá)的閾值
- only.pos = TRUE表示只求高表達(dá)的基因
- test.use表示檢測差異基因所用的方法。
SingleR注釋
細(xì)胞亞群注釋
library(SingleR)
library(celldex)
refdata <- MouseRNAseqData()
refdata$label.ont
testdata <- GetAssayData(scRNA, slot="data")
clusters <- scRNA@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.fine,
# label.finea耗時比較長一點(diǎn)
method = "cluster", clusters = clusters,
assay.type.test = "logcounts", assay.type.ref = "logcounts")
rm(refdata, testdata) #珍惜內(nèi)存
table(cellpred$labels)
給scRNA增添celltype注釋信息
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = F)
table(celltype$ClusterID,celltype$celltype)
為singleR的細(xì)胞cluster鑒定結(jié)果痊乾。
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]}
DimPlot(scRNA, group.by="celltype", label=F , reduction='tsne')
