seurat提取表达矩阵_Seurat
library(dplyr)
library(Seurat)
# 10X的数据可以使用Read10X这个函数,会返回一个UMI count矩阵,其中的每个值表示每个基因(行)在每个细胞(列)的分子数量
pbmc.data
# class(pbmc.data)
# str(pbmc.data)
# 利用这个UMI count矩阵构建Seurat对象,(使用非标准化的原始数据先构建一个对象);
# Seurat对象是一个容器,里面装了数据(比如表达矩阵)和分析结果(比如PCA、聚类)。关于seurat对象的详细解释,见: https://github.com/satijalab/seurat/wiki/Seurat#slots
# 在CreateSeuratObject时,会自动计算unique基因数量和总分子数,可以在object meta data中找到。
pbmc
pbmc #结果有个active assay,它存的就是表达矩阵,用pbmc[["RNA"]]@counts就可以访问
str(pbmc) #可以利用str来查看seurat对象的信息。
pbmc[[]] #使用[[ ]]是对assay的快速访问。
pbmc@meta.data #通过object@meta.data也可以访问。
pbmc@version #如果要查看其他的信息,可以用@
pbmc[["RNA"]]@counts #
#关于seurat对象还需要进一步学习!
# ------------------------------------------------------------------------
# 比较稀疏矩阵和矩阵占空间大小(可无!) ----------------------------------------------------------
# # What does data in a count matrix look like?
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30] #这个UMI count矩阵长什么样?和我们平常见到的bulk转录组矩阵一样吗?找几个基因,看看前30个细胞的情况
#从结果可以发现它是sparse Matrix,很杂乱。用 . 表示空着,数值为0,即没有检测到该分子。因为单细胞数据和bulk转录组数据还不太一样,scRNA的大部分数值都是0,原因一是由于建库时细胞的起始反应模板浓度就很低;二是测序存在dropout现象(本来基因在这个细胞有表达量,但没有测到)
# Seurat采用了一种聪明的再现方式,比原来用0表示的矩阵大大减小了空间占用,可以对比一下:
dense.size
dense.size
sparse.size
sparse.size
dense.size/sparse.size #比较稀疏矩阵跟矩阵之间占内存的大小,发现稀疏矩阵更省空间;
# -----------------------------------------------------------------------
# Standard pre-processing workflow ----------------------------------------
#以下为Seurat中scRNA-seq数据的标准预处理流程。包括:基于QC指标选择和过滤细胞;数据的normalization 和 scaling;寻找高异质性基因。
# * * * 质控(QC)及筛选细胞 -----------------------------------------------------
# Seurat常用的QC指标:
# 1)* 每个细胞中检测到的unique基因的数量。(因为劣质细胞或空droplets基因数很少;一个droplet中有多个细胞基因数会很多;)
# 2)* Similarly, the total number of molecules detected within a cell (correlates strongly with unique genes) ???
# 3)* map到线粒体基因组上的reads比例。 (因为低质量或死亡的细胞会有大量线粒体污染;)
# 这篇文章 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4758103/ 介绍了一些QC的标准
# Where are QC metrics stored in Seurat? ----------------------------------
# 质控筛选,通过[[ ]]向对象添加列,这里[[ ]]符号在pbmc的metadata这一子集里面添加新列,利用这个方法还可以挑出其他的feature;把表达矩阵里的线粒体QC单独存放在Seurat对象里,等于添加一列到metadata的对象里,
head(pbmc@meta.data, 5) #这时候看一下发现只有三列
pbmc[["percent.mt"]]
# # 检查下metadata, 查看QC指标存储的地方。
head(pbmc@meta.data, 5) #如果没有上一步, pbmc[["percent.mt"]]则就不会有percent.mt这一列!meta.data里面多了一列内容。
# 对QC指标可视化,并依QC指标过滤细胞;# unique基因数大于2500或者小于200的都去掉,mitochondrial counts >5%的也去掉。
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# 还可以用FeatureScatter函数画散点图来可视化两组feature信息的相关性,可以使用Seurat对象里的任何东西,如对象中的列、PC分数等。
plot1
plot2
CombinePlots(plots = list(plot1, plot2)) #这个图的结果说明有表达量的reads和在线粒体的reads数量一般成反比,线粒体reads占比很高的情况一般有表达量的很少;相反,如果是真正表达的基因reads,很少会来自线粒体基因组;(右边????)就是刚才所说的比对结果中的基因数量(feature)一般会和测序得到的reads count值成正比
# 细胞过滤 --------------------------------------------------------------------
pbmc 200 & nFeature_RNA < 2500 & percent.mt < 5) # 过滤掉不想要的细胞,
# 数据归一化(normalize) --------------------------------------------------------
# 为何normalize要叫“归一化”?
# 不必纠结名称,原理很好理解:因为我们关心的是差异表达基因,即一个基因在不同样本(不同细胞)的差异,但这种差异往往不能直接比较,因为不同样本的测序深度和文库大小可能不同(比如单细胞建库时每个EP管加的引物,酶之类的可能不同),那么基因表达量差异就存在一部分因素来自这里。
#为了更真实反映基因差异的生物学意义,就要将数据进行一个转换(例如RPKM、TPM等),可以让同一基因在不同样本中具有可比性。 “归一”归的就是这个起跑线。另外,原始表达量是一个离散程度很高的值,
#有的高表达基因表达量可能成千上万,可有的却只有几十。即使采用了一些归一化的算法消除了一些技术性误差,但真实存在的表达量极差往往会影响之后绘制热图、箱线图的美观性,#于是可以另外采用log 进行一个区间缩放(却不会影响真实值),
#比如原来表达量为1的,用log2(1+1)转换后结果依然为1;原来表达量为1000的,log2(1000+1)转换后约为10。那么原来相差1000倍的变化,现在只差10倍,在不破坏原有数据差异的同时,使数据更加集中。
#关于Seurat归一化,可以看这一篇:https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf
pbmc
# 得到的结果存储在:pbmc[["RNA"]]@data
# 鉴定差异基因 ------------------------------------------------------------------
#在细胞与细胞间进行比较,选择表达量差别最大的(即同一个基因在有的细胞中表达量很高,同时在部分细胞中表达很少),利用FindVariableFeatures函数,会计算一个mean-variance结果,也就是给出表达量均值和方差的关系并且得到top variable features。
# 计算top variable features有三种算法的,分别是:vst(默认)、mean.var.plot、dispersion
pbmc
top10
# plot variable features with and without labels
plot1
plot2
CombinePlots(plots = list(plot1, plot2))
# 数据标准化(scale) ------------------------------------------------------------
# scale的目的是实现数据的线性转换(scaling),这也是降维处理(如PCA)之前的标准预处理。主要利用了ScaleData函数。前面的归一化(normalize)是log处理,它是对所有基因的表达量统一对待的,最后放在了一个起跑线上。但是为了真正找到差异,我们还要基于这个起跑线,考虑不同样本对表达量的影响。
# scale使得每个基因在所有细胞中的表达量均值为0,方差为1
# 这一步主要目的就是下一步的降维做准备,使用全部基因是比较慢。如果想提高速度,可以只对鉴定的HVGs(高异质性基因, 之前FindVariableFeatures设置的是2000个)进行操作,其实这个函数默认就是对部分差异基因进行操作。
# 这样操作的结果中降维和聚类不会被影响,但是只对HVGs基因进行标准化,下面画的热图依然会受到全部基因中某些极值的影响。所以为了热图的结果,还是对所有基因进行标准化比较好。
# ScaleData这个函数有两个默认参数:do.scale = TRUE和do.center = TRUE,然后需要输入进行scale/center的基因名(默认是取前面FindVariableFeatures的结果)。center就是对每个基因在不同细胞的表达量都减去均值;scale就是对每个进行center后的值再除以标准差(就是进行了一个z-score的操作).
# * * * 使用全部基因进行scale ----------------------------------------------------
all.genes
length(all.genes)
pbmc
# * * * 只对鉴定的HVGs进行scale -------------------------------------------------
pbmc
# * * * 移除不想要的差异来源 -------------------------------------------------------
# 怎样移除不想要的差异来源?比如说不想线粒体污染导致的差异影响整体差异分布,可以按下列方式进行scale
pbmc
# 线性降维 --------------------------------------------------------------------
# 降维的真实目的是尽可能减少具有相关性的变量数目,例如原来有700个样本,也就是700个维度/变量,但实际上根据样本中的基因表达量来归类,它们或多或少都会有一些关联,
# 这样有关联的一些样本就会聚成一小撮,表示它们的”特征距离“相近。最后直接分析这些”小撮“,就用最少的变量代表了实际最多的特征
# 使用scale后的数据进行 PCA降维。它默认选择之前鉴定的差异基因(2000个)作为input,但可以使用features进行设置;默认分析的主成分数量是50个
pbmc
# 检查下PCA结果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5) #或者用 print(pbmc@reductions$pca, dims = 1:3, nfeatures = 5)
#通过 VizDimReduction, DimPlot和 DimHeatmap 对PCA结果进行可视化;
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca") #提取两个主成分(默认前两个,可以修改dims选项)绘制散点图。 这个图中每个点就是一个样本,它根据PCA的结果坐标进行画图,这个坐标保存在了cell.embeddings中:
# head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
# #还支持定义分组:根据pbmc@meta.data中的ident分组
# table(pbmc@meta.data$orig.ident) #当然这里只有一个分组
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE) #每个细胞和基因都根据PCA结果得分进行了排序,默认画前30个基因(nfeatures设置),1个主成分(dims设置),细胞数没有默认值(cells设置);其中balanced表示正负得分的基因各占一半
#也可以 DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
# 确定数据集的维数 ----------------------------------------------------------------
# 既然每个主成分都表示具有相关性的一小撮变量(称之为”metafeature“,元特征集),那么问题来了:降维后怎么选择合适数量的主成分来代表整个数据集?
# 可以通过JackStraw或者碎石图来判断;# JackStrawPlot更耗时,碎石图更快。
# * * * JackStraw法 -------------------------------------------------------
pbmc
pbmc
JackStrawPlot(pbmc, dims = 1:15)
# * * * 碎石图法 -------------------------------------------------------------
ElbowPlot(pbmc) #它是根据每个主成分对总体变异水平的贡献百分比排序得到的图,我们主要关注“肘部”的PC,它是一个转折点,一般认为拐点前的 PC(主成分)可以比较好地代表总体变化。
# 开始不确定时要多选一些主成分,选太少后续分析可能会出错。
# 细胞聚类 --------------------------------------------------------------------
pbmc
pbmc
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5) ## 结果聚成几类,用Idents查看
table(Idents(sce)) #也可以查看分成了几类。
#-----------------------------------------------------------------------
# 线性降维PCA的一个特点就是:它将高维中不相似的数据点放在较低维度时,会尽可能多地保留原始数据的差异,因此导致这些不相似的数据相距甚远,大多的数据重叠放置
# tSNE兼顾了高维降维+低维展示,降维后的那些不相似的数据点依然可以放得靠近一些,保证了点既在局部分散,又在全局聚集
# 后来开发的UMAP相比tSNE又能展示更多的维度,但是tSNE的图更好看一些(参考文献:https://sci-hub.tw/https://doi.org/10.1038/nbt.4314)
# 另一种降维方法—非线性降维(UMAP/tSNE) ------------------------------------------------
pbmc
DimPlot(pbmc, reduction = "umap") # 还可以设置label = TRUE让数字显示在每个cluster上
saveRDS(pbmc, file = "C:/Users/HASEE/Desktop/Seurat/pbmc_tutorial.rds") #对计算过程很费时的结果保存一下
# 找差异表达基因(cluster biomarkers) ---------------------------------------------
# 找差异表达基因的目的是根据表达量不同这个特征而分出不同细胞群的基因们(就是找具有生物学意义的HVGs,即biomarkers),marker可以理解成标记,就是一看到有这些基因,就说明是这个细胞群体。
# 利用FindAllMarkers()函数对单独的一个细胞群与其他几群细胞进行比较,找到正、负表达biomarker(这里的正负有点上调、下调基因的意思;正marker表示在我这个cluster中表达量高,而在其他的cluster中低)
cluster1.markers
head(cluster1.markers, n = 5)
cluster5.markers
head(cluster5.markers, n = 5)
pbmc.markers
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC) ## 这一步过滤好好理解!(进行了分类、排序、挑前2个)
cluster1.markers
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# 多种可视化方法 -----------------------------------------------------------------
# VlnPlot、FeaturePlot、DoHeatmap
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE) #用小提琴图对某些基因在全部cluster中表达量进行绘制
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A")) #最常用的可视化方法,将基因表达量投射到降维聚类结果中
top10 % group_by(cluster) %>% top_n(n = 10, wt = avg_logFC) #对给定细胞和基因绘制热图,例如对每个cluster绘制前20个marker
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
# 赋予每个cluster细胞类型 ---------------------------------------------------------
# 重点在于背景知识的理解,能够将marker与细胞类型对应起来,例如这里从各个cluster中找到的marker对应到了不同的细胞(这一过程是比较考验课题理解程度的)
new.cluster.ids
names(new.cluster.ids)
pbmc
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
saveRDS(pbmc, file = "C:/Users/HASEE/Desktop/Seurat/pbmc3k_final.rds")
seurat提取表达矩阵_Seurat相关推荐
- seurat提取表达矩阵_10X scRNA免疫治疗学习笔记-3-走Seurat标准流程
刘小泽写于19.10.15 笔记目的:根据生信技能树的单细胞转录组课程探索10X Genomics技术相关的分析 课程链接在:http://jm.grazy.cn/index/mulitcourse/ ...
- seurat提取表达矩阵_如何改造Seurat包的DoHeatmap函数?
刘小泽写于19.12.4 分析过单细胞数据的小伙伴应该都使用过Seurat包,其中有个函数叫DoHeatmap,具体操作可以看: 单细胞转录组学习笔记-17-用Seurat包分析文章数据 前言 走完S ...
- seurat提取表达矩阵_Hemberg-lab单细胞转录组数据分析
单细胞RNA-seq简介 混合RNA-seq2000年末的重大技术突破,取代微阵列表达芯片被广泛使用 通过混合大量细胞获取足够RNA用于建库测序,来定量每个基因的平均表达水平 用于比较转录组,例如比较 ...
- seurat提取表达矩阵_本周最新文献速递20200719
本周最新文献速递 一 (将不完全相似的表型进行荟萃分析,找出新的遗传位点,这是我第二次见到文献进行这种操作,一个新思路,值得借鉴) 文献题目: Multivariate genomic scan im ...
- seurat提取表达矩阵_GPL17586、GPL19251和GPL16686平台芯片ID转换
芯片分析中经常会遇到Affymetrix Human Transcriptome Array 2.0芯片,由于目前还没有现成的R包可以用,因此分析方法也不统一.见生信技能树Jimmy老师HTA2.0芯 ...
- seurat提取表达矩阵_单细胞分析实录(5): Seurat标准流程
前面我们已经学习了单细胞转录组分析的:使用Cell Ranger得到表达矩阵和doublet检测,今天我们开始Seurat标准流程的学习.这一部分的内容,网上有很多帖子,基本上都是把Seurat官网P ...
- seurat提取表达矩阵_什么,你一定要基于FPKM标准化表达矩阵做单细胞差异分析...
学徒和学员已经陆续出师,是时候把生信技能树的舞台交给后辈了!下面是<生信入门第7期>学员的分享最近看到有一个简单的数据挖掘文章,就做了一个单细胞表达矩阵的差异分析,然后强行解释一波.而且使 ...
- seurat提取表达矩阵_单细胞数据分析神器——Seurat
近年来,单细胞技术日益火热,并且有着愈演愈烈的趋势.在2015年至2017年,甚至对某细胞群体或组织进行单细胞测序,解析其细胞成分就能发一篇CNS级别的文章.近两三年,单细胞技术从最开始的基因组,转录 ...
- R语言 | GEO数据库下载GSE基因芯片 以及表达矩阵和临床信息的提取
目录 1.载入R包 2.利用AnnoProbe下载GEO数据库中的数据 3.提取表达矩阵和临床信息 4.输出文件 1.获得GEO数据库中的数据 下面以GSE14520数据系为例: 获得GEO数据库中的 ...
最新文章
- ubuntu系统无法ssh登录--安装openssh
- 尝试为文件附加自动命名的数据库,但失败。已存在同名的数据库,或指定的文件无法打开或位于 UNC 共享目录中...
- ASP无法上传大文件的解决方法
- 一个线性几何不等式猜想
- LeetCode - 121. 买卖股票的最佳时机
- C++实现邻接表存储的图及bfs遍历
- 简单的故事品味生活,
- python matplotlib 绘图操作
- 嵌入式工程师是青春饭吗?越老越吃香吗?
- 奥迪A8的L3级自动驾驶方案---奥迪A8的zFAS
- 小学听力测试英语软件,小学英语听力测试
- 分享电脑中截图的五种方法(包括截长图)
- 专科计算机专业软件工程,软件工程专业专科学校都有哪些?
- 凌晨3点不回家:成年人的世界不是他们说的那样
- Qt编写安防视频监控系统46-视频存储
- SpringCloud学习笔记(八)Gateway 网关
- Linux之dos2unix和unix2dos
- 《爱上Pandas》系列-你还在用VLookup吗?
- 最新C语言/C++语言培训项目实战(完整)
- c语言中vector的用法,c中vector的用法
热门文章
- 非u盘怎么安装linux系统文件,如何不用U盘和光盘装linux?
- css实现右上角角标
- 亿发软件:中小企业定制一体化管理解决方案,全面提升数据价值
- 第八章软件构造的性能——构造性能的度量、原则与方法(java中的垃圾回收机制及算法)
- R语言贝叶斯广义线性混合(多层次/水平/嵌套)模型GLMM、逻辑回归分析教育留级影响因素数据...
- 【17】AMOLED屏幕子像素定位
- HLSL中mul()函数的解释
- java斗图表情_java程序员斗图表情包 为何总是输
- macOS Catalina 10.15.6(19G73)原版镜像 by OpenCore-0.6.0-07-15编译版
- 研究生计算机专业知识复试面试常见问题