专栏首页单细胞天地聚类算法之PCA与tSNE

聚类算法之PCA与tSNE

前 · 言

第二单元第六讲:聚类算法之PCA与tSNE

还是之前文章附件的图片,其中b图是选取两个主成分做的PCA图,c图是tSNE图:

几个常用函数的转置t(transpose),傻傻分不清?: 计算距离介绍过dist()函数,它是按行为操作对象,而聚类是要对样本聚类,因此要先将我们平时见到的表达矩阵(行为基因,列为样本)转置;同样PCA也是对行/样本进行操作,也是需要先转置;另外归一化的scale()函数虽然是对列进行操作,但它的对象是基因,因此也需要转置

关于PCA的学习,之前写过:

  • StatQuest-PCA学习:https://www.jianshu.com/p/b83ac8f7f5a7
  • StatQuest--在R中拆解PCA:https://www.jianshu.com/p/8a74508c3737

先构建一个非常随机的测试数据

# 设置随机种子,可以重复别人使用的随机数
set.seed(123456789)
library(pheatmap)
library(Rtsne)
library(ggfortify)
library(mvtnorm)
# 设置两个正态分布的随机矩阵(500*20)
ng=500 
nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc)
a2=rnorm(ng*nc);dim(a2)=c(ng,nc)
a3=cbind(a1,a2)
> dim(a3)
[1] 500  40
# 添加列名
colnames(a3)=c(paste0('cell_01_',1:nc),
               paste0('cell_02_',1:nc))
# 添加行名
rownames(a3)=paste('gene_',1:ng,sep = '')
# 先做个热图
pheatmap(a3)

没有体现任何的基因差异或者样本聚类(热图中的聚类是自然层次聚类),可以看到样本名都是无规律的交叉显示

如果做PCA呢?
# 先转置一下,让行为样本
>  a3=t(a3);dim(a3)
[1]  40 500

# prcomp()主成分分析
pca_dat <- prcomp(a3, scale. = TRUE)
p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
print(p)

可以看到每组的20个细胞都分不开,但每组具体有哪些样本还是看不出来,因此这里为每组加上颜色来表示

# 先在原来数据的基础上添加样本分组信息(别忘了a3是一个矩阵,先转换成数据框)
df=cbind(as.data.frame(a3),group=c(rep('b1',20),rep('b2',20)))
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()
另外看下tsne

利用了一个核心函数Rtsne()

set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0)
# 结果得到一个列表
> str(tsne_out)
List of 14
 $ N                  : int 40
 $ Y                  : num [1:40, 1:2] -36.7 -28 -168 -33.4 22.4 ...
 $ costs              : num [1:40] 0.00348 -0.00252 0.01496 0.01646 0.00951 ...
# 其中在Y中存储了画图坐标
> head(tsne_out$Y,3)
           [,1]      [,2]
[1,]  -36.72621 -78.03709
[2,]  -28.00151  33.30229
[3,] -167.98560 -80.26850

tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2") #为坐标添加列名
# 基础作图代码
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()
# 在此基础上添加颜色分组信息,首先还是将tsnes这个矩阵变成数据框,然后增加一列group信息,最后映射在geom_point中
tsnes=as.data.frame(tsnes)
group=c(rep('b1',20),rep('b2',20))
tsnes$group=group
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point(aes(col=group))

构建一个有规律的测试数据

ng=500
nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc)
# 和之前的区别就在a2这里,都加了3
a2=rnorm(ng*nc)+3;dim(a2)=c(ng,nc)
a3=cbind(a1,a2)
colnames(a3)=c(paste0('cell_01_',1:nc),paste0('cell_02_',1:nc))
rownames(a3)=paste('gene_',1:ng,sep = '')
pheatmap(a3)

热图已经能看出来差异了,再看看PCA

a3=t(a3);dim(a3)
df=cbind(as.data.frame(a3),group=c(rep('b1',20),rep('b2',20)))
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()

tsne也是如此

set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0)
tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2")
tsnes=as.data.frame(tsnes)
group=c(rep('b1',20),rep('b2',20))
tsnes$group=group
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point(aes(col=group))

真实数据演练

载入RPKM数据
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
# 表达量信息
> dat[1:2,1:3]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610007P14Rik              0              0       74.95064
0610009B22Rik              0              0        0.00000
# 样本属性
> head(metadata,3)
               g plate  n_g all
SS2_15_0048_A3 1  0048 3065 all
SS2_15_0048_A6 2  0048 3036 all
SS2_15_0048_A5 1  0048 3742 all
#所有数据的聚类分组信息
group_list=metadata$g 
#批次信息
plate=metadata$plate 
> table(plate)
plate
0048 0049 
 384  384 
对数据进行PCA
# 操作前先备份
dat_back=dat
# 先对表达矩阵进行转置,然后转换成数据框,就可以添加批次信息了
dat=dat_back
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,plate )

> dim(dat_back)
[1] 12689   768
> dim(dat)
[1]   768 12690

library("FactoMineR")
library("factoextra")
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca, # repel =T,
             geom.ind = "point", # 只显示点,不显示文字
             col.ind = dat$plate, # 按分组上色
             #palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # 添加晕环
             legend.title = "Groups"

可以看到两个批次之间分不开,说明没有批次效应

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:刘小泽

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-07-16

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • scRNA-seq原始数据的处理

    好的书籍是人类进步的阶梯,但有些人却找不到优秀的阶梯,为此我们开设了书籍翻译这个栏目,作为你学习之路的指路明灯;分享国内外优秀书籍,弘扬分享精神,做一个知识的传...

    生信技能树jimmy
  • scHLAcount || 单细胞转录组HLA等位基因分析

    scHLAcount允许我们使用个性化的参考基因组计算HLA I类基因HLA-A、B和C的单细胞转录组序列数据中的分子数;和HLA II类基因DPA1, DPB...

    生信技能树jimmy
  • 综述:高维单细胞RNA测序数据分析工具(中)

    当你的才华还撑不起你的野心时,请潜下心来,脚踏实地,跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了,通过文献速递这个栏目很幸运聚集了一些小伙伴...

    生信技能树jimmy
  • 【Android菜鸟学习之路】环境搭建问题-修改AVD Path

    悟空聊架构
  • Python 随机数生成

    以上实例我们使用了 random 模块的 randint() 函数来生成随机数,你每次执行后都返回不同的数字(0 到 9),该函数的语法为:

    于小勇
  • 『高级篇』docker之kubernetes基础集群部署(上)(33)

      ETCD保证了数据的存储,保证了数据的高可用,还有数据的一致性,它跟zookeeper类似。kubernetes需要存储很多东西,像它本身的节点信息,组件信...

    IT故事会
  • 分布式环境下对部分热数据(如redis热key,热请求)进行探测,并对探测结果及时同步到各个client实例的JVM内存的方案简述

    譬如突然大量请求都命中了redis的某个分片,造成该redis卡顿,影响其他请求。热key特性如 goodsId=100,突发1万请求该key。

    天涯泪小武
  • 安利3个Python数据分析EDA神器!

    EDA是数据分析必须的过程,用来查看变量统计特征,可以此为基础尝试做特征工程。东哥这次分享3个EDA神器,其实之前每一个都分享过,这次把这三个工具包汇总到一起来...

    用户2769421
  • HLA Dictionary-Allel和抗原之间的对应关系

    HLA Dictionary , 中文名叫做HLA 字典,记录了HLA-A, -B, -C, -DRB1/3/4/5 和 -DQB1 等HLA基因的不同Alle...

    生信修炼手册
  • 周期串

    Vincent-yuan

扫码关注云+社区

领取腾讯云代金券