4. 差异分析 - 数据

rm(list = ls()) 
load(file = "step2output.Rdata")

#差异分析，用limma包来做
#需要表达矩阵（exp）和分组信息（group_list），不需要改
library(limma)
design=model.matrix(~group_list)              #分组信息为二分类。
fit=lmFit(exp,design)                         #线性拟合。
fit=eBayes(fit)                               #统计计算。
deg=topTable(fit,coef=2,number = Inf)         #得到每个探针的logFC和P.Value，并按P.Value从小到大排序。


#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))     #给数据框新增一列，列名为probe_id，值为数据框的行名。
                                                #使用tidyverse的函数处理数据框会丢失行名。
              #deg$probe_id=rownames(deg)     #此行代码新增列的功能与上一行相同，但不会丢失行名。
                                              #以上两种代码运行任意一种均可。
head(deg)                                     #查看前6行。
#2.加symbol列，火山图要用
deg <- inner_join(deg,ids,by="probe_id")      #按probe_id匹配deg和ids，保留匹配的行，并连接deg和ids的列。
                                                #inner_join不会改变行的顺序，行仍按deg的顺序排列。
head(deg)                                     #查看前6行。
  #按照symbol列去重复
deg <- deg[!duplicated(deg$symbol),]          #消除多个探针（probe_id）对应同一个基因（symbol）的情况。
                                                #重复的基因（symbol）只保留其第一次出现的探针（probe_id）。
#3.加change列,标记上下调基因
logFC_t=1                                                 #logFC的预定值，可修改。
P.Value_t = 0.01                                          #P.Value的预定值，可修改。
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)     #寻找下调基因。
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)      #寻找上调基因。
change = ifelse(k1,"down",ifelse(k2,"up","stable"))       #标记上下调基因。上调基因为up，下调基因为down，其余基因为stable。
deg <- mutate(deg,change)                                 #将标记添加到deg中成为新的一列。列名为change，值为up，down，stable。
#4.加ENTREZID列，用于富集分析
  #（symbol转entrezid，然后inner_join）
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)         #人类物种的R包（.db）
                              #其他物种：http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
s2e <- bitr(deg$symbol,                  #需要转换的基因id（symbol）列。
            fromType = "SYMBOL",         #需要转换的原基因类型，需大写。
            toType = "ENTREZID",         #需要转换的目标基因类型，需大写。
            OrgDb = org.Hs.eg.db)        #人类。
                                         #s2e为两列的数据框，列名为SYMBOL和ENTREZID。
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))    #按symbol和SYMBOL匹配deg和s2e，保留匹配的行，并连接deg和s2e的列。


save(group_list,deg,logFC_t,P.Value_t,file = "step4output.Rdata")