admin管理员组

文章数量:1487745

GEO数据挖掘补充(一)——TinyArray简化流程

GEO数据挖掘补充—TinyArray简化流程

来自生信技能树课程

打破下载时间限制

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

R包的加载和数据的下载

代码语言:txt复制
library(tinyarray)
packageVersion("tinyarray")#查看数据包版本信息
library(stringr)
gse = "GSE56649"#填写所要处理的GSE号
geo = geo_download(gse)
exp=geo$exp#从名为 geo 的对象中提取名为exp的组件,并将提取的组件赋值给一个新的变量exp。
exp[1:4,1:4]#自行判断是否需要进行log
geo$exp = log2(geo$exp+1)
boxplot(geo$exp)#检查数据是否异常

对数据进行分组

代码语言:txt复制
Group=ifelse(str_detect(geo$pd$source_name_ch1,"control"),
             "control",
             "RA")
Group = factor(Group,levels = c("control","RA"))

探针注释

代码语言:txt复制
find_anno(geo$gpl)#输出下步操作代码
#> [1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"
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

差异分析及其可视化

代码语言:txt复制
dcp = get_deg_all(geo$exp,Group,ids)
table(dcp$deg$change)
#> 
#>   down stable     up 
#>    828  19339    657
head(dcp$deg)
#>       logFC   AveExpr         t      P.Value    adj.P.Val        B    probe_id
#> 1 -1.633518 11.204696 -16.14895 8.628803e-14 2.358899e-09 20.92048 208629_s_at
#> 2 -1.713884 11.965832 -15.06767 3.572073e-13 4.882577e-09 19.68466 201861_s_at
#> 3 -2.274611  7.651265 -13.43343 3.618635e-12 3.297481e-08 17.61994 203479_s_at
#> 4 -1.935657  8.628444 -13.04331 6.501205e-12 3.554534e-08 17.08858   203391_at
#> 5 -1.893808  8.896532 -12.38173 1.812689e-11 6.607251e-08 16.15073   205457_at
#> 6 -1.972646  8.437584 -11.83695 4.353131e-11 1.328288e-07 15.34203 202354_s_at
#>    symbol change ENTREZID
#> 1   HADHA   down     3030
#> 2 LRRFIP1   down     9208
#> 3   OTUD4   down    54726
#> 4   FKBP2   down     2286
#> 5   ILRUN   down    64771
#> 6  GTF2F1   down     2962
dcp$plots
library(ggplot2)
ggsave("deg.png",width = 15,height = 5)

富集分析

代码语言:txt复制
genes=dcp$deg$ENTREZID[dcp$deg$change!="stable"]
head(genes)
#> [1] "3030"  "9208"  "54726" "2286"  "64771" "2962"

g=quick_enrich(genes,destdir=tempdir())
#quick_enrich()快速的富集分析,此步可能因为网络问题报错

names(g)
#> [1] "kk"     "go"     "kk.dot" "go.dot"
g[[1]][1:4,1:4]
#>                ID                  Description GeneRatio  BgRatio
#> hsa05322 hsa05322 Systemic lupus erythematosus    24/568 136/8217
#> hsa05034 hsa05034                   Alcoholism    27/568 187/8217
#> NA           <NA>                         <NA>      <NA>     <NA>
#> NA.1         <NA>                         <NA>      <NA>     <NA>
library(patchwork)
g[[3]]+g[[4]]
ggsave("enrich.png",width = 12,height = 7)

本文标签: GEO数据挖掘补充(一)TinyArray简化流程