查看原文
其他

由表达矩阵看内部异质性

刘小泽 单细胞天地 2022-06-06

课程笔记




粉丝:有单细胞线上课程吗?

小编:什么? 我们的单细胞转录组分析线上课程已经上线好久了,你们竟然都不知道吗,每篇推文后面的课程推荐没人看的吗,小编已哭晕在厕所

好了,戏演完了,下面郑重介绍下我们的单细胞线上课程:(详情戳下方链接) 

全网第一个单细胞转录组数据分析实战视频教程(预售第二波通知)

课程链接在: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加上去

最后设置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结果




一篇文章带你走进单细胞的天地

关于人类细胞图谱项目详细再现性的案例研究

单细胞测序助力人体肠道类器官培养

第一届生物信息人才发展峰会诚邀演讲嘉宾及参会者

不同物种肾脏常驻巨噬细胞候选基因表达标签的确定

基质微环境对胰腺癌有何影响?

听说我们生信技能树论坛搜索功能失效?

单细胞测序揭示毛囊真皮微环境的分子分化轨迹

表达矩阵处理—表达QC(reads)

单细胞测序探究小鼠心脏的动态变化



如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程

一年一度的生信技能树单细胞线下培训班火热招生

生信技能树(爆款入门培训课)全国巡讲约你

全国巡讲第11站-港珠澳专场(生信技能树爆款入门课)






您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存