非模式生物富集背景库构建

返回

非模式生物富集背景库构建

发布于: 更新于:
统计加载中...

1 非模式生物富集背景库构建

1.1 一、蛋白质序列文件获取

1.1.1 从 NCBI 获取蛋白序列 (以拟南芥为例)

  1. 打开 NCBI 首页 National Center for Biotechnology Information，选择 Genome，输入目标种属拉丁文，最后点击 Search。(以拟南芥【Arabidopsis thaliana】为例)

  2. 选择基因组组装版本 Genome - NCBI - NLM

  3. 进入所选数据页面，选择下载 FTP 文件 Arabidopsis thaliana genome assembly TAIR10.1 - NCBI - NLM

  4. 在 GTF 下载界面，选择蛋白质序列文件，以.ffa.gz 为后缀 Index of /genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1

1.1.2 从 ensembl 获取蛋白序列 (以拟南芥为例)

  1. 打开 ensembl 官网链接

    1. 植物：Ensembl Plants
    2. 动物：Ensembl genome browser 115

  2. 搜索种属 (以拟南芥【Arabidopsis thaliana】为例)

  3. 进入目标种属数据页面，点击 Gene annotation 栏目中的 FASTAArabidopsis_thaliana - Ensembl Genomes 62

  4. 进入 FASTA 数据链接页面，下载.pep.all.fa.gz 后缀文件 Index of /pub/plants/release-62/fasta/arabidopsis_thaliana

1.2 二、使用 eggNOG-mapper 获取注释文件

  1. 打开 eggNOG-mappereggNOG-mapper
  2. 选择 Protein，上传获取的目标种属的蛋白质序列文件，输入可用的邮箱，用于接收进程信息以及处理后的注释文件
  3. 打开收到的邮件，点击 Click to manage your job
  4. 进入跳转页面后，点击 Start job 开始运行【本示例已运行完毕】
  5. 任务完成后，会发送新邮件提示任务已完成，打开邮件，点击 Access your job files here
  6. 点击下载 out.emapper.annotations 文件

1.3 三、使用 R 语言对注释文件进行预处理

1.3.1 加载所需包

setwd('~/Projects/BioInfo-CodeCollection/eggNOG_mapper_enrich/')

########## 1. 加载所需包 ##########
library(tidyverse)
library(clusterProfiler)
library(ontologyIndex)

1.3.2 读取数据

########## 2. 读取数据 ##########
degs <- read_tsv("./data/A-vs-B.DESeq2.total.DEG.tsv")
genelist <- degs$GeneID

1.3.3 eggNOG_mapper 注释结果整理

########## 3. eggNOG_mapper注释结果整理 ##########
# 读 eggNOG-mapper 输出，跳过注释行；自动修列名、去空格

emapper <- read_tsv(
  file = './data/out.emapper.annotations',
  comment = '##', # 忽略所有以 ## 开头的注释/日志行
  name_repair = 'universal', # 把非法列名（空格、-、/等）→下划线
  trim_ws = TRUE # 去除所有空格
) |>
  rename(GeneID = '.query') |> # 列名“.query”→GeneID（反引号必加）
  filter(str_detect(GeneID, '\\.1$')) |> #存在.01等，则只保留主转录本.1.不运行则保留所有，后续会去掉冗余.
  mutate(GeneID = str_remove(GeneID, '\\..*$')) #去掉.之后的所有，与基因组ID一致(若想做转录本富集，则不去掉)
head(emapper)

1.3.4 GO_term2gene[GO 编号与基因对应表]

########## 4. GO_term2gene[GO编号与基因对应表] ##########
GO_term2gene <- emapper |>
  select(GeneID, GOs) |>
  filter(GOs != "-") |>
  separate_rows(GOs, sep = ",") |>
  mutate(GOs = str_trim(string = GOs, side = 'both')) |> #去除字符串首尾空白字符
  rename(GO_ID = GOs) |>
  distinct(GO_ID, GeneID) # 去重
write_tsv(GO_term2gene, "./data/GO_term2gene.tsv")

1.3.5 GO_term2name[GO 编号与名称对应表]

########## 5. GO_term2name[GO编号与名称对应表] ##########--------------
if (!file.exists('./data/go_basic.obo')) {
  download.file(
    url = "http://purl.obolibrary.org/obo/go/go-basic.obo",
    destfile = './data/go_basic.obo'
  )
}

# 解析go_basic.obo
go.obo <- get_ontology(
  file = './data/go_basic.obo',
  propagate_relationships = NULL,
  extract_tags = "everything",
  merge_equivalent_terms = TRUE

)

summary(go.obo)
head(go.obo$name)
head(names(go.obo$name))
str(go.obo$name)

GO_term2name <- data.frame(
  GO_ID = names(go.obo$name), # GO编号
  GO_name = as.character(go.obo$name), # GO词条名称
  stringsAsFactors = FALSE
) 
write_tsv(GO_term2name, './data/GO_term2name.tsv')

1.3.6 KEGG_term2gene[KEGG 通路与基因对应表]

########## 6. KEGG_term2gene[KEGG通路与基因对应表] ##########
KEGG_term2gene <- emapper |>
  select(KEGG_Pathway, GeneID) |>
  filter(KEGG_Pathway != "-") |>
  separate_rows(KEGG_Pathway, sep = ",") |>
  filter(!grepl(pattern = "^map", x = KEGG_Pathway)) |> # 去掉map开头（非ko通路）
  rename(GeneID = GeneID, Pathway_ID = KEGG_Pathway)
write_tsv(KEGG_term2gene, './data/KEGG_term2gene.tsv')

1.3.7 KEGG_term2name[KEGG 通路与名称对应表]

########## 7. KEGG_term2name[KEGG通路与名称对应表] ##########
KEGG_term2name <- read_tsv(
  file = "https://rest.kegg.jp/list/pathway",
  col_names = c("Pathway_ID", "Pathway_Name")
) |>
  mutate(
    # Pathway_ID = str_remove(Pathway_ID, "path:"), # 去掉 path:
    Pathway_ID = str_replace(Pathway_ID, "^map", "ko") # map→ko
  ) 
write_tsv(KEGG_term2name, './data/KEGG_term2name.tsv')

1.3.8 GO 富集分析

########## 8. GO富集分析 ##########
GO_term2gene <- read_tsv('./data/GO_term2gene.tsv')
GO_term2name <- read_tsv('./data/GO_term2name.tsv')
GO_filtered_genes <- genelist[genelist %in% GO_term2gene$GeneID]  

# 用 clusterProfiler 的 enricher 做超几何检验
GO_enrich_result <- enricher(
  gene = GO_filtered_genes,
  pvalueCutoff = 0.05, # 原始 p 值阈值
  qvalueCutoff = 0.2, # BH 校正后 q 值阈值
  pAdjustMethod = "BH", # 多重检验校正方法
  TERM2GENE = GO_term2gene,
  TERM2NAME = GO_term2name
)

# 解析obo文件
go.obo <- get_ontology(
  file = "./data/go_basic.obo",
  extract_tags = "everything"
)

# 生成 GO_ID→ONTOLOGY（BP/CC/MF）映射表
go2ns <- tibble(
  ID = names(go.obo$namespace),
  ONTOLOGY = as.character(go.obo$namespace)
) |>
  mutate(
    ONTOLOGY = if_else(is.na(ONTOLOGY) | ONTOLOGY == "", "obsolete", ONTOLOGY)
  ) |>
  distinct() # 去重
  
# 将命名空间并入富集结果，便于按 BP/CC/MF 分面绘图
GO_enrich_result@result <- GO_enrich_result@result |>
  left_join(go2ns, by = c("ID" = "ID")) %>% # 按 GO ID 连接
  filter(!is.na(ONTOLOGY)) %>% # 去掉未匹配上的条目
  relocate(ONTOLOGY, .before = everything()) # 把 ONTOLOGY 放第一列，查看方便
  
# 转成 data.frame 并写出，方便 Excel 打开或后续 ggplot 自用
go_res_df <- as.data.frame(GO_enrich_result)
write_tsv(go_res_df, "./data/GO_enrich.tsv") # 结果文件：GO_enrich.tsv

1.3.9 KEGG 富集分析

########## 9. KEGG富集分析 ##########
KEGG_term2gene <- read_tsv("./data/KEGG_term2gene.tsv")
KEGG_term2name <- read_tsv("./data/KEGG_term2name.tsv")
  
KEGG_filtered_genes <- genelist[genelist %in% KEGG_term2gene$GeneID]
  
# 用 enricher 做 KEGG 超几何检验
KEGG_enrich_result <- enricher(
  gene = KEGG_filtered_genes,
  pvalueCutoff = 0.05, # 原始 p 值阈值
  qvalueCutoff = 0.2, # BH 校正后 q 值阈值
  pAdjustMethod = "BH", # 多重检验校正
  TERM2GENE = KEGG_term2gene,
  TERM2NAME = KEGG_term2name
)

1.3.10 KEGG 通路分类

########## 10. KEGG通路分类 ##########
# 从 KEGG 官网下载通路分类表
br_raw <- readLines("https://rest.kegg.jp/get/br:br08901")
ko_tbl <- tibble(line = br_raw) %>%
  mutate(A = str_match(line, "^A([A-Za-z\\s]+)$")[, 2] %>% str_trim()) %>%
  fill(A) %>%
  mutate(B = str_match(line, "^B\\s+(.+)$")[, 2] %>% str_trim()) %>%
  fill(B) %>%
  mutate(
    map = str_match(line, "^C\\s{2,}(\\d{5})\\s{2,}")[, 2],
    descr = str_match(line, "^C\\s{2,}\\d{5}\\s{2,}(.+)$")[, 2]
  ) %>%
  filter(!is.na(map)) %>%
  mutate(ko = str_c("ko", map)) %>%
  select(ko, category = A, subcategory = B, description = descr)
  
KEGG_enrich_result@result <- KEGG_enrich_result@result %>%
  left_join(ko_tbl, by = c("ID" = "ko")) %>%
  filter(!is.na(category)) %>% # 去掉未匹配上的通路
  filter("category" != "Human Diseases") %>% #过滤掉人类疾病
  filter("category" != "Drug Development") %>% #过滤掉药物开发
  relocate(category, subcategory, .before = 1)

kegg_res_df <- as.data.frame(KEGG_enrich_result)
write_tsv(kegg_res_df, "./data/KEGG_enrich.tsv")

1.3.11 绘图

########## 11. 绘图 ##########

# 设置全局主题
theme_set(
  theme_bw(base_family = "sans") + # 白底网格主题 + 无衬线字体
    theme(
      text = element_text(face = "bold"), # 全部文字加粗
      axis.text.y = element_text(size = 12), # Y 轴刻度字号 12
      axis.text.x = element_text(size = 10), # X 轴刻度字号 10
      axis.title.x = element_text(size = 14), # X 轴标题字号 14
      plot.title = element_text(size = 18, hjust = 0.5), # 图标题 18 号并水平居中
      legend.text = element_text(size = 10) # 图例文字 10 号
    )
)
  
# GO 条形图
p1 <- barplot(
  GO_enrich_result,
  x = "Count",
  showCategory = 6, # 每类显示 6 条
  split = "ONTOLOGY", # 按命名空间拆分
  color = "p.adjust", # 颜色对应校正后 p 值
  label_format = 50, # 标签最长 50 字符
  title = "GO Enrichment Analysis"

) +
  facet_grid(ONTOLOGY ~ ., scale = 'free') # 三分面，坐标轴自由缩放

# GO 气泡图
p2 <- dotplot(
  GO_enrich_result,
  showCategory = 6, # 每类（BP/CC/MF）各显 6 条
  split = "ONTOLOGY", # 依据 ONTOLOGY 列拆分成 3 张图
  color = "p.adjust", # 颜色映射 BH 校正后 p 值
  label_format = 50, # 通路名字符上限 50，超长自动换行
  title = "GO Enrichment Analysis"
) +
  facet_grid(ONTOLOGY ~ ., scale = 'free') # 三分面，Y 轴长度随条数自适应
print(p1)
print(p2)

# 保存
ggsave("GO_barplot.pdf", p1, width = 8, height = 8) # 条形图（前面已生成 p1）
ggsave("GO_dotplot.pdf", p2, width = 8, height = 8) # 气泡图  

# KEGG 条形图
p3 <- barplot(
  KEGG_enrich_result,
  x = "Count", # 横轴用基因数
  showCategory = 20, # 显示校正后最显著的 20 条通路
  color = "p.adjust", # 颜色按 p.adjust 渐变
  label_format = 50, # 通路名字符上限 50
  title = "KEGG Enrichment Analysis"
)

# KEGG 气泡图
p4 <- dotplot(
  KEGG_enrich_result,
  showCategory = 20, # 同样 top18
  color = "p.adjust", # 气泡颜色映射 p.adjust
  label_format = 50, # 通路名长度控制
  title = "KEGG Enrichment Analysis"
)

# 保存
ggsave("KEGG_barplot.pdf", p3, width = 8, height = 8) # 条形图
ggsave("KEGG_dotplot.pdf", p4, width = 8, height = 8) # 气泡图

评论