非模式生物自建库后进行GSEA分析

最新推荐文章于 2024-09-14 17:22:51 发布

刘融晨

最新推荐文章于 2024-09-14 17:22:51 发布

阅读量871

点赞数 7

文章标签： r语言

本文链接：https://blog.csdn.net/m0_73361516/article/details/135431704

版权

library(readr)
annotations <- read_csv("~/Program/crab/advantage/data/annotations.csv") %>% 
  mutate(Gene_Name = if_else(Gene_Name != "-",
                             Gene_Name,
                             gene_id))

gene_info = annotations %>% dplyr::select(GID = gene_id,Gene_Name) %>% unique()

# 将基因与GO的对应关系整理出来
gene2go <- dplyr::select(annotations,GID = gene_id,GO = GOs) %>%
  separate_rows(GO, sep = ',', convert = F) %>%
  filter(GO!="-",!is.na(GO)) %>% 
  mutate(EVIDENCE = 'A') %>% 
  unique()
gene2go %>% head()

# 将基因与pathway的对应关系整理出来
gene2pathway <- dplyr::select(annotations,GID = gene_id,Pathway = KEGG_Pathway) %>%
  separate_rows(Pathway, sep = ',', convert = F) %>%
  filter(str_detect(Pathway, 'map')) %>% 
  unique()%>%na.omit()
pathway2gene %>% head()

## 准备GO注释内容
library(GO.db)
goterms <- Term(GOTERM)
GOlist=as.data.frame(goterms)
go2name <- rownames_to_column(GOlist,var = 'GO_id')
go2name %>% head()

## 这里提供两种获得KEGG注释的方法，kEGG数据库经常抽风(这里用不到)
library(magrittr)
get_path2name <- function(){
  keggpathid2name.df <- clusterProfiler:::kegg_list("pathway")
  keggpathid2name.df[,1] %<>% gsub("path:map", "", .)
  colnames(keggpathid2name.df) <- c("path_id","path_name")
  return(keggpathid2name.df)
}
pathway2name <- get_path2name()

AnnotationForge::makeOrgPackage(gene_info=gene_info,
                                go=gene2go,
                                ko = gene2pathway,
                                maintainer='rongchen.liu <rongchen.liu@foxmail.com>',
                                author='rongchen.liu',
                                version="0.1" ,
                                outputDir="./", 
                                tax_id="95602",
                                genus="C",
                                species="rab",
                                goTable = 'go')
install.packages('./org.Crab.eg.db/',repos = NULL, type="source")

library(org.Crab.eg.db)
library(clusterProfiler)

GSEA_input <- de_result$log2FoldChange
names(GSEA_input) = de_result$Row.names
GSEA_input = sort(GSEA_input, decreasing = TRUE)


### 2.1 比较 受精卵大规格（7公5母）与 受精卵小规格（3公2母）-1
################################################################################
GSEA_input1 <- de_result1$log2FoldChange
names(GSEA_input1) = de_result1$Row.names
GSEA_input1 = sort(GSEA_input1, decreasing = TRUE)

Go_gseresult_bp1 <- gseGO(GSEA_input1,
                         'org.Crab.eg.db',
                         keyType = "GID",
                         ont="BP",
                         nPerm = 1000,
                         minGSSize = 10,
                         maxGSSize = 1000,
                         pvalueCutoff= 1)
Go_gseresult_mf1 <- gseGO(GSEA_input1,
                         'org.Crab.eg.db',
                         keyType = "GID",
                         ont="MF",
                         nPerm = 1000,
                         minGSSize = 10,
                         maxGSSize = 1000,
                         pvalueCutoff= 1)

Go_gseresult_cc1 <- gseGO(GSEA_input1,
                         'org.Crab.eg.db',
                         keyType = "GID",
                         ont="CC",
                         nPerm = 1000,
                         minGSSize = 10,
                         maxGSSize = 1000,
                         pvalueCutoff= 1)

Go_gseresult_all1 <- gseGO(GSEA_input1,
                          'org.Crab.eg.db',
                          keyType = "GID",
                          ont="ALL",
                          nPerm = 1000,
                          minGSSize = 10,
                          maxGSSize = 1000,
                          pvalueCutoff= 1)
tpm = Go_gseresult_all1 %>% as.data.frame() %>% filter(pvalue <= 0.05 & qvalue <= 0.25) 

library(enrichplot)

gseaplot2(Go_gseresult_all1, 1:5, pvalue_table = F)