ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

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

2022-09-10 01:02:35  阅读:260  来源: 互联网

标签:分析 fit list 流程 基因 down design GEO deg


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

 

标签:分析,fit,list,流程,基因,down,design,GEO,deg
来源: https://www.cnblogs.com/xiaogaobugao/p/16675823.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有