GEO数据挖掘系列，第14篇学习笔记：进行差异分析，得到deg差异基因，并为后续富集分析而补充些必要的数据
为了无缝衔接上一篇学习笔记，该文中的序号将接着上一篇来标注
6.获取差异基因deg，准备富集分析所需数据
本篇代码运行后将最终得到如图1所示的表格：①前六列是limma差异分析的输出结果；②后四列是为了画热图、火山图、做富集分析而补充上去的
图1
6.1前六列：差异分析
以下是二分组（实验组vs对照组）差异分析的固定步骤，且无需改动任何代码，只关注
deg
即可，不用管中间变量
design
和
fit
rm(list = ls())
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)# 得到模型矩阵
fit = lmFit(exp,design)# 线性拟合
fit = eBayes(fit)# 贝叶斯检验
deg = topTable(fit,coef = 2,number = Inf)# 提取结果
如果不想看中间变量，可以把上面四行代码写成一个函数，在运行自定义函数后，将只显示deg结果，如下
# chayi = function(exp,Group){
# design = model.matrix(~Group)
# fit = lmFit(exp,design)
# fit = eBayes(fit)
# deg = topTable(fit,coef = 2,number = Inf)
# }
以上这几行代码是可以根据
limma
R包的帮助文档写出来的，或者在网上搜索其他人分享的代码
此外，如果是多分组，则模型矩阵不一样，以及还要构建一个比较矩阵，这种情况后续再讨论
6.2后四列：补充数据
为deg数据框新添加四列
6.2.1加探针列
将行名即探针名，作为单独一列，加到数据框中，将这列名为probe_id
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#另外一个方法：tibble::rownames_to_column()
6.2.2
加symbol
列
为数据框加上探针注释，将symbol列也连接进去
【注意】
上一篇提到了“一个探针对应多个基因”，这里我们又要注意“多个探针对应到同一个基因”如图2这三个探针都对应到AAGAB基因了，这种情况会发生在很多芯片数据中，遇到了必须要处理
图2
解决方法有3种，但没有哪个是正确答案，请选取自己认为较为合适的一种来用
①随机去重（本文代码所用方法）：从这三个中随机选一个保留，如选了第一个
图3
②保留行和/行平均值最大的探针：行和max、行平均值max，这两个其实是一个意思
图4
③取多个探针的平均值
图5
# 探针的排序本身就是随机的，所以只要能够做到每种探针只剩一个，那就算是完成随机去重了
ids = distinct(ids,symbol,.keep_all = T)# 所以用distinct()函数按照symbol列去重复即可
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")# 按照两个共有的probe_id列来进行合并
nrow(deg)
6.2.3加change列
用来标记上下调基因，这样才能在画火山图的时候，让up的为红色，down的为蓝色，stable的为灰色
logFC_t = 1 #设置logFC的阈值
p_t = 0.05 #设置p的阈值
k1 = (deg$P.Value $logFC k2 = (deg$P.Value $logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
6.2.4加ENTREZID列（这一步的代码在作图完成后执行）
用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db) # 如果物种不是人类，而假如是小鼠/大鼠的，则这里要改成对应的R包
# Hs人类、Mn小鼠、Rn大鼠
s2e = bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#从SYMBOL转成ENTREZID，并且用org.Hs.eg.db这个R包的数据作参考
#一部分基因没匹配上是正常的,<30%甚至50%的失败都没事，只要不是只剩三五千行那么少就行
#该链接告诉我们哪个物种对应哪个org.xxx.db，如果找不到则无法使用bitr()这个转换函数了
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
nrow(deg)
deg = inner_join(deg,s2e,by=c("symbol"
