刘小泽写于19.6.18、23-第二单元第四讲：得到表达矩阵后看内部异质性
笔记目的：根据生信技能树的单细胞转录组课程探索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个基因呢？很显然，利用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加上去
# 先构建一个分组的数据框，列是原来的分组信息
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"
)
