scRNA-单细胞转录组学习笔记-6

刘小泽写于19.6.18 笔记目的:根据生信技能树的单细胞转录组课程探索smart-seq2技术相关的分析技术 课程链接在:http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53 第二单元第四讲:得到表达矩阵后看内部异质性

前言

依然是主要代码跟着流程走,这里只写一写我认为比较重要的内容

上次我们得到了原始表达矩阵a,过滤后的表达矩阵dat ,样本信息meata

简单看下:

> dim(a)
[1] 24582   768
> dim(dat)
[1] 12198   768
> head(meta)
               g plate  n_g all
SS2_15_0048_A3 1  0048 2624 all
SS2_15_0048_A6 1  0048 2664 all
SS2_15_0048_A5 2  0048 3319 all
SS2_15_0048_A4 3  0048 4447 all
SS2_15_0048_A1 2  0048 4725 all
SS2_15_0048_A2 3  0048 5263 all
# 样本信息中的g表示cutree分的4大聚类结果;plate为细胞板批次;n_g是每个细胞样本中有表达的基因数量;all暂时用不到

另外,注意最好每次运行代码之前,都要清空一下变量,然后设置不要将字符型变成因子型向量

rm(list = ls())  
options(stringsAsFactors = F)

想要做个热图

先构建分组信息,也就是提取出层次聚类信息

需要注意一点,count的表达矩阵和rpkm表达矩阵的聚类分组结果是不一样的

# 我们这里是count表达矩阵分组
> grp=meta$g
> table(grp)
grp
  1   2   3   4 
312 300 121  35 

然后想想,热图需要什么信息?

主要就是行、列,行是基因,列是样本。那么先对基因(行)进行设置:

因为dat矩阵相对于a虽然过滤掉了一万多基因,但是依然还剩一万多,然后我们有700多样本,那么可以算一下,这样的结果是10000*700的图,相当大,并且看不出什么含义。我们可以将基因数设置小一点,可以设置成1000,先看一下

好了,基因数有了(1000个),**但是取哪1000个基因呢?**很显然,利用headtail直接取前/后1000个基因是不能使人信服的,这里可以用sd 进行筛选,也就是取表达量标准差最大的1000个基因(也即是说,这1000个基因在所有的样本中表达差异最大,这样更像差异表达基因)

tail(sort(apply(dat,1,sd)),1000)
# 解释下代码:从里向外看=》apply对dat矩阵的每一行求sd值,然后用sort排序,默认从小到大,然后用tail从后到前,也即是从大到小取1000个
# 最后取出基因名
top_g=names(tail(sort(apply(dat,1,sd)),100))
> head(top_g)
[1] "Comt"    "Mrgprf"  "Stard13" "Cdipt"   "Ifnar1"  "Pdcd6ip"

画第一张热图

1000基因 X 768个细胞 = 70多万的点

这个热图反映了log2(cpm+1)的表达量,范围是0-15

library(pheatmap)
pheatmap(dat[top_g,],show_colnames =F,show_rownames = F,
           filename = 'all_cells_top_1000_sd.png')

热图基础上增加归一化

上面👆那个图,可以看出基因绝对表达量 ,颜色越偏红色表示绝对表达量越高,比如顶部那些基因的表达量就是要比底部那些基因的高

但是,有个问题,这样会受到某些特高表达基因的影响,导致其他基因的差异就不明显;另外,我们真正关心的是一个基因在不同样本中的差异,是一种相对的表达量。

可以这么理解:有的基因本身就是表达量小,但不能因为小就认为它在每个样本都是一样的。虽然小,也是有差异的小。但往往这种差异会由于"强者"的存在而被忽视。 因此,这里我们要做的,就是将这种"弱小"的差异给拎出来

主要利用scale() ,先理解一下:

# 构建一个测试数据
x=1:10;plot(x)
scale(x);plot(scale(x))# 结果就是变成从-1.4到+1.4的范围

可以看到,scale后并不改变数据分布,只是修改了坐标,让结果取值更加集中

注意:scale是对列进行操作,而我们是想对基因(也就是按行操作),这个函数有两个主要的选项:centerscale ,其中center是将每列的元素减去这一列的均值(这个选项是默认TRUE的);scale 是在center操作后,再将处理过的元素除以标准差(同样是默认TRUE的)。另外,处理完别忘了再转换回来

n=t(scale(t(dat[top_g,])))

画个新的热图

可以看到之前热图显示的坐标范围是0- 15,当然这里我们可以设置上限、下限,比如可以将上限设为2,下限设为-2

n[n>2]=2
n[n< -2]= -2

范围设置好以后,可以再将分组信息grp加上去

# 先构建一个分组的数据框,列是原来的分组信息
anno_col=data.frame(g=grp)
# 行名是样本名,也就是归一化后的n矩阵的列名
rownames(anno_col)=colnames(n)
# 看一下结果
> head(anno_col)
               g
SS2_15_0048_A3 1
SS2_15_0048_A6 1
SS2_15_0048_A5 2
SS2_15_0048_A4 3
SS2_15_0048_A1 2
SS2_15_0048_A2 3

最后设置pheatmap的选项:

    pheatmap(n,
             show_colnames =F, #不显示列名
             show_rownames = F, #不显示行名
             annotation_col=anno_col, # 在列上加注释(也就是分组信息)
             filename = 'all_cells_top_1000_new.png')

可以看到这张图和画的第一张图的区别是:出现了一些红色,说明新出现了一些基因在某些样本中高表达,并且很明显。但是仍然很有可能它们的实际表达量并不高,仅仅是玩了一个"样本排位赛“(即使数值再小,也有甲乙丙丁)

关于分组有一点奇怪

可以看到这里的分组信息有点散乱,想到:这里使用的anno_col 是利用grp得到的,而grp是基于一万多基因的dat矩阵得到的(回忆下第5篇内容)

因此这里的分组信息可以更新一下,基于我们这里的top1000基因,只需要将原来的dat换成现在的n矩阵就好,依然选取前4个聚类分群

# 将原来dat换为n
hc=hclust(dist(t(n))) 
clus = cutree(hc, 4)
top1000_grp=as.factor(clus)
> table(top1000_grp)
top1000_grp
  1   2   3   4 
462 166  42  98 

看一下当前基于1000个基因的前4组和原来基于所有基因的前4组之间重合度,虽然他们总和一样,都是1000,而且也都是按照1-4的顺序排列,但实际上背后的意义千差万别

> table(top1000_grp,grp)
           grp
top1000_grp   1   2   3   4
          1 167 295   0   0
          2   7   3 121  35
          3  42   0   0   0
          4  96   2   0   0

举个例子,有462个基因属于新分组的1号组,但其中有295个属于原来分组的2号组(这个数量超过了原来分组的1号组),可以看出**新分组和原分组的重合度并不高,**因此更加说明我们重新分组的重要性

new_anno_col=data.frame(g=top1000_grp)
rownames(new_anno_col)=colnames(n)
 pheatmap(n,
             show_colnames =F, #不显示列名
             show_rownames = F, #不显示行名
             annotation_col=new_anno_col, # 在列上加注释(也就是分组信息)
             filename = 'all_cells_top_1000_new_2.png')

做个PCA

之前好不容易过滤得到的dat矩阵,不能因为下面分析的失误被"污染”,因此再进行下一个分析之前先做一个数据备份是个好习惯

dat_bk=dat
# 然后我们就能放心对dat进行操作了
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,grp)

PCA分析需要行是样本,列是基因表达量的数据框(和聚类一样,是对行/样本进行操作,最后做的图中一个点就表示一个样本/细胞)

最后用PCA进行计算分析,用fviz_pca_ind函数进行可视化

这里用到的分组还是之前基于全部基因进行聚类的cutree结果

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$grp, # 按组分颜色
               # palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, 
               legend.title = "Groups"
  )

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related