RNA 2. SCI文章中基于GEO的差异表达基因之 limma
01. 前言
目前随着测序成本的降低,越来越多的文章中都将使用基因的表达数据,那么就会涉及到差异基因,最基础的思路就是获得表达数据,根据临床信息进行分组,利用各种软件计算出差异表达基因,这里主要基于GEO的数据介绍 limma 软件包,这个非常常用,尤其在做数据库的表达分析中。
02. 软件包安装
limma, GEOquery安装软件包也非常简单,无非就是那么几种方式,我们看下这两个软件怎么安装,如下:
#BiocManager::install("GEOquery")
#BiocManager::install("limma")
03. GEO 数据读取
GEO (Gene Expression Omnibus),是由美国国立生物技术信息中心 NCBI 创建并维护的基因表达数据库,创建于2000年,收录了世界各国研究机构提交的高通量基因表达数据,也就是说只要是目前已经发表的论文,论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。GEOquery 是专门用于读取 GEO 数据库的数据,用起来很方便,我们这里就是利用其有缺的特点成功获得了 prob 和样本的表达矩阵,但是这并不是目的,还需要我们将prob转换为基因symbol,这才是我们想要的结果,如下:
魔幻操作,一键清空~当前环境中对象全部删除
rm(list = ls())
在调用as.data.frame的时,将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
options(stringsAsFactors = F)
把GSE41258_eSet.Rdata赋值给f,方便后面流程化处理
f='GSE41258_eSet.Rdata'
这是一个函数,利用包将数据集的表达信息下载下来,赋值给了gset,而不下载注释信息和平台信息,并保存到本地,文件名为f。
if(!file.exists(f)){gset <- getGEO('GSE41258', destdir=".",AnnotGPL = F, ## 注释文件getGPL = F) ## 平台文件save(gset,file=f) ## 保存到本地
}
载入数据
load('GSE41258_eSet.Rdata')
查看数据类型
class(gset)
[1] “list”
看一下有几个元素
length(gset)
[1] 1
取第一个元素
gset[[1]]
查看元素的数据类型
class(gset[[1]])
[1] “ExpressionSet”
attr(,“package”)
[1] “Biobase”
因为这个GEO数据集只有一个GPL平台,所以下载到的是一个含有一个元素的list 取出第一个元素赋值给一个对象a
a=gset[[1]]
a现在是一个对象,取a这个对象通过看说明书知道要用exprs这个函数,该函数得到表达矩阵,现在 得到的dat就是一个表达矩阵,只不过基因的ID是探针名
dat=exprs(a)
看一下dat这个矩阵的维度
dim(dat)
[1] 22283 390
查看dat这个矩阵的1至4行和1至4列,逗号前为行,逗号后为列。这个表达矩阵是已经log之后的,表达量一般是0-10左右,如果是原始芯片表达的信号值一般是几千到一万,则需要log处理。
dat[1:5,1:5]## GSM1012278 GSM1012279 GSM1012280 GSM1012281 GSM1012282
## 1007_s_at 705.0 825.0 755.0 803.0 541.0
## 1053_at 22.0 17.1 17.4 11.6 47.1
## 117_at 34.3 61.8 59.8 61.3 51.3
## 121_at 188.0 263.0 190.0 197.0 241.0
## 1255_g_at 29.9 39.2 32.9 35.4 34.1
画个图看一下各样本之间有没有批次效应,一般中位数都差不多,las是将横坐标样本信息竖着排列
boxplot(dat,las=2)
通过查看说明书知道取对象a里的临床信息用pData
pd=pData(a)
查看一下,挑选一些感兴趣的临床表型,这里我们将得到其分组title信息。
View(pd)
运行一个字符分割包
library(stringr)
临床数据分组,获得复发组和未复发组,用于后续比较分析,如下
RE=rownames(pd[grepl('recurrence: Yes',as.character(pd$characteristics_ch1.13)),]) #复发组
PRI=rownames(pd[grepl('recurrence: No',as.character(pd$characteristics_ch1.13)),]) #未复发组
dat=dat[,c(RE,PRI)]
group_list=c(rep('RE',length(RE)),rep('PRI',length(PRI))) #分组信息
看一下两个分组各有几个
table(group_list)## group_list
## PRI RE
## 134 55
我们下载到的表达数据是芯片探针表示,需要我们把探针对应到基因symbol号上,经过下面操作即可获得,但是这个需要我们找对配置文件,一般来说都会有对应的数据包,如果没有需要我们自己弄,这个有时间会专题讲一下。首先我们采用的数据集 GSE41258 对应数据包为hgu133a.db,安装并加载数据包,对这个包进行探索,看一下有多少元素,找到有SYMBOL的那个元素就是我们需要的对应关系。例如:[34] “hgu133aSYMBOL”
ls("package:hgu133a.db")## [1] "hgu133a" "hgu133a.db" "hgu133a_dbconn"
## [4] "hgu133a_dbfile" "hgu133a_dbInfo" "hgu133a_dbschema"
## [7] "hgu133aACCNUM" "hgu133aALIAS2PROBE" "hgu133aCHR"
## [10] "hgu133aCHRLENGTHS" "hgu133aCHRLOC" "hgu133aCHRLOCEND"
## [13] "hgu133aENSEMBL" "hgu133aENSEMBL2PROBE" "hgu133aENTREZID"
## [16] "hgu133aENZYME" "hgu133aENZYME2PROBE" "hgu133aGENENAME"
## [19] "hgu133aGO" "hgu133aGO2ALLPROBES" "hgu133aGO2PROBE"
## [22] "hgu133aMAP" "hgu133aMAPCOUNTS" "hgu133aOMIM"
## [25] "hgu133aORGANISM" "hgu133aORGPKG" "hgu133aPATH"
## [28] "hgu133aPATH2PROBE" "hgu133aPFAM" "hgu133aPMID"
## [31] "hgu133aPMID2PROBE" "hgu133aPROSITE" "hgu133aREFSEQ"
## [34] "hgu133aSYMBOL" "hgu133aUNIPROT"
toTable这个函数:通过看hgu133plus2.db这个包的说明书知道提取probe_id(探针名)和symbol(基因名)的对应关系的表达矩阵的函数为toTable。并赋值给ids 这时候我们得到19825个探针对应的基因名。刚才我们的表达矩阵中是22283个基因探针,这就意味着刚才的表达矩阵中可能存在多个探针重复对应一个基因名。这就需要我们对数据进行进一步筛选、处理。
ids=toTable(hgu133aSYMBOL)
将列名统一改为’probe_id’,'symbol’方便后续统一操作。
colnames(ids)=c('probe_id','symbol')
12645个独特的基因探针,意味着本来19825个里面有一部分是重复的
length(unique(ids$symbol))## [1] 12645
每个对象出现的个数
tail(sort(table(ids$symbol)))
##
## PDLIM5 FGFR2 PTGER3 CFLAR TCF3 HFE
## 9 10 10 11 11 13
table(sort(table(ids$symbol)))##
## 1 2 3 4 5 6 7 8 9 10 11 13
## 8066 2704 1184 446 148 62 16 9 5 2 2 1
```
画图观察
plot(table(sort(table(ids$symbol))))
%in%用于判断是否匹配,然后取匹配的几行,去掉无法匹配的信息。
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in% rownames(dat),]
取表达矩阵中可以与探针名匹配的那些,去掉无法匹配的表达数据,这时只剩下19825个探针及表达信息,其余已被剔除。
dat[1:5,1:5]
dat=dat[ids$probe_id,]
```
ids新建median这一列,列名为median,同时对dat这个矩阵按照行操作,取每一行的中位数,将结果给到median这一列的每一行
ids$median=apply(dat,1,median)
对ids按照median中位数从大到小排列的顺序排序 ##即先按symbol排序,相同的symbol再按照中位数从大到小排列,方便后续保留第一个值。将对应的行赋值为一个新的ids,这样order()就相当于sort()
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
将symbol这一列取取出重复项,’!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果.最后得到18832个基因。
ids=ids[!duplicated(ids$symbol),]
新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat
dat=dat[ids$probe_id,]
把ids的symbol这一列中的每一行给dat作为dat的行名
rownames(dat)=ids$symbol
查看最后表达数据前五行五列,如下:
dat[1:5,1:5]
04. limma计算差异表达基因
首先加载软件包,如下:
limma包的具体用法参考 limma Users Guide 构建分组信息,构建好比较矩阵是关键 注意这里的表达矩阵信息dat是经过处理后的,为转置,行为gene,列为sample
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)## PRI RE
## 1 0 1
## 2 0 1
## 3 0 1
## 4 0 1
## 5 0 1
## 6 0 1
比较信息,contrast.matrix##查看比较矩阵的信息,这里我们设置的是RE Vs. PRI
contrast.matrix<-makeContrasts("RE-PRI",levels = design)
拟合模型
fit <- lmFit(dat,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
差异表达基因计算,其中coef比较分组 n基因数
DEG<-topTable(fit2, coef=1, n=Inf) %>% na.omit()
head(DEG)
查看个数,如下:
dim(DEG)
## [1] 12645 6
保存文件,如下:
save(DEG,file = "DEG_all.Rdata")
05. 差异表达基因展示方式
绘制火山图
借助于有人开发的更高级的包,用于完成某些特殊的功能,或者更美观。安装软件包并加载,如下:
绘制火山图利用EnhancedVolcano,这个包功能强大,可以任意修改图上的标签、字体、增加主副标题等,作图效果如下:
EnhancedVolcano(DEG,lab = rownames(DEG),labSize = 2,x = "logFC",y = "P.Value",# selectLab = rownames(DEG)[1:4],xlab = bquote(~Log[2]~ "fold change"),ylab = bquote(~-Log[10]~italic(P)),pCutoff = 0.01,## pvalue阈值FCcutoff = 2,## FC cutoffxlim = c(-5,5),ylim = c(0,5),colAlpha = 0.6,legendLabels =c("NS","Log2 FC"," p-value"," p-value & Log2 FC"),legendPosition = "top",legendLabSize = 10,legendIconSize = 3.0,pointSize = 1.5,title ="limma results",subtitle = 'Differential expression',)
- 绘制热图
热图的使用比较频繁,得到差异基因后可以直接绘制热图相对简单好用的要属pheatmap包了管道中的常规提取需要加上特殊的占位符。首先提取出想要画的数据,提取FC高的top 50,管道提取
up_50<-DEG %>% as_tibble() %>% mutate(genename=rownames(DEG)) %>% dplyr::arrange(desc(logFC)) %>% .$genename %>% .[1:50]
提取FC低的前50个基因
down_50<-DEG %>% as_tibble() %>% mutate(genename=rownames(DEG)) %>% dplyr::arrange(logFC) %>% .$genename %>% .[1:50]
index<-c(up_50,down_50)
对于表达数据来说,在做热图或者聚类分析时,需要我们进行归一化,这样数据才不会过于分散,如下:
index_matrix<-t(scale(t(dat[index,])))
index_matrix[index_matrix>1]=1
index_matrix[index_matrix<-1]=-1
index_matrix[1:5,1:5]
注释信息的数据框
anno=data.frame(group=group_list)
rownames(anno)=colnames(index_matrix)
head(anno)## group
## GSM1012314 RE
## GSM1012315 RE
## GSM1012317 RE
## GSM1012318 RE
## GSM1012319 RE
## GSM1012320 RE
绘制热图,如下:
pheatmap(index_matrix,show_colnames =F,show_rownames = F,cluster_cols = T, annotation_col=anno)
第一次利用 R markdown 生产 md 文件之后在导入公众号,特别好用,果然明白高手们都是怎样炼成的了,虽然稍微好用些,但是还是有很多弊端,需要之后调整,总的来说,学来的都是自己的,蛮有意思的!
这期主要就是基于GEO 芯片数据 利用 limma 做差异表达分析,我觉得算是全网最完成的版本,下期再说几个软件。
Reference:
Viechtbauer, W. (2010). Conducting Meta-Analyses in R with the metafor Package. Journal of Statistical Software, 36(3), 1–48.
Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., & Mosteller, F. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. Journal of the American Medical Association, 271(9), 698–702.
Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A random-effects regression model for meta-analysis. Statistics in Medicine, 14(4), 395–411.
van Houwelingen, H. C., Arends, L. R., & Stijnen, T. (2002). Advanced methods in meta-analysis: Multivariate approach and meta-regression. Statistics in Medicine, 21(4), 589–624.
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48.
关注公众号,扫码进群,免费解答,后期会有免费直播教程,敬请期待!
桓峰基因
生物信息分析,SCI文章撰写及生物信息基础知识学习:R语言学习,perl基础编程,linux系统命令,Python遇见更好的你
35篇原创内容
公众号
RNA 2. SCI文章中基于GEO的差异表达基因之 limma相关推荐
- RNA 24. SCI文章中基于TCGA的免疫浸润细胞分析的在线小工具——TIMER
点击关注,桓峰基因 桓峰基因 生物信息分析,SCI文章撰写及生物信息基础知识学习:R语言学习,perl基础编程,linux系统命令,Python遇见更好的你 135篇原创内容 公众号 今天来介绍一个使 ...
- RNA 30. SCI文章中基于TCGA和GTEx数据挖掘神器(GEPIA2)
这期介绍一个基于TCGA和GTEx数据挖掘神器(GEPIA2),个人觉得如果没有编程基础的可以直接利用这个在线小工具分析自己的研究的单个基因或者多个基因,效果还是蛮好的! 桓峰基因公众号推出转录组分析 ...
- RNA 29. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER2.0)
桓峰基因公众号推出转录组分析教程,有需要生信的老师可以联系我们!转录分析教程整理如下: RNA 1. 基因表达那些事--基于 GEO RNA 2. SCI文章中基于GEO的差异表达基因之 limma ...
- RNA 4. SCI 文章中基于TCGA 差异表达之 edgeR
01. 前言 在不久的将来,预计新兴的数字基因表达(DGE) 技术在许多功能基因组学应用方面将超越微阵列技术.其中的一个基础数据分析任务,尤其是基因表达研究,涉及确定是否有证据表明对于转录本或外显子在 ...
- RNA 1. SCI 文章中读取 GEO 数据
依稀记得10年前一个样本的RNA-SEQ费用还是蛮高的,而现在就是洒洒水啦,所以这种转录组的数据也已经成为文章的主流,基本上就是WES结合RNA,但是个人感觉,真的能关联到突变基因和表达的水平上去的, ...
- RNA 27 SCI文章中转录因子结合motif富集到调控网络 (RcisTarget)
点击关注,桓峰基因 桓峰基因公众号推出转录组分析和临床预测模型教程,有需要生信的老师可以联系我们!首选看下转录分析教程整理如下: RNA 1. 基因表达那些事–基于 GEO RNA 2. SCI文章中 ...
- RNA 25. SCI文章中只有生信没有实验该怎么办?
今天来介绍一个使用非常方便的在线免疫组化分析工具--PHA (The Human atlas),免疫组化时单基因干湿结合生信中最常见的补充实验的方法之一,性价比较高.但是如果没有条件进行自己样本的免疫 ...
- RNA 25. SCI文章中估计组织浸润免疫细胞和基质细胞群的群体丰度(MCP-counter)
点击关注,桓峰基因 今天来介绍一个利用基因表达估计组织浸润免疫细胞和基质细胞群的群体丰度的软件包--MCP-counter,亲试,非常好用. 桓峰基因的教程不但教您怎么使用,还会定期分析一些相关的文章 ...
- RNA 15. SCI 文章中的融合基因之 FusionGDB2
基于 RNA 数据分析1-13期基本介绍完成,而基因融合同样也是转录组测序中能够获得的对于临床上非常有意义的数据,这期就看看融合基因该怎么分析,增添文章的内容. 一. 融合基因 融合基因就是两个基因& ...
最新文章
- java basic认证_Basic认证
- 【机器学习】通俗易懂的条件概率公式
- 程序员食品营养(1)-面包基础
- linux 多线程 多进程同步
- 将图片序列化和反序列化
- poj3981 字符串替换-字符串的基本操作
- UI登录表单使用模板素材
- Java常见设计模式总结
- 网易云音乐APP(基于APICloud平台)
- 制造业管理软件如何帮助企业解决仓库管理难题?
- gcc开启C99或C11标准支持
- @Valid和@Validated注解校验List<Object>
- Ubuntu环境下moos-Ivp编译
- mac网页java无法加载,chrome浏览器mac版无法加载怎么办_chrome浏览器mac版打不开网页解决方法-win7之家...
- 莫烦---Pytorch学习
- 打造H5动感影集的爱恨情仇(动画性能篇)
- 51单片机8*8点阵显示“中国”
- Themleaf中html向controller传递参数 th:onclick传参
- 如何动态获得view的大小
- java萍方字体_html苹方字体 - osc_wv1mxwu2的个人空间 - OSCHINA - 中文开源技术交流社区...