-
Notifications
You must be signed in to change notification settings - Fork 2
/
Limma包差异分析+heapmap+vol+GO+KEGG+GSEA
481 lines (404 loc) · 21.1 KB
/
Limma包差异分析+heapmap+vol+GO+KEGG+GSEA
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
######################################################################
####### Limma包差异分析+热图+火山图+GO分析+KEGG分析+GSEA分析 #########
######################################################################
### 1.读入探针矩阵数据
### 2.数据预处理：NA去除+log2转换+样本芯片矫正
### 3.探针ID和基因名symbol转换
### 4.基因Symbol表达矩阵获取：探针矩阵转换为基因Symbol矩阵+同名symbol保留表达量的最大值。
### 5.使用limma包来做芯片的差异分析(非配对)------选一个就行
### 5.使用limma包来做芯片的差异分析(配对)------选一个就行
### 6.输出差异基因
### 7. 作图，差异表达图、热图、火山图
### 8. 差异基因集的GO分析
### 9. 差异基因集的KEGG分析
### 10. 差异基因集的GSEA富集分析
### 400多行代码，你只要改18个地方，这个代码就能为你所用。
##################################################################################################
##################################################################################################
### 0. 镜像
##################################################################################################
options(download.file.method = 'libcurl')
options(url.method='libcurl')
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
##################################################################################################
### 1.读入探针矩阵数据
##################################################################################################
exprSet <- read.table("1_input.txt",
comment.char="!",
stringsAsFactors=F,
header=T) ### comment.char="!" 意思是！开头的行作为注释
rownames(exprSet) <- exprSet[,1]
exprSet <- exprSet[,-1]### 第一列变成行名
##################################################################################################
### 2.数据预处理：NA去除+log2转换+样本芯片矫正
##################################################################################################
## NA的去除
ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
## 自动判断是否需要log2转换，需要则进行log2转换
if (LogC) {
ex[which(ex <= 0)] <- NaN
## 取log2
exprSet <- log2(ex)
print("log2 transform finished")
}else{
print("log2 transform not needed")
}
## 样本芯片矫正
library(limma)
boxplot(exprSet,outline=FALSE, notch=T, las=2)#矫正之前单个样本值，此时平不平都无所谓。
exprSet=normalizeBetweenArrays(exprSet)### 该函数默认使用quntile 矫正样本间差异。
boxplot(exprSet,outline=FALSE, notch=T, las=2)#矫正之后单个样本值，此时是平的。
exprSet <- as.data.frame(exprSet)
#代码输出图片方式
pdf("Row data plot.pdf",width = 12,height = 9)
boxplot(exprSet,outline=FALSE, notch=T, las=2)
dev.off()
##################################################################################################
### 3.探针ID和基因名symbol转换
##################################################################################################
platformMap <- data.table::fread("1_platformMap.txt",data.table = F) ## 读入platformMap.txt文件，platformMap 中有常见的GPL平台和ID转换需要的R包（一一对应）。
index <- "GPL6102" ## 第11111111111处需要改的地方，改成你自己数据的GPL平台,输入GPL号，找到该平台对应的ID转
## 换R包，这个在GEO检索的页面有
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db") ## 我们要看这一步结果，看这个GPL平台是通过什么R包来进行ID转换的的。
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("illuminaHumanv2.db") ## 第22222222222处需要改的地方，将illuminaHumanv2.db改为你需要安装的包，##安装这包。
library(illuminaHumanv2.db) ## 第3333333333333333处需要改的地方， ## 加载该R包。
probe2symbol_df <- toTable(get("illuminaHumanv2SYMBOL")) ## 第44444444444处需要改的地方，用该R包里面的illuminaHumanv2SYMBOL函数将探针ID转换成
## 基因symbol，得到ID-symbol对应表
##################################################################################################
### 4.基因Symbol表达矩阵获取：探针矩阵转换为基因Symbol矩阵+同名symbol保留表达量的最大值。
##################################################################################################
library(dplyr)
library(tibble)
exprSet <- exprSet %>% ##万能代码，不管用##
rownames_to_column("probe_id") %>%
inner_join(probe2symbol_df,by="probe_id") %>%
dplyr::select(-probe_id)%>%
dplyr::select(symbol,everything()) %>%
mutate(rowMean =rowMeans(.[,-1])) %>%
dplyr::arrange(desc(rowMean)) %>%
distinct(symbol,.keep_all = T) %>%
dplyr::select(-rowMean) %>%
column_to_rownames("symbol")
save(exprSet,file = "2_exprSet_rmdup.Rdata") #将基因表达矩阵保存起来
##################################################################################################
### 5.使用limma包来做芯片的差异分析---非配对差异分析
##################################################################################################
rm(list = ls())#清除
## 0. 加载表达数据
load(file = "2_exprSet_rmdup.Rdata")
####################在TCGA数据库中创建分组########################
group=sapply(strsplit(colnames(exprSet),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)
sum(group==1)
sum(group==0)
####################在TCGA数据库中创建分组########################
## 1.创建分组
group <- c(rep("con",18),rep("treat",18))##第555555555555处需要改的地方，意思是在矩阵里面前面18个对照组，后面18个实验组
group <- factor(group,levels = c("con","treat"))
## 样本PCA主成分分析：提前预测结果，看两组样本是否分开，分开则说明差异挺大的，可以进行差异分析。
res.pca <- prcomp(t(exprSet), scale = TRUE)
library(factoextra)
fviz_pca_ind(res.pca,col.ind = group)
# #
# 如果有个别样本分的不好，应该在input文件里面将这个样本删除才行#
# #
## 构建比较矩阵
design <- model.matrix(~group)
## 比较矩阵命名
colnames(design) <- levels(group)
design
## 2.线性模型拟合
fit <- lmFit(exprSet,design)
## 3.贝叶斯检验
fit2 <- eBayes(fit)
## 4.输出差异分析结果,其中coef的数目不能超过design的列数
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf, p.value=0.05)
## 这个数据很重要需要保存一下
save(allDiff,file = "2_allDiff.Rdata")
##################################################################################################
### 5.使用limma包来做芯片的差异分析---配对差异分析
##################################################################################################
rm(list = ls())#清除
## 0. 加载表达数据
load(file = "2_exprSet_rmdup.Rdata")
## 1.创建分组
group <- c(rep("con",18),rep("treat",18)) ##第5555555处需要改的地方，意思是在矩阵里面前面18个对照组，后面18个实验组
group <- factor(group,levels = c("con","treat"))
## 2.样本PCA主成分分析：提前预测结果，看两组样本是否分开，分开则说明差异挺大的，可以进行差异分析。
res.pca <- prcomp(t(exprSet), scale = TRUE)
#install.packages("FactoMineR")
#install.packages("factoextra")
library(factoextra)
fviz_pca_ind(res.pca,col.ind = group)
design <- model.matrix(~group)
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
### 3. 在此基础上进行配对分析
group <- c(rep("con",18),rep("treat",18)) ##第5555555555555处需要改的地方，意思是在矩阵里面前面18个对照组，后面18个实验组
group <- factor(group,levels = c("con","treat"))
pairinfo = factor(rep(1:18,2))##第555555555555处需要改的地方，加一个配对信息，几个样本就是1：几
### 4. 配对分析
design=model.matrix(~ pairinfo+group)
fit <- lmFit(exprSet,design)
fit2 <- eBayes(fit)
allDiff=topTable(fit2,adjust='fdr',coef=ncol(design),number=Inf, p.value=0.05)
## 这个数据很重要需要保存一下
save(allDiff,file = "2_allDiff.Rdata")
write.csv(allDiff,file= "2_limma_result.csv",row.names = T)
##################################################################################################
### 6.输出差异基因
##################################################################################################
## 定义差异基因：logFC为0.5，矫正后的p值小于0.05
library(dplyr)
diffgene <- allDiff %>%
filter(adj.P.Val < 0.05) %>% #第6666666处需要改的地方
filter(abs(logFC) >0.5) #第777777777处需要改的地方
## 差异基因和矩阵合并
diffgene_count <- merge(as.data.frame(diffgene),as.data.frame(as.data.frame(exprSet)),by="row.names",sort=FALSE)
## 保存差异基因的表达矩阵
write.csv(diffgene_count,file= "2_diffgene_count.csv",row.names = T)
##################################################################################################
### 7. 作图，差异表达图，热图，火山图。
##################################################################################################
## 1. 数据处理成作图的格式
exprSet <- as.data.frame(t(exprSet))#矩阵转置
dd <- cbind(group=group,exprSet)#添加分组信息
test = dd[,1:10]
## 2.作图展示一个差异基因,我这里以CD36为例
my_comparisons <- list(c("treat", "con"))
library(ggpubr)
ggboxplot(
dd, x = "group", y = "CD36", #第888888888888处需要改的地方，把CD36改成自己的基因。
color = "group", palette = c("#00AFBB", "#E7B800"),
add = "jitter")+
stat_compare_means(comparisons = my_comparisons, method = "t.test")
## 改写成函数
diffplot <- function(gene){
my_comparisons <- list(
c("treat", "con")
)
library(ggpubr)
ggboxplot(
dd, x = "group", y = gene,
color = "group", palette = c("#00AFBB", "#E7B800"),
add = "jitter"
)+
stat_compare_means(comparisons = my_comparisons, method = "t.test")
}
diffplot("CD36") ##第999999999999处需要改的地方，把CD36改成自己的基因
diffplot("MOXD1") ##第10101010101010处需要改的地方，把CD36改成自己的基因
### 4.画热图，选取一群adj.P.Val < 0.05和abs(logFC) >1的差异基因做热图
rm(list = ls())
## 1.加载表达数据
load(file = "2_exprSet_rmdup.Rdata")
## 2.加载差异列表
load(file = "2_allDiff.Rdata")
library(dplyr)
diffgene <- allDiff %>%
filter(adj.P.Val < 0.05) %>% ##第111111111111处需要改的地方
filter(abs(logFC) >1) ##第12121212121212处需要改的地方，注意看这个阈值有没有结果
## 3.用名称提取部分数据用作热图绘制
heatdata <- exprSet[rownames(diffgene),]
## 4.制作一个分组信息用于注释
group <- c(rep("con",18),rep("treat",18)) ##第131313131313处需要改的地方，意思是在矩阵里面前面18个对照组，后面18个实验组
annotation_col <- data.frame(group)
rownames(annotation_col) <- colnames(heatdata)
## 加载热图的R包
library(pheatmap)
library(viridisLite)
### 经过修饰的图
pheatmap(heatdata, #热图的数据
cluster_rows = TRUE,#行聚类
cluster_cols = FALSE,#列聚类，可以看出样本之间的区分度
annotation_col =annotation_col, #标注样本分类
annotation_legend=TRUE, # 显示注释
show_rownames = T,# 显示行名
scale = "row", #以行来标准化，这个功能很不错
color = viridis(10, alpha = 1, begin = 0.5, end = 1, direction = 1),#调色
#filename = "heatmap_F.pdf",#是否保存
cellwidth = 9, # 格子宽度
cellheight = 9,# 格子高度
fontsize = 9 # 字体大小
)
### 5. 画火山图，即展示所有基因差异变化
rm(list = ls())
library(ggplot2)
library(ggrepel)
library(dplyr)
load(file = "2_allDiff.Rdata")
### 无论是什么名称，都改为data
data <- allDiff
data$gene <- rownames(data)
logFCfilter = 0.5 ### 0.5处画一条阈值线， ##第1414141414处需要改的地方，
logFCcolor = 1 ### 1以上的基因标记出来为绿色， ##第1515151515处需要改的地方，
### 标记上下调
index = data$adj.P.Val <0.05 & abs(data$logFC) > logFCfilter
data$group <- 0
data$group[index & data$logFC>0] = 1
data$group[index & data$logFC<0] = -1
data$group <- factor(data$group,levels = c(1,0,-1),labels =c("Up","NS","Down") )
### 正式画图
ggplot(data=data, aes(x=logFC, y =-log10(adj.P.Val),color=group)) +
geom_point(alpha=0.8, size=1.2)+
scale_color_manual(values = c("darkred", "black", "royalblue4"))+
labs(x="log2 (fold change)",y="-log10 (adj.P.Val)")+
theme(plot.title = element_text(hjust = 0.4))+
geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(-logFCfilter,logFCfilter),lty=4,lwd=0.6,alpha=0.8)+
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
theme(legend.position="top")+
geom_point(data=subset(data, abs(logFC) >= logFCcolor & adj.P.Val <0.05),alpha=0.8, size=3,col="green2")+
geom_text_repel(data=subset(data, abs(logFC) >= logFCcolor & adj.P.Val <0.05),
aes(label=gene),col="black",alpha = 0.8)
##################################################################################################
### 8. 差异基因集的GO富集分析。
##################################################################################################
rm(list = ls())
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
#ID转换，Symbol转换为entrezID
rt <-read.table(file="2_diffgene_count.csv", header=TRUE, sep=",")
rt <- rt[,-1]### 删除第一列
genes=as.vector(rt[,1])
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA) #找出基因对应的entrezID
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
#GO富集分析
rt=out[is.na(out[,"entrezID"])==F,] #去除基因id为NA的基因
gene=rt$entrezID
geneFC=rt$logFC
#geneFC=2^geneFC
names(geneFC)=gene
kk=enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =1, qvalueCutoff = 1, ont="all", readable =T)#人的富集分析
GO=as.data.frame(kk)
GO=GO[(GO$pvalue<0.05),] #将P小于0.05的结果保留下来。
write.table(GO,file="2_GO_P0.05.txt",sep="\t",quote=F,row.names = F) #保存富集结果
#定义显示Term数目
showNum=10
if(nrow(GO)<30){
showNum=nrow(GO)
}
#画柱状图,label_format=100参数让字体不重叠
pdf(file="2_GO_P0.05_Top10_barplot.pdf",width = 9,height =7)
bar=barplot(kk, drop = TRUE, showCategory =showNum,split="ONTOLOGY",color = "pvalue", label_format=100) + facet_grid(ONTOLOGY~., scale='free')
print(bar)
dev.off()
#画气泡图
pdf(file="2_GO_P0.05_Top10_bubble.pdf",width = 9,height =7)
bub=dotplot(kk,showCategory = showNum, orderBy = "GeneRatio",split="ONTOLOGY", color = "pvalue", label_format=100) + facet_grid(ONTOLOGY~., scale='free')
print(bub)
dev.off()
#画圈图
pdf(file="2_GO_P0.05-Top5_circos.pdf",width = 9,height = 6.5)
cnet=cnetplot(kk, foldChange=geneFC, showCategory = 5, circular = TRUE, colorEdge = TRUE)
print(cnet)
dev.off()
##################################################################################################
### 9. 差异基因集的KEGG富集分析。
##################################################################################################
rm(list = ls())
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
library("pathview")
#ID转换，Symbol转换为entrezID
rt <-read.table(file="2_diffgene_count.csv", header=TRUE, sep=",")
rt <- rt[,-1]### 删除第一列
genes=as.vector(rt[,1])
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA) #找出基因对应的entrezID
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
#kegg富集分析-2022年7月7日更新此处--公众号：科研部
rt=out[is.na(out[,"entrezID"])==F,] #去除基因id为NA的基因
gene=rt$entrezID
geneFC=rt$logFC
#geneFC=2^geneFC
names(geneFC)=gene
# install.packages("R.utils")
library(R.utils)
R.utils::setOption("clusterProfiler.download.method","auto")
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =1, qvalueCutoff =1)
library(DOSE)#entrez gene id转换成gene symbol
kk<-setReadable(kk, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
KEGG=as.data.frame(kk)
KEGG=KEGG[(KEGG$pvalue<0.05),]
write.table(KEGG,file="2_KEGG_P0.05.txt",sep="\t",quote=F,row.names = F)
#定义显示Term数目
showNum=30
if(nrow(KEGG)<showNum){
showNum=nrow(KEGG)
}
#柱状图
pdf(file="2_KEGG_P0.05_Top30_barplot.pdf",width = 9,height = 7)
barplot(kk, drop = TRUE, showCategory = showNum, color = "pvalue", label_format=100)
dev.off()
#气泡图
pdf(file="2_KEGG_P0.05_Top30_bubble.pdf",width = 9,height = 7)
dotplot(kk, showCategory = showNum, orderBy = "GeneRatio",color = "pvalue", label_format=100)
dev.off()
#圈图
pdf(file="2_KEGG_P0.05_Top5_circos.pdf",width = 11,height = 7)
kkx=setReadable(kk, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(kkx, foldChange=geneFC,showCategory = 5, circular = TRUE, colorEdge = TRUE,node_label="all")
dev.off()
#通路图---会map多很多信号通路，并下载下来
keggxls=read.table("2_KEGG_P0.05.txt",sep="\t",header=T)
dir.create("2_KEGG_P0.05_Mapplot") #创建名为2_KEGG_P0.05_Mapplot的文件夹
setwd("C:\\Users\\A\\Desktop\\GEO_limma（配对+非配对）\\2_KEGG_P0.05_Mapplot") ###第161616161616处需要改的地方，进入这个文件夹，将后面的生成的文件保存在这个文件夹下。
for(i in keggxls$ID){
pv.out <- pathview(gene.data = gene, pathway.id = i, species = "hsa", out.suffix = "pathview")
}
##################################################################################################
### 10. 所有基因集的GSEA富集分析。
##################################################################################################
rm(list = ls())
library(clusterProfiler)
setwd("C:\\Users\\A\\Desktop\\GEO_limma（配对+非配对）") ##第17171717171717处需要改的地方，回到原来的文件夹
load(file = "2_allDiff.Rdata")
library(dplyr)
df <- select(allDiff,-c(2:6))#删除后面不要的2-6列
df$SYMBOL <- rownames(df)#行名变为第一列
head(df)#查看前面几行
dim(df)#数据总共几行几列
#用clusterProfiler包里面的bitr转换ENTREZID和SYMBOL。
df.id<-bitr(df$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db") #这是人的，小鼠就把org.Hs.eg.db改为org.Mm.eg.db就行
#让基因名、ENTREZID、logFC对应起来
easy.df<-merge(df,df.id,by="SYMBOL",all=F)
#按照logFC排序
sortdf<-easy.df[order(easy.df$logFC, decreasing = T),]#降序排序
#把logFC按照从大到小提取出来
gene.expr = sortdf$logFC
#给上面提取的flogFC对应上ENTREZID
names(gene.expr) <- sortdf$ENTREZID
####现在进行GSEA分析####
library(R.utils)
R.utils::setOption("clusterProfiler.download.method","auto")
kk <- gseKEGG(gene.expr, organism = "hsa", pvalueCutoff = 1) ####如果报错，在human和hsa，rno和rat之间来回切换，这跟您的网络相关，时间可能需要比较长。
library(DOSE)#entrez gene id转换成gene symbol
kk<-setReadable(kk, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
dim(kk)
#按照enrichment score从高到低排序，便于查看富集通路
sortkk<-kk[order(kk$enrichmentScore, decreasing = T),]
head(sortkk)
dim(sortkk)
write.table(sortkk,"2_Allgene_gsea_output.txt",sep = "\t",quote = F,col.names = T,row.names = F) #最后我们ID还要换成名字。
#画自己想画的，这些hsa开头的都是对应刚刚生成的文件里面的，选择性的画就行。 #第18181818181818处需要改的地方，
paths <- c("hsa03320", "hsa03010", "hsa04350", "hsa00062", "hsa04510", "hsa04022", "hsa03030", "hsa04020", "hsa04151", "hsa04152")
library(enrichplot)
gseaplot2(kk, paths)
gseaplot2(kk, paths, color = colorspace::rainbow_hcl(10), pvalue_table = TRUE) ####把p值标上。