背景
之所以要和大家分享“认识Seurat数据结构”这样的文章,往小了说是因为以后大家在使用Seurat从V4过渡到V5的时候,一定会遇到因为数据结构的细微变化或者数据的指定、索取方式的改变,从而出现各种奇奇怪怪的问题。甚至如果大家是现在才入手学习的话,默认是安装的V5版本,而网络上比较成熟的Seurat单细胞分析的中文教程还是以V4为主的,直接去套用他们现成的代码一定会遇到很多麻烦。
前两篇笔记中,我们详细介绍了Seurat V5的标准流程还有一些可视化的方法。第三篇笔记我们初步学习了R中的数据类型,并利用公开的一个PBMC数据集看了未经加工处理的数据集原始结构。
手把手教学|单细胞分析基础流程(二):Seurat数据可视化
前面三篇内容共同构成了我们今天学习内容的基础。
这一次我们就一起来看看对这个数据集创建Seurat以后,那些数据分析的标准流程会对数据结构产生什么样的影响。
01 创建Seurat对象
pbmc.data
可以看到,我们创建的这个Seurat对象中也有类似于pbmc.data的数据结构的信息,存储在pbmc@assays$RNA@layers$counts中,我们可以看到主要的区别在于两个维度(Dim)数量的变化,基因的数量从32738个变成了13714个,细胞数没发生变化,这是因为我们在创建对象的时候简单地对基因数和细胞数进行了过滤(min.cells = 3, min.features = 200)。
02 标准预处理
之前说过标准预处理主要是设置质控标准、进行质控以及选择细胞进行后续分析。
我们通常使用PercentageFeatureSet()函数来计算线粒体基因相关的reads的百分比。代码如下:
我们就不再对线粒体基因进行可视化了,也不再去摸索质控的阈值,这部分之前已经提到过了。因此我们直接按照官网给的阈值进行过滤
可以看到细胞数量和有表达量的基因的数量都发生了相应的变化。
因为我们设置的过滤条件是只保留那些基因数量在200到2500之间且线粒体基因比例不超过5%的细胞。
因此基因数量没有变,依旧是13714,细胞数量从2700变成了2638,所有细胞当中有表达量的基因的数量同样变成了2238732。
03 数据标准化
这一步默认使用的是LogNormalize
在这一步之后可以看到明显的变化。pbmc@assays$RNA@layers原来是只有一个list的,那就是counts,但现在在数据标准化以后,多了一个list,即data。并且data当中的基因表达量已经经过了标准化,不再是整数。
这里面精确地记录了标准化的时候使用的参数,可以根据这些信息回溯、优化分析过程。
04 鉴定高可变基因
在鉴定高可变基因以后,多了一个pbmc@assays$RNA@meta.data(注意不是下面的那个meta.data),这里面储存着各个基因的统计值。
同时,鉴定高可变基因这一步的参数信息也被储存在pbmc@commands$FindVariableFeatures.RNA中。
05 数据缩放
这次我们还是和上次一样,对全部基因都进行缩放。
06 PCA降维
PCA的结果放在pbmc@reductions$pca中,这一步使用的参数信息存储在pbmc@commands$RunPCA.RNA中。
07 确定PCA维度
08 细胞聚类及其可视化
参数信息储存在pbmc@commands$FindNeighbors.RNA.pca中。但因为environment中pbmc@commands内容太多,所以参数信息直接被折叠隐藏了。
我们可以通过View(pbmc)进行查看
也可以直接输入pbmc@commands$FindNeighbors.RNA.pca命令,把详细信息调取出来。
09 寻找差异表达基因
差异表达基因的可视化就不展示了,不会影响到数据结构的。
10 细胞注释
我们还是根据官网给出的注释结果来做
提取信息
上面的部分介绍了Seurat数据结构的变化过程,从创建Seurat对象到标准预处理、数据标准化、鉴定高可变基因、数据缩放、PCA降维、确定PCA维度、细胞聚类及其可视化、寻找差异表达基因和细胞注释。在每一步操作中,数据结构都会发生相应的变化。
在一步一步认识了每一层数据所代表的含义以后,我们就可以尝试从pbmc这个大的数据集当中提取出想要的信息来进行个性化的分析。
提取特定的表达矩阵
例如,首先取未经标准化counts的1,3行,然后取2,4列。行是每一个细胞,列是每一个基因。
在R中可以使用方括号来选择行、列、元素。[i,]表示第i行,[,j]表示第j列,[i, j]就表示第i行第j列的那个元素。
上面是因为如果我想展示第1,3行的全部元素,就意味着要展示第1个和第三个细胞的全部基因,那样太占篇幅了,所以只展示了前100个。
然后我们取第1,3行的第2,4列
我们再换一个例子,直接取标准化之后data的前15行及其前20列
提取细胞分类信息
提取细胞子集
提取细胞子集是非常常用的一个功能。比如我们对某个器官组织进行了单细胞测序以后,想要在整个组织的众多类型细胞当中注释出来某一种数量很少、没有特异性marker的细胞是基本不可能的。最好的办法就是先进行大类的注释,然后提取细胞子集进行更加详细的注释。
比如这篇2023年的Cell Metabolism上发表的文献,就是从整群细胞→免疫细胞→肿瘤相关髓系细胞逐层细分,最后将研究中心放在了某一两群TAM的亚群上。
从上面的分析流程来看,非常多的关于筛选细胞的信息都可以在pbmc@meta.data找到,因此我们可以用这里的信息作为筛选标准来提取细胞子集。
但是这里面居然没有非常重要的我们注释后的celltype信息,那假如我们想根据注释后的celltype来提取子集该怎么办呢?
可以看到在图中有两行信息:
active.assay指的是当前活跃的测定,即正在被操作或分析的数据集。对于大多数操作,这个数据集将会被默认使用。
active.ident则是当前活跃的身份类别,通常是指细胞的分群结果。在我们进行细胞注释以后,这里自然就变成了注释后的结果。因此我们可以直接调用idents,它默认是注释之后的celltype。
所以要根据注释后的celltype来提取子集有两种办法,第一种是根据“当前活跃的身份类别”直接提取
第二种方法是把celltype信息给添加到meta.data当中去,这种方法能成立的原因其实和第一种方法是一样的
这里有个细节,用pbmc@meta.data$celltypes的信息提取时,要用==而不是之前的=。
以上就是我们本次分享的全部内容了。希望大家一路跟下来能够理解数据结构变化的思路,在单细胞的分析上能越来越熟练。
本文作者是"Algernon"同学,在获得授权后,实验老司机将本文发表于公众号。
文稿:Algernon
校对:煲仔饭
素材:
https://images.app.goo.gl/h1G593XFqEQoDg636
经验总结|生物信息学的学习方法
手把手教学|做完蛋白质谱后的KEGG、GO分析
经验总结|实验新手必读:提高效率与质量的实用技巧