第 32 课- GSVA 分析

#1. 读取标准化后的表达矩阵
input <- read.table("tpms_mrna_tumor_normal.txt", sep="\\t",
                    header=T, row.names = 1, check.names=F)
exp <- as.matrix(input)
dim(exp)

#2. R包与文件准备
##2.1 加载R包
library(GSEABase)
library(GSVA)

##2.2 读取基因列表文件
msigdb_GMTs <- "msigdb_v7.0_GMTs"
msigdb <- "msigdb.v7.0.symbols.gmt"

# 获取gmt文件内容
geneset <- getGmt(file.path(msigdb_GMTs, msigdb))

#3. 进行GSVA分析（使用1线程）
t1=proc.time()
es.max <- gsva(exp,
               geneset,
               mx.diff=FALSE, verbose=TRUE,
               parallel.sz=1) # CPU功率稳定在18W左右，还有时会在20W
## 平静时，我电脑的CPU功率是7W
t2=proc.time()
t=t2-t1
print(paste0("执行时间：",t[3]相关主题：1,"秒"))
# [1] "执行时间：3146.29000000001秒",也就是52mins,说明还是要慢很多的，如果是单线程的话~

#3. 进行GSVA分析 （使用16线程）
t1=proc.time()
es.max <- gsva(exp,
               geneset,
               mx.diff=FALSE, verbose=TRUE,
               parallel.sz=1) # CPU功率稳定在14W左右
## 平静时，我电脑的CPU功率是7W
t2=proc.time()
t=t2-t1
print(paste0("执行时间：",t[3]相关主题：1,"秒"))
# [1] "执行时间：2028.75秒",也就是半个多小时

##保存结果
write.table(es.max, file="TCGA.GSVA_MsigDb_result.txt",sep="\\t",
            quote=F,row.names = T)
save(es.max, file = "TCGA.GSVA_MsigDb_result.Rda")
# load(file = "TCGA.GSVA_MsigDb_result.Rda")
#4. GSVA差异分析
##4.1 分组信息
group_list <- ifelse(substr(colnames(es.max), 14, 16) == "01A", "Tumor",
                     "Normal")
group_list <- factor(group_list, levels = c("Normal", "Tumor"))
table(group_list)

##4.2 差异分析
library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(es.max)

contrast.matrix<-makeContrasts("Tumor-Normal",
                               levels = design)
contrast.matrix

###拟合limma线性模型
fit <- lmFit(es.max,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
###查看是否差异显著
res <- decideTests(fit2, p.value=0.05)
summary(res)

###提取差异分析结果
tmp <- topTable(fit2)
DEG <- na.omit(tmp)
head(DEG[1:4,])

###将差异分析结果保存到文件
write.table(DEG, file="TCGA.GSVA_DEG_limma_result.txt",sep="\\t",
            quote=F,row.names = T)
save(DEG, group_list, file = "TCGA.GSVA_DEG_limma.Rda")

#5. 可视化
# 绘制差异分析热图
rm(list = ls())
load(file = "TCGA.GSVA_DEG_limma.Rda")
load(file = "TCGA.GSVA_MsigDb_result.Rda")

DEG_sig <- DEG[DEG$P.Value<0.01 & abs(DEG$logFC) > 0.5,]
dat <- es.max[match(rownames(DEG_sig),rownames(es.max)),]

annotation_col <- data.frame(group_list)
rownames(annotation_col) <- colnames(dat)
library(pheatmap)
pheatmap(dat,
         annotation_col = annotation_col,
         show_colnames = F,
         scale = "row",
         cluster_cols = F,
         )