第一，同样的，清除变量+加载之前的数据
> rm(list = ls())
> load(file = "step2output.Rdata")
> group_list
[1] control control control treat treat treat
Levels: control treat
差异分析，用limma包来做
需要表达矩阵exp和group_list分组对象，不需要改
{
> library(limma)#加载limma包
#需要自己做好三个数据（表达矩阵，分组矩阵，差异比较矩阵）
#总共三个步骤（lmFit,eBayes,topTable）
+ design=model.matrix(~group_list)#把group_list设置成一个model.matrix
#> design
# (Intercept) group_listtreat
#1 1 0
#2 1 0
#3 1 0
#4 1 1
#5 1 1
#6 1 1
#attr(,"assign")
#[1] 0 1
#attr(,"contrasts")
#attr(,"contrasts")$group_list
#[1] "contr.treatment"
+ fit=lmFit(exp,design)#lmFit用于线性拟合；至少两个输入，一个是表达矩阵，一个是分组对象。
#表达矩阵必须是matrix类数据结构，每一列都是存放一个样本，每一行是一个探针信息或者是注释后的基因名
+ fit=eBayes(fit)#根据lmFit的拟合结果进行统计推断
+ deg=topTable(fit,coef=2,number = Inf)#得出差异比较矩阵
}
#view(deg)
差异表达矩阵deg.png
image.png
image.png
为deg数据框添加几列
{
#1.加probe_id列：把deg的行名（探针名）赋值给probe_id，并加一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))#mutate增加一列
head(deg)
#2.加symbol列，火山图要用
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
#按照symbol列去重复#因为存在一个探针映射多个基因的情况
deg <- deg[!duplicated(deg$symbol),]
}
#3.加change列,下面两行是阈值，可修改
logFC_t=1
logP_t = 0.01
test1 = deg$P.Value < logP_t
test2 = deg$logFC < -logFC_t
test3 = deg$logFC > logFC_t
{
change = ifelse(test1 & test2 ,
"down" ,
ifelse(test1 & test3 ,
"up",
"stable"))#不满足down和up时，就是stable;分成三种
deg <- mutate(deg,change)#加一列有down和up和stable
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类基因ID转换库
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
}
save(group_list,deg,file = "step4output.Rdata")
image.png
