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.framecolMeans(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 + plot2

5. 归一化表达值

### 归一化表达值
#对基因在细胞间的表达量进行线性变换，使平均值为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

Seurat包分析单细胞转录组数据代码相关推荐

monocle3包分析单细胞转录组数据
1. 构建new_cell_data_set对象 Usage new_cell_data_set(expression_data, cell_metadata = NULL, gene_metadat ...
scater分析单细胞转录组数据代码
1. 演示数据和包的加载 # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages(&q ...
代码分析 | 单细胞转录组数据整合详解
两种整合方法详解 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅 ...
单细胞转录组数据整合分析专题研讨会(2019.11)
2019年10月9日,单细胞转录组再等Nature.题为Decoding human fetal liver haematopoiesis的研究,对受孕后4周至17周的人胚胎肝脏.卵黄囊.肾脏和皮肤组 ...
代码分析 | 单细胞转录组质控详解
前言 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅综述:三万字 ...
Nat. Commun. | 从单细胞转录组数据中学习可解释的细胞和基因签名嵌入
本文介绍由加拿大麦吉尔大学与蒙特利尔高等商学院.北京大学.复旦大学的研究人员联合发表在Nature Communications的研究成果:本文作者提出了单细胞嵌入式主题模型scETM(single- ...
SCS【3】单细胞转录组数据 GEO 下载及读取
点击关注,桓峰基因今天来介绍一下GEO单细胞转录组下载数据以及整理,单细胞测序的原理以及数据结果都与bulk测序的方式有一定的差距,所以我们单独说一下. 桓峰基因的教程不但教您怎么使用,还会定期分析 ...
scDML：单细胞转录组数据的批次对齐
scRNA-seq支持在单细胞分辨率呈现基因表达谱,这提高了对已知细胞类型的检测,以及异质细胞与疾病失调的理解.然而,scRNA-seq技术的广泛应用产生了许多庞大而复杂的数据集,这给集成来自不同批次 ...
使用ArchR分析单细胞ATAC-seq数据(第十四章)
本文首发于我的个人博客, http://xuzhougeng.top/ 往期回顾: 使用ArchR分析单细胞ATAC-seq数据(第一章) 使用ArchR分析单细胞ATAC-seq数据(第二章) 使用 ...

Seurat包分析单细胞转录组数据代码

1. 构建Seurat对象

2. 预处理与质量控制

3. 校正表达值

4. 鉴定高变异基因

5. 归一化表达值

6. PCA降维

7. UMAP/tSNE降维

8. 差异表达分析

9. 鉴定细胞类型

参考

Seurat包分析单细胞转录组数据代码相关推荐

最新文章

热门文章