生信技能树的生信入门课程，第四周转录组测序数据分析给大家讲解了如何在KEGG数据库中标记显著差异表达基因，使用的是网页版操作，详细过程可以看看我们的优秀学员笔记：优秀学员笔记：转录组课堂笔记Day5

PPT课件中的对应的内容：

标记效果如下：

本次代码版本使用的是一个R包 ggkegg ：https://github.com/noriakis/ggkegg

完整版的说明文档在这里：https://noriakis./software/ggkegg/

0.安装R包

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors./bioconductor")
options("repos"=c(CRAN="https://mirrors./CRAN/"))

## 去github上下载下来压缩包
devtools::install_local("ggkegg-main.zip")
install.packages("ggfx")

rm(list=ls())
## 加载包
library(ggfx)
library(cols4all)
library(tidygraph)
library(dplyr)
library(igraph)
library(ggkegg)

绘图

这个包有很多出色的绘图功能，下面来看看！

1.在某通路中标记一个基因

如在 hsa04110（Cell cycle – Homo sapiens (human)，细胞周期通路） 通路中突出标记基因 CDKN2A：

fill_color：指定标记基因的颜色；

directory=”./”：设置下载的通路xml文件路径，默认为临时路径

## Highlight genes in the pathway, and overlay raw map.
p <- highlight_entities("hsa04110", "CDKN2A",directory="./",fill_color = "tomato",use_cache = T)
p
ggsave("hsa04110.png",plot = p, bg="white",width = 9, height = 8)

结果如下：

2.同时标记多个基因：

给一个基因名字向量即可：

# 同时标记多个基因
p <- highlight_entities("hsa04110", c("CDKN2A", "CDC45"),directory="./",
                        fill_color = "tomato",use_cache = T)
p

3.还可以标记连续值

输入一个带名字的数值向量，可以用于展示差异基因的FC值，pvalue值，相关性等：

# 使用连续值标注
## Highlight genes in the pathway, and overlay raw map.
vecs <- c(-5,-2,2,4,5,7,8)
names(vecs) <- c("CDKN2A", "CDC45","PDPK1","PIK3CA","GSK3B","PIK3R6","THEM4")
vecs

cs <- highlight_entities("hsa04110", vecs) + 
  scale_fill_viridis(alpha = 0.5)
cs

结果如下：

4.保存图片

可以使用 ggkeggsave 保存图片，原始分辨率，我前面用了ggsave，来看看这个函数：

# 保存
ggkeggsave(filename="test.pdf", cs, dpi=300)
knitr::include_graphics("test.pdf")
kegg pathway通路高亮标记基因（代码版）

5.映射差异分析结果

可以使用 assign_deseq2 函数来指定要反映在图中的数值（例如，log2FoldChange）作为列参数，可以将该值分配给节点。如果多个基因被击中，numeric_combine参数指定如何组合多个值（默认为平均值）。

示例数据来自https:///10.1038/s41388-022-02235-8，该数据集分析了感染BK多瘤病毒的人类尿路上皮细胞的转录组变化（Baker等人，2022年）。

准备好的示例数据 uro.deseq.res.rda 在这里下载：https://github.com/noriakis/ggkegg/issues/8

# 映射差异分析结果
library(ggkegg)
library(DESeq2)
library(org.Hs.eg.db)
library(dplyr)

## The file stores DESeq() result on transcriptomic dataset deposited by Baker et al. 2022.
load("uro.deseq.res.rda") 
res

vinf <- results(res, contrast=c("viral_infection","BKPyV (Dunlop) MOI=1","No infection"))
head(vinf)


## LFC
g <- pathway("hsa04110") |> mutate(deseq2=assign_deseq2(vinf, org_db=org.Hs.eg.db),
                                   padj=assign_deseq2(vinf, column="padj", org_db=org.Hs.eg.db),
                                   converted_name=convert_id("hsa"))

g

ggraph(g, layout="manual", x=x, y=y) + 
  geom_edge_parallel(width=0.5, arrow = arrow(length = unit(1, 'mm')), 
                     start_cap = square(1, 'cm'),
                     end_cap = square(1.5, 'cm'), aes(color=subtype_name)) +
  geom_node_rect(aes(fill=deseq2, filter=type=="gene"), color="black") +
  ggfx::with_outer_glow(geom_node_text(aes(label=converted_name, 
                                           filter=type!="group"), size=2.5), colour="white", expand=1)+
  scale_fill_gradient(low="lightblue",high="red", name="LFC") +
  theme_void()

也可以是pvalue：

## Adjusted p-values
ggraph(g, layout="manual", x=x, y=y) + 
  geom_edge_parallel(width=0.5, arrow = arrow(length = unit(1, 'mm')), 
                 start_cap = square(1, 'cm'),
                 end_cap = square(1.5, 'cm'), aes(color=subtype_name))+
  geom_node_rect(aes(fill=padj, filter=type=="gene"), color="black")+
  ggfx::with_outer_glow(geom_node_text(aes(label=converted_name, filter=type!="group"), size=2.5), colour="white", expand=1)+
  scale_fill_gradient(name="padj")+
  theme_void()

今天分享到这~

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：