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")