一个R包，得到差异分析结果列表，差异基因，三张图(热图，PCA，火山图)

---

   大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

   就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

   这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

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