柚子快報邀請碼778899分享：R GO和kEGG富集分析

http://yzkb.51969.com/

文章目錄

前言一、GO和KEGG1. **GO 富集分析：**2. KEGG 富集分析：

二、使用步驟1.數(shù)據(jù)處理2. GO分析3.KEGG富集

總結(jié)

前言

GO（Gene Ontology，基因本體）富集和 KEGG（Kyoto Encyclopedia of Genes and Genomes，京都基因與基因組百科全書）富集分析能夠從不同角度揭示基因的功能和生物學(xué)意義

一、GO和KEGG

1. GO 富集分析：

說明基因在分子功能（Molecular Function）、生物過程（Biological Process）和細(xì)胞組成（Cellular Component）三個方面的特征和傾向。幫助了解基因參與的具體生物學(xué)活動，例如基因是具有催化活性、結(jié)合能力，還是參與細(xì)胞分裂、信號轉(zhuǎn)導(dǎo)等過程，以及在細(xì)胞的哪個部位發(fā)揮作用。

2. KEGG 富集分析：

反映基因所參與的代謝通路、信號轉(zhuǎn)導(dǎo)通路、疾病相關(guān)通路等。揭示基因在細(xì)胞整體的生化反應(yīng)和生理過程中的協(xié)同作用和調(diào)控關(guān)系。有助于理解基因在疾病發(fā)生發(fā)展、藥物作用機(jī)制等方面的作用。

綜合來看，GO 富集更側(cè)重于基因功能的分類描述，而 KEGG 富集則更側(cè)重于基因在生物系統(tǒng)中的通路級別的相互作用和調(diào)控。兩者結(jié)合可以更全面深入地理解基因的功能和在生物學(xué)系統(tǒng)中的角色

二、使用步驟

1.數(shù)據(jù)處理

library(tidyverse)

library(ggplot2)

library(ggrepel)

library(readxl)

setwd('G:/R/TCGA/venn')

# 這里的excel文檔是之前TCGA差異分析得到的數(shù)據(jù),主要用longFC,p值,FDR和基因名

COVIN_24_wpi_DEGs <- read_xlsx('long_CONVID_19_deg.xlsx',sheet='24 wpi DEGs')

COVIN_24_wpi_DEGs <- data.frame(COVIN_24_wpi_DEGs$logFC,COVIN_24_wpi_DEGs$PValue,COVIN_24_wpi_DEGs$FDR,COVIN_24_wpi_DEGs$external_gene_name,COVIN_24_wpi_DEGs$ensembl_gene_id,COVIN_24_wpi_DEGs$description,row.names = COVIN_24_wpi_DEGs$external_gene_name)

colnames(COVIN_24_wpi_DEGs) <- c('logFC','FDR','PValue','genename','gene_id','description')

COVIN_24_wpi_DEGs$time <- c('24 w.p.i.')

common_difference_gene_24_wpi <- read.csv('common_difference_gene_24_wpi.csv',row.names = 1)

selected_data_24_wpi <- COVIN_24_wpi_DEGs[COVIN_24_wpi_DEGs$genename %in% t(common_difference_gene_24_wpi), ]

threshold <- 0.3

selected_data_24_wpi$regulation <- ifelse(selected_data_24_wpi$logFC > threshold, "Up-Regulated", ifelse(selected_data_24_wpi$logFC < -threshold, "Down-Regulated", ifelse(abs(selected_data_24_wpi$logFC) < threshold, "NotSig", "Unknown")))

# 獲取上下調(diào)基因

top <- 100

top_genes_24_wpi <- bind_rows(

selected_data_24_wpi %>%

filter(regulation == 'Up-Regulated') %>%

arrange(desc(abs(logFC))) %>%

head(top)

)

bottom_genes_24_wpi <- bind_rows(

selected_data_24_wpi %>%

filter(regulation == 'Down-Regulated') %>%

arrange(desc(abs(logFC))) %>%

head(top)

)

# 保存數(shù)據(jù)之后方便之后的通路分析

library(openxlsx)

dir.create('differential gene/result/24_wpi/data_100')

write.xlsx(selected_data_24_wpi,'./differential gene/result/24_wpi/data_100/data_24_wpi.xlsx',rowNames = TRUE)

write.xlsx(top_genes_24_wpi,'./differential gene/result/24_wpi/data_100/top_genes_24_wpi.xlsx',rowNames = TRUE)

write.xlsx(bottom_genes_24_wpi,'./differential gene/result/24_wpi/data_100/bottom_genes_24_wpi.xlsx',rowNames = TRUE)

2. GO分析

setwd('G:/R/TCGA/venn/differential gene')

library(clusterProfiler)

# 這里org.Hs.eg.db需要通過BiocManager包進(jìn)行下載,如果多次嘗試下載不成功的惡化,可以先更改鏡像源,之后再嘗試下載

# options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")

# BiocManager::install('org.Hs.eg.db')

library(org.Hs.eg.db)

library(ggplot2)

library(readxl)

# 讀取數(shù)據(jù)

upgene <- read_xlsx('./result/24_wpi/data_100/up_genes_24_wpi.xlsx')

rownames(upgene) <- upgene$...1

Genes <- bitr(rownames(upgene),

fromType = 'SYMBOL',

toType = c('ENTREZID'),

OrgDb = org.Hs.eg.db)

# GO富集分析和結(jié)果存儲

GO <- enrichGO(gene = Genes$ENTREZID, # 輸入基因的ENTREZID

OrgDb = org.Hs.eg.db, # 注釋信息

keyType = 'ENTREZID',

ont = 'ALL',# 可選條目BP/CC/MF

pAdjustMethod = 'BH',# p值的校正方式

qvalueCutoff = 1,# q的閾值

pvalueCutoff = 1, # p的閾值

minGSSize = 5,

maxGSSize = 500,

readable = TRUE # 是否將ENTREZID轉(zhuǎn)換為symbol

)

dir.create('./result/24_wpi/up_100')

write.csv(data.frame(ID=rownames(GO@result),GO@result),'./result/24_wpi/up_100/up_genes_24_wpi_GO.csv')

# 柱狀圖

pdf(file = './result/24_wpi/up_100/up_genes_24_wpi_GO_barplot.pdf',width = 10,height = 8)

barplot(GO,drop=TRUE,

showCategory = 6,

split='ONTOLOGY') +

facet_grid(ONTOLOGY~.,scales = 'free')

dev.off()

# 點(diǎn)圖

pdf(file = './result/24_wpi/up_100/up_genes_24_wpi_GO_dotplot.pdf',width = 10,height = 8)

dotplot(GO,showCategory=6, split='ONTOLOGY')+

facet_grid(ONTOLOGY~.,scales = 'free')

dev.off()

# 不分面

pdf(file = './result/24_wpi/up_100/up_genes_24_wpi_GO_dotplot_nosplit.pdf',width = 10,height = 8)

dotplot(GO,

x='GeneRatio',

color='p.adjust',

showCategory=15,

size=NULL,

split=NULL,

font.size=13,

title='',

orderBy='x',

label_format=30)

dev.off()

3.KEGG富集

KEGG <- enrichKEGG(Genes$ENTREZID,

organism = 'hsa',

keyType = 'kegg',# kegg數(shù)據(jù)庫

pAdjustMethod = 'BH',

pvalueCutoff = 1,

qvalueCutoff = 1)

pdf(file = './result/24_wpi/up_100/up_genes_24_wpi_KEGG.pdf',width = 10,height = 8)

barplot(KEGG,

x='GeneRatio',

color='p.adjust',

showCategory=15)

dev.off()

pdf(file = './result/24_wpi/up_100/up_genes_24_wpi_KEGG_dotplot.pdf',width = 10,height = 8)

dotplot(KEGG)

dev.off()

# View(as.data.frame(KEGG)) # 以數(shù)據(jù)框的形式展示KEGG的結(jié)果

# browseKEGG(KEGG,'hsa03010') # 打開KEGG官網(wǎng)查看相關(guān)通路的信息

# 結(jié)果存儲

KEGG_results <- DOSE::setReadable(KEGG,

OrgDb = org.Hs.eg.db,

keyType = 'ENTREZID')

write.csv(data.frame(ID=rownames(GO@result),GO@result),'./result/24_wpi/up_100/up_genes_24_wpi_KEGG.csv')

總結(jié)

圖片很多,這里就不再一一展示了,有興趣的朋友可以自己常以跑一下,了解一下它們之間的區(qū)別,需要數(shù)據(jù)的話私我領(lǐng)取!!!

柚子快報邀請碼778899分享：R GO和kEGG富集分析

http://yzkb.51969.com/

參考鏈接