GEO数据挖掘补充(二)——多分组数据

拉布拉多--奶芙

已于 2024-10-26 14:34:53 修改

阅读量684

点赞数 2

分类专栏： GEO数据挖掘文章标签： windows r语言数据分析数据挖掘

于 2024-10-22 17:35:43 首次发布

本文链接：https://blog.youkuaiyun.com/MJ0417/article/details/143117297

版权

GEO数据挖掘专栏收录该内容

6 篇文章

订阅专栏

来自--生信技能树课程

加速下载

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

加载R包,数据获取及整理

library(stringr)
library(tinyarray)
gse = "GSE474"#加入GSE号
geo = geo_download(gse)
exp=geo$exp
range(exp)#变量的变化范围，判断数据是否要取log
exp=log2(exp+1)
boxplot(exp,las=2)#检查数据有无异常

生成分组向量

#对照组放在levels的第一个位置，三分组一般是两两比较，后面差异分析时看清楚谁比谁就可以了。
Group=str_split(geo$pd$title,"-",simplify = T)[,3]
Group=factor(Group,levels = c("NonObese","Obese","MObese"))

探针注释

geo$gpl#输出GPL(平台)号，寻找相应R包进行探针注释
find_anno(geo$gpl)#输出下步操作代码
#> [1]library(hgu133a.db);ids <- toTable(hgu133aSYMBOL)` and `ids <- AnnoProbe::idmap('GPL96')` are both avaliable
library(hgu133a.db)
ids <- toTable(hgu133aSYMBOL)
head(ids)
#>    probe_id symbol
#> 1 1007_s_at   DDR1
#> 2   1053_at   RFC2
#> 3    117_at  HSPA6
#> 4    121_at   PAX8
#> 5 1255_g_at GUCA1A
#> 6   1294_at   UBA7

差异分析

#此处应用limma分析，如果使用tinyarray简化，这段代码就不需要了
library(limma)
design=model.matrix(~0+Group)
colnames(design)=levels(Group)
px = levels(Group)
px
#> [1] "NonObese" "Obese"    "MObese"
contrast.matrix <- makeContrasts(paste0(px[2],"-",px[1]), 
                                 paste0(px[3],"-",px[2]), 
                                 paste0(px[3],"-",px[1]), 
                                 levels=design)
fit <- lmFit(exp, design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
#用lapply分别提取出差异分析结果表格
deg = lapply(1:ncol(design), function(i){
  topTable(fit2,coef = i,number = Inf)
})
names(deg) = colnames(contrast.matrix)
names(deg)
#> [1] "Obese-NonObese"  "MObese-Obese"    "MObese-NonObese"
#Obese-NonObese的意思就是Obese为处理，NoObese为对照的差异分析结果。
head(deg[[1]])
#>                              logFC  AveExpr         t      P.Value adj.P.Val
#> 219601_s_at              0.6071528 4.655855  4.904042 5.554551e-05 0.8001469
#> AFFX-M27830_M_at         1.2562538 9.251799  4.751155 8.162170e-05 0.8001469
#> AFFX-HUMRGE/M10098_3_at  1.4103245 6.962377  4.570331 1.287731e-04 0.8001469
#> 216546_s_at              1.3149112 7.869546  4.433599 1.818218e-04 0.8001469
#> 205605_at                1.2603083 5.866266  4.294402 2.582900e-04 0.8001469
#> 214208_at               -0.7614676 5.430652 -4.146651 3.747303e-04 0.8001469
#>                                 B
#> 219601_s_at             -1.140167
#> AFFX-M27830_M_at        -1.280408
#> AFFX-HUMRGE/M10098_3_at -1.449781
#> 216546_s_at             -1.580239
#> 205605_at               -1.715011
#> 214208_at               -1.860056
#此时得到的deg就是三次差异分析的结果表格组成的列表。
#之后就分别提取出结果画图即可。

TinyArray简化流程

#多分组的数据，get_deg_all仍然可以简化操作，目前是三分组就两两差异分析，四个或五个分组的数据是后面几个组与第一个组差异分析
#暂不支持其他的做法和更多的分组。
dcp = get_deg_all(exp,Group,ids,symmetry = T,logFC_cutoff = 0.585,cluster_cols = F,entriz = F)
dcp$plots
ggplot2::ggsave("deg.png",width = 15,height = 10)

富集分析

#富集分析
#富集分析的输入数据是差异基因名字，只要你提供一组基因名就能做，要分别做三次还是取交集或合集一起做都可以，没有标准答案。
deg=get_deg(exp,Group,ids,logFC_cutoff = 0.585)
Group
g1=deg$'Obese-NonObese'$ENTREZID[deg$'Obese-NonObese'$ENTREZID!="stable"]
g2=deg$'MObese-Obese'$ENTREZID[deg$'MObese-Obese'$ENTREZID!="stable"]
g3=deg$'MObese-NonObese'$ENTREZID[deg$'MObese-NonObese'$ENTREZID!="stable"]

if(F){
  genes = intersect_all(g1,g2,g3)#取交集
}else{
  genes = union_all(g1,g2,g3)#取合集
}

head(genes)

#有可能因为网络问题报错
g = quick_enrich(genes,destdir = tempdir())
names(g)
g[[1]][1:4,1:4]
library(patchwork)
g[[3]]+g[[4]]
ggplot2::ggsave("enrich.png",width = 12,height = 7)