Seurat包基本分析实战—文献图表复现
时间:2022-07-27
本文章向大家介绍Seurat包基本分析实战—文献图表复现,主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。
分享是一种态度
一、背景知识
- 文献:https://www.aging-us.com/article/103695/text
- GSE:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465
- 目的是复现文献的第一幅总图
文献
二、大致思路
(1)根据GEO号下载表达矩阵,以及meta信息;
(2)根据meta信息筛选tumor cell 的表达矩阵;
(3)根据表达矩阵构建Seurat对象,以及添加分组信息、质控;
(4)归一化→降维→聚类;
(5)组图
三、具体实现
1、下载数据
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465
1-1
1-2
1-3
2、导入数据
a <- read.table("GSE84465_GBM_All_data.csv.gz")
a[1:4,1:4]
dim(a)
#[1] 23465 3589
#原始矩阵有3598个cell,23465个gene
b <- read.table("SraRunTable.txt",sep = ",",header = T)
dim(b)
b[1:4,1:4]
#[1] 3589 33
- 通过观察,表达矩阵a:3589个cell(name感觉是不同信息的组合),23465个gene(symbol ID)
- 通过观察meta b了解到cell name是主要根据每个细胞的plate_id与well组合而成;并且有tissue、patient_id两种分类方式在后面分析会用到。
3、筛选目标类型cell
new.b <- data.frame(names=paste0("X",b$plate_id,".",b$Well),
#names即为a的colname
tissue=b$Tissue, #tumor二分类
group=b$Patient_ID) #样本分类
new.b <- new.b[match(colnames(a),new.b$names),]
#使顺序一致,以便后期的筛选
table(new.b$names==colnames(a))
rownames(new.b) <- colnames(a)
index <- new.b$tissue=="Tumor"
a_filt <- a[,index]
dim(a_filt)
#挑选到2343个tumor cell 表达矩阵
b$Patient_ID有四类,说明这些细胞是从四个病人的病灶组织中取的。
4、构建Seurat对象,以及质控、可视化
group=data.frame(group=new.b[index,3],
row.names = colnames(a_filt))
library("Seurat")
scRNA = CreateSeuratObject(counts=a_filt,
meta.data = group,
min.cells = 3, min.features = 1500)
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
#标记线粒体基因
head(scRNA@meta.data)
summary(scRNA@meta.data)
#结果显示没有线粒体基因,因此感觉过滤也没有意义
pctMT=5
scRNA <- subset(scRNA, subset = percent.mt < pctMT)
- 通过结果可以看到仅过滤掉一个细胞;
- 但有趣的是我是按照文献的三个标准过滤的,而文献是过滤掉了194个cell。
- 从这里开始,分析结果就与文献产生了差异。但目前也不知道是为什么。
4-1
可视化三个图
#图A:观察不同组cell的counts、feature分布
col.num <- length(levels(scRNA@meta.data$group))
p1_1 <- VlnPlot(scRNA,
features = c("nFeature_RNA","nCount_RNA"),
group.by = "group",
cols =rainbow(col.num)) +
theme(legend.position = "none") +
labs(tag = "A")
#图B:nCount_RNA与对应的nFeature_RNA关系
library(ggplot2)
p1_2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",
group.by = "group",pt.size = 1.3) +
labs(tag = "B")
5、挑选hvg高变基因
#标准化
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
#归一化
scRNA <- ScaleData(scRNA, features = (rownames(scRNA)))
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 1500)
top10 <- head(VariableFeatures(scRNA), 10)
plot1 <- VariableFeaturePlot(scRNA)
p1_3 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) +
theme(legend.position = c(0.1,0.8)) +
labs(tag = "C")
#标记top10 hvg
library(patchwork)
p1_1 | p1_2 | p1_3
上三个子图
6、降维
#默认使用前2000个hvg降维处理
scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA))
p2_1 <- DimPlot(scRNA, reduction = "pca", group.by="group")+
labs(tag = "D")
#挑选主成分
scRNA <- JackStraw(scRNA,reduction = "pca")
scRNA <- ScoreJackStraw(scRNA,dims = 1:20)
p2_2 <- JackStrawPlot(scRNA,dims = 1:20, reduction = "pca") +
theme(legend.position="bottom") +
labs(tag = "E")
p2_3 <- ElbowPlot(scRNA, ndims=20, reduction="pca")
#结果显示可挑选前20个pc
p2_1| (p2_2 | p2_3)
中三个子图
7、聚类
pc.num=1:20
scRNA <- FindNeighbors(scRNA, dims = pc.num)
# dims参数,需要指定哪些pc轴用于分析;这里利用上面的分析,选择20
scRNA <- FindClusters(scRNA, resolution = 0.4)
table(scRNA@meta.data$seurat_clusters)
# 0 1 2 3 4 5 6 7 8 9 10 11 12
#393 380 281 271 177 140 122 108 72 51 41 35 18
scRNA = RunTSNE(scRNA, dims = pc.num)
p3_1 <- DimPlot(scRNA, reduction = "tsne") +
labs(tag = "E")
- 如上结果,可知分为13个clust;
- 寻找每个clust的marker gene,这里使用的是
FindAllMarkers()
函数
diff.wilcox = FindAllMarkers(scRNA)
#大概2-3min
head(diff.wilcox)
library(tidyverse)
all.markers = diff.wilcox %>% select(gene, everything()) %>%
subset(p_val<0.05 & abs(diff.wilcox$avg_logFC) > 0.5)
top20 = all.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
top20
#0-9共10个cluster,每个cluster选取top10高变基因
top20 = CaseMatch(search = as.vector(top20$gene), match = rownames(scRNA))
length(top20)
length(unique(sort(top20)))
#这里选取的是wilcox方法挑选的差异基因
p3_2 <- DoHeatmap(scRNA, features = top20, group.by = "seurat_clusters")
此处的结果也是与原文差别比较大的地方。均是对每个clust寻找top20 marker gene。但是原文使用的limma包识别,去重后仅有96个gene,而我自己尝试的或还有227个,相差比较大。The differential analysis identified 8,025 marker genes. The top 20 marker genes of each cell cluster are displayed in the heatmap. A total of 96 genes are listed beside of the heatmap after omitting the same top marker genes among clusters. The colors from purple to yellow indicate the gene expression levels from low to high.
p3_1 | p3_2
下二个子图
8、总图合并
p <- (p1_1 | (p1_2 | p1_3) ) /
((p2_1| p2_2 | p2_3) /
(p3_1 | p3_2))
ggsave("my_try.pdf", plot = p, width = 15, height = 18)
- SEO技巧:Shell脚本自动提交网站404死链到搜索引擎
- Nginx发布1.9.0版本,新增支持TCP代理和负载均衡的stream模块
- WordPress4.2升级修复补丁:解决大量404请求以及评论表情路径及尺寸异常问题
- Linux系统编译安装Redis以及主从复制配置小记
- Go-Maps
- 为WordPress开启Nginx缩略图功能,七牛从此陌路
- 为网站开启Nginx缓存加速,支持html伪静态页面
- 解决WordPress升级4.2后调用国外图片导致大量404请求的问题
- JS代码实现浏览器网页标题的动态切换,略微提高网站粘性
- Go-List
- 分享张戈博客自用的php网址在线转换二维码的API源码
- zabbix agentd客户端插件Shell一键自动安装脚本
- SendCloud邮件队列状态和已使用额度的Python监控脚本
- linux/scp命令报“bash: scp: command not found lost connection”错误的解决办法
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法