基于基因集的样品队列分组之层次聚类
随着ngs价格的持续走低，转录组测序项目早就走入了大样品时代，当然了，早在芯片价格亲民的时候就有这样的趋势，目前单细胞转录组价格也是在走这个老路。
那么，对于大样品队列的转录组，很多时候是没有已知的合理的分组， 这个时候会人为的去分组后看队列异质性，比如根据免疫高低进行分组。
那么这个根据免疫高低进行分组就有多种实现方式，我们这里简单的演示一下PCA和热图的层次聚类以及gsea或者gsva这样的打分的分组，看看是否有区别。
首先看看热图的层次聚类分组
需要载入 step1-output.Rdata 这个文件里面的表达量矩阵哦，如果你不知道 step1-output.Rdata 如果得到，看文末的代码。
首先，挑选目标基因集的表达量矩阵，进行热图并且顺便层次聚类，然后简单的暴力分组；
load(file = 'step1-output.Rdata')
cg=c('CD3D','CD3G CD247','IFNG','IL2RG','IRF1','IRF4','LCK','OAS2,STAT1')
cg
cg=cg[cg %in% rownames(dat)]
library(pheatmap)
p=pheatmap(dat[cg,])
hc=cutree(p$tree_col,5)
ac=data.frame(hc=as.character(hc))
rownames(ac)=colnames(dat)
pheatmap(dat[cg,],annotation_col = ac)
得到如下所示：
层次聚类暴力分组
可以看到， 1和2在热图的左右两边，而3，4，5在中间，其中5个分组里面居然就一个样品。所以我们需要把暴力分组调整为合理的免疫基因高低分组，代码如下所示：
group_list=ifelse(hc <3 ,'low','high')
table(group_list)
ac=data.frame(group_list)
rownames(ac)=colnames(dat)
pheatmap(dat[cg,],annotation_col = ac)
这个时候可以看到样品很清晰的分成了免疫高低两个组：
层次聚类合理分组
不过，这样的分组，数量并不是均等的哦！
> table(group_list)
group_list
high low
38 69
值得一提的是 这样的免疫基因的高低分组是一个数据集内部的高低概念哦，并不能跨越数据集去合并哦。
我们也可以随便把这个免疫基因集的表达量矩阵进行PCA看看，高低分组后全局表达量矩阵其实很难在PCA上面区分开来，但是在这个免疫基因集的小表达量矩阵是没有问题，如下所示：
library("FactoMineR") #画主成分分析图需要加载这两个包
library("factoextra")
dat.pca <- PCA(as.data.frame(t(dat[cg,])) )
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
addEllipses = T,
legend.title = "Groups"
)
清晰可见的两个分组的界限：
免疫高低两个组的清晰差异界限
关于 step1-output.Rdata 这个文件
上面的代码载入 step1-output.Rdata 这个文件，下面给出来这个文件的制作方式，代码如下所示：
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
gset <- geoChina("GSE58812")
gset
gset[[1]]
a=gset[[1]] #
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
# GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
boxplot(dat[,1:4],las=2)
dat=log2(dat)
boxplot(dat[,1:4],las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat[,1:4],las=2)
pd=pData(a)
#通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
library(stringr)
table(pd$`meta:ch1`)
group_list= ifelse(pd$`meta:ch1` == '0' ,'stable','meta')
table(group_list)
dat[1:4,1:4]
dim(dat)
a
ids=idmap('GPL570','soft')
head(ids)
ids=ids[ids$symbol != '',]
dat=dat[rownames(dat) %in% ids$ID,]
ids=ids[match(rownames(dat),ids$ID),]
head(ids)
colnames(ids)=c('probe_id','symbol')
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
dat[1:4,1:4]
ids=ids[ids$probe_id %in% rownames(dat),]
dat[1:4,1:4]
dat=dat[ids$probe_id,]
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
dat['ACTB',]
dat['GAPDH',]
save(dat,group_list,
file = 'step1-output.Rdata')
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
腾讯云开发者
