单细胞实战2: Seurat+SingleR--从矩阵到细胞类型

单细胞数据标准分析Seurat流程：

Seurat官网：https://satijalab.org/seurat/articles/get_started.html（官网教程极为详细，建议仔细反复阅读）

1. 建立Seurat对象

scRNA1 <- CreateSeuratObject(counts = all.matrix)
save.image(file = "scRNA1.RData") #保存

2.cutoff值的确定

library(limma)
library(Seurat)
library(dplyr)
library(magrittr)

table(Idents(scRNA1)) #查看每个样多少个细胞
#  con    FM 
#21963 17682 

#使用PercentageFeatureSet函数计算线粒体基因的百分比
head(scRNA1@meta.data)
scRNA1[["percent.mt"]] <- PercentageFeatureSet(object = scRNA1, pattern = "^mt-")
#保存基因特征小提琴图
pdf(file="02.featureViolin.pdf",width=10,height=6)    
VlnPlot(object = scRNA1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0)
dev.off() #不带点
pdf(file="02.featureViolin_withPoint.pdf",width=10,height=6)           
VlnPlot(object = scRNA1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off() #带点

pdf(file="02.featureViolin_log.pdf",width=10,height=6)     
VlnPlot(object = scRNA1, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0,log = T)
dev.off()  ##不带点且y轴取log（以这个图确定cutoff值比上图更好）

#测序深度的相关性绘图
pdf(file="02.featureCor.pdf",width=10,height=6)         
plot1 <- FeatureScatter(object = scRNA1, feature1 = "nCount_RNA", feature2 = "percent.mt",pt.size=1.5)
plot2 <- FeatureScatter(object = scRNA1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",,pt.size=1.5)
CombinePlots(plots = list(plot1, plot2)) 
dev.off()   #保存基因特征相关性图
#只画一个样的细胞
FeatureScatter(scRNA1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", cells = WhichCells(scRNA1, expression = orig.ident == "con"))
FeatureScatter(scRNA1, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", cells = WhichCells(scRNA1, expression = orig.ident == "FM"))

#画频率分布条形图
pdf(file="02.featurebar.pdf",width=10,height=6)     
hist(scRNA1$nCount_RNA,breaks = 100,xlim = range(0:10000))+abline(v=(1000+100*c(0:10)),col=rainbow(11),lwd=1)
dev.off()
pdf(file="02.mtbar.pdf",width=10,height=6)
hist(scRNA1$percent.mt,xlim = range(0:40),freq=T)
dev.off()

3. 过滤细胞

selected <- WhichCells(scRNA1, expression = nFeature_RNA > 400  & nFeature_RNA < 3000 & percent.mt < 15 & nCount_RNA > 1500 & nCount_RNA < 10000) 
>selected #查看符合标准的细胞有多少
#过滤不合条件的细胞
scRNA <- subset(x = scRNA1, subset = nFeature_RNA > 400  & nFeature_RNA < 3000 & percent.mt < 15 & nCount_RNA > 1500 & nCount_RNA < 10000)

#过滤后的feature图
pdf(file="03.featureViolin_withPoint.pdf",width=10,height=6)         
VlnPlot(object = scRNA, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()
pdf(file="03.featureViolin_log.pdf",width=10,height=6)     #(带点且y轴取log)
VlnPlot(object = scRNA, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,log = T)
dev.off()
#在过滤前后检查每个样本的细胞数
table(Idents(scRNA))
#  con    FM 
#18704 13420 

4. 对数据进行标准化

scRNA <- NormalizeData(object = scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
scRNA <- FindVariableFeatures(object = scRNA, selection.method = "vst", nfeatures = 2000)

#输出特征方差图
top10 <- head(x = VariableFeatures(object = scRNA), 10)
pdf(file="04.featureVar.pdf",width=10,height=6)              #保存基因特征方差图
plot1 <- VariableFeaturePlot(object = scRNA)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2))
dev.off()

05.归一化

all.genes <- rownames(scRNA)
scRNA <- ScaleData(scRNA, features = all.genes)                   

6. PCA

scRNA=RunPCA(object= scRNA,npcs = 50,pc.genes=VariableFeatures(object = scRNA))     

#PCA可视化1: 一维可视化，对每个维度中，权重较大的基因进行可视化
pdf(file="06.pcaGene.pdf",width=10,height=8)
VizDimLoadings(object = scRNA, dims = 1:4, reduction = "pca",nfeatures = 20)
dev.off()

#PCA可视化2: 二维可视化
pdf(file="06.PCA.pdf",width=6.5,height=6)
DimPlot(object = scRNA, reduction = "pca")
dev.off()

#PCA可视化3:主成分分析热图 确定选多少个pc
pdf(file="06.pcaHeatmap.pdf",width=20,height=20)
DimHeatmap(object = scRNA, dims = 1:15, cells = 500, balanced = TRUE,nfeatures = 30,ncol=3)
dev.off()

#PCA可视化4: JackStrawPlot 确定选多少个pc
scRNA <- JackStraw(object = scRNA, num.replicate = 100)
scRNA <- ScoreJackStraw(object = scRNA, dims = 1:20)
pdf(file="06.pcaJackStraw.pdf",width=8,height=6)
JackStrawPlot(object = scRNA, dims = 1:20)
dev.off()

#PCA可视化5: ElbowPlot 确定选多少个pc
ElbowPlot(scRNA)
pdf(file="06.elbowplot.pdf",width=10,height=8)
ElbowPlot(scRNA)
dev.off()

07.细胞聚类

pcSelect=20
scRNA <- FindNeighbors(object = scRNA, dims = 1:pcSelect,k.param = 20)    #计算邻接距离（细胞间邻接节点的距离）

#看不同resolution对分群的影响 
scRNA <- FindClusters(object = scRNA, resolution = c(seq(0,1.6,.2))) #计算resolution为从0到1.6间隔为0.2的各个值时的分群情况
#以上操作结果储存在scRNA@metadata,每个分辨率都有单独一列，使用active.ident也可以查看，active.ident的level就是各个分群
library(clustree)
pdf(file="06.clustree.pdf",width=10,height=14)  #可视化不同resolution对分群的影响
clustree(scRNA@meta.data, prefix = "RNA_snn_res.")
dev.off()
#对非线性降维的结果可视化可以通过idents()函数来指定分辨率
Idents(object = scRNA) <- "RNA_snn_res.0.6"
levels(scRNA)

#也可直接指定resolution
#scRNA <- FindClusters(object = scRNA, resolution = 0.5)   

8. 降维后用聚类的类别可视化（tSNE/uMAP)

#TSNE聚类
scRNA <- RunTSNE(object = scRNA, dims = 1:pcSelect)                      
pdf(file="08.TSNE.pdf",width=6.5,height=6)
TSNEPlot(object = scRNA, pt.size = 0.5, label = TRUE)   
dev.off()
write.table(scRNA$seurat_clusters,file="08.tsneCluster.txt",quote=F,sep="\t",col.names=F)

#umap聚类 
scRNA <- RunUMAP(scRNA, dims = 1:pcSelect)
pdf(file="08.uMAP.pdf",width=6.5,height=6)
DimPlot(scRNA, reduction = "umap",pt.size=0.5,label = TRUE)
dev.off()

9. 寻找差异表达的特征

logFCfilter=0.5
adjPvalFilter=0.05
scRNA.markers <- FindAllMarkers(object = scRNA,
                               only.pos = FALSE,
                               min.pct = 0.25,
                               logfc.threshold = logFCfilter)
sig.markers=scRNA.markers[(abs(as.numeric(as.vector(scRNA.markers$avg_log2FC)))>logFCfilter & as.numeric(as.vector(scRNA.markers$p_val_adj))<adjPvalFilter),]
write.table(sig.markers,file="09.markers.xls",sep="\t",row.names=F,quote=F)
top10 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

#绘制marker在各个cluster的热图
pdf(file="09.tsneHeatmap.pdf",width=48,height=36)
DoHeatmap(object = scRNA, features = top10$gene) + NoLegend()
dev.off()  #绘制速度较慢

#绘制marker的小提琴图
pdf(file="02.markerViolin.pdf",width=10,height=6)
VlnPlot(object = scRNA, features = c("Cd83", "Cxcl2"))
dev.off()

#绘制marker在各个cluster的散点图
pdf(file="02.markerScatter.pdf",width=10,height=6)
FeaturePlot(object = scRNA, features = c("Cd83", "Cxcl2"))
dev.off()

#绘制marker在各个cluster的气泡图
pdf(file="06.markerBubble.pdf",width=12,height=6)
cluster0Marker=c("Cbr2", "Lyve1", "Selenop", "Folr2", "Ednrb", "F13a1", "Mrc1", "Igf1", "Slc40a1
", "Cd163")
DotPlot(object = scRNA, features = cluster0Marker)
dev.off()

10. SingleR注释细胞类型

library(celldex)
ref <- ImmGenData()
library(SingleR)
library(BiocParallel)

pred.scRNA <- SingleR(test = scRNA@assays$RNA@data, ref = ref,labels = ref$label.main, clusters = scRNA@active.ident, fine.tune = TRUE, BPPARAM = MulticoreParam(40))
#鉴定各聚类的细胞类型，可选择method 参数为“single”，返回每个细胞的鉴定结果。clusters参数指定的各聚类的名称
pred.scRNA$pruned.labels

#查看注释准确性 
pdf(file="10.celllabel.pdf",width=10,height=10)   
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)
dev.off()

#绘制带cell label的tsne和umap图
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
levels(scRNA)
# 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
scRNA <- RenameIdents(scRNA,new.cluster.ids)
levels(scRNA)

pdf(file="10.TSNE_lable.pdf",width=6.5,height=6)
TSNEPlot(object = scRNA1, pt.size = 0.5, label = TRUE)  
dev.off()

pdf(file="10.uMAP_label.pdf",width=6.5,height=6)
DimPlot(scRNA1, reduction = "umap",pt.size=0.5,label = TRUE)
dev.off()