GEO代码分析流程 - 4. 差异分析 - 数据

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列，火山图要用

1 2	`deg <-` `inner_join(deg,ids,by="probe_id")` `head(deg)`

1. 按probe_id匹配deg和ids，保留匹配的行，并连接deg和ids的列，生成一个新数据框。（inner_join不会改变行的顺序，行仍按deg的顺序排列。）
2. 查看前6行。

按照symbol列去重复

1	`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。