0. 准备
setwd("D:/R/CHOL")
rm(list = ls())
load(file = "step2output.Rdata")
1. 差异分析
差异分析，用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)
2. 分组信息为group_list，为二分类。
3. 线性拟合。
4. 统计计算。
5. 得到每个探针的logFC和P.Value，并按P.Value从小到大排序。
2. 为deg数据框添加列
加probe_id列，把行名变成一列
library(dplyr) deg <- mutate(deg,probe_id=rownames(deg)) #deg$probe_id=rownames(deg) head(deg)
2. 给数据框新增一列，列名为probe_id，值为数据框的行名。
3. 此行代码新增列的功能与上一行相同，但不会丢失行名。以上两种代码运行任意一种均可。
4. 查看前6行。
加symbol列，火山图要用
deg <- inner_join(deg,ids,by="probe_id") head(deg)
1. 按probe_id匹配deg和ids，保留匹配的行，并连接deg和ids的列，生成一个新数据框。 （inner_join不会改变行的顺序，行仍按deg的顺序排列。）
2. 查看前6行。
按照symbol列去重复
deg <- deg[!duplicated(deg$symbol),]
消除多个探针（probe_id）对应同一个基因（symbol）的情况。 重复的基因（symbol）只保留其第一次出现的探针（probe_id）。
加change列，标记上下调基因
logFC_t=1 P.Value_t = 0.01 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")) table(change) deg <- mutate(deg,change)
1. logFC的阈值，可修改。
2. P.Value的阈值，可修改。
3. 寻找下调基因。（若要按其他矫正后的p值，则修改代码中的“P.Value”为其他相应的矫正p值的列名）
4. 寻找上调基因。
5. 标记上下调基因。上调基因为up，下调基因为down，其余基因为stable。
6. 查看上调、下调、其余基因的数目。
7. 将标记添加到deg中成为新的一列。列名为change，值为up，down，stable。
加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)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#deg <- merge(deg,s2e,by.x = "symbol",by.y = "SYMBOL")
3. 人类物种的R包（.db）
其他物种：http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
4. 需要转换的基因id（symbol）列。
5. 需要转换的原基因类型，需大写。
6. 需要转换的目标基因类型，需大写。
7. 人类物种的R包。若研究对象为小鼠或大鼠等其他物种，则需修改。
得到的s2e为包含两列向量的数据框，列名分别为SYMBOL和ENTREZID。
8. 按symbol和SYMBOL匹配deg和s2e，保留匹配的行，并连接deg和s2e的列。
9. 此行代码功能同上，但inner_join不会改变deg本身的顺序。
保存数据
save(group_list,deg,logFC_t,P.Value_t,file = "step4output.Rdata")
