使用tinyarray简化geo常规芯片数据分析流程

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

library(tinyarray)
packageVersion("tinyarray")
#> [1] '2.3.2'
library(stringr)
geo = geo_download("GSE7305")
exp = geo$exp
exp = log2(exp+1)
boxplot(exp,las = 2)

pd = geo$pd
gpl_number = geo$gpl
# 分组信息
k = str_detect(pd$title,"Normal");table(k)
#> k
#> FALSE  TRUE 
#>    10    10
Group = ifelse(k,"Normal","Disease")
Group = factor(Group,levels = c("Normal","Disease"))
# 探针注释
find_anno(geo$gpl)
library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)
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
#差异分析和它的可视化
dcp = get_deg_all(exp,Group,ids,entriz = F)
table(dcp$deg$change)
#> 
#>   down stable     up 
#>    579  19621    624
head(dcp$deg)
#>      logFC  AveExpr        t                          P.Value
#> 1 6.270309 8.436140 45.39552 0.000000000000000000000009106509
#> 2 3.943359 7.351799 35.25755 0.000000000000000000002600407155
#> 3 2.318498 6.631187 32.33367 0.000000000000000000017855505829
#> 4 4.905540 8.140399 30.78154 0.000000000000000000053206731115
#> 5 4.878195 6.815838 29.02740 0.000000000000000000195062651758
#> 6 4.106051 9.045949 28.82714 0.000000000000000000227319306208
#>                     adj.P.Val        B    probe_id    symbol change
#> 1 0.0000000000000000002489492 41.58809   202992_at        C7     up
#> 2 0.0000000000000000473924204 37.34483   204971_at      CSTA     up
#> 3 0.0000000000000001952499562 35.77275   228564_at LINC01116     up
#> 4 0.0000000000000004155825748 34.85700 208131_s_at     PTGIS     up
#> 5 0.0000000000000013331313106 33.74579   210002_at     GATA6     up
#> 6 0.0000000000000013809647852 33.61341   212190_at  SERPINE2     up
dcp$plots
library(ggplot2)
ggsave("deg.png",width = 15,height = 5)

#富集分析
deg = get_deg(exp,Group,ids)
genes = deg$ENTREZID[deg$change!="stable"]
head(genes)
#> [1] "730"    "1475"   "375295" "5740"   "2627"   "5270"
#有可能因为网络问题报错
g = quick_enrich(genes,destdir = tempdir())
names(g)
#> [1] "kk"     "go"     "kk.dot" "go.dot"
g[[1]][1:4,1:4]
#>                                      category
#> hsa04610                   Organismal Systems
#> hsa04514 Environmental Information Processing
#> hsa05150                       Human Diseases
#> hsa04110                   Cellular Processes
#>                                  subcategory       ID
#> hsa04610                       Immune system hsa04610
#> hsa04514 Signaling molecules and interaction hsa04514
#> hsa05150       Infectious disease: bacterial hsa05150
#> hsa04110               Cell growth and death hsa04110
#>                                  Description
#> hsa04610 Complement and coagulation cascades
#> hsa04514             Cell adhesion molecules
#> hsa05150     Staphylococcus aureus infection
#> hsa04110                          Cell cycle
library(patchwork)
g[[3]]+g[[4]]

ggsave("enrich.png",width = 12,height = 7)