load('GPRscore_and_group.Rdata')
load('KIRC_tpm.Rdata')

rt = exprSet_tcga_mRNA[,rownames(data)]
group_list=data$group
group_list=factor(group_list,levels = c('Low','High'))

library(limma)
design=model.matrix(~ group_list)

fit=lmFit(rt,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)

gene <- rownames(allDiff)
## 转换
library(clusterProfiler)
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(logFC=allDiff$logFC,
                      SYMBOL = rownames(allDiff))
gene_df <- merge(gene_df,gene,by="SYMBOL")

## geneList 三部曲
## 1.获取基因logFC
geneList <- gene_df$logFC
## 2.命名
names(geneList) = gene_df$ENTREZID
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)

library(clusterProfiler)

## 读入hallmarks gene set，从哪来？
hallmarks <- read.gmt("c2.cp.kegg.v7.5.1.entrez.gmt")
# 需要网络
y <- GSEA(geneList,TERM2GENE =hallmarks)
yd <- data.frame(y)

library(enrichplot)
gseaplot2(y,"KEGG_RENAL_CELL_CARCINOMA",color = "#0072b5",pvalue_table = T)


######################tme
load('TMEscore_and_group.Rdata')
load('KIRC_tpm.Rdata')

rt = exprSet_tcga_mRNA[,rownames(data_immune)]
group_list=data_immune$group
group_list=factor(group_list,levels = c('Low','High'))

library(limma)
design=model.matrix(~ group_list)

fit=lmFit(rt,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)

gene <- rownames(allDiff)
## 转换
library(clusterProfiler)
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(logFC=allDiff$logFC,
                      SYMBOL = rownames(allDiff))
gene_df <- merge(gene_df,gene,by="SYMBOL")

## geneList 三部曲
## 1.获取基因logFC
geneList <- gene_df$logFC
## 2.命名
names(geneList) = gene_df$ENTREZID
## 3.排序很重要
geneList = sort(geneList, decreasing = TRUE)

library(clusterProfiler)

## 读入hallmarks gene set，从哪来？
hallmarks <- read.gmt("c2.cp.kegg.v7.5.1.entrez.gmt")
# 需要网络
y <- GSEA(geneList,TERM2GENE =hallmarks)
yd <- data.frame(y)

library(enrichplot)
gseaplot2(y,"KEGG_RENAL_CELL_CARCINOMA",color = "#0072b5",pvalue_table = T)


###########WGCNA
load('KIRC_tpm.Rdata')
load('GPRscore_and_group.Rdata')
load('TMEscore_and_group.Rdata')
library(WGCNA)

data$GPR_group=ifelse(data$group=='High','GPR_high','GPR_low')
data$TME_group=ifelse(data_immune$group=='High','TME_high','TME_low')
data$group=paste0(data$GPR_group,'+',data$TME_group)
table(data$group)

rt=exprSet_tcga_mRNA[,rownames(data)]
WGCNA_matrix = t(rt[order(apply(rt,1,mad), decreasing = T)[1:5000],])
datExpr0 <- WGCNA_matrix  ## top mad genes
datExpr0 <- as.data.frame(datExpr0)

gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK # 返回TRUE则继续

if (!gsg$allOK){
  
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
 
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

## 样本过滤前
dev.off()
sampleTree = hclust(dist(datExpr0), method = "average")
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)

## 根据图片挑选cutheight,我们不丢样本了!!!!!
clust = cutreeStatic(sampleTree, cutHeight = 145, minSize = 10)
table(clust) # 0代表切除的，1代表保留的
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]

## 更新anno
anno=data[,'group',drop=F]
anno=anno[rownames(datExpr),,drop=F]

## 样本过滤后
sampleTree = hclust(dist(datExpr), method = "average")
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers (after)", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)

## 制作更适合WGCNA的银屑病表型矩阵
anno$GPR_high_TME_low=ifelse(anno$group=='GPR_high+TME_low',1,0)
anno$GPR_high_TME_high=ifelse(anno$group=='GPR_high+TME_high',1,0)

anno$GPR_low_TME_low=ifelse(anno$group=='GPR_low+TME_low',1,0)
anno$GPR_low_TME_high=ifelse(anno$group=='GPR_low+TME_high',1,0)


datTraits =anno[,-1]
datExpr=datExpr[rownames(datTraits),]
sampleNames = rownames(datExpr)
# 能全部对上
traitRows = match(sampleNames, rownames(datTraits))  

###power值散点图
enableWGCNAThreads()   #多线程工作
powers = c(1:30)       #幂指数范围1:20
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5,blockSize = 100000)

dev.off()
par(mfrow = c(1,2))
cex1 = 0.9
###拟合指数与power值散点图
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
abline(h=0.9,col="red") #可以修改
###平均连通性与power值散点图
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

dev.off()
###邻接矩阵转换
sft #查看最佳power值
softPower =sft$powerEstimate #最佳power值
# 发现不合适就自定义，此处我自定义!!!!!
softPower=22
adjacency = adjacency(datExpr, power = softPower)

net = blockwiseModules(datExpr, power = softPower,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "PRTOM",
                       verbose = 3)
# 显示模块数量以及各自包含的基因数目
# 0表示未分入任何模块的基因
# 1是最大的模块，往后依次降序排列，分别对应各自模块的基因
table(net$colors)


mergedColors = labels2colors(net$colors)
mergedColors
##手动保存
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()


moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
# 用color labels重新计算MEs（Module Eigengenes:模块的第一主成分）
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p") #（这是重点）计算ME和表型相关性
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)


# 设置热图上的文字（两行数字：第一行是模块与各种表型的相关系数；
# 第二行是p值）
# signif 取有效数字
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# 然后对moduleTraitCor画热图
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))


##把turquoise挑出来
# 选择导出模块
module = c('turquoise')
# 选择模块中基因/探针
probes = names(datExpr)
inModule = (moduleColors %in% module)
modProbes_POS = probes[inModule]
#modprobes可以后续分析
modProbes_POS
write.table(modProbes_POS,file ='turquoise-high_low_WGCNA.txt',row.names = F,col.names = F,quote=F)

##把green,yellow,red,blue,turquoise负相关的出来

# 选择导出brown模块
module = c('brown')
# 选择模块中基因/探针
probes = names(datExpr)
inModule = (moduleColors %in% module)
modProbes_NEG = probes[inModule]
#modprobes可以后续分析
modProbes_NEG
write.table(modProbes_NEG,file ='low_high_WGCNA.txt',row.names = F,col.names = F,quote=F)

#######################################
#FGSEA#####
gc()
load('KIRC_tpm.Rdata')
load('GPRscore_and_group.Rdata')
load('TMEscore_and_group.Rdata')

data$GPR_group=ifelse(data$group=='High','GPR_high','GPR_low')
data$TME_group=ifelse(data_immune$group=='High','TME_high','TME_low')
data$group=paste0(data$GPR_group,'+',data$TME_group)
table(data$group)

data$group=stringr::str_replace(data$group,pattern = 'GPR_high\\+TME_high',replacement = 'Mixed')
data$group=stringr::str_replace(data$group,pattern = 'GPR_low\\+TME_low',replacement ='Mixed')

data$group
save(data,file ='GPR_TME_combined_group.Rdata')
rt=exprSet_tcga_mRNA[,rownames(data)]

## 导入基因集
#BiocManager::install('fgsea')
library(fgsea)
library(msigdbr)
msigdbr_species()
a=msigdbr_collections()
# 假设做鼠，人就Homo sapiens
m_df<- msigdbr(species = "Homo sapiens", category = "C5",subcategory = 'BP')
# 变list
BP <- m_df %>% split(x = .$gene_symbol, f = .$gs_name)

head(BP)

# 构造预制函数，内置流程是先排序差异基因，针对GO_BP的fgsea
preranked_BP <- function(x) {
  ranks <- x %>% 
    na.omit()%>%
    mutate(ranking=logFC)
  ranks <- ranks$ranking
  names(ranks) <- rownames(x)
  head(ranks, 10)
  set.seed(123456)
  BP_x <- fgsea(pathways = BP, 
                stats = ranks,
                minSize=10,
                maxSize=500,
                nperm=1000)
  
  BP_x$pathway<-gsub("GOBP_","",BP_x$pathway)
  BP_x$pathway<-gsub("_"," ",BP_x$pathway)
  return(BP_x)
}

# 获得各亚型特征的通路
## high_low的marker基因和fgsea
rt = exprSet_tcga_mRNA[,rownames(data)]
group_list=data$group
group_list=ifelse(group_list=='GPR_high+TME_low','GPR_high+TME_low','other')
group_list=factor(group_list,levels = c('other','GPR_high+TME_low'))
library(limma)
design=model.matrix(~ group_list)

fit=lmFit(rt,design)
fit=eBayes(fit) 
hl_allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

# 使用预置函数直接进行fgsea
library(dplyr)
BP_hl <- preranked_BP(hl_allDiff)
sig_BP_hl <- BP_hl %>% filter(abs(NES)>1 & padj<0.05)
sig_BP_hl <- sig_BP_hl[order(sig_BP_hl$NES,decreasing = T),]

## mixed类细胞的marker基因和fgsea
rt = exprSet_tcga_mRNA[,rownames(data)]
group_list=data$group
group_list=ifelse(group_list=='Mixed','Mixed','other')
group_list=factor(group_list,levels = c('other','Mixed'))
library(limma)
design=model.matrix(~ group_list)

fit=lmFit(rt,design)
fit=eBayes(fit) 
mixed_allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

# 使用预置函数直接进行fgsea
library(dplyr)
BP_mixed <- preranked_BP(mixed_allDiff)
sig_BP_mixed <- BP_mixed %>% filter(abs(NES)>1 & padj<0.1)
sig_BP_mixed <- sig_BP_mixed[order(sig_BP_mixed$NES,decreasing = T),]

## lh类细胞的marker基因和fgsea
rt = exprSet_tcga_mRNA[,rownames(data)]
group_list=data$group
group_list=ifelse(group_list=='GPR_low+TME_high','GPR_low+TME_high','other')
group_list=factor(group_list,levels = c('other','GPR_low+TME_high'))
library(limma)
design=model.matrix(~ group_list)

fit=lmFit(rt,design)
fit=eBayes(fit) 
lh_allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

# 使用预置函数直接进行fgsea
library(dplyr)
BP_lh <- preranked_BP(lh_allDiff)
sig_BP_lh <- BP_lh %>% filter(abs(NES)>1 & padj<0.05)
sig_BP_lh <- sig_BP_lh[order(sig_BP_lh$NES,decreasing = T),]


# 合并作图
##A和B合并 suffixes后缀，一个是A，一个是B

merged <- merge(sig_BP_hl[sample(1:nrow(sig_BP_hl),20),c(1,5)], sig_BP_mixed[sample(1:nrow(sig_BP_mixed),20),c(1,5)], by = "pathway" , all = T,
                suffixes = c(".hl",".mixed"))
## 重复合并，合并C
merged2 <- merge(merged,sig_BP_lh[sample(1:nrow(sig_BP_lh),20),c(1,5)], by = "pathway" , all = T,
                 suffixes = c(".hl",".mixed"))

merged3=merged2
colnames(merged3) <- c("pathway","GPR_high+TME_low", "Mixed", "GPR_low+TME_high")
merged3[is.na(merged3)] <- 0
merged3=as.data.frame(merged3)
rownames(merged3) <- merged3$pathway
rownames(merged3)<- tolower(rownames(merged3))
merged3[,1] <- NULL

library(pheatmap)
pheatmap(merged3, cluster_cols = F, cluster_rows = T, border_color=NA, 
         cellwidth =30,  color = colorRampPalette(c(rep("Darkblue",1), "white", rep("red",1)))(1000) ,
         breaks = seq(-2,2,length.out = 1000), main = "Biological process enrichment", 
         angle_col = 45, fontsize=14, fontsize_row = 10)


####TIP数据库
## http://www.360doc.com/content/18/1230/16/19913717_805477264.shtml
load('KIRC_tpm.Rdata')
load('GPR_TME_combined_group.Rdata')
data_hl=data[data$group=='GPR_high+TME_low',]
data_lh=data[data$group=='GPR_low+TME_high',]
data_mix=data[data$group=='Mixed',]

set.seed(123456)
a1=sample(1:nrow(data_hl),10)
set.seed(123456)
a2=sample(1:nrow(data_mix),10)
set.seed(123456)
a3=sample(1:nrow(data_lh),10)

data_hl=data_hl[a1,]
data_mix=data_mix[a2,]
data_lh=data_lh[a3,]

data_tip=rbind(data_hl,data_mix,data_lh)
rt=exprSet_tcga_mRNA[,rownames(data_tip)]
rt=2^(rt)-1
write.table(rt,file ='KIRC_for_TIP.txt',quote = F,sep = '\t',col.names = NA)


### 保存结果到TIP_result
tip=read.table('TIP_result.txt',sep = '\t',header = T,row.names = 1,check.names = F)
identical(colnames(tip),colnames(rt))

library(pheatmap)
anno=data_tip[,'group',drop=F]
pheatmap::pheatmap(tip,annotation_col = anno,
                   cluster_rows = F,cluster_cols = F,show_colnames = F,
                   border_color = NA,color = colorRampPalette(c(rep("Darkblue",1), "white", rep("red",1)))(1000))

pheatmap::pheatmap(tip,annotation_col = anno,
                   cluster_rows = T,cluster_cols = F,show_colnames = F,
                   border_color = NA,color = colorRampPalette(c(rep("Darkblue",1), "white", rep("red",1)))(1000))
