单细胞转录组数据处理之降维聚类分群

前面我们一起学习了单细胞转录组数据的上游分析,而且了解了自己的项目的样本数量和测序量,还过滤了不合格的细胞和基因, 系列教程目录如下:
理论上我们已经足够认识表达矩阵了,现在可以开始单细胞转录组数据的主角:降维聚类分群。如果你的背景知识不足,也可以先读一下综述,我们单细胞天地有中文指引
降维聚类分群是一条龙分析
我们并不是开发单细胞数据处理算法,所以大概率上,大家其实会把降维聚类分群一起做了,在seurat3里面的代码是:
sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000) 
# 步骤 ScaleData 的耗时取决于电脑系统配置(保守估计大于一分钟)
sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.2)
# 这个 resolution 可以调整,值越大,分出来的细胞亚群越多,默认是 0.8 
table(sce@meta.data$RNA_snn_res.0.2)

首先看两种降维

简单解释一下,这代码里面的FindVariableFeaturesRunPCA函数,是两种不同策略的降维。
  • 首先FindVariableFeatures是硬过滤,根据一些统计指标,比如sd,mad,vst等等来判断你输入的单细胞表达矩阵里面的2万多个基因里面,最重要的2000个基因,其余的1.8万个基因下游分析就不考虑了。
  • 然后RunPCA函数其实跑完之后2000个基因会转变为2000个维度,但是我们通常看前面的十几个维度就ok了,所以也是一个效率非常高的降维方式。
这个时候,你一定有疑问,为什么FindVariableFeatures是挑选2000个基因而不是其它数量呢?RunPCA函数跑完后我们应该是挑选前多少个维度呢?10个还是15个,还是20个还是50呢。

然后看聚类分群

聚类分群是紧密连接的,细胞可以看做是空间的不同点,如果是二维平面空间,点与点之间的距离很方便计算,距离的远近就决定着细胞是否属于一个类群。
计算距离的公式很多,我们就不一一展开,但是需要注意的是二维平面空间,三维球体空间的细胞距离很方便计算,但是如果是50个维度的空间,计算几万个细胞之间的距离就很可怕了,如果是2000个维度,甚至是2万个维度,基本上个人计算机就可以放弃了。这就是为什么我们前面通常是需要降维的。
下面是一个典型的单细胞转录组项目数据处理的描述:
降维聚类分群
可以看到他们的第一步降维是,选取top 5000的表达量离散度大的基因,第二步降维是选取top20的主成分。使用KNN-graph的聚类,最终定下来了10个细胞亚群。
一般来说,如果单细胞转录组数据仅仅是文章生物学故事的一个环节,就会采取标准的seurat流程,如下所示:
标准的seurat流程
如果你看的文献足够多,还会发现,在降维聚类分群之后,通常是有一个细胞在二维平面的散点图展示,如下所示:
细胞在二维平面的散点图
如果你足够心细,也会发现其实细胞的空间距离排布坐标通常是tSNE和umap来展现。
那我们说的tSNE和umap是怎么回事呢
tSNE(t-distributed stochastic neighbor embedding) 是一种非线性降维的算法,是流形学习的一种,当然,你大概率上是无需理解流形学习的,认识一下其它流形学习的经典方法即可,包括:Isomap,LLE,LE和diffusion maps等。我发现这篇推文介绍的非常好:单细胞中的流形(一):理解 tSNE中的perplexity,看完你需要记住的是;
  • 困惑度(perplexity)可以表示细胞的邻近个数,在tSNE图上的直观反映是细胞点的分布是否紧凑。perplexity设置越大,细胞分布越紧凑。
  • tSNE的参数设置:perplexity < (细胞数-1)/3,建议perplexity = 细胞数 / 50;
  • tSNE倾向于保留数据的局部结构。
本来我还期待着他们继续讲解umap,但是看起来似乎公众号已经停更了,唉,现在搞分享,写原创教程也挺不容易的, 很容易就quit了。
我给大家的策略是,反正你得多尝试,umap啥都尝试一遍,最后选择效果好的图表去展现即可。
可以修改的参数
第一次降维,表达矩阵里面通常是2~3万的基因数量,需要筛选到千这个数量级,可以是1000-5000,都是可以

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:
(0)

相关推荐

  • 我从Science中,偷学到这个聚类分析技能!真舍不得分享

    大家好,我是风.欢迎来到风风的从零开始单细胞系列.前面我们已经学习了数据下载.构建分析对象和数据质控.如果你的scater出现了一些警告内容,提示函数被替代,那也不用着急.正如其他内容一样,scate ...

  • 单细胞转录组分析肿瘤异质性

    欢迎来到医科研,这里是白介素2的读书笔记,跟我一起聊临床与科研的故事, 生物医学数据挖掘,R语言,TCGA.GEO, SEER数据挖掘. 肿瘤异质性指的是不同肿瘤细胞可以表现出不同的遗传和表型特征,包 ...

  • 明码标价之单细胞转录组的质控降维聚类分群和生物学注释

    一般来说,公共数据集都会给出表达量矩阵和具体不同细胞亚群特异性基因,比如 GSE122083 数据集背后的文献,就给出来了这些分群: NK (NKG7 and GNLY), NKT (CD3D and ...

  • 单细胞转录组的质控降维聚类分群和注释哪个步骤最关键

    我们非常强调进入一个领域需要读综述来获取基本认知,尤其是单细胞,我们在<单细胞天地>公众号给大家精选了2017-2020的4篇综述: 2017年7月的 Identifying cell p ...

  • 四行代码完成单细胞转录组的降维聚类分群

    昨天在<生信技能树>公众号的一个教程:这也能画?,我提到了一个很无聊的R包,名字是:scRNAstat ,它可以4行代码进行单细胞转录组的降维聚类分群,其实完全没有技术含量, 就是把 Se ...

  • 表达量矩阵全部更改为0-1矩阵会影响降维聚类分群吗?

    常规的读入10x的3个文件,需要自己根据下面的网址去下载 pbmc3k_filtered_gene_bc_matrices.tar.gz 文件,并且解压哦,然后  Read10X 函数读入解压后的文件 ...

  • 单细胞转录组数据处理之细胞亚群继续分群

    前面我们一起学习到了单细胞转录组数据的降维聚类分群,而且拿到的亚群,也可以解释其生物学意义,见:单细胞转录组数据处理之细胞亚群注释 其实理论上细胞亚群是可以无限划分的,因为世界上没有两个一模一样的细胞 ...

  • 人人都能学会的单细胞聚类分群注释

    作为生物信息学教学队伍的财务一名,我旁观了大量代码实战技巧,也勉强是学会了一下R语言,恰好看到朋友圈单细胞比较火爆,而且群主的CNS图表复现超级容易理解,我也跟着学习了一下,目录如下: CNS图表复现 ...

  • CNS图表复现04—单细胞聚类分群的resolution参数问题

    回顾 我们的CNS图表复现之旅已经开始,前面3讲是: CNS图表复现01-读入csv文件的表达矩阵构建Seurat对象 CNS图表复现02-Seurat标准流程之聚类分群 CNS图表复现03-单细胞区 ...

  • 单细胞转录组数据处理综述

    很久以前无意中翻译过一篇单细胞的新闻,单细胞测序 也关注过这方面进展,北大谢晓亮组又更新了他们的单细胞全基因组扩展方法 正好我们博士阶段有一门课是写一个综述系统性的总结一个研究领域,我就很自然的选择了 ...

  • cytofWorkflow之聚类分群(四)

    前面我们公布了<cytof数据资源介绍(文末有交流群)>,现在就开始正式手把手教学. 上一讲我们构造好了SingleCellExperiment对象,后续全部的分析都会以这个SingleC ...