GEO数据挖掘-基于芯片

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor") 
cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'survminer',
                   "tinyarray") 
Biocductor_packages <- c('GEOquery',
                         'GO.db',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot")

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


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T,quietly = 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) 
}
#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。
#library报错，就单独安装。

require(ggplot2)  # 加载ggplot2包

package_name <- "ggplot2"
require(package_name, character.only = TRUE)  # 加载ggplot2包

rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000) 
options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
#研究一下这个eSet
class(eSet)
length(eSet)

eSet = eSet[[1]] 
class(eSet)

#(1)提取表达矩阵exp
exp <- exprs(eSet)
#⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数，
#如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析
dim(exp)
#⭐二个要检查的地方👇
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#⭐可能要修改的地方👇
exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句
#⭐第三个要检查的地方👇
boxplot(exp,las = 2) #看是否有异常样本

#(2)提取临床信息
pd <- pData(eSet)
#⭐多分组中提取两分组的代码示例，二分组不需要
if(F){
  #因为现在这个例子不是多分组，所以编造一列做示例。
  pd$fake = paste0(rep(c("a","b","c","d"),each = 5),1:5)
  k1 = str_detect(pd$fake,"b");table(k1)
  k2 = str_detect(pd$fake,"c");table(k2)
  pd = pd[k1|k2,]
}
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

# 原始数据处理的代码，按需学习
# https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw

if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
#⭐要修改的地方：分组信息,必须学会ifelse和str_detect
k = str_detect(pd$title,"Normal");table(k) #不在title就在pd的其他列
Group = ifelse(k,"Normal","Disease")

# 需要把Group转换成因子，并设置参考水平，指定levels
#⭐要修改的地方，对照组在前，处理组在后
Group = factor(Group,levels = c("Normal","Disease"))
Group

#⭐检查自己得到的分组是否正确
data.frame(pd$title,Group)
#2.探针注释的获取-----------------
#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
#⭐要操作的地方
library(tinyarray)
gpl_number #首先看看编号是多少
#View(pkg_all) 
#然后在pkg_all里搜索gpl编号，找到对应的R包前缀(第二列)，没搜到就是没有R包,再看方法2。
#也可以用代码直接得到对应的R包前缀：
pkg_all[pkg_all$gpl==gpl_number,2]
#⭐要操作的地方
#用上面找到的前缀替换下面所有的hgu133plus2，共5处
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db",ask = F,update = F)
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出R包里都有啥
ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
# 方法2 下载并读取GPL网页的表格文件，按列取子集
#⭐要操作的地方
library(tinyarray)
get_gpl_txt(gpl_number) #获取表格文件的下载链接
# 接下来是复制网址去浏览器下载、放在工作目录下、读取、提取探针id和基因symbol(没有现成的需要拆分和转换)，不同文件代码不统一，等看同学们的例子。
# 注意:最终的数据ids只能有两列，第一列列名是probe_id,第二列列名是symbol,且都是字符型，否则后面代码要报错咯。
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,file = "step2output.Rdata")

#比较复杂的探针注释参考资料
#资料1：拆分取列https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5
#资料2：多种id的转换https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/pn0s1mmsaxocfynb?singleDoc# 《又一个有点难的探针注释（多种id的转换）》
#资料3：其他id转换 https://www.jianshu.com/p/f4e799f06b52

rm(list = ls())  
#⭐不用改代码，会看结果就行
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)
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"
)
#如果少于4个点就不会画出圈，是正常的。因为圈是置信区间，样本太少无法计算，不是必须的。
# 2.top 1000 sd 热图---- 
g = names(tail(sort(apply(exp,1,sd)),1000)) #day7-apply的思考题
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
                            Group = Group)
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row", #按行标准化，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
         breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，超出此范围的数字显示极限颜色
         ) 

# 关于scale的进一步学习：zz.scale.R

rm(list = ls()) 
# 加载之前保存的数据，包括 exp（表达矩阵）、Group（样本分组）和 ids（探针注释）。
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#⭐检查
nrow(deg) #如果行数为0就是你找的探针注释是错的。

#3.加change列,标记上下调基因
#⭐阈值，可按需修改
logFC_t = 1
p_t = 0.05
#⭐思考，如何使用padj而非p值
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#火山图
library(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()

# 差异基因热图----
# 表达矩阵行名替换为基因名
exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,]
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n) 
pheatmap(n,show_colnames =F,
         show_rownames = F,
         scale = "row",
         #cluster_cols = F, 
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
) 

#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol, 
           fromType = "SYMBOL",
           toType = "ENTREZID",
           OrgDb = org.Hs.eg.db)#⭐Hs是人类,注意物种
#一部分基因没匹配上是正常的。只要物种正确，<30%的失败都没事。
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
#⭐检查行数，减少太多说明不正常
nrow(deg) 
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#多了一些行少了一些行都正常，SYMBOL与ENTREZID不是一对一的。
nrow(deg)
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

#(1)输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"] 

#(2)富集
#⭐下面的三句都要注意物种
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
#其他物种https://www.genome.jp/kegg/catalog/org_list.html
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")
#如果ekk是空的，这句就会报错，因为没富集到任何通路。
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
                ont = "ALL",readable = TRUE)
#setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)

#(3)可视化
#barplot可以换成dotplot
barplot(ego, split = "ONTOLOGY") + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
barplot(ekk)
#如果ekk中没有padj<0.05的通路，就会报错，因为只画padj<0.05,没有参数


# 更多资料---
# GSEA：https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
# Y叔的书：http://yulab-smu.top/clusterProfiler-book/index.html
# GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
# 网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~

ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")