批次效应差异分析，文氏图富集

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("GEOquery")
rm(list=ls())  #清空环境内变量
options(stringsAsFactors = F)  #避免自动将字符串转换为R语言因子
suppressMessages(library(GEOquery))
library(stringr)
library(ggplot2)
library(reshape2)
library(limma)
gset = getGEO('GSE116959', destdir=".", AnnotGPL = F, getGPL = F)  #设置后两个为T比较耗时，而且作用不大
exp<-exprs(gset[[1]])  #exp即为表达矩阵
#使用pData()函数提取临床信息
pdata<-pData(gset[[1]]) 
#临床信息中哪一列提供了分组信息需要自己去鉴别
#采用字符串分隔函数按空格进行分隔取第三位
group_list<-str_split(pdata$characteristics_ch1.4,' ',simplify = T)[,3]

#设置样本分组
group_list = factor(group_list,levels = c("Normal","Tumor"))  #设置为因子变量
table(group_list)
library(limma) 
exp=normalizeBetweenArrays(exp)  #limma中的校正函数

exp_LL=melt(exp)     #melt表示拆分数据
colnames(exp_LL)=c('probe','sample','value') #设置拆分后列名
exp_LL$group=rep(group_list,each=nrow(exp))  #获得其分组
head(exp_LL)
pp=ggplot(exp_LL,aes(x=sample,y=value,fill=group))+  
  geom_boxplot()+
  ggtitle("校正后数据分布")+  #标题
  theme(plot.title = element_text(hjust = 0.5)) +  #标题居中
  #theme(axis.text.x = element_blank()  #隐去x轴文字
  theme(axis.text.x=element_text(angle=90,hjust = 1,size=7)  #调整x轴名称文字属性
  )#绘制箱线图
print(pp)

range(exp) #查看表达矩阵的取值范围，一般而言，范围在20以内的表达值基本已经经过了log对数转换。
library(FactoMineR)
library(factoextra)
dat=as.data.frame(t(exp))
dat.pca <- PCA(dat, graph = FALSE)
fviz_pca_ind(dat.pca,
             geom.ind = "point",
             col.ind = group_list,
             addEllipses = TRUE,
             palette = "jco", 
             legend.title = "Groups"
)
gpl <- getGEO('GPL17077', destdir=".")
colnames(Table(gpl))
head(Table(gpl)[,c(1,7)]) ## 选择自己需要的列
probe2gene=Table(gpl)[,c(1,7)]
probe2gene<-subset(probe2gene,GENE_SYMBOL!='')
ids<-probe2gene
head(probe2gene)
library(devtools)
#install_github("jmzeng1314/idmap3") #安装方法
library(idmap3)
ids=idmap3::get_pipe_IDs('GPL17077')
head(ids) 
library(tidyverse)
exp <- as.data.frame(exp) #也可以直接加到管道操作中
#id匹配
library(dplyr)
exp <- exp %>% 
  mutate(probe_id=rownames(exp)) %>% 
  inner_join(ids,by="probe_id") %>%
  select(probe_id, symbol, everything())

exp <- exp[!duplicated(exp$symbol),]  #去除重复
rownames(exp) <- exp$symbol
exp <- exp[,-(1:2)] #去掉第一列和第二列
exp[1:5,1:4]  #这里得到处理完成后最终的表达矩阵
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)  #提取所有基因的差异分析结果
colnames(deg)
##设置logFC值和P值来标记上下调基因
logFC=1
P.Value = 0.05
#type1 = (deg$P.Value < P.Value)&(deg$logFC < -logFC)
#type2 = (deg$P.Value < P.Value)&(deg$logFC > logFC)
#deg$Group = ifelse(type1,"Down",ifelse(type2,"Up","Not-Sig"))
#table(deg$Group)
# 直接索引
deg2 <- deg[deg$P.Value < P.Value&(abs(deg$logFC) > logFC), ]
write.csv(deg2[0], "deg2.csv")

library(org.Hs.eg.db)
library(AnnotationDbi) 
library(DOSE)
library(clusterProfiler)
gene_symbol=read.csv("deg1.1.csv",header=FALSE)
#install.packages("VennDiagram") 
library(VennDiagram)
venn.plot<-venn.diagram(x=list(前=gene_symbol$V1,后=gene_symbol$V2),
                        filename ="venntu.png", #输出的图片名称，可以选不同的图片格式，如png,tiff等
                        fill=c("cornflowerblue","red"),##圆圈颜色
                        col = "transparent",alpha = 0.50,cex = 1.0,fontface = "bold", ##透明度调整与字体调整等
                        cat.col = c("darkblue", "darkgreen"), ##字体颜色调整
                        cat.cex = 1.5,cat.dist = 0.07,margin = 0.2) ##字体大小调整

entrez_id = mapIds(x = org.Hs.eg.db, keys = gene_symbol$V1, keytype = "SYMBOL", column = "ENTREZID")
result2=enrichDO(gene = entrez_id,ont = "DO",pvalueCutoff = 0.05,qvalueCutoff = 0.05)
barplot(result2)
cnetplot(result2,categorySize="pvalue")