由表达矩阵看内部异质性
- 2020 年 3 月 31 日
- 筆記
依然是主要代码跟着流程走，这里只写一写我认为比较重要的内容
上次我们得到了原始表达矩阵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个基因呢？很显然，利用head或tail直接取前/后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是对列进行操作，而我们是想对基因(也就是按行操作)，这个函数有两个主要的选项：center和scale ，其中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结果
