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]]## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 22283 features, 390 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: GSM1012278 GSM1012279 ... GSM1012667 (390 total)
## varLabels: title geo_accession ... tissue:ch1 (65 total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## pubMedIds: 19359472
## 28958617
## Annotation: GPL96查看元素的数据类型
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 13table(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]## GSM1012314 GSM1012315 GSM1012317 GSM1012318 GSM1012319
## 1007_s_at 406.0 341.0 340.0 272.0 390.0
## 1053_at 51.7 57.8 30.2 52.8 46.9
## 117_at 48.3 53.7 66.0 49.0 52.1
## 121_at 193.0 173.0 214.0 183.0 193.0
## 1255_g_at 30.2 29.1 29.4 37.7 31.0dat=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]## GSM1012314 GSM1012315 GSM1012317 GSM1012318 GSM1012319
## ZZZ3 70.8 68.9 51.8 84.1 105.00
## ZZEF1 111.0 108.0 112.0 85.6 94.90
## ZYX 395.0 190.0 445.0 295.0 200.00
## ZXDC 48.1 45.2 40.8 41.3 72.00
## ZXDB 13.7 16.4 13.2 14.5 9.6904. 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)## logFC AveExpr t P.Value adj.P.Val B
## TMEM135 14.465834 61.894709 4.632299 6.742306e-06 0.04883784 2.891058
## PIEZO2 6.138996 25.060212 4.535583 1.022129e-05 0.04883784 2.565965
## USP5 -18.161940 128.376720 -4.473758 1.329130e-05 0.04883784 2.360903
## PDE5A 6.223826 19.803704 4.438066 1.544890e-05 0.04883784 2.243508
## A2M 160.977775 477.118519 4.303340 2.704258e-05 0.05789707 1.806981
## SENP7 1.736058 8.290053 4.277159 3.010643e-05 0.05789707 1.723380查看个数，如下：
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 cutoff
xlim = 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]## GSM1012314 GSM1012315 GSM1012317 GSM1012318 GSM1012319
## PLG -1.0000000 -0.1079513 -0.1078853 -0.1079805 -0.1078910
## KNG1 -0.1081490 -0.1080782 -0.1079174 -0.1080653 -0.1079914
## HAND2-AS1 -0.1149248 -0.1147553 -0.1147930 -0.1143299 -0.1147591
## HPD -0.1032319 -0.1030374 -0.1031344 -0.1031754 -0.1032063## SERPINA10 -0.1016934 -0.1013278 -0.1008606 -0.1013278 -0.1011856注释信息的数据框
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.
关注公众号，扫码进群，免费解答，后期会有免费直播教程，敬请期待！
本文使用 文章同步助手 同步
网友评论