Java教程

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

本文主要是介绍GEO代码分析流程 - 4. 差异分析 - 数据,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

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

 

这篇关于GEO代码分析流程 - 4. 差异分析 - 数据的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!