【原】基于基因集的样品队列分组之gsea等打分

随着ngs价格的持续走低，转录组测序项目早就走入了大样品时代，当然了，早在芯片价格亲民的时候就有这样的趋势，目前单细胞转录组价格也是在走这个老路。

那么，对于大样品队列的转录组，很多时候是没有已知的合理的分组， 这个时候会人为的去分组后看队列异质性，比如根据免疫高低进行分组。

那么这个根据免疫高低进行分组就有多种实现方式，我们这里简单的演示一下PCA和热图的层次聚类以及gsea或者gsva这样的打分的分组，看看是否有区别。

gsea等打分后对样品队列的高低分组

前面我们已经分享了：基于基因集的样品队列分组之层次聚类，以及 基于基因集的样品队列分组之PCA，还剩下看gsea等打分后对样品队列的高低分组。

同样的需要载入 step1-output.Rdata 这个文件里面的表达量矩阵哦，如果你不知道 step1-output.Rdata 如果得到，看文末的代码。

针对指定基因集，去构建gsva函数所需要的 GeneSetCollection对象， 代码如下所示：

load(file = 'step1-output.Rdata')
cg=c('CD3D','CD3G CD247','IFNG','IL2RG','IRF1','IRF4','LCK','OAS2,STAT1')
cg
cg=cg[cg %in% rownames(dat)]
cgl = list(cg=cg)
library(GSVA)
library(GSEABase)  
geneSet <- GeneSetCollection(mapply(function(geneIds, keggId) {
  GeneSet(geneIds, geneIdType=EntrezIdentifier(),
          collectionType=KEGGCollection(keggId),
          setName=keggId)
}, cgl, names(cgl)))
geneSet

有了表达量矩阵， 以及指定基因集，运行gsva函数就是一句话代码的事情：

dat.gsva <- gsva( dat , geneSet, 
                 mx.diff=TRUE, verbose=FALSE,
                 parallel.sz= 8) 
dat.gsva =dat.gsva[1,]
group.gsva <- ifelse( dat.gsva>median(dat.gsva) , 'group1', 'group2')
table(group.gsva) 

因为我们这次就输入了一个基因集，所以就是每个样品对这个基因集打分即可。然后根据不同样品的打分进行高低分组后可视化。

group.gsva
group1 group2 
    53     54 
ac=data.frame(group.gsva)
rownames(ac)=colnames(dat)
library(pheatmap)
pheatmap(dat[cg ,],annotation_col = ac)

可以看到样品比较好的分成了数量比较一致的两个组：

热图如下所示：

热图

然后看主成分：

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.gsva, # color by groups 
             addEllipses = T,
             legend.title = "Groups"
) 

基本上也是类似的：

主成分

也可以自行去和已经分享了：基于基因集的样品队列分组之层次聚类，以及 基于基因集的样品队列分组之PCA，对比看看，加深你的理解哦。

关于  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')

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 [e-mail]

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。