單細(xì)胞實(shí)戰(zhàn)2: Seurat+SingleR--從矩陣到細(xì)胞類型注釋

單細(xì)胞數(shù)據(jù)標(biāo)準(zhǔn)分析Seurat流程：

Seurat官網(wǎng)：https://satijalab.org/seurat/articles/get_started.html（官網(wǎng)教程極為詳細(xì)，建議仔細(xì)反復(fù)閱讀）

  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)) #查看每個(gè)樣多少個(gè)細(xì)胞
#  con    FM 
#21963 17682 

#使用PercentageFeatureSet函數(shù)計(jì)算線粒體基因的百分比
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() #不帶點(diǎn)
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() #帶點(diǎn)

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()  ##不帶點(diǎn)且y軸取log（以這個(gè)圖確定cutoff值比上圖更好）

#測序深度的相關(guān)性繪圖
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()   #保存基因特征相關(guān)性圖
#只畫一個(gè)樣的細(xì)胞
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()
  1. 過濾細(xì)胞
selected <- WhichCells(scRNA1, expression = nFeature_RNA > 400  & nFeature_RNA < 3000 & percent.mt < 15 & nCount_RNA > 1500 & nCount_RNA < 10000) 
>selected #查看符合標(biāo)準(zhǔn)的細(xì)胞有多少
#過濾不合條件的細(xì)胞
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)     #(帶點(diǎn)且y軸取log)
VlnPlot(object = scRNA, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,log = T)
dev.off()
#在過濾前后檢查每個(gè)樣本的細(xì)胞數(shù)
table(Idents(scRNA))
#  con    FM 
#18704 13420
  1. 對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化
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)
  1. PCA
scRNA=RunPCA(object= scRNA,npcs = 50,pc.genes=VariableFeatures(object = scRNA))     

#PCA可視化1: 一維可視化，對每個(gè)維度中，權(quán)重較大的基因進(jìn)行可視化
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:主成分分析熱圖 確定選多少個(gè)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 確定選多少個(gè)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 確定選多少個(gè)pc
ElbowPlot(scRNA)
pdf(file="06.elbowplot.pdf",width=10,height=8)
ElbowPlot(scRNA)
dev.off()

07.細(xì)胞聚類

pcSelect=20
scRNA <- FindNeighbors(object = scRNA, dims = 1:pcSelect,k.param = 20)    #計(jì)算鄰接距離（細(xì)胞間鄰接節(jié)點(diǎn)的距離）

#看不同resolution對分群的影響 
scRNA <- FindClusters(object = scRNA, resolution = c(seq(0,1.6,.2))) #計(jì)算resolution為從0到1.6間隔為0.2的各個(gè)值時(shí)的分群情況
#以上操作結(jié)果儲(chǔ)存在scRNA@metadata,每個(gè)分辨率都有單獨(dú)一列，使用active.ident也可以查看，active.ident的level就是各個(gè)分群
library(clustree)
pdf(file="06.clustree.pdf",width=10,height=14)  #可視化不同resolution對分群的影響
clustree(scRNA@meta.data, prefix = "RNA_snn_res.")
dev.off()
#對非線性降維的結(jié)果可視化可以通過idents()函數(shù)來指定分辨率
Idents(object = scRNA) <- "RNA_snn_res.0.6"
levels(scRNA)

#也可直接指定resolution
#scRNA <- FindClusters(object = scRNA, resolution = 0.5)
  1. 降維后用聚類的類別可視化（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()
  1. 尋找差異表達(dá)的特征
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在各個(gè)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在各個(gè)cluster的散點(diǎn)圖
pdf(file="02.markerScatter.pdf",width=10,height=6)
FeaturePlot(object = scRNA, features = c("Cd83", "Cxcl2"))
dev.off()

#繪制marker在各個(gè)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()
  1. SingleR注釋細(xì)胞類型
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))
#鑒定各聚類的細(xì)胞類型，可選擇method 參數(shù)為“single”，返回每個(gè)細(xì)胞的鑒定結(jié)果。clusters參數(shù)指定的各聚類的名稱
pred.scRNA$pruned.labels

#查看注釋準(zhǔn)確性 
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()
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容