我写了一个R包，简化芯片的差异分析

1.由来

近期由于工作需要，有大批的芯片数据等着我分析。我想着简化一下代码，一步到位出来差异分析结果。配合练习写R包，今天算是搞定了大头，分享一下给有缘人使用~

目前差异分析仅支持二分组数据，多分组的后面再说~

2.R包安装和准备

我的包托管在Github上，并且依赖了曾老板写的AnnoProbe包，他的包也在github，所以没办法自动帮你安装，需要先安装好哦。

/分割的是用户名和包名，知道了用户名，你就可以在github上搜索到包对应的页面啦。

if(!require(AnnoProbe))devtools::install_github("jmzeng1314/AnnoProbe")
if(!require(tinyarray))devtools::install_github("xjsun1221/tinyarray")

github的包安装需要多折腾，实在折腾不了就放弃吧。如果仅仅是网络问题，可以使用本地安装，到github上下载这个包，下来放在工作目录下，然后：

#devtools::install_local("tinyarray-master.zip",upgrade = F)

只要有一些R包安装的基础知识和解决报错能力，就可以搞定啦。

library(stringr)
library(AnnoProbe)
library(tinyarray)

3.GEO数据下载

rm(list = ls())
options(stringsAsFactors = F)
gse = "GSE4107"
geo = geo_download(gse)

写成了一个函数geo_download，默认使用AnnoProbe下载，如果报错就说明这个GSE没有被AnnoProbe收录，可以

geo = geo_download(gse,by_annopbrobe = F)

就会使用GEOquery下载了，如果网速实在堪忧，可以去GEO下载matrix.gz文件，放在工作目录下。

4.得到表达矩阵和分组信息

geo$exp = log(geo$exp+1)
group_list = ifelse(str_detect(geo$pd$title,"control"),"control","treat")
group_list = factor(group_list,levels = c("control","treat"))

5.找探针注释信息

结合两个来源，一个是bioconductor的包，一个是idmap。find_anno会返回可用的命令，复制下来运行即可。

find_anno(geo$gpl)

## [1] "`ids <- toTable(hgu133plus2SYMBOL)` or `ids <- AnnoProbe::idmap('GPL570')` is both avaliable"

ids <- AnnoProbe::idmap('GPL570')

6. 完成差异分析及可视化

把很多代码集成到了一起，得到的dcp是一个列表里面包括了差异分析结果表格，差异基因以及三张图。

dcp = get_deg_all(geo$exp,
                  group_list,
                  ids,
                  logFC_cutoff = 1,
                  scale_before = T,
                  cluster_cols = F)

## [1] "362 down genes,1361 up genes"

head(dcp[[1]])

##      logFC  AveExpr         t      P.Value    adj.P.Val        B    probe_id
## 1 4.755775 5.433371 15.686305 9.489229e-14 2.091036e-09 19.58011   202768_at
## 2 3.617970 7.889055 15.545324 1.147345e-13 2.091036e-09 19.43798   209189_at
## 3 2.808211 7.686054 14.835515 3.051878e-13 3.595391e-09 18.69553 201694_s_at
## 4 1.700875 9.251381 11.538846 5.019678e-11 3.920727e-07 14.57411 201041_s_at
## 5 3.146005 5.175579 10.617313 2.533831e-10 1.731715e-06 13.19062   223316_at
## 6 3.482923 4.763760  9.961547 8.508265e-10 4.651894e-06 12.13533   220276_at
##   symbol change ENTREZID
## 1   FOSB     up     2354
## 2    FOS     up     2353
## 3   EGR1     up     1958
## 4  DUSP1     up     1843
## 5  CCDC3     up    83643
## 6  RERGL     up    79785

str(dcp[[2]])

## List of 2
##  $ upgenes  : chr [1:1361] "FOSB" "FOS" "EGR1" "DUSP1" ...
##  $ downgenes: chr [1:362] "NR1H4" "SLC51A" "NETO2" "ETNK1" ...

dcp[[3]]