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

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

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就是三次差异分析的结果表格组成的列表。
#之后就分别提取出结果画图即可。