單個基因集的GSEA分析
http://software.broadinstitute.org/gsea/msigdb/index.jsp
GSEA官網(wǎng)列出關(guān)于msigdb基因集的詳細(xì)信息
C1是通過染色體位置進(jìn)行富集分析
load('./CORO2A/CORO2A_up_and_down_gene.Rdata')
load(file = './Rdata/DEGs_results_CORO2A_Si2_vs_NC.Rdata')
gene_diff<- c(up_gene,Down_gene)
# DEG<- DEG[rownames(DEG)%in%up_gene,]
# DEG<- DEG[rownames(DEG)%in%Down_gene,]
geneList=DEG$log2FoldChange
names(geneList)=rownames(DEG)
geneList=sort(geneList,decreasing = T)
#選擇gmt文件（MigDB中的全部基因集）
d='D:\\GSVA_and_GSEA_Gene_Sets/MsigDB/symbols/' ### 特別注意這里是兩杠，第一根\的意思就是轉(zhuǎn)字符？
gmts <- list.files(d,pattern = 'all')
gmts
#### 第二：關(guān)于GSEA的理解
library(GSEABase) # BiocManager::install('GSEABase')
## 下面使用lapply循環(huán)讀取每個gmt文件，并且進(jìn)行GSEA分析
gmtfile=gmts[1]
geneset <- read.gmt(file.path(d,gmtfile))
print(paste0('Now process the ',gmtfile))
egmt <- GSEA(geneList,
TERM2GENE=geneset,
pvalueCutoff = 0.9,
verbose=FALSE)
head(egmt)
gseaplot(egmt, geneSetID = rownames(egmt[1,]))
多個基因集和GSEA和GSVA分析
第一：要下載相應(yīng)的GSEA基因集
第二：關(guān)于GSEA 理解輸入輸出
1.需要輸入一個geneList選取差異表達(dá)的基因
2.需要輸入gmts數(shù)據(jù)集（從GSEA官網(wǎng)下載）：
d='./MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
gmts3.lapply函數(shù)是對gmts基因集的循環(huán)奈籽，如果下載全基因集或者只關(guān)注一個基因集，那就不需要循環(huán)：如注釋信息 # gmtfile=gmts[2]:相當(dāng)于只取第二個基因集
- geneset <- read.gmt(file.path(d,gmtfile)) ：用read.gmt 這個函數(shù)讀取gmt文檔信息
- egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
? head(egmt)：相當(dāng)于還是主函數(shù)胞得。主要的數(shù)據(jù)分析gseaplot(egmt, geneSetID = rownames(egmt[1,])) ：選取了一條基因集畫圖：選取的基因集的名字是rownames(egmt[1,])
7： gsea_results_list<- lapply(gsea_results, function(x){
? cat(paste(dim(x@result)),'\n')
? x@result
}) ##########這個的意思應(yīng)該是吧gsea_results進(jìn)行了操作列出 gsea_results_list8：gsea_results_df <- do.call(rbind, gsea_results_list) ###########把gsea_results_list做轉(zhuǎn)化成data.frame慢洋。do.call神奇操作
9：gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE') 垮抗。選擇有差異的基因集進(jìn)行畫圖刚陡。此處要注意gsea_results[[2]]其中的數(shù)字2 要和KEGG_CELL_CYCLE是對應(yīng)的惩妇。所以畫圖前一定要用notepad++查看下載下來的基因集和數(shù)字的對應(yīng)關(guān)系。
load(file = 'anno_DEG.Rdata')
geneList=DEG$logFC
names(geneList)=DEG$symbol
geneList=sort(geneList,decreasing = T)
#選擇gmt文件（MigDB中的全部基因集）
d='./MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
gmts
#### 第二：關(guān)于GSEA的理解
library(GSEABase) # BiocManager::install('GSEABase')
## 下面使用lapply循環(huán)讀取每個gmt文件筐乳，并且進(jìn)行GSEA分析
gsea_results <- lapply(gmts, function(gmtfile){
# gmtfile=gmts[2]
geneset <- read.gmt(file.path(d,gmtfile))
print(paste0('Now process the ',gmtfile))
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
head(egmt)
# gseaplot(egmt, geneSetID = rownames(egmt[1,]))
return(egmt)
})
# 上面的代碼耗時屿附，所以保存結(jié)果到本地文件
save(gsea_results,file = 'gsea_results.Rdata')
#提取gsea結(jié)果，熟悉這個對象
gsea_results_list<- lapply(gsea_results, function(x){
cat(paste(dim(x@result)),'\n')
x@result
})
gsea_results_df <- do.call(rbind, gsea_results_list)
gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE')
gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6')
gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE')
}
第三：GSVA的輸入輸出
計算每個樣本在基因集中的富集程度哥童，然后計算Fold change 和P value。做類似于熱圖的圖褒翰。所以后面的很多算法和制作熱圖一致贮懈。
關(guān)注input：函數(shù)需要輸入X=dat這個expression data，row為基因名优训，col為sample 朵你。另外還需要gmt這個數(shù)據(jù)集。
d='./MSigDB/symbols/'
gmts=list.files(d,pattern = 'all') 通過list.files這個可以列出當(dāng)前路徑下的所有g(shù)mt文件es_max <- lapply(gmts, function(gmtfile) 揣非；使用lapply函數(shù)進(jìn)行循環(huán)抡医，對內(nèi)一個下載的基因集進(jìn)行操作；如果我們只關(guān)心一個基因集早敬，就不需要這個循環(huán)了比如 #gmtfile=gmts[8]
geneset <- getGmt(file.path(d,gmtfile))
? es.max <- gsva(X, geneset,
? mx.diff=FALSE, verbose=FALSE,
? parallel.sz=1)#######這一步很重要忌傻，是函數(shù)的最重要步驟。geneset <- getGmt(file.path(d,gmtfile)) 通過這個getGmt這個函數(shù)可以從gmt 文件中獲得geneset搞监；gsva這個主函數(shù)通過輸入X這個表達(dá)矩陣和geneset數(shù)據(jù)集進(jìn)行計算Fold change 和P value水孩。
5adjPvalueCutoff <- 0.001
? logFCcutoff <- log2(2)####進(jìn)行閾值調(diào)整6.帥選差異表達(dá)的基因集##類似做熱圖和火山圖（因為實驗設(shè)計是本身存在不同的分組信息）；
6.1
es_deg <- lapply(es_max, function(es.max)這個函數(shù)對所有的基因集得到的es.max進(jìn)行循環(huán)操作
6.2
分組內(nèi)容如下：design <- model.matrix(~0+factor(group_list))
? colnames(design)=levels(factor(group_list))
? rownames(design)=colnames(es.max)6.3 deg = function(es.max,design,contrast.matrix) 這個函數(shù)對es.max按照分組信息進(jìn)行基因集富集情況進(jìn)行差異分析（類似于尋找差異表達(dá)基因）
7.最終會輸出ex.max 和ex.deg這兩個數(shù)據(jù)
8.最后作圖琐驴，做差異富集的圖俘种，（類似差異表達(dá)的基因）作圖的時候回報錯秤标，我標(biāo)注在代碼上：詳見代碼注釋
### 對 MigDB中的全部基因集 做GSVA分析。
## 還有ssGSEA, PGSEA
{
load(file = 'step1-output.Rdata')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dat[1:4,1:4]
ids$median=apply(dat,1,median)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,]
rownames(dat)=ids$symbol
dat[1:4,1:4]
X=dat
table(group_list)
## Molecular Signatures Database (MSigDb)
d='./MSigDB/symbols/'
gmts=list.files(d,pattern = 'all')
gmts
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSVA) # BiocManager::install('GSVA')
if(T){
es_max <- lapply(gmts, function(gmtfile){
#gmtfile=gmts[8];gmtfile
geneset <- getGmt(file.path(d,gmtfile))
es.max <- gsva(X, geneset,
mx.diff=FALSE, verbose=FALSE,
parallel.sz=1)
return(es.max)
})
adjPvalueCutoff <- 0.001
logFCcutoff <- log2(2)
es_deg <- lapply(es_max, function(es.max){
table(group_list)
dim(es.max)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(es.max)
design
library(limma)
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),
levels = design)
contrast.matrix<-makeContrasts("TNBC-noTNBC",
levels = design)
contrast.matrix ##這個矩陣聲明宙刘，我們要把progres.組跟stable進(jìn)行差異分析比較
deg = function(es.max,design,contrast.matrix){
##step1
fit <- lmFit(es.max,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
##這一步很重要苍姜，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
res <- decideTests(fit2, p.value=adjPvalueCutoff)
summary(res)
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
re = deg(es.max,design,contrast.matrix)
nrDEG=re
head(nrDEG)
return(nrDEG)
})
}
gmts
save(es_max,es_deg,file='gsva_msigdb.Rdata')
load(file='gsva_msigdb.Rdata')
library(pheatmap)
lapply(1:length(es_deg), function(i){
# i=2
print(i)
dat=es_max[[i]]
df=es_deg[[i]]
df=df[df$P.Value<0.01 & abs(df$logFC) > 0.3,] #######在這里調(diào)閾值，保證輸出
print(dim(df))
if(nrow(df)>5){ #######nrow(df)>5:只有差異富集的基因集大于五個以上才畫圖
n=rownames(df)
dat=dat[match(n,rownames(dat)),] ######n=rownames(df),篩選表達(dá)矩陣數(shù)據(jù)
ac=data.frame(g=group_list) ########對應(yīng)后面的圖中的annotation_col
rownames(ac)=colnames(dat)
rownames(dat)=substring(rownames(dat),1,50)
pheatmap::pheatmap(dat,
fontsize_row = 8,height = 11,
annotation_col = ac,show_colnames = F,
filename = paste0('gsva_',strsplit(gmts[i],'[.]')[[1]][1],'.pdf'))
}
})
adjPvalueCutoff <- 0.001
logFCcutoff <- log2(2)
df=do.call(rbind ,es_deg)
es_matrix=do.call(rbind ,es_max)
df=df[df$P.Value<0.01 & abs(df$logFC) > 0.5,]
write.csv(df,file = 'GSVA_DEG.csv')
}
