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]]
