1 画热图,针对top100的sd的基因集的表达矩阵,没有聚类分组
cg=names(tail(sort(apply(dat,1,sd)),100)) ##取表达量标准差最大的100行的行名
pheatmap(dat[cg,],show_colnames =F,show_rownames = F,
filename = 'all_cells_top_100_sd.png')
2 针对top100的sd的基因集的表达矩阵，归一化，分组画图
cg=names(tail(sort(apply(dat,1,sd)),100)) ##取表达量标准差最大的100行的行名
n=t(scale(t(dat[cg,]))) #scale()函数去中心化和标准化
#对每个探针的表达量进行去中心化和标准化
n[n>2]=2 #矩阵n中归一化后，大于2的项，赋值使之等于2（相当于设置了一个上限）
n[n< -2]= -2 #小于-2的项，赋值使之等于-2（相当于设置了一个下限）
n[1:4,1:4]
ac=data.frame(g=group_list) #制作细胞（样本）分组矩阵
rownames(ac)=colnames(n) ##ac的行名（样本名）等于n的列名（样本名）
##判断分组矩阵的行（样本数）和表达矩阵的列（样本数）是否相等
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree1.png')
3 针对top100的sd的基因集的表达矩阵 进行重新 聚类并且分组。
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
##这个聚类分组只是对top100的sd的基因集
hc=hclust(dist(t(n)))
clus = cutree(hc, 4)
group_list=as.factor(clus)
table(group_list) ##这个聚类分组信息是针对top100的sd的基因集的，和针对全部基因集的分组结果不一样
table(group_list,df$g) ## 其中 df$g 是前面步骤针对全部表达矩阵的层次聚类结果。
## 下面针对本次挑选100个基因的表达矩阵的层次聚类结果进行热图展示。
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree_2.png')
dev.off() ##关闭画板
