以复现图表的方式来学习一篇文章
1.文章简介
本次学习的2019年7月发表于Circulation.文章题目为:Single-Cell Analysis of the Normal Mouse Aorta Reveals Functionally Distinct Endothelial Cell Populations ”
2. 文章思路很直接套路
- 做Normal Aorta的 Single-Cell,对其进行分群,其中有一个亚群叫EC
- 重点关注亚群EC,对EC群再分群,做gene set enrichment 分析,做marker gene - 表达的小提琴图,heatmap图
- 实验验证EC群的特性。比如不同的diet之后(心血管疾病肯定和饮食之类的因素有关),做不同condition 下的single-Cell,看EC群基因变化,然后就和疾病关系靠一靠
作业复现的图表是Figure1 的1B 和1C
Figure1_C.png
3. 文件下载和解读
主要用Seurat V3. Broad Institute的single cell portal上面存放了他们发表的single-cell文章的数据 这里一共是3个文件
- Meta_DATA_Chow_12PCs_outfile.gz
- Cluster_File_Chow_12PCs_outfile.gz
- Seurat_Chow_12PCs_outfile.mtx 通过命名可以看出,文件均为跑过了PCA,tSNE分群后的输出数据,也就是说,这次的任务是非常下游的可视化过程。
读取文件1和2
安装加载包,我用的是Seurat v3
library(data.table)
library(R.utils)
library(Seurat)
library(ggplot2)
读取文件1和文件2
cluster_file <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/Cluster_File_Chow_12PCs_outfile.gz',
header = T,sep = 't',data.table = F)
meta_data <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/META_DATA_Chow_12PCs_outfile.gz',
header = T, sep = 't',data.table = F)
这里看一下两个文件。cluster_file显示的是tSNE_x和tSNE_y的坐标。这也是我们画fig1B tSNE图需要的。这个nametype列名是样品名 meta_data显示了分群信息。在Jimmy大神前面单细胞的讲座中,了解到需要将这些数据的行名(rownames)改为样品名。
head(cluster_file[1:3,1:3])
## NAME X Y
## 1 TYPE numeric numeric
## 2 AAACCTGAGACGCTTT-1 -19.0225914077646 0.924488198164653
## 3 AAACCTGCACTGTGTA-1 -26.7488891822262 -0.784745741413262
head(meta_data[1:3,1:3])
## NAME Cluster Sub-Cluster
## 1 TYPE group group
## 2 AAACCTGAGACGCTTT-1 VSMC VSMC
## 3 AAACCTGCACTGTGTA-1 VSMC VSMC
读取文件3 我在读取.mtx的文件上,入了好大的坑。看到后缀是.mtx文件,我的第一个反应是用readMM()读。但是,持续报错,error msg说这个文件不是稀疏矩阵。后来还是和前面一样用fread()读的,居然读出来了。然后我看了一下,是每个基因在每个样品中的表达,在数值上也有稀疏矩阵的.这个符号。读入后,我用了typeof()看这个.mxt,它显示的是list.这一部分我至今困惑。言归正传,为了后续的分析,需要将读入的.mtx文件的行名改成基因名。
mxt <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/Seurat_Chow_12PCs_outfile.mtx',header = T, sep = 't',data.table = F)
head(mxt[1:3,1:3])
## GENE AAACCTGAGACGCTTT-1 AAACCTGCACTGTGTA-1
## 1 0610007P14Rik 0.7264328 0.0000000
## 2 0610009B22Rik 0.0000000 0.6094115
## 3 0610009L18Rik 0.7264328 0.0000000
4. 重现fig1B
由于作者已经给出了tSNE的x,y轴坐标,这一步就特别简单。直接把有坐标信息的Cluster_File_Chow_12PCs_outfile.gz文件和有分群信息的Meta_DATA_Chow_12PCs_outfile.gz一整合,生成一个既有坐标信息,又有分群信息的文件(该文件叫dat),然后用ggplot画图就可以。
# 整合坐标信息和分群信息的文件
phe=data.frame(cell=meta_data$NAME[-1],
cluster =meta_data$Cluster[-1])
cluster_file_tSNE <- data.frame(cell = cluster_file$NAME[-1],
tSNE_1 =cluster_file$X[-1],
tSNE_2 = cluster_file$Y[-1],
avg_intensity =meta_data$`Average Intensity`[-1])
dat=merge(cluster_file_tSNE,phe,by.x = "cell",by.y="cell")
dat$tSNE_1 <- as.numeric(dat$tSNE_1)
dat$tSNE_2 <- as.numeric(dat$tSNE_2)
#画出图1B
dat_tsne =ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
Figure 1B
5.重现fig1C
fig1C的主体思路是找到每个群的marker gene,然后用Seurat v3自带Dotplot()可视化。marker gene的找法用Seurat自带的FindMarkers().
首先FindMarkers()里面的一个输入对象是object,这里需要使用CreateSeuratObject用.mtx这个文件创建一个对象。
# 创建Seurat对象
mxt_obj <- CreateSeuratObject(counts = mxt, project = "mxt", min.cells = 3, min.features = 200)
找marker gene需要分群信息。这里可以深刻体会AddMetaData()
的意义。通过AddMetaData()
,给CreateSeuratObject创建的对象提供更多的信息。作者已经在metadata文件里面提供了分群,那么通过AddMetaData()
,就直接可以把分群信息赋予创建的对象,不用再通过单细胞分析全流程FindCluster()
等函数计算。
另外,需要注意的是,如果想用FindMarkers()
这个函数找差异基因,它认定的分群需要用数字表示,不能识别字符串比如’VSMC’这种分群,所以得把作者提供的meta-data里面的分群信息改成数字1,2表示
# 将分群信息改成数字形式
meta_data_3 <- meta_data
meta_data_3$Cluster <- as.numeric(meta_data_3$Cluster)
head(meta_data_3)
## Cluster
## AAACCTGAGACGCTTT-1 6
## AAACCTGCACTGTGTA-1 6
## AAACGGGCACGGATAG-1 6
## AAACGGGGTGTTCTTT-1 6
## AAAGATGAGGCTCATT-1 6
## AAAGATGCAAGACACG-1 6
# 创建的对象中使用AddMetaData()添加分群信息
mxt_obj <- AddMetaData(object = mxt_obj,
metadata = meta_data_3,
col.name = 'Cluster')
进入FindMarkers()
这个函数。注意两个参数,group.by
和ident.1
-
-
ident.1
参数指明你要找的marker gene 所在的群。前面说了,ident.1
参数只认数字代表的分群,不认识字符串。所以,在meta-data里面,将分群信息改为数字
-
-
-
group.by
是划分的标准,比如,你找不同分群的marker gene,分群信息就是在meta_data里面的哪一列保存,列名就是group.by
的参数值。在这里,我们这个分群信息保存在列名为"Cluster"的列下面
-
# find marker gene for clsuter VSMC。VSMC对应的群的数字是6
markers_VSMC <- FindMarkers(mxt_obj,
group.by = 'Cluster',
ident.1 = 6,
test.use = "wilcox",
min.pct = 0.25)
文中提到,marker gene 的定义是 genes with the highest differential expression relative to all other genes. 所以得到的结果就按照avg_logFC排序,取前5。
注明:我找的marker gene和文章的做对比。就EC群的Ccl21a对不上。作者找的叫Ccl21a,我找了个,叫Cldn5.两不是一个基因
全部marker gene得到后,用seurat 自带的dotplot()
画图
DotPlot(object = mxt_meta2, features = Marker_gene,group.by = 'Cluster')+RotatedAxis()
Figure 1C
- ASP.NET MVC 4 - 测试驱动 ASP.NET MVC
- LVS+Keepalived高可用环境部署梳理(主主和主从模式)
- 随着区块链的火爆,相关顶级域名“矿池”KC.com已建站
- Flash/Flex学习笔记(50):3D线条与填充
- LVM常规操作记录梳理(扩容/缩容/快照等)
- Flash/Flex学习笔记(55):背面剔除与 3D 灯光
- 资源等待类型sys.dm_os_wait_stats
- NVIDIA不再允许数据中心用GeForce驱动,提供区块链服务除外
- 非常强悍并实用的双机热备+负载均衡线上方案
- Apache 压力测试工具ab
- SQL之收集SQL Server线程等待信息
- 聚合索引(clustered index) / 非聚合索引(nonclustered index)
- 域名资讯:单词域名can.com以15.5万美金成功交易
- jQuery无缝图片横向(水平)/竖向(垂直)滚动
- 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 数组属性和方法
- 打卡群刷题总结1005——跳跃游戏
- 真是活久见,在 Minecraft 的虚拟游戏里竟然还能管理 Kubernetes!
- 打卡群2刷题总结1007——反转链表
- 打卡群2刷题总结1001——两数之和 II - 输入有序数组
- 复杂一点的SQL语句
- PL/SQL Developer连接本地Oracle 11g 64位数据库
- 打卡群刷题总结1007——买卖股票的最佳时机 II
- 事务Transaction
- 打卡群2刷题总结1006—— 删除链表的倒数第N个节点
- 打卡群刷题总结1006——跳跃游戏 II
- 面试官常问的Spring依赖注入和Bean的装配问题,今天给大家讲清楚!
- 打卡群刷题总结1003——分割等和子集
- 打卡群2刷题总结1005——有效的括号
- 腾讯云服务器操作系统TencentOS安装与体验
- 打卡群2刷题总结1004——无重复字符的最长子串