Gene set and pathway

Differential expression (DE) analysis typically yields a list of genes or proteins. Our intention is to use such lists to gain novel insights about genes and proteins that may have roles in a given phenomenon, phenotype or disease progression. However, in many cases, gene lists generated from DE analysis are difficult to interpret due to their large size and lack of useful annotations. Hence, pathway analysis (also known as gene set analysis or over-representation analysis), aims to reduce the complexity of interpreting gene lists via mapping the listed genes to known (i.e. annotated) biological pathways, processes and functions. This learning submodule introduces common curated biological databases including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG).

Learning Objectives:

1. Introduction to Ontology and Gene Ontology.

2. Introduction to KEGG Pathway Database.

3. Download terms, pathway gene set from GO and KEGG.

4. Save results to GMT file format.

Ontology and Gene Ontology

Overview

In this section we will learn about the concept of gene ontology in Bioinformatics. Ontology is set of concepts and categories defined by a shared vocabulary to denote properties of the concepts, as well as the relationships between the concepts. Ontology plays an important role in the field of bioinformatics. Ontology enables unambiguous communication e.g., a way to understand different groups’ annotations of various genomes. Also, it allows the knowledge to be structured to perform automated analyses by computer programs.

The Gene Ontology (GO) database defines a structured, common, and controlled vocabulary to describe attributes of genes and gene products across organisms. Collaboration is key to build a consensus vocabulary. But the term gene ontology, or GO, is commonly used to refer to both the terms as well as the associations between genes, which is sometimes a source of potential confusion. In order to avoid this, here we will use the term “GO” to describe the set of terms and their hierarchical structure and “GO annotations” to describe the set of associations between genes and GO terms. The GO is divided into three categories to describe the genes and gene products from three different angles: Molecular Function, Biological Process, and Cellular Component.

The structure of GO can be described in terms of into directed acyclic graphs (DAGs), where each GO term is a node, and the relationships between the terms are edges between the nodes. GO is loosely hierarchical, with ‘child’ terms being more specialized than their ‘parent’ terms, but unlike a strict hierarchy, a term may have more than one parent term (note that the parent/child model does not hold true for all types of relations). The structure of the controlled vocabularies are intended to reflect true, biological relationships. In contrast to strict hierarchies, DAGS allow multiple relationships between a more granular (child) term and a more general parent term. The relationship between terms affects how queries are made. For example,a query for all genes with binding activity would include transcription factors as well as genes with other types of binding activity (such as protein binding, ligand binding). The illustration of category and structure of GO is shown in the figure below:

(Source: https://www.ebi.ac.uk/, http://geneontology.org/)

Gene ontology relationship

In DAGs graph, terms are represented as nodes and relations (also known as object properties) between the terms are edges. There are commonly used relationships in GO such as is a (is a subtype of), part of, has part, regulates, negatively regulates and positively regulates. All terms (except from the root terms representing each aspect) have a sub-class relationship to another terms.

Examples:

Example: **GO:1904659:glucose transport** *is* a **GO:0015749:monosaccharide transport**.

The is a relation forms the basic structure of GO. If we say A is a B, we mean that node A is a subtype of node B

Example: **GO:0031966:mitochondrial membrane** *is part of* **GO:0005740:mitochondrial envelope**

The part of relation is used to represent part-whole relationships. A part of relation would only be added between A and B if B is necessarily part of A: wherever B exists, it is as part of A, and the presence of the B implies the presence of A. However, given the occurrence of A, we cannot say for certain that B exists.

Example: **GO:0098689:latency-replication decision** *regulates* **GO:0019046:release from viral latency**

A relation that describes case in which one process directly affects the manifestation of another process or quality, i.e. the former regulates the latter.

A more specific case with more nodes and edges can be seen at the figure below:
(Source: https://advaitabio.com/)
For more technical information about relations and their properties used in GO and other ontologies see the OBO Relations Ontology (RO)

GO storage file formats

GO terms are updated monthly in the following formats:

• OBO 1.4 files are human-readable (in addition to machine-readable) and can be opened in any text editor.

• OWL files can be read by Protégé text editor.

In this learning submodule, we will only use “.OBO” to obtain GO terms.The OBO file format is for representing ontologies and controlled vocabularies. The format itself attempts to achieve the following goals:

• Human readability

• Ease of parsing

• Extensibility

• Minimal redundancy

The file structure can is shown in the following figure.

The OBO file has a header, which is an unlabeled section at the beginning of the document. The header ends when the first term is encountered. Next, term is represented in labeled section with the tag [Term]. Under each term, we can find other information such as term ID, official name, category (namespace), term definition, synonym and relation to other GO terms.

At this step, we still don’t know what genes are related to which GO terms. In order to retrieve custom sets of gene ontology annotations for any list of genes from organisms, NCBI has published a Gene2GO database that obtain GO terms and the entrez gene ids related to those go terms. The database can be retrieve from here text editor. The Gene2GO database can be viewed using text editor, the file structure is presented in the figure below:

The OBO and Gene2GO databases will be used in combination to obtain GO term and related genes for enrichment analysis.

Retrieving GO terms from DE gene list

This section focuses on downloading related GO terms based on the DE genelist obtained from the DE analysis in the previous section. Here, we will use topGO and hgu133plus2.db R packages to obtain GO terms. The topGO package has built-in functions that use Gene2GO databases to retrieve GO terms from the gene ID give by DE analysis. Since the dataset we used in submodule 02 was generated for human, we will use hgu133plus2.db database to map probeIDs to gene symbols. The installation process of two package can be done by the script below:

# Installation of topGO and hgu133plus2.db package
suppressMessages({if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
  suppressWarnings(BiocManager::install("topGO", update = F))
  suppressWarnings(BiocManager::install("hgu133plus2.db", update = F))
  suppressWarnings(BiocManager::install("AnnotationDbi", update = F))
})

# Loading the library
suppressPackageStartupMessages({
  suppressWarnings(library("topGO"))
  suppressWarnings(library("hgu133plus2.db"))
  suppressWarnings(library("AnnotationDbi"))
})

groupGOTerms: 	GOBPTerm, GOMFTerm, GOCCTerm environments built.

Load the DE genelist generated from the DE analysis using limma.

# Loading the DE result
data = readRDS("./data/DE_genes.rds")

By default, the DE analysis performed by limma contains multiple features. However, adjusted p-value and gene ID are the most important features for enrichment analysis. We can use the following code to list of gene IDs and their p-value.

# Get p-value from DE results
genelist <- data$adj.P.Val
# Assign gene IDs to associated p-values
names(genelist) <- data$PROBEID

After successfully obtaining the genelist, we need to map the gene IDs to the gene symbols using hgu133plus2.db.

# Map gene IDs to gene symbols
gene <- suppressMessages(AnnotationDbi::select(hgu133plus2.db, names(genelist), "SYMBOL"))

# Remove duplicated gene IDs
gene <- gene[!duplicated(gene[,1]),]

# Assign result to a new genlist with gene symbols
geneList2 <- genelist
names(geneList2) <- gene[,2]

Now, we can search for related GO terms based on the new gene list using topGO package. First, we need to create a topGOdata object.

# Retrieve all the GO terms related to the genelist obtained from the expression matrix
GOdata <- new("topGOdata", description = "",ontology = "BP",
                    allGenes = geneList2, geneSel = function(x)x, nodeSize = 10,
                    annot = annFUN.org, ID = "alias", mapping = "org.Hs.eg")

Building most specific GOs .....

	( 12457 GO terms found. )

Build GO DAG topology ..........

	( 15840 GO terms and 35901 relations. )

Annotating nodes ...............

	( 16332 genes annotated to the GO terms. )

We can search for related GO terms using geneInTerm function and view the term with associated genes.

# Obtain a list of genes for each GO term
allGO = genesInTerm(GOdata)
allGO[1:5]

$`GO:0000002`

1. 'AKT3'
2. 'DNA2'
3. 'DNAJA3'
4. 'ENDOG'
5. 'FLCN'
6. 'LIG3'
7. 'LONP1'
8. 'MDP1'
9. 'MEF2A'
10. 'METTL4'
11. 'MGME1'
12. 'MPV17'
13. 'OPA1'
14. 'PARP1'
15. 'PIF1'
16. 'PIM1'
17. 'POLB'
18. 'POLG'
19. 'POLG2'
20. 'PPARGC1A'
21. 'PRIMPOL'
22. 'RRM1'
23. 'RRM2B'
24. 'SESN2'
25. 'SLC25A33'
26. 'SLC25A36'
27. 'SLC25A4'
28. 'SSBP1'
29. 'STOML2'
30. 'STOX1'
31. 'TOP3A'
32. 'TP53'
33. 'TWNK'
34. 'TYMP'

$`GO:0000003`

1. 'A1CF'
2. 'A2M'
3. 'AAAS'
4. 'ABAT'
5. 'ABCC8'
6. 'ABHD2'
7. 'ACE'
8. 'ACE2'
9. 'ACOD1'
10. 'ACOX1'
11. 'ACR'
12. 'ACRBP'
13. 'ACRV1'
14. 'ACSBG1'
15. 'ACSBG2'
16. 'ACSL4'
17. 'ACTL7A'
18. 'ACTL9'
19. 'ACTR2'
20. 'ACTR3'
21. 'ACVR1'
22. 'ACVR1B'
23. 'ACVR1C'
24. 'ACVR2A'
25. 'ADA'
26. 'ADAD1'
27. 'ADAD2'
28. 'ADAM15'
29. 'ADAM18'
30. 'ADAM2'
31. 'ADAM20'
32. 'ADAM21'
33. 'ADAM28'
34. 'ADAM29'
35. 'ADAM32'
36. 'ADAM5'
37. 'ADAMTS1'
38. 'ADAMTS16'
39. 'ADAMTS2'
40. 'ADCY10'
41. 'ADCY3'
42. 'ADCY7'
43. 'ADCYAP1'
44. 'ADCYAP1R1'
45. 'ADGRG1'
46. 'ADGRG2'
47. 'ADIG'
48. 'ADM'
49. 'ADNP'
50. 'ADRA2A'
51. 'ADRA2B'
52. 'AFF4'
53. 'AFP'
54. 'AGFG1'
55. 'AGO2'
56. 'AGO4'
57. 'AGRP'
58. 'AGT'
59. 'AKAP3'
60. 'AKAP4'
61. 'AKR1C3'
62. 'AKT1'
63. 'ALDOA'
64. 'ALKBH5'
65. 'ALOX15B'
66. 'ALPL'
67. 'AMBP'
68. 'AMD1'
69. 'AMH'
70. 'AMHR2'
71. 'ANAPC1'
72. 'ANAPC10'
73. 'ANAPC11'
74. 'ANAPC13'
75. 'ANAPC15'
76. 'ANAPC16'
77. 'ANAPC2'
78. 'ANAPC4'
79. 'ANAPC5'
80. 'ANAPC7'
81. 'ANG'
82. 'ANGPT2'
83. 'ANKLE1'
84. 'ANKRD31'
85. 'ANKRD49'
86. 'ANTXR1'
87. 'ANXA1'
88. 'AP3B1'
89. 'APC2'
90. 'APELA'
91. 'APLF'
92. 'APOB'
93. 'APOL2'
94. 'APOL3'
95. 'APP'
96. 'AQP4'
97. 'AR'
98. 'AREG'
99. 'ARG1'
100. 'ARHGDIB'
101. 'ARID4A'
102. 'ARID4B'
103. 'ARID5B'
104. 'ARMC12'
105. 'ARMC2'
106. 'ARRB2'
107. 'ASB1'
108. 'ASF1B'
109. 'ASH1L'
110. 'ASPM'
111. 'ASZ1'
112. 'ATAT1'
113. 'ATM'
114. 'ATN1'
115. 'ATP1A1'
116. 'ATP1A4'
117. 'ATP2B2'
118. 'ATP2B4'
119. 'ATP7A'
120. 'ATP8B3'
121. 'ATR'
122. 'ATRX'
123. 'AURKA'
124. 'AURKC'
125. 'AVP'
126. 'AVPR1A'
127. 'AXL'
128. 'AZI2'
129. 'AZIN2'
130. 'B4GALNT1'
131. 'B4GALT1'
132. 'BACH1'
133. 'BAD'
134. 'BAG6'
135. 'BAK1'
136. 'BAP1'
137. 'BASP1'
138. 'BAX'
139. 'BBS1'
140. 'BBS2'
141. 'BBS4'
142. 'BCAP31'
143. 'BCL2'
144. 'BCL2L1'
145. 'BCL2L10'
146. 'BCL2L11'
147. 'BCL2L2'
148. 'BCL6'
149. 'BIK'
150. 'BIRC3'
151. 'BMAL1'
152. 'BMP15'
153. 'BMP4'
154. 'BMP5'
155. 'BMP6'
156. 'BMP7'
157. 'BMPR1A'
158. 'BMPR1B'
159. 'BMPR2'
160. 'BNC1'
161. 'BOK'
162. 'BOLL'
163. 'BPY2'
164. 'BRCA2'
165. 'BRD2'
166. 'BRDT'
167. 'BRINP1'
168. 'BRIP1'
169. 'BRME1'
170. 'BSG'
171. 'BTBD18'
172. 'BUB1'
173. 'BUB1B'
174. 'BUB3'
175. 'C11orf80'
176. 'C14orf39'
177. 'C16orf92'
178. 'C1QBP'
179. 'C2CD6'
180. 'C3'
181. 'C3orf62'
182. 'C9orf24'
183. 'CABS1'
184. 'CABYR'
185. 'CACNA1H'
186. 'CACNA1I'
187. 'CAD'
188. 'CADM1'
189. 'CALCA'
190. 'CALR'
191. 'CALR3'
192. 'CAPN2'
193. 'CASP2'
194. 'CASP3'
195. 'CASP8'
196. 'CATSPER1'
197. 'CATSPER2'
198. 'CATSPER3'
199. 'CATSPERB'
200. 'CATSPERD'
201. ⋯
202. 'TBC1D21'
203. 'TBPL1'
204. 'TBX3'
205. 'TCF21'
206. 'TCFL5'
207. 'TCP1'
208. 'TCP11'
209. 'TCTE1'
210. 'TDRD1'
211. 'TDRD10'
212. 'TDRD12'
213. 'TDRD5'
214. 'TDRD6'
215. 'TDRD7'
216. 'TDRD9'
217. 'TDRKH'
218. 'TDRP'
219. 'TEAD3'
220. 'TEAD4'
221. 'TEKT2'
222. 'TEKT3'
223. 'TEP1'
224. 'TERB1'
225. 'TERB2'
226. 'TERF1'
227. 'TESC'
228. 'TESK1'
229. 'TESK2'
230. 'TESMIN'
231. 'TEX101'
232. 'TEX11'
233. 'TEX12'
234. 'TEX14'
235. 'TEX15'
236. 'TEX19'
237. 'TEX43'
238. 'TFAP2C'
239. 'TFPT'
240. 'TGFB2'
241. 'TGFB3'
242. 'TGFBR1'
243. 'TGFBR2'
244. 'TGS1'
245. 'TH'
246. 'THBD'
247. 'THEG'
248. 'THRA'
249. 'THRB'
250. 'TIAL1'
251. 'TIFAB'
252. 'TIMP1'
253. 'TIMP4'
254. 'TIPARP'
255. 'TLE6'
256. 'TLR3'
257. 'TLR9'
258. 'TMED2'
259. 'TMEM119'
260. 'TMEM203'
261. 'TMEM95'
262. 'TMF1'
263. 'TMPRSS12'
264. 'TNC'
265. 'TNFAIP6'
266. 'TNFSF10'
267. 'TNP1'
268. 'TNP2'
269. 'TOB2'
270. 'TOP2A'
271. 'TOP2B'
272. 'TOP3A'
273. 'TP63'
274. 'TPGS1'
275. 'TPPP2'
276. 'TPPP3'
277. 'TRAC'
278. 'TRIM27'
279. 'TRIM28'
280. 'TRIM36'
281. 'TRIP13'
282. 'TRO'
283. 'TRPC3'
284. 'TRPC6'
285. 'TRPC7'
286. 'TRPM2'
287. 'TSGA10'
288. 'TSNAX'
289. 'TSNAXIP1'
290. 'TSPAN8'
291. 'TSPY1'
292. 'TSSK1B'
293. 'TSSK2'
294. 'TSSK3'
295. 'TSSK4'
296. 'TSSK6'
297. 'TTC12'
298. 'TTC21A'
299. 'TTC26'
300. 'TTF1'
301. 'TTK'
302. 'TTLL1'
303. 'TTLL3'
304. 'TTLL9'
305. 'TUBG1'
306. 'TUBG2'
307. 'TUBGCP2'
308. 'TUBGCP3'
309. 'TUBGCP4'
310. 'TUBGCP5'
311. 'TUBGCP6'
312. 'TUT4'
313. 'TUT7'
314. 'TXNDC2'
315. 'TXNDC8'
316. 'TXNRD3'
317. 'TYRO3'
318. 'UBAP2L'
319. 'UBB'
320. 'UBE2B'
321. 'UBE2J1'
322. 'UBE2Q1'
323. 'UBE3A'
324. 'UBR2'
325. 'UBXN8'
326. 'UCHL1'
327. 'UCN'
328. 'UCP2'
329. 'UMODL1'
330. 'UMPS'
331. 'UNC13B'
332. 'UNC5C'
333. 'USP42'
334. 'USP9X'
335. 'USP9Y'
336. 'UTF1'
337. 'UTP14C'
338. 'VDAC2'
339. 'VDR'
340. 'VEGFA'
341. 'VGF'
342. 'VIP'
343. 'VIPAS39'
344. 'VMP1'
345. 'VPS13A'
346. 'VPS13B'
347. 'VPS54'
348. 'WASHC5'
349. 'WBP2NL'
350. 'WDR19'
351. 'WDR33'
352. 'WDR48'
353. 'WDR77'
354. 'WFDC2'
355. 'WIPF3'
356. 'WNT2B'
357. 'WNT3'
358. 'WNT4'
359. 'WNT5A'
360. 'WNT7A'
361. 'WNT9B'
362. 'WT1'
363. 'XDH'
364. 'XKRY'
365. 'XRCC2'
366. 'XRN2'
367. 'YBX2'
368. 'YBX3'
369. 'YTHDC1'
370. 'YTHDC2'
371. 'YTHDF2'
372. 'YTHDF3'
373. 'YY1'
374. 'ZAN'
375. 'ZBTB16'
376. 'ZCWPW1'
377. 'ZFP41'
378. 'ZFP42'
379. 'ZFPM2'
380. 'ZFY'
381. 'ZGLP1'
382. 'ZMIZ1'
383. 'ZMYND15'
384. 'ZNF148'
385. 'ZNF225'
386. 'ZNF296'
387. 'ZNF318'
388. 'ZNF32'
389. 'ZNF35'
390. 'ZNF449'
391. 'ZNF541'
392. 'ZNF628'
393. 'ZNF830'
394. 'ZP1'
395. 'ZP2'
396. 'ZP4'
397. 'ZPBP'
398. 'ZPBP2'
399. 'ZSCAN2'
400. 'ZSCAN21'
401. 'ZW10'

$`GO:0000012`

1. 'APLF'
2. 'APTX'
3. 'ERCC6'
4. 'ERCC8'
5. 'LIG4'
6. 'PARP1'
7. 'SIRT1'
8. 'TDP1'
9. 'TERF2'
10. 'TNP1'
11. 'XNDC1N'
12. 'XRCC1'

$`GO:0000018`

1. 'ACTB'
2. 'ACTL6A'
3. 'ACTR2'
4. 'ALYREF'
5. 'ANKLE1'
6. 'APLF'
7. 'ARID2'
8. 'ATAD5'
9. 'BCL6'
10. 'BLM'
11. 'BRD8'
12. 'CD28'
13. 'CD40'
14. 'CFL1'
15. 'CGAS'
16. 'CHEK1'
17. 'CLC'
18. 'CLCF1'
19. 'DMAP1'
20. 'EAF2'
21. 'EP400'
22. 'EPC1'
23. 'EPC2'
24. 'ERCC2'
25. 'ERCC6'
26. 'EXOSC3'
27. 'EXOSC6'
28. 'FANCB'
29. 'FBH1'
30. 'FIGNL1'
31. 'FOXP3'
32. 'FUS'
33. 'H1-0'
34. 'H1-1'
35. 'H1-10'
36. 'H1-2'
37. 'H1-3'
38. 'H1-4'
39. 'H1-5'
40. 'H1-6'
41. 'H1-8'
42. 'H1-9P'
43. 'HDGFL2'
44. 'HELB'
45. 'HELQ'
46. 'HMCES'
47. 'IL10'
48. 'IL2'
49. 'IL27RA'
50. 'IL4'
51. 'IL7R'
52. 'ING2'
53. 'ING3'
54. 'KAT5'
55. 'KDM1A'
56. 'KHDC3L'
57. 'KLHL15'
58. 'KMT5B'
59. 'KMT5C'
60. 'KPNA1'
61. 'KPNA2'
62. 'MAD2L2'
63. 'MAGEF1'
64. 'MBTD1'
65. 'MEAF6'
66. 'MLH1'
67. 'MMS19'
68. 'MORF4L1'
69. 'MORF4L2'
70. 'MRE11'
71. 'MRGBP'
72. 'MRNIP'
73. 'MSH2'
74. 'MSH3'
75. 'MSH6'
76. 'NDFIP1'
77. 'NSD2'
78. 'OOEP'
79. 'PARP1'
80. 'PARP3'
81. 'PARPBP'
82. 'PAXIP1'
83. 'PIAS4'
84. 'POGZ'
85. 'POLQ'
86. 'PPP4C'
87. 'PPP4R2'
88. 'PRDM7'
89. 'PRDM9'
90. 'PTPRC'
91. 'RAD50'
92. 'RAD51'
93. 'RAD51AP1'
94. 'RADX'
95. 'RBBP8'
96. 'RECQL5'
97. 'RIF1'
98. 'RMI2'
99. 'RPA2'
100. 'RTEL1'
101. 'RUVBL1'
102. 'RUVBL2'
103. 'SETD2'
104. 'SHLD1'
105. 'SHLD2'
106. 'SIRT6'
107. 'SLC15A4'
108. 'SMAP2'
109. 'SMARCAD1'
110. 'SMCHD1'
111. 'SPIDR'
112. 'STAT6'
113. 'SUPT6H'
114. 'TBX21'
115. 'TERF2'
116. 'TERF2IP'
117. 'TEX15'
118. 'TFRC'
119. 'TGFB1'
120. 'THOC1'
121. 'TIMELESS'
122. 'TNFSF13'
123. 'TNFSF4'
124. 'TP53BP1'
125. 'TRRAP'
126. 'UBE2B'
127. 'UBQLN4'
128. 'USP51'
129. 'VPS72'
130. 'WAS'
131. 'WDR48'
132. 'WRAP53'
133. 'YEATS4'
134. 'ZCWPW1'
135. 'ZNF365'
136. 'ZRANB3'
137. 'ZSCAN4'

$`GO:0000022`

1. 'AURKB'
2. 'AURKC'
3. 'BIRC5'
4. 'CDCA8'
5. 'INCENP'
6. 'KIF23'
7. 'KIF4A'
8. 'MAP10'
9. 'NUMA1'
10. 'PRC1'
11. 'RACGAP1'

Now, we already had GO terms with genes. However, we still do not know the meaning of GO terms related to biological process. We can use GO.db database to get a set of annotation maps describing the entire Gene Ontology assembled using data from GO. We can use the following code to install the GO.db R package.

suppressMessages({if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
  suppressWarnings(BiocManager::install("GO.db", update = F))
})
library(GO.db)

Then, we can use the following command to obtain the GO terms description.

# Getting the name of each GO term
terms <- names(allGO)
# Getting the description of each GO term
descriptions <-lapply(Term(terms), `[[`, 1)

In order to perform enrichment analysis in later submodules, we need to save the GO terms and genesets to the standard output. One commonly used format is Gene Matrix Transposed file format (*.gmt). The GMT file format is a tab delimited file format that describes gene sets. In the GMT format, each row represents a gene set; in the GMX format, each column represents a gene set. Here, we can save GO terms and genesets to the *.gmt using the following function:

# A function to save the GO terms with geneset to the local repository
writeGMT <- function(genesets, descriptions, outfile) {

  if (file.exists(outfile)) {
    file.remove(outfile)
  }
  for (gs in names(genesets)) {
    write(c(gs, gsub("\t", " ", descriptions[[gs]]), genesets[[gs]]), file=outfile, sep="\t", append=TRUE, ncolumns=length(genesets[[gs]]) + 2)
  }
}
outfile <- "./data/GO_terms.gmt"
writeGMT(allGO, descriptions, outfile)

Tip: Saving data to the Google Cloud Bucket

gsutil cp ./data/GO_terms.gmt gs://cpa-output

Kyoto Encyclopedia of Genes and Genomes (KEGG)

Overview

KEGG is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances. KEGG is utilized for bioinformatics research and education, including data analysis in genomics, metagenomics, metabolomics and other omics studies, modeling and simulation in systems biology, and translational research in drug development. The KEGG database project was initiated in 1995 by Minoru Kanehisa, professor at the Institute for Chemical Research, Kyoto University, under the then ongoing Japanese Human Genome Program. Foreseeing the need for a computerized resource that can be used for biological interpretation of genome sequence data, he started developing the KEGG PATHWAY database. It is a collection of manually drawn KEGG pathway maps representing experimental knowledge on metabolism and various other functions of the cell and the organism. Each pathway map contains a network of molecular interactions and reactions and is designed to link genes in the genome to gene products (mostly proteins) in the pathway. This has enabled the analysis called KEGG pathway mapping, whereby the gene content in the genome is compared with the KEGG PATHWAY database to examine which pathways and associated functions are likely to be encoded in the genome. KEGG is a “computer representation” of the biological system. It integrates building blocks and wiring diagrams of the system—more specifically, genetic building blocks of genes and proteins, chemical building blocks of small molecules and reactions, and wiring diagrams of molecular interaction and reaction networks. The illustrative structure of KEGG is presented as figure below.

Retrieving pathways from KEGG databases

In this section, we will retrieve pathways and related genesets from the KEGG database using R command line. Here we will use KEGGREST R package that provides a client interface to the KEGG REST server. KEGGREST can be installed from the Bioconductor using following command.

suppressMessages({if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
  suppressWarnings(BiocManager::install("KEGGREST", update = F))
})
suppressPackageStartupMessages({
  library(KEGGREST)
})

KEGG contains a number of databases. To get an idea of what is available, run listDatabases():

KEGGREST::listDatabases()

1. 'pathway'
2. 'brite'
3. 'module'
4. 'ko'
5. 'genome'
6. 'vg'
7. 'ag'
8. 'compound'
9. 'glycan'
10. 'reaction'
11. 'rclass'
12. 'enzyme'
13. 'disease'
14. 'drug'
15. 'dgroup'
16. 'environ'
17. 'genes'
18. 'ligand'
19. 'kegg'

We can use these databases in further queries. Note that in many cases you can also use a three-letter KEGG organism code or a “T number” (genome identifier) in the same place you would use one of these database names.

We can obtain the list of organisms available in KEGG with the keggList() function:

organism <- keggList("organism")

print(paste0("KEGG supports ",dim(organism)[1]," organisms"))

[1] "KEGG supports 8889 organisms"

To view the supported organisms we can use the following command:

# View several supported organism
head(organism)

A matrix: 6 × 4 of type chr

T.number	organism	species	phylogeny
T01001	hsa	Homo sapiens (human)	Eukaryotes;Animals;Vertebrates;Mammals
T01005	ptr	Pan troglodytes (chimpanzee)	Eukaryotes;Animals;Vertebrates;Mammals
T02283	pps	Pan paniscus (bonobo)	Eukaryotes;Animals;Vertebrates;Mammals
T02442	ggo	Gorilla gorilla gorilla (western lowland gorilla)	Eukaryotes;Animals;Vertebrates;Mammals
T01416	pon	Pongo abelii (Sumatran orangutan)	Eukaryotes;Animals;Vertebrates;Mammals
T03265	nle	Nomascus leucogenys (northern white-cheeked gibbon)	Eukaryotes;Animals;Vertebrates;Mammals

In submodule 02, we performed DE analysis on a human dataset. Therefore, we need to download pathways for humans. The abbreviation of human pathway in KEGG is hsa and we can use keggList function to get the pathway list.

# Obtain the pathways belong to human
pathways.list <- keggList("pathway", "hsa")

The pathway list contains pathway description and pathway code in a single line of text. To see the first five pathways, we can use the following command:

#list the specific pathways to view
pathway_ids <- c("hsa00010", "hsa00020", "hsa00030", "hsa00040", "hsa00051")
# View the first five pathways
pathways.list[pathway_ids]

hsa00010

'Glycolysis / Gluconeogenesis - Homo sapiens (human)'

hsa00020

'Citrate cycle (TCA cycle) - Homo sapiens (human)'

hsa00030

'Pentose phosphate pathway - Homo sapiens (human)'

hsa00040

'Pentose and glucuronate interconversions - Homo sapiens (human)'

hsa00051

'Fructose and mannose metabolism - Homo sapiens (human)'

We can see that, in each line, the text in the quotation mark contains pathway information while the later part contains pathway code leading by a prefix path:. To get pathway codes from the pathway list, we can use the following command:

# Retrieve all the pathway IDs belong to human
pathway.codes <- sub("path:", "", names(pathways.list))
pathway.codes

1. 'hsa01100'
2. 'hsa01200'
3. 'hsa01210'
4. 'hsa01212'
5. 'hsa01230'
6. 'hsa01232'
7. 'hsa01250'
8. 'hsa01240'
9. 'hsa00010'
10. 'hsa00020'
11. 'hsa00030'
12. 'hsa00040'
13. 'hsa00051'
14. 'hsa00052'
15. 'hsa00053'
16. 'hsa00500'
17. 'hsa00520'
18. 'hsa00620'
19. 'hsa00630'
20. 'hsa00640'
21. 'hsa00650'
22. 'hsa00562'
23. 'hsa00190'
24. 'hsa00910'
25. 'hsa00920'
26. 'hsa00061'
27. 'hsa00062'
28. 'hsa00071'
29. 'hsa00100'
30. 'hsa00120'
31. 'hsa00140'
32. 'hsa00561'
33. 'hsa00564'
34. 'hsa00565'
35. 'hsa00600'
36. 'hsa00590'
37. 'hsa00591'
38. 'hsa00592'
39. 'hsa01040'
40. 'hsa00230'
41. 'hsa00240'
42. 'hsa00250'
43. 'hsa00260'
44. 'hsa00270'
45. 'hsa00280'
46. 'hsa00290'
47. 'hsa00310'
48. 'hsa00220'
49. 'hsa00330'
50. 'hsa00340'
51. 'hsa00350'
52. 'hsa00360'
53. 'hsa00380'
54. 'hsa00400'
55. 'hsa00410'
56. 'hsa00430'
57. 'hsa00440'
58. 'hsa00450'
59. 'hsa00470'
60. 'hsa00480'
61. 'hsa00510'
62. 'hsa00513'
63. 'hsa00512'
64. 'hsa00515'
65. 'hsa00514'
66. 'hsa00532'
67. 'hsa00534'
68. 'hsa00533'
69. 'hsa00531'
70. 'hsa00563'
71. 'hsa00601'
72. 'hsa00603'
73. 'hsa00604'
74. 'hsa00511'
75. 'hsa00730'
76. 'hsa00740'
77. 'hsa00750'
78. 'hsa00760'
79. 'hsa00770'
80. 'hsa00780'
81. 'hsa00785'
82. 'hsa00790'
83. 'hsa00670'
84. 'hsa00830'
85. 'hsa00860'
86. 'hsa00130'
87. 'hsa00900'
88. 'hsa00232'
89. 'hsa00524'
90. 'hsa00980'
91. 'hsa00982'
92. 'hsa00983'
93. 'hsa03020'
94. 'hsa03022'
95. 'hsa03040'
96. 'hsa03010'
97. 'hsa00970'
98. 'hsa03013'
99. 'hsa03015'
100. 'hsa03008'
101. 'hsa03060'
102. 'hsa04141'
103. 'hsa04130'
104. 'hsa04120'
105. 'hsa04122'
106. 'hsa03050'
107. 'hsa03018'
108. 'hsa03030'
109. 'hsa03410'
110. 'hsa03420'
111. 'hsa03430'
112. 'hsa03440'
113. 'hsa03450'
114. 'hsa03460'
115. 'hsa03250'
116. 'hsa03260'
117. 'hsa03264'
118. 'hsa03265'
119. 'hsa03266'
120. 'hsa03267'
121. 'hsa02010'
122. 'hsa04010'
123. 'hsa04012'
124. 'hsa04014'
125. 'hsa04015'
126. 'hsa04310'
127. 'hsa04330'
128. 'hsa04340'
129. 'hsa04350'
130. 'hsa04390'
131. 'hsa04392'
132. 'hsa04370'
133. 'hsa04371'
134. 'hsa04630'
135. 'hsa04064'
136. 'hsa04668'
137. 'hsa04066'
138. 'hsa04068'
139. 'hsa04020'
140. 'hsa04070'
141. 'hsa04072'
142. 'hsa04071'
143. 'hsa04024'
144. 'hsa04022'
145. 'hsa04151'
146. 'hsa04152'
147. 'hsa04150'
148. 'hsa04080'
149. 'hsa04060'
150. 'hsa04061'
151. 'hsa04512'
152. 'hsa04514'
153. 'hsa04144'
154. 'hsa04145'
155. 'hsa04142'
156. 'hsa04146'
157. 'hsa04140'
158. 'hsa04136'
159. 'hsa04137'
160. 'hsa04110'
161. 'hsa04114'
162. 'hsa04210'
163. 'hsa04215'
164. 'hsa04216'
165. 'hsa04217'
166. 'hsa04115'
167. 'hsa04218'
168. 'hsa04510'
169. 'hsa04520'
170. 'hsa04530'
171. 'hsa04540'
172. 'hsa04550'
173. 'hsa04814'
174. 'hsa04810'
175. 'hsa04640'
176. 'hsa04610'
177. 'hsa04611'
178. 'hsa04613'
179. 'hsa04620'
180. 'hsa04621'
181. 'hsa04622'
182. 'hsa04623'
183. 'hsa04625'
184. 'hsa04650'
185. 'hsa04612'
186. 'hsa04660'
187. 'hsa04658'
188. 'hsa04659'
189. 'hsa04657'
190. 'hsa04662'
191. 'hsa04664'
192. 'hsa04666'
193. 'hsa04670'
194. 'hsa04672'
195. 'hsa04062'
196. 'hsa04911'
197. 'hsa04910'
198. 'hsa04922'
199. 'hsa04923'
200. 'hsa04920'
201. 'hsa03320'
202. 'hsa04929'
203. 'hsa04912'
204. 'hsa04913'
205. 'hsa04915'
206. 'hsa04914'
207. 'hsa04917'
208. 'hsa04921'
209. 'hsa04926'
210. 'hsa04935'
211. 'hsa04918'
212. 'hsa04919'
213. 'hsa04928'
214. 'hsa04916'
215. 'hsa04924'
216. 'hsa04614'
217. 'hsa04925'
218. 'hsa04927'
219. 'hsa04260'
220. 'hsa04261'
221. 'hsa04270'
222. 'hsa04970'
223. 'hsa04971'
224. 'hsa04972'
225. 'hsa04976'
226. 'hsa04973'
227. 'hsa04974'
228. 'hsa04975'
229. 'hsa04979'
230. 'hsa04977'
231. 'hsa04978'
232. 'hsa04962'
233. 'hsa04960'
234. 'hsa04961'
235. 'hsa04964'
236. 'hsa04966'
237. 'hsa04724'
238. 'hsa04727'
239. 'hsa04725'
240. 'hsa04728'
241. 'hsa04726'
242. 'hsa04720'
243. 'hsa04730'
244. 'hsa04723'
245. 'hsa04721'
246. 'hsa04722'
247. 'hsa04744'
248. 'hsa04740'
249. 'hsa04742'
250. 'hsa04750'
251. 'hsa04360'
252. 'hsa04380'
253. 'hsa04211'
254. 'hsa04213'
255. 'hsa04710'
256. 'hsa04713'
257. 'hsa04714'
258. 'hsa05200'
259. 'hsa05202'
260. 'hsa05206'
261. 'hsa05205'
262. 'hsa05204'
263. 'hsa05207'
264. 'hsa05208'
265. 'hsa05203'
266. 'hsa05230'
267. 'hsa05231'
268. 'hsa05235'
269. 'hsa05210'
270. 'hsa05212'
271. 'hsa05225'
272. 'hsa05226'
273. 'hsa05214'
274. 'hsa05216'
275. 'hsa05221'
276. 'hsa05220'
277. 'hsa05217'
278. 'hsa05218'
279. 'hsa05211'
280. 'hsa05219'
281. 'hsa05215'
282. 'hsa05213'
283. 'hsa05224'
284. 'hsa05222'
285. 'hsa05223'
286. 'hsa05166'
287. 'hsa05170'
288. 'hsa05161'
289. 'hsa05160'
290. 'hsa05171'
291. 'hsa05164'
292. 'hsa05162'
293. 'hsa05168'
294. 'hsa05163'
295. 'hsa05167'
296. 'hsa05169'
297. 'hsa05165'
298. 'hsa05110'
299. 'hsa05120'
300. 'hsa05130'
301. 'hsa05132'
302. 'hsa05131'
303. 'hsa05135'
304. 'hsa05133'
305. 'hsa05134'
306. 'hsa05150'
307. 'hsa05152'
308. 'hsa05100'
309. 'hsa05146'
310. 'hsa05144'
311. 'hsa05145'
312. 'hsa05140'
313. 'hsa05142'
314. 'hsa05143'
315. 'hsa05310'
316. 'hsa05322'
317. 'hsa05323'
318. 'hsa05320'
319. 'hsa05321'
320. 'hsa05330'
321. 'hsa05332'
322. 'hsa05340'
323. 'hsa05010'
324. 'hsa05012'
325. 'hsa05014'
326. 'hsa05016'
327. 'hsa05017'
328. 'hsa05020'
329. 'hsa05022'
330. 'hsa05030'
331. 'hsa05031'
332. 'hsa05032'
333. 'hsa05033'
334. 'hsa05034'
335. 'hsa05417'
336. 'hsa05418'
337. 'hsa05410'
338. 'hsa05412'
339. 'hsa05414'
340. 'hsa05415'
341. 'hsa05416'
342. 'hsa04930'
343. 'hsa04940'
344. 'hsa04950'
345. 'hsa04936'
346. 'hsa04932'
347. 'hsa04931'
348. 'hsa04933'
349. 'hsa04934'
350. 'hsa01521'
351. 'hsa01524'
352. 'hsa01523'
353. 'hsa01522'

We can use the following command to check how many pathways are available for human

print(paste0("Number of available pathways for human are: ", length(pathway.codes)))

[1] "Number of available pathways for human are: 353"

The following code will help to get list of genes and pathway’s description for all pathways available in human.

# Function to get all the gene names for each pathway
genes.by.pathway <- sapply(pathway.codes,
                           function(pwid){
                             pw <- keggGet(pwid)
                             if (is.null(pw[[1]]$GENE)) return(NA)
                             pw2 <- pw[[1]]$GENE[c(FALSE,TRUE)]
                             pw2 <- unlist(lapply(strsplit(pw2, split = ";", fixed = T), function(x)x[1]))
                             return(pw2)
                           }
)
# Function to get description for each pathway
description.by.pathway <- sapply(pathway.codes,
                                 function(pwid){
                                   pw <- keggGet(pwid)
                                   if (is.null(pw[[1]]$NAME)) return(NA)
                                   pw2 <- pw[[1]]$NAME
                                   return(pw2)
                                 }
)
# Convert the pathway description to a list
description.by.pathway <- as.list(description.by.pathway)

We can view the first five pathways with their genesets using the following command

# View the five pathways with the genesets
genes.by.pathway <- genes.by.pathway[!is.na(genes.by.pathway)]
genes.by.pathway[1:5]

$hsa00010

1. 'HK3'
2. 'HK1'
3. 'HK2'
4. 'HKDC1'
5. 'GCK'
6. 'GPI'
7. 'PFKM'
8. 'PFKP'
9. 'PFKL'
10. 'FBP1'
11. 'FBP2'
12. 'ALDOC'
13. 'ALDOA'
14. 'ALDOB'
15. 'TPI1'
16. 'GAPDH'
17. 'GAPDHS'
18. 'PGK2'
19. 'PGK1'
20. 'PGAM1'
21. 'PGAM2'
22. 'PGAM4'
23. 'ENO3'
24. 'ENO2'
25. 'ENO1'
26. 'ENO4'
27. 'PKM'
28. 'PKLR'
29. 'PDHA2'
30. 'PDHA1'
31. 'PDHB'
32. 'DLAT'
33. 'DLD'
34. 'LDHAL6A'
35. 'LDHAL6B'
36. 'LDHA'
37. 'LDHB'
38. 'LDHC'
39. 'ADH1A'
40. 'ADH1B'
41. 'ADH1C'
42. 'ADH7'
43. 'ADH4'
44. 'ADH5'
45. 'ADH6'
46. 'AKR1A1'
47. 'ALDH2'
48. 'ALDH3A2'
49. 'ALDH1B1'
50. 'ALDH7A1'
51. 'ALDH9A1'
52. 'ALDH3B1'
53. 'ALDH3B2'
54. 'ALDH3A1'
55. 'ACSS1'
56. 'ACSS2'
57. 'GALM'
58. 'PGM1'
59. 'PGM2'
60. 'G6PC1'
61. 'G6PC2'
62. 'G6PC3'
63. 'ADPGK'
64. 'BPGM'
65. 'MINPP1'
66. 'PCK1'
67. 'PCK2'

$hsa00020

1. 'CS'
2. 'ACLY'
3. 'ACO2'
4. 'ACO1'
5. 'IDH1'
6. 'IDH2'
7. 'IDH3B'
8. 'IDH3G'
9. 'IDH3A'
10. 'OGDHL'
11. 'OGDH'
12. 'DLST'
13. 'DLD'
14. 'SUCLG1'
15. 'SUCLG2'
16. 'SUCLA2'
17. 'SDHA'
18. 'SDHB'
19. 'SDHC'
20. 'SDHD'
21. 'FH'
22. 'MDH1'
23. 'MDH2'
24. 'PC'
25. 'PCK1'
26. 'PCK2'
27. 'PDHA2'
28. 'PDHA1'
29. 'PDHB'
30. 'DLAT'

$hsa00030

1. 'GPI'
2. 'G6PD'
3. 'PGLS'
4. 'H6PD'
5. 'PGD'
6. 'RPE'
7. 'RPEL1'
8. 'TKT'
9. 'TKTL2'
10. 'TKTL1'
11. 'TALDO1'
12. 'RPIA'
13. 'DERA'
14. 'RBKS'
15. 'PGM1'
16. 'PGM2'
17. 'PRPS1L1'
18. 'PRPS2'
19. 'PRPS1'
20. 'RGN'
21. 'IDNK'
22. 'GLYCTK'
23. 'ALDOC'
24. 'ALDOA'
25. 'ALDOB'
26. 'FBP1'
27. 'FBP2'
28. 'PFKM'
29. 'PFKP'
30. 'PFKL'

$hsa00040

1. 'GUSB'
2. 'KL'
3. 'UGT2A1'
4. 'UGT2A3'
5. 'UGT2B17'
6. 'UGT2B11'
7. 'UGT2B28'
8. 'UGT1A6'
9. 'UGT1A4'
10. 'UGT1A1'
11. 'UGT1A3'
12. 'UGT2B10'
13. 'UGT1A9'
14. 'UGT2B7'
15. 'UGT1A10'
16. 'UGT1A8'
17. 'UGT1A5'
18. 'UGT2B15'
19. 'UGT1A7'
20. 'UGT2B4'
21. 'UGT2A2'
22. 'UGDH'
23. 'UGP2'
24. 'AKR1A1'
25. 'CRYL1'
26. 'RPE'
27. 'RPEL1'
28. 'XYLB'
29. 'AKR1B1'
30. 'AKR1B10'
31. 'DCXR'
32. 'SORD'
33. 'DHDH'
34. 'FGGY'
35. 'CRPPA'

$hsa00051

1. 'MPI'
2. 'PMM2'
3. 'PMM1'
4. 'GMPPB'
5. 'GMPPA'
6. 'GMDS'
7. 'GFUS'
8. 'FPGT'
9. 'FCSK'
10. 'ENOSF1'
11. 'HK3'
12. 'HK1'
13. 'HK2'
14. 'HKDC1'
15. 'PFKM'
16. 'PFKP'
17. 'PFKL'
18. 'FBP1'
19. 'FBP2'
20. 'PFKFB1'
21. 'PFKFB2'
22. 'PFKFB3'
23. 'PFKFB4'
24. 'TIGAR'
25. 'KHK'
26. 'SORD'
27. 'AKR1B1'
28. 'AKR1B10'
29. 'ALDOC'
30. 'ALDOA'
31. 'ALDOB'
32. 'TPI1'
33. 'TKFC'

Use the following command to see the description of the first five pathways

# View the description of the first five pathways
description.by.pathway <- description.by.pathway[!is.na(description.by.pathway)]
description.by.pathway[names(genes.by.pathway[1:5])]

$hsa00010

'Glycolysis / Gluconeogenesis - Homo sapiens (human)'

$hsa00020

'Citrate cycle (TCA cycle) - Homo sapiens (human)'

$hsa00030

'Pentose phosphate pathway - Homo sapiens (human)'

$hsa00040

'Pentose and glucuronate interconversions - Homo sapiens (human)'

$hsa00051

'Fructose and mannose metabolism - Homo sapiens (human)'

Then we can save the output to *.gmt file using the following commands

# Saving the pathway information to the local repository
outfile <- "./data/KEGG_pathways.gmt"
writeGMT(genes.by.pathway, description.by.pathway, outfile)

gsutil cp ./data/KEGG_pathways.gmt gs://cpa-output

In the next submodule, we will do Pathway Analysis.