0.准备工作-下载R包（待修改）

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F，options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'patchwork') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "KEGG.db",
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot",
                         "ggplotify")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#没有error就是成功！

#哪个报错，就回去安装哪个。如果你没有安装xx包，却提示你xx包不存在，这也正常，是因为复杂的依赖关系，缺啥补啥。
if(!require(AnnoProbe))devtools::install_local("./AnnoProbe-master.zip",upgrade = F)
library(AnnoProbe)

1.数据下载

rm(list = ls())
library(GEOquery)
gse_number = "GSE56649"
eSet <- getGEO(gse_number, 
               destdir = '.', #下载后的数据放在工作目录下(不用改),".."代表工作目录的上级文件夹
               getGPL = F)#一起下载数据时，把GPL的大表格也下载下来，时间巨长（不用改）
class(eSet)
length(eSet)
eSet = eSet[[1]]
#1.提取表达矩阵exp
exp <- exprs(eSet)#提取表达矩阵
exp[1:4,1:4]#结果都是几百，说明没取对数
exp = log2(exp+1)#如果结果是1-15之间的数值，不用取对数，如果数值很大，那就取；为了防止exp为0，所以+1
boxplot(exp)#箱线图查看样本的分布情况
#2.提取临床信息
pd <- pData(eSet)#通常情况下，表达矩阵的行名和PD的列名是一一对应的
dim(exp)#查看表达矩阵exp的具体内容，核对并确认样本数量及样本顺序都是相同的
#3.调整pd的行名顺序与exp列名完全一致#（都是GSM的编号）
p = identical(rownames(pd),colnames(exp));p#判断pd的rownames与exp的colnames是否一致
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]#如果不一致的话变得一致
#4.提取芯片平台编号
eSet@annotation#对象除了用$之外，还可以用@取子集
gpl_number <- eSet@annotation
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")#将数据保存至Rdata

2.Group(实验分组)和ids(探针注释)

rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
View(pd)#查看pd的内容，挑选可以作为分组信息的列
# 1.Group----表示实验分组的因子及探针注释（对应那些基因）的数据框
# 第一类，有现成的可以用来分组的列
if(F) Group = pd$`disease state:ch1` 

#第二类，自己生成
if(F){
  Group=c(rep("RA",times=13),
          rep("control",times=9))
}
rep(c("RA","control"),times = c(13,9))#有风险，13和9都是自己数出来的

#*第三类，匹配关键词，自行分类*
Group=ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA")
##找到了只显示分组情况的一列，提取关键字并命名对照组和实验组，注意一定让control在前！

#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,
               levels = c("control","RA"))
Group

# 注意levels与因子内容必须对应一致
# Group = pd$`disease state:ch1`
# Group = factor(Group,
#                 levels = c("healthy control","rheumatoid arthritis"))

#2.ids-----------------
#方法1 BioconductorR包(最常用)
gpl_number 
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# 方法2 读取GPL平台的soft文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：soft文件列名不统一，活学活用，有的GPL平台没有提供注释，如GPL16956
  a = getGEO(gpl_number,destdir = ".")
  b = a@dataTable@table#选取soft文件中带有Gene symbol的表格命名为b
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")] #提取b中"ID","Gene Symbol"两列组成"data.frame"
  colnames(ids2) = c("probe_id","symbol")# 将ids2的列名更改为"probe_id","symbol"
  ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]
}#反选了ids2中symbol一列无内容即有内容的元素，以及不含///的元素，清楚无效数据

# 方法3 官网下载,文件读取
##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus

# 方法4 自主注释 
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

3.PCA图和热图

rm(list = ls())  
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")#输入数据：exp和Group
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

# 1.PCA 图----
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra) 
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = Group, # color by groups
                         palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot
ggsave(plot = pca_plot,filename = paste0(gse_number,"_PCA.png"))
save(pca_plot,file = "pca_plot.Rdata")

# 2.top 1000 sd 热图---- 
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]

# 直接画热图，对比不鲜明
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col
)

# 用标准化的数据画热图，两种方法的比较：https://mp.weixin.qq.com/s/jW59ujbmsKcZ2_CM5qRuAg

## 1.使用热图参数
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",#由于我们需要比较的是一个基因在不同样本中的表达变化，选择按行标准化，比如标准化后数据范围由0-15缩至-4到4
         breaks = seq(-3,3,length.out = 100)#设置范围，比如当其中某几个数值太大或太小导致颜色极化，设置范围后超过范围的数将以极限值的颜色展示出来。比如1000用3的颜色，即3的颜色代表了大于等于3的数值。
         ) #breaks 参数解读在上面链接
dev.off()

## 2.自行标准化再画热图
n2 = t(scale(t(n)))#scale强化列与列之间的比较，弱化行与行之间的比较
pheatmap(n2,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
)
dev.off()

# 关于scale的进一步探索：zz.scale.R

# 3.相关性热图---- 
#直接用全部的表达矩阵exp找两个基因之间的相关性
pheatmap::pheatmap(cor(exp),
                   annotation_col = annotation_col)

#用差异较大的top1000个基因n找相关性
pheatmap::pheatmap(cor(n),
                   annotation_col = annotation_col)

#用top1000个基因再标准化后的n2找相关性
pheatmap::pheatmap(cor(n2),
                   annotation_col = annotation_col
                   )

dev.off()

# 关于相关性背后的故事：https://mp.weixin.qq.com/s/IqMW6Qjf64dn30F4RQg5kQ

04.提取差异基因和ID转换

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)#提取数据fit形成差异分析后的表达矩阵

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))#新增一列：把DEG的行名变成一列，增加probe_id
head(deg)
#2.加上探针注释
#当一个探针对应多个基因时称为非特异性探针，Bioconductor中自动去除了
#而当多个探针对应一个基因的时候，按照基因去重复，一个基因保留一个探针即可。
table(!duplicated(ids$probe_id))
table(!duplicated(ids$symbol))#去掉了第二次至睇多次出现的同一探针对应不同基因的基因，即只保留第一个
#按symbol列去重，常见标准有3个：最大值/平均值/随机去重
#随机去重，另两个见 zz.去重方式.R
ids = ids[!duplicated(ids$symbol),]

deg <- inner_join(deg,ids,by="probe_id")#把deg和ids连接在一起
head(deg)
nrow(deg)

#3.加change列,标记上下调基因
logFC_t=1
P.Value_t = 0.01
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)#下调基因的条件
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)#上调基因的条件
change = ifelse(k1,"down",ifelse(k2,"up","stable"))
deg <- mutate(deg,change)#将判断结果加入表格中
table(deg$change)#查看具体差异基因个数

#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)#ID转换的R包，Hs是人的意思，其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类

dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
length(unique(deg$symbol))
save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

05.富集分析

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

# 1.GO 富集分析----

#(1)输入数据
#head(deg)
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']

#(2)富集
#以下步骤耗时很长，设置了存在即跳过
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)#gene ID自动转换成gene symbol
  #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,file = paste0(gse_number,"_GO.Rdata"))
}
load(paste0(gse_number,"_GO.Rdata"))
#class(ego)
#z=ego@result;z

#(3)可视化
#条带图
barplot(ego)
#气泡图
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + facet_grid(ONTOLOGY ~ ., scale = "free") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))#当x轴太长时设置了折叠

#geneList 用于设置下面图的颜色
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)

#(3)展示top通路的共同基因，要放大看。
#Gene-Concept Network
cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map,这个函数最近更新过，版本不同代码会不同

Biobase::package.version("enrichplot")

if(F){
  emapplot(pairwise_termsim(ego)) #新版本
}else{
  emapplot(ego)#老版本
}
#(4)展示通路关系 https://zhuanlan.zhihu.com/p/99789859
#goplot(ego)

#(5)Heatmap-like functional classification
heatplot(ego,foldChange = geneList,showCategory = 8)

# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#（1）输入数据
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#（2）对上调/下调/所有差异基因进行富集分析

if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa')
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')
  save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
}
load(paste0(gse_number,"_KEGG.Rdata"))

#(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)#结果为FAUSE即没有富集的结果，不要怀疑自己
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)

#(4)p.adjust不适用可以按照pvalue筛选通路
down_kegg <- kk.down@result %>%
  filter(pvalue<0.05) %>% #筛选行
  mutate(group=-1) #新增列

up_kegg <- kk.up@result %>%
  filter(pvalue<0.05) %>%
  mutate(group=1)

#(5)可视化
source("kegg_plot_function.R")#source是不打开脚本的情况下运行代码
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4))#改横坐标轴，应全为正值
ggsave(g_kegg,filename = 'kegg_up_down.png')

# 3.gsea作kegg和GO富集分析----
## https://www.jianshu.com/p/c5b7b7dbf29b
#GSEA是把全部基因进行富集，在GO和KEGG无法富集到的时候可以选择
#(1)查看示例数据
data(geneList, package="DOSE")
#(2)将我们的数据转换成示例数据的格式
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
#(3)富集分析
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  verbose      = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可视化
g2 = kegg_plot(up_kegg,down_kegg)
g2

# 4.能看懂的资料越来越多----
# GSEA学习更多：https://www.jianshu.com/p/baf85b51752e
# 富集分析学习更多：http://yulab-smu.top/clusterProfiler-book/index.html
# 弦图：https://www.jianshu.com/p/e4bb41865b7f

06.

07.