Seurat包分析单细胞转录组数据代码
作者：互联网
1. 构建Seurat对象
# install.packages('Seurat')
rm(list=ls())
# 表达矩阵来自CellRanger分析结果
data_dir <- "/run_count_1kpbmcs/outs/filtered_feature_bc_matrix"
setwd("project_path")
library(Seurat)
library(dplyr)
### 构建Seurat对象
expr_dgCMatrix <- Read10X(data.dir = data_dir)
#class(expr_dgCMatrix)
seurat_object <- CreateSeuratObject(counts = expr_dgCMatrix)
# counts:Either a matrix-like object with unnormalized data
# with cells as columns and features as rows
# or an Assay-derived object
seurat_object
# sparse Matrix of class "dgCMatrix" 行:基因，列:细胞barcode
seurat_object@assays$RNA@counts
seurat_object@assays$RNA@data[1:3,1:2]
seurat_object@meta.data # data.frame
colMeans(seurat_object)
rowMeans(seurat_object)
2. 预处理与质量控制
### 预处理与质量控制
#计算线粒体RNA比例
seurat_object[["percent.mt"]] <- PercentageFeatureSet(seurat_object,
pattern = "^MT-")
head(seurat_object@meta.data, 5)
#使用violin plot查看QC指标分布，如图1所示
VlnPlot(seurat_object, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3)
#作图显示feature之间的皮尔逊相关性
plot1 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
#筛选数据根据之前violin plot来确定。
seurat_object <- subset(seurat_object, subset = nFeature_RNA > 500 &
nFeature_RNA < 6000 & percent.mt < 25)
dim(seurat_object)
saveRDS(seurat_object, file = "subset.rds")
#readRDS("subset.rds")3. 校正表达值
### 校正表达值
#校正细胞的表达值，将count数量根据表达总量矫正之后，乘以10000，并进行对数转换
seurat_object_norm <- NormalizeData(seurat_object)
# # 查看矫正前后的数据
# head(colMeans(seurat_object))
# head(colMeans(seurat_object_norm))4. 鉴定高变异基因
### 鉴定高变异基因
seurat_object_norm <- FindVariableFeatures(seurat_object_norm, nfeatures = 2000)
#查看基因的平均表达与标准方差之间的关系，如图3所示
plot1 <- VariableFeaturePlot(seurat_object_norm)
top10_genes <- head(VariableFeatures(seurat_object_norm), 10)
# repel参数用了后作图出错？
# plot2 <- LabelPoints(plot = plot1, points = top10_genes,
# repel = TRUE,xnudge = 0,ynudge = 0)
plot2 <- LabelPoints(plot = plot1, points = top10_genes, repel = FALSE)
plot1 + plot25. 归一化表达值
### 归一化表达值
#对基因在细胞间的表达量进行线性变换，使平均值为0，方差为1，去除高表达基因的影响
all.genes <- rownames(seurat_object_norm)
seurat_object_scaled <- ScaleData(seurat_object_norm, features = all.genes)
#以上代码的计算原理如下:
#pbmc@assays$RNA@scale.data <- t(scale(t(pbmc@assays$RNA@data)))
#scale是个十分耗时的过程，为缩短计算时间，ScaleData默认仅对高变异基因进行scale；
#不论是对全部基因或仅对高变异基因进行scale，并不影响下游的PCA和聚类过程
#ScaleData同样可以用来去除其他因素（如细胞周期状态或线粒体RNA比例）带来的影响
seurat_object_scaled <- ScaleData(seurat_object_scaled, vars.to.regress = "percent.mt")
dim(seurat_object_scaled)6. PCA降维
### 线性降维
#在scale.data上运行PCA，默认使用的基因为高变异基因
seurat_object_scaled <- RunPCA(seurat_object_scaled,
features = VariableFeatures
(object = seurat_object_scaled))
#观察PC上重要基因
print(seurat_object_scaled[["pca"]], dims = 1:5, nfeatures = 5)
#在PC1与PC2上表示基因，如图4所示，这里展示的实际上是细胞进行PCA降维时每个基因的系数，
#即PCA feature loading，根据PCA算法，有pbmc@reductions$pca@cell.embeddings ＝
#t(pbmc@assays$RNA@scale.data[VariableFeatures(pbmc),])
#%*% pbmc@reductions$pca@feature.loadings
VizDimLoadings(seurat_object_scaled, dims = 1:2, reduction = "pca")
#graph2ppt(file="4.PC1nPC2_value.pptx", width=9, aspectr=1.5)
#细胞的PCA降维作图
DimPlot(seurat_object_scaled, reduction = "pca")
#graph2ppt(file="5.PCA降维.pptx", width=9, aspectr=1.5)
#查看每个细胞中单个基因在不同PC维度上的得分
DimHeatmap(seurat_object_scaled, dims = 1, cells = 500, balanced = TRUE)
#graph2ppt(file="6.单个PC维度_heatmap.pptx", width=9, aspectr=1.5)
DimHeatmap(seurat_object_scaled, dims = 1:9, cells = 500, balanced = TRUE)
#graph2ppt(file="7-1.1到9个PC维度_heatmap.pptx", width=12, aspectr=1.5)
DimHeatmap(seurat_object_scaled, dims = 10:18, cells = 500, balanced = TRUE)
#graph2ppt(file="7-2.9到18个PC维度_heatmap.pptx", width=12, aspectr=1.5)
DimHeatmap(seurat_object_scaled, dims = 19:27, cells = 500, balanced = TRUE)
#graph2ppt(file="7-3.18到27个PC维度_heatmap.pptx", width=12, aspectr=1.5)
###保存数据
saveRDS(seurat_object_scaled, file = "PCA.rds")
pbmc <- seurat_object_scaled
### 确定PCA维度数量
#根据JackStrawPlot来确定下一步降维作用的PC数量，如图７所示。P值越小说明PC越重要，应该纳入下一步分析
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)#这里dims = 1:20表示看前20个PC的P值，因为第一次算出来的结果真的惊人。而且这个dim最大只能是20.
JackStrawPlot(pbmc, dims = 1:20) #用JackStrawPlot绘图，绘制前10个PC的P值分布曲线
#graph2ppt(file="8.PCA维度_JackStraw.pptx", width=12, aspectr=1.5)
#由图可以看出，在PC10-PC12之后，PC的P值出现显著降低
#或根据ElbowPlot确定PC数量，如图所示。标准差较大的PC所解释的数据变异性更大，应当纳入下一步分析
ElbowPlot(pbmc)
#graph2ppt(file="9.PCA维度_ElbowPlot.pptx", width=12, aspectr=1.5)
#由图可看出，在PC9-PC10后出现标准差数据拐点。综上，下游分析中选择了前10个PC
### 细胞聚类
#选取前10个主成分来分类细胞。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#这里的resolution参数可以修改细胞类型的多少。resolution越大细胞的类型越多，
#所以可以尝试。0.01-1群细胞，0.05-2群细胞，0.1-4群细胞，0.5-8群细胞。
#pbmc <- FindClusters(pbmc, resolution = 0.05)
##或者上述的resolution还可以用另外一种方式将其划分，便于找到最优的分群数目。
pbmc <- FindClusters(pbmc, resolution = c(0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1))
head(Idents(pbmc), 5)
##查看上述不同的resolution的时候，聚类的情况。
#install.packages('magrittr')
library(magrittr)
#pbmc@meta.data %>% View()
#install.packages('clustree')
library(clustree)
#Creates a plot of a clustering tree showing the relationship
#between clusterings at different resolutions.
clustree(pbmc@meta.data, prefix = "RNA_snn_res.")
#graph2ppt(file="10.clustree.pptx", width=9, aspectr=1.5)
#根据上述的结果，找到最佳的分群resolution的值为0.3。
pbmc <- FindClusters(pbmc, resolution = 0.3)7. UMAP/tSNE降维
### 非线性降维（UMAP/tSNE）
#UMAP降维
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap")
#graph2ppt(file="11.umap.pptx", width=9, aspectr=1.5)
#tSNE降维,利用pc1-pc20。pbmc[["pca"]] Number of dimensions: 50
pbmc <- RunTSNE(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "tsne")
#graph2ppt(file="12.tsne.pptx", width=9, aspectr=1.5)
###这里保存数据，用于鉴定singlelets和doublelets
saveRDS(pbmc, file = "ForDoublets.rds")8. 差异表达分析
### 差异表达分析
# Idents(pbmc)
# find all markers of cluster 1# 找到cluster1当中的所有marker
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
cluster3.markers <- FindMarkers(pbmc, ident.1 = 3, min.pct = 0.25)
head(cluster3.markers, n = 5)
#使用不同的假设检验方法（test.use = "roc"）来进行差异表达分析
# 默认 test.use = "wilcox"
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25,
test.use = "roc", only.pos = TRUE)
# find all markers distinguishing cluster 5 from clusters 0 and 3#找到cluster5当中和clusters0和clusters3能够区分的marker。
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# 找到每一个cluster当中的marker，并且只展示阳性的marker。
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE,
min.pct = 0.25, logfc.threshold = 0.25)
##上述步骤的原理如下：
#pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
head(pbmc.markers, n = 10)
##将这个数据保存一下，用于亚群鉴定得时候用。
write.table(pbmc.markers,file = 'pbmc.markers.txt',sep = '\t')
#差异基因可视化，此外还可以通过RidgePlot, CellScatter, DotPlot等进行展示，这里可以每个亚群筛选一个，也可以根据需要。
##选择每个cluster的前三个基因。此外，这里的这个名字识别的是最后一列gene的名字。
VlnPlot(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D"))
#graph2ppt(file="13.clusterMarker_VlnPlot.pptx", width=20, aspectr=1.5)
# 使用原始count绘制
VlnPlot(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D"),
slot = "counts", log = TRUE)
#graph2ppt(file="14.clusterMarker_VlnPlot_count.pptx", width=20, aspectr=1.5)
#差异基因在细胞分区图上的展示。
FeaturePlot(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D"))
#graph2ppt(file="15.clusterMarker_FeaturePlot.pptx", width=15, aspectr=1.5)
#差异基因可视化RidgePlot
RidgePlot(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D"))
#graph2ppt(file="16.clusterMarker_RidgePlot.pptx", width=18, aspectr=1.5)
#差异基因可视化CellScatter,参数还没搞懂。
#CellScatter(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D")
#差异基因可视化DotPlot
DotPlot(pbmc, features = c("IL7R", "IL32", "TRAC", "CD3D"))
#graph2ppt(file="17.clusterMarker_DotPlot.pptx", width=12, aspectr=1.5)
#热图显示。每个cluster按照avg_log2FC筛选排名前10的，进行热图展示。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
#graph2ppt(file="18.markersHeatmap.pptx", width=9, aspectr=1.5)
saveRDS(pbmc, file = "ClustersPlot_beyond.rds")9. 鉴定细胞类型
### 鉴定细胞类型
#根据已知细胞类型的marker gene对各个cluster的细胞进行命名
new.cluster.ids <- paste0("group", 1:11)
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
#graph2ppt(file="鉴定细胞类型.pptx", width=9, aspectr=1.5)
### 保存结果
saveRDS(pbmc, file = "final.rds")
参考
https://zhuanlan.zhihu.com/p/90575225?from=singlemessage
https://www.singlecellcourse.org/single-cell-rna-seq-analysis-using-seurat.html
标签：pptx,Seurat,pbmc,单细胞,转录,file,seurat,###,aspectr 来源： https://blog.csdn.net/qq_27390023/article/details/121440469