不用Python不用R？一行命令搞定它——获取KEGG通路基因列表

原创西瓜站长生信益站 2023-08-13 08:00

生信益站，一点就有益！祝友友们天天开心，早日发 CNS~

咱们接着上一讲：

你的 GSEA 做对了吗？MSigDB 的 C2 通路只有 186 个，而 KEGG 收录的通路却有 355 个？

https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp#C2

俗话说，授人以鱼不如授人以渔。

今天站长就给站友们解密——如何 DIY 自己的 KEGG 基因集/GMT 格式。DIY的意义：

GSVA分析；
GSEA分析；
研究特定的通路，比如拿你感兴趣的基因与特定通路取交集等。

获取通路列表

这里以人/hsa为例，提取人类所有的 KEGG通路id和通路name，可以看到一共有355个：

curl -s https://rest.kegg.jp/list/pathway/hsa -o hsa.ko2pathway

head -5 hsa.ko2pathway
# hsa01100  Metabolic pathways - Homo sapiens (human)
# hsa01200  Carbon metabolism - Homo sapiens (human)
# hsa01210  2-Oxocarboxylic acid metabolism - Homo sapiens (human)
# hsa01212  Fatty acid metabolism - Homo sapiens (human)
# hsa01230  Biosynthesis of amino acids - Homo sapiens (human)

less hsa.ko2pathway |wc -l
# 355

sed -i 's# - Homo sapiens (human)##' hsa.ko2pathway
head -5 hsa.ko2pathway
# hsa01100   Metabolic pathways
# hsa01200   Carbon metabolism
# hsa01210   2-Oxocarboxylic acid metabolism
# hsa01212   Fatty acid metabolism
# hsa01230   Biosynthesis of amino acids

提取通路内的基因列表

下面代码展示了如何提取hsa00010通路内的基因列表：

# 下载`hsa00010`通路的基因/KO信息
# 返回文件为`json`格式
curl -s http://togows.dbcls.jp/entry/pathway/hsa00010/genes.json \
    | awk 'BEGIN{OFS="\t"}$0~/KO:/{
        match($0,/"([^"]+)": "([^;[]+);?(.+?) \[(KO:[^\]]+)/,a);
        a[3]=a[3]==""?"-":a[3];
        print a[1],a[2],a[3],a[4]
        }' \
    | head -10

# 修修改改就可以整理成GSEA富集分析所需的基因集文件/gmt格式
# entreze-ID  Gene-Symbol  Description  KO-ID
3101    HK3      hexokinase 3                      KO:K00844
3098    HK1      hexokinase 1                      KO:K00844
3099    HK2      hexokinase 2                      KO:K00844
80201   HKDC1    hexokinase domain containing 1    KO:K00844
2645    GCK      glucokinase                       KO:K12407
2821    GPI      glucose-6-phosphate isomerase     KO:K01810
5213    PFKM     phosphofructokinase, muscle       KO:K00850
5214    PFKP     phosphofructokinase, platelet     KO:K00850
5211    PFKL     phosphofructokinase, liver type   KO:K00850
2203    FBP1     fructose-bisphosphatase 1         KO:K03841

在浏览器上打开这个链接：

http://togows.dbcls.jp/entry/pathway/hsa00010/

可以发现hsa00010除了gene列表，还可以提取module、drug，cpd等列表：

提取modules的代码如下（大家自行对比区别，其实很简单）：

curl http://togows.dbcls.jp/entry/pathway/hsa00010/modules.json
[
  {
    "hsa_M00001": "Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate [PATH:hsa00010]",
    "hsa_M00002": "Glycolysis, core module involving three-carbon compounds [PATH:hsa00010]",
    "hsa_M00003": "Gluconeogenesis, oxaloacetate => fructose-6P [PATH:hsa00010]",
    "hsa_M00307": "Pyruvate oxidation, pyruvate => acetyl-CoA [PATH:hsa00010]"
  }
]

提取compounds的代码如下（compounds：代谢物）：

curl http://togows.dbcls.jp/entry/pathway/hsa00010/compounds.json
[
  {
    "C00022": "Pyruvate",
    "C00024": "Acetyl-CoA",
    "C00031": "D-Glucose",
    "C00033": "Acetate",
    "C00036": "Oxaloacetate",
    "C00068": "Thiamin diphosphate",
    "C00074": "Phosphoenolpyruvate",
    "C00084": "Acetaldehyde",
    "C00103": "D-Glucose 1-phosphate",
    "C00111": "Glycerone phosphate",
    "C00118": "D-Glyceraldehyde 3-phosphate",
    "C00186": "(S)-Lactate",
    "C00197": "3-Phospho-D-glycerate",
    "C00221": "beta-D-Glucose",
    "C00236": "3-Phospho-D-glyceroyl phosphate",
    "C00267": "alpha-D-Glucose",
    "C00469": "Ethanol",
    "C00631": "2-Phospho-D-glycerate",
    "C00668": "alpha-D-Glucose 6-phosphate",
    "C01159": "2,3-Bisphospho-D-glycerate",
    "C01172": "beta-D-Glucose 6-phosphate",
    "C01451": "Salicin",
    "C05125": "2-(alpha-Hydroxyethyl)thiamine diphosphate",
    "C05345": "beta-D-Fructose 6-phosphate",
    "C05378": "beta-D-Fructose 1,6-bisphosphate",
    "C06186": "Arbutin",
    "C06187": "Arbutin 6-phosphate",
    "C06188": "Salicin 6-phosphate",
    "C15972": "Enzyme N6-(lipoyl)lysine",
    "C15973": "Enzyme N6-(dihydrolipoyl)lysine",
    "C16255": "[Dihydrolipoyllysine-residue acetyltransferase] S-acetyldihydrolipoyllysine"
  }
]

知识拓展

KEGG MODULE 数据库由M号标识的 KEGG 模块和 RM 号标识的 KEGG 反应模块组成，分别是手动定义的基因集和反应集的功能单元。KEGG 模块进一步分为通路模块和特征模块，如下所示：

途径模块——代谢途径中基因组的功能单元，包括分子复合物
特征模块——表征表型特征的基因集功能单元
反应模块——代谢途径中连续反应步骤的功能单元
https://www.genome.jp/kegg/module.html

KEGG DRUG是日本、美国和欧洲批准药物的综合药物信息资源，基于化学结构和/或活性成分的化学成分进行统一。每个 KEGG DRUG 条目均由D号标识，并与 KEGG 原始注释相关联，包括治疗靶点、药物代谢和其他分子相互作用网络信息。

https://www.genome.jp/kegg/drug/

KEGG COMPOUND是小分子、生物聚合物和其他与生物系统相关的化学物质的集合。每个条目均由 C 编号标识，例如 L-赖氨酸为 C00047，并包含化学结构和相关信息，以及指向其他 KEGG 数据库和外部数据库的各种链接。一些 COMPOUND 条目还表示为带有“相同”链接的 GLYCAN 和 DRUG 条目。KEGG COMPOUND 中代表性条目的分类在以下 BRITE 层次结构文件中给出。

https://www.genome.jp/kegg/compound/

批量下载通路基因集

❝
因为我们要制作KEGG基因集进行GSEA富集，因此只下载基因列表即可。其实 y 叔的clusterProfiler也提供了针对 KEGG 模块 module 的富集分析方法enrichMKEGG/gseMKEGG。
❞

批量下载所有通路的基因集：

cat hsa.ko2pathway| while IFS=$'\t' read -r id name
do
   gene_set=`curl -s http://togows.dbcls.jp/entry/pathway/${id}/genes.json \
      | awk 'BEGIN{OFS="\t"}$0~/KO:/{match($0,/"([^"]+)": "([^;[]+);?(.+?) \[(KO:[^\]]+)/,a); print a[1]}' \
      | paste -d"\t" -s`
   echo -e "${id}\t${name}\t${gene_set}" >> kegg.hsa.gmt
done

可以看出前 8 个通路只有 2 列，没有基因列表。分析的时候，可以把他们删除，这些都是正经的纯代谢通路，可以丢掉（剩 355-8=347 个）。

head -10 kegg.hsa.gmt
# hsa01100   Metabolic pathways
# hsa01200   Carbon metabolism
# hsa01210   2-Oxocarboxylic acid metabolism
# hsa01212   Fatty acid metabolism
# hsa01230   Biosynthesis of amino acids
# hsa01232   Nucleotide metabolism
# hsa01250   Biosynthesis of nucleotide sugars
# hsa01240   Biosynthesis of cofactors
# hsa00010   Glycolysis / Gluconeogenesis 3101 3098 3099 80201 2645 2821 5213 5214 5211 2203 8789 230 226 229 7167 2597 26330 5232 5230 5223 5224 441531 2027 2026 2023 387712 5315 5313 5161 5160 5162 1737 1738 160287 92483 3939 3945 3948 124 125 126 131 127 128 130 10327 217 224 219 501 223 221 222 218 84532 55902 130589 5236 55276 2538 57818 92579 83440 669 9562 5105 5106
# hsa00020   Citrate cycle (TCA cycle) 1431 47 50 48 3417 3418 3420 3421 3419 55753 4967 1743 1738 8802 8801 8803 6389 6390 6391 6392 2271 4190 4191 5091 5105 5106 5161 5160 5162 1737

友情提示

如果纯粹做 KEGG 的 GSEA 富集分析，使用clusterProfiler::gseKEGG()即可（基因集是自动下载的，所以理论上是完整滴），所以GSEA大部分情况是不需要自己整理GMT基因集滴。如果你有自制的GMT基因集，可以使用clusterProfiler::GSEA(geneList, TERM2GENE = your_gmt_object)，也可以使用博德研究所java版本的GESA软件。教程可以参考这篇文章：

7 篇文章带你学会 GSEA 富集分析！

怎么样，学会了吗？是不是很简单...

OK，今天的分享到此为止，希望能对您有所帮助。

您的关注、点赞、在看、转发是对益站最大的鼓励和支持哈。

联系站长

❝
对本篇文章有疑问，可以在益站发消息留言，也欢迎各位童鞋扫码加入我们的 QQ 交流群。
❞

喜欢此内容的人还喜欢

微信扫一扫
关注该公众号