首先感谢生信技能树大神jmzeng1314提供的github包，由于我这边访问github比较困难，因此我已经导入到我的 gitee 托管平台 https://gitee.com/swcyo，
特别声明：正版托管地址是https://github.com/jmzeng1314
由于GEO官网下载网址的服务器均位于美国，在国内访问下载麻烦，因此，国内大神jmzeng1314开发了“GEOmirror和AnnoProbe两个神包，二者结合起来可迅速下载GEO数据，并且可以直接注释基因。极力推荐
-需要的包为GEOquery、GEOmirror和AnnoProbe
#安装bioconductor包GEOquery
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GEOquery")
#安装github包GEOmirror和AnnoProbe
remotes::install_github("jmzeng1314/GEOmirror")
remotes::install_github("jmzeng1314/AnnoProbe")
#三个包同时加载
library(AnnoProbe)
library(GEOmirror)
library(GEOquery)
- 操作演示，如GSE13507的下载，矩阵、临床信息获得、并基因注释
#下载获取GSE13507数据
gset=AnnoProbe::geoChina('GSE13507')
gset
数据
gset
#检查表达式集
eSet=gset[[1]] # 提取表达矩阵
probes_expr <- exprs(eSet)
dim(probes_expr)
probes_expr=log2(probes_expr+1) #表达矩阵行log2+1归一化处理
# 提取表型数据信息
phenoDat <- pData(eSet)
表达矩阵
log2+1处理
临床表型数据
- 得到的表达矩阵并不是
symbol，因此需要芯片注释，也就是ID转换，平台注释是在线加载的，但也仅需数秒即可完成
## 对表达芯片的探针进行基因注释
(gpl=eSet@annotation)
checkGPL(gpl)
printGPLInfo(gpl)
probe2gene=idmap(gpl)
genes_expr <- filterEM(probes_expr,probe2gene )
![]()
ID转换
ID转换
以为到这就结束了吗？并不是，还可以继续
- limma的经典2组差异分析走一下
## define the group
group_list=factor(c(rep('Normal',68),rep('Cancer',188))) #自己定义分组和数量
##按title自定义，如包含‘bladder’定义为‘cancer’
library(stringr)
group_list=ifelse(str_detect(phenoDat$title,"bladder"),"cancer","normal")
#设置参考水平，对照在前，处理在后
group_list = factor(group_list,
levels = c("normal","cancer"))
table(group_list)
library(limma)
design=model.matrix(~factor(group_list))
design
fit=lmFit(genes_expr,design)
fit=eBayes(fit)
DEG=topTable(fit,coef=2,n=Inf)
head(DEG)
差异基因
- 对差异分析结果进行一些检验
need_deg=data.frame(symbols=rownames(DEG), logFC=DEG$logFC, p=DEG$P.Value)
image.png
火山图走一下，由于没有明显差异，因此没有上下调
deg_volcano(need_deg,1) #第一种图
deg_volcano(need_deg,1,logFC_thred = 0) #设置logFC范围
deg_volcano(need_deg,2) # 第二种图
image.png
image.png
热图走一下
deg_heatmap(DEG,genes_expr,group_list)
deg_heatmap(DEG,genes_expr,group_list,5) #显示前5对，数字可以自己定义
image.png
显示前5对
- boxplot也来一下，一看就是基于ggpurb，当然可以自己后续DIY，加显著性标记
check_diff_genes('RAC3',genes_expr,group_list)
check_diff_genes('RAC3',genes_expr,group_list)+stat_compare_means(method = "t.test") #加个P值，手到擒来
d r
image.png
有了差异基因和logFC值，结合前面Y叔的神包，GO和KEGG富集分析速速的
