隨著ngs價格的持續(xù)走低，轉(zhuǎn)錄組測序項目早就走入了大樣品時代，當(dāng)然了，早在芯片價格親民的時候就有這樣的趨勢，目前單細(xì)胞轉(zhuǎn)錄組價格也是在走這個老路。

那么，對于大樣品隊列的轉(zhuǎn)錄組，很多時候是沒有已知的合理的分組， 這個時候會人為的去分組后看隊列異質(zhì)性，比如根據(jù)免疫高低進(jìn)行分組。

那么這個根據(jù)免疫高低進(jìn)行分組就有多種實(shí)現(xiàn)方式，我們這里簡單的演示一下PCA和熱圖的層次聚類以及gsea或者gsva這樣的打分的分組，看看是否有區(qū)別。

gsea等打分后對樣品隊列的高低分組

前面我們已經(jīng)分享了：基于基因集的樣品隊列分組之層次聚類，以及 基于基因集的樣品隊列分組之PCA，還剩下看gsea等打分后對樣品隊列的高低分組。

同樣的需要載入 step1-output.Rdata 這個文件里面的表達(dá)量矩陣哦，如果你不知道 step1-output.Rdata 如果得到，看文末的代碼。

針對指定基因集，去構(gòu)建gsva函數(shù)所需要的 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

有了表達(dá)量矩陣， 以及指定基因集，運(yùn)行g(shù)sva函數(shù)就是一句話代碼的事情：

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) 

因?yàn)槲覀冞@次就輸入了一個基因集，所以就是每個樣品對這個基因集打分即可。然后根據(jù)不同樣品的打分進(jìn)行高低分組后可視化。

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

可以看到樣品比較好的分成了數(shù)量比較一致的兩個組：

熱圖如下所示：

然后看主成分：

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"
) 

基本上也是類似的：

也可以自行去和已經(jīng)分享了：基于基因集的樣品隊列分組之層次聚類，以及 基于基因集的樣品隊列分組之PCA，對比看看，加深你的理解哦。

關(guān)于  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現(xiàn)在是一個對象，取a這個對象通過看說明書知道要用exprs這個函數(shù)
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這個矩陣按行操作，取每一行的中位數(shù)，將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序，將對應(yīng)的行賦值為一個新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項，'!'為否，即取出不重復(fù)的項，去除重復(fù)的gene ，保留每個基因最大表達(dá)量結(jié)果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第一次出現(xiàn)的信息

dat['ACTB',]
dat['GAPDH',]
save(dat,group_list,
     file = 'step1-output.Rdata')

寫在文末

我在《生信技能樹》，《生信菜鳥團(tuán)》，《單細(xì)胞天地》的大量推文教程里面共享的代碼都是復(fù)制粘貼即可使用的， 有任何疑問歡迎留言討論，也可以發(fā)郵件給我，詳細(xì)描述你遇到的困難的前因后果給我，我的郵箱地址是 jmzeng1314@163.com

如果你確實(shí)覺得我的教程對你的科研課題有幫助，讓你茅塞頓開，或者說你的課題大量使用我的技能，煩請日后在發(fā)表自己的成果的時候，加上一個簡短的致謝，如下所示：

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

十年后我環(huán)游世界各地的高校以及科研院所（當(dāng)然包括中國大陸）的時候，如果有這樣的情誼，我會優(yōu)先見你。