GOpro: Determine groups of genes and find their most characteristic GO term
30 October 2018
Contents
1 Installation
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("GOpro", dependencies = TRUE)2 Loading
After the package is installed, it can be loaded into R workspace typing
library(GOpro)3 Overview
This document presents an overview of the GOpro package. This package is for determining groups of genes and finding characteristic functions for these groups. It allows for interpreting groups of genes by their most characteristic biological function. It provides one function findGO which is based on the set of methods. One of these methods allows for determining significantly different genes between at least two distinct groups (i.e. patients with different medical condition) - the ANOVA test with correction for multiple testing. It also provides two methods for grouping genes. One of them is so-called all pairwise comparisons utilizing Tukey’s method. By this method profiles of genes are determined, i.e. in terms of the gene expression genes are grouped according to the differences in the expressions between given cohorts. Another method of grouping is hierarchical clustering. This package provides a method for finding the most characteristic gene ontology terms for anteriorly obtained groups using the one-sided Fisher’s test for overrepresentation of the gene ontology term. If genes were grouped by the hierarchical clustering, then the most characteristic function is found for all possible groups (for each node in the dendrogram).
4 Details
Genes must be named with the gene aliases and they must be arranged in the same order for each cohort.
4.1 Determining significantly different genes based on their expressions
Genes which are statistically differently expressed are selected for the
further analysis by ANOVA test.
The topAOV parameter denotes the maximum number of significantly
different genes
to be selected. The significance level of ANOVA test is specified
by the sig.levelAOV parameter.
This threshold is used as the significance level in the BH correction
for multiple testing.
In the case of equal p-values of the test (below the given threshold),
all genes for which the p-value of the test is the same as for the gene
numbered with the topAOV value are included in the result.
4.2 Grouping genes based on their similarity
There are two methods provided for grouping genes. They are specified by the grouped parameter. The first one using Tukey’s test is called when grouped equals ‘tukey’ and the second one can be called by using the ‘clustering’ value.
4.2.1 All pairwise comparisons by Tukey’s test
The Tukey’s test is applied to group genes based on their profiles. The sig.levelTUK parameter denotes the significance level of Tukey’s test. For each gene two-sided Tukey’s test is conducted among cohorts. The mean expressions in the cohorts are arranged in ascending order and the result of the test is adapted. All genes with the same order of means and the same result of the test are grouped together. I.e. notation colon=bladder<leukemia denotes that the mean expression calculated for a particular gene in the colon cancer cohort is statistically the same as for the bladder cancer cohort. Both mean values determined for aforementioned cohorts are statistically lower than the mean expression measured for the leukemia cohort.
4.2.2 Hierarchical clustering
Hierarchical clustering is utilized for grouping genes based on dissimilarity measures. In this case all clusters are subjected to the further analysis. The clust.metric parameter is a method to calculate a distance measure between genes, the clust.method is the agglomeration method used to cluster genes, and the dist.matrix is an optional parameter for distance matrix if available clust.metric methods are not sufficient for the user.
4.3 Finding characteristic gene ontology terms
For each specified group the org.Hs.eg.db is searched for all relevant GO terms. The number of counts of each GO term is calculated for each group. Then the Fisher’s test is applied in order to find the most characteristic GO terms for each group of genes. The minGO and maxGO parameters denote the range of counts of genes annotated with each GO term. All GO terms with counts above or below this range are omitted in the analysis. It enables for the exclusion of very rare or very common gene ontology terms. Gene ontology domains to be searched for GO terms can be specified by the onto parameter. Possible domains are: ‘MF’ (molecular function), ‘BP’ (biological process), and ‘CC’ (cellular component). The sig.levelGO parameter specifies the significance level of the Fisher’s test (correction for multiple testing is included).
5 Data
Data used in this example comes from The Cancer Genome Atlas. They were downloaded via RTCGA.PANCAN12 package. The data represents expressions of 300 genes randomly chosen from 16115 genes determined for each patient (sample). Three cohorts are included: acute myleoid leukemia, colon cancer, and bladder cancer. The data is stored in this GOpro package as a MultiAssayExperiment object.
An example of the data structure:
exrtcga## A MultiAssayExperiment object of 3 listed
## experiments with user-defined names and respective classes.
## Containing an ExperimentList class object of length 3:
## [1] leukemia: matrix with 300 rows and 173 columns
## [2] colon: matrix with 300 rows and 190 columns
## [3] bladder: matrix with 300 rows and 122 columns
## Features:
## experiments() - obtain the ExperimentList instance
## colData() - the primary/phenotype DataFrame
## sampleMap() - the sample availability DataFrame
## `$`, `[`, `[[` - extract colData columns, subset, or experiment
## *Format() - convert into a long or wide DataFrame
## assays() - convert ExperimentList to a SimpleList of matrices6 Example
To run the analysis with default parameters on the exrtcga object call:
findGO(exrtcga)## DataFrame with 7 rows and 4 columns
## profile GOs
## <list> <list>
## bladder=colon<leukemia bladder=colon<leukemia GO:0044877
## colon<bladder<leukemia colon<bladder<leukemia c("GO:0000981", "GO:0043565")
## colon=bladder<leukemia colon=bladder<leukemia
## leukemia<bladder<colon leukemia<bladder<colon
## leukemia<bladder=colon leukemia<bladder=colon c("GO:0005783", "GO:0042802")
## leukemia<colon<bladder leukemia<colon<bladder GO:0046982
## leukemia<colon=bladder leukemia<colon=bladder
## p.values
## <list>
## bladder=colon<leukemia c(`GO:0044877` = 0.011)
## colon<bladder<leukemia c(`GO:0000981` = 0, `GO:0043565` = 0.042)
## colon=bladder<leukemia c(NA)
## leukemia<bladder<colon c(NA)
## leukemia<bladder=colon c(`GO:0005783` = 0.023, `GO:0042802` = 0.023)
## leukemia<colon<bladder c(`GO:0046982` = 0.047)
## leukemia<colon=bladder c(NA)
## GENES
## <list>
## bladder=colon<leukemia FAM46A WASF2
## colon<bladder<leukemia ARID2 BPTF CCDC88A CEP110 CLCN6 ELF2 FOSB KIAA0226 LUC7L3 RNF24 UBR2 USP15 ZCCHC11 ZNF592
## colon=bladder<leukemia DENND4B RCOR3
## leukemia<bladder<colon AK1 BRI3 CCDC51 CD24 COX8A DCI FDPS FKBP9 HSD3B7 MRPL54 SDHAF1 SMAGP STAU1 SYAP1 TBRG4 TPRN WDR34 ZNF511
## leukemia<bladder=colon NSDHL PFKM
## leukemia<colon<bladder BOLA1 EDF1 NDUFA1 NEDD8 PAF1 PIPSL PLXNA1 PPAP2A SAE1 SERTAD1 SLC39A1
## leukemia<colon=bladder NA NA NAThe results of the analysis can be presented in a more descriptive way or in a concise one. To get more descriptive results use extend=TRUE option. Additionally, the TERM, DEFINITION, and ONTOLOGY for each ontology term are returned.
findGO(exrtcga, extend = TRUE)## DataFrame with 1 row and 5 columns
## GROUP GOID
## <IntegerList> <CharacterList>
## 1 1,2,2,... GO:0044877,GO:0000981,GO:0043565,...
## TERM
## <CharacterList>
## 1 protein-containing complex binding,DNA-binding transcription factor activity, RNA polymerase II-specific,sequence-specific DNA binding,...
## DEFINITION
## <CharacterList>
## 1 Interacting selectively and non-covalently with a macromolecular complex.,A protein or a member of a complex that interacts selectively and non-covalently with a specific DNA sequence (sometimes referred to as a motif) within the regulatory region of a RNA polymerase II-transcribed gene to modulate transcription. Regulatory regions include promoters (proximal and distal) and enhancers. Genes are transcriptional units.,Interacting selectively and non-covalently with DNA of a specific nucleotide composition, e.g. GC-rich DNA binding, or with a specific sequence motif or type of DNA e.g. promotor binding or rDNA binding.,...
## ONTOLOGY
## <CharacterList>
## 1 MF,MF,MF,...In order to find top 2 GO terms for genes grouped by hierarchical clustering run the following call. The result of clustering is presented on the plot.
findGO(exrtcga, topGO = 2, grouped = 'clustering')## DataFrame with 99 rows and 4 columns
## profile GOs
## <list> <list>
## G1 G1
## G2 G2
## G3 G3
## G4 G4 NA
## G5 G5
## ... ... ...
## G95 G95 c("GO:0000981", "GO:0006357")
## G96 G96
## G97 G97 GO:0003872
## G98 G98 GO:0000981
## G99 G99 GO:0003872
## p.values
## <list>
## G1 c(NA)
## G2 c(NA)
## G3 c(NA)
## G4 c(`NA` = NA)
## G5 c(NA)
## ... ...
## G95 c(`GO:0000981` = 0.011, `GO:0006357` = 0.038)
## G96 c(NA)
## G97 c(`GO:0003872` = 0.001)
## G98 c(`GO:0000981` = 0)
## G99 c(`GO:0003872` = 0.011)
## GENES
## <list>
## G1 SMAGP
## G2 NA NA NA
## G3 HSD3B7
## G4 PLXNA1
## G5 AK1
## ... ...
## G95 CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24
## G96 CD24 FKBP9
## G97 SMAGP PPAP2A HSD3B7 PLXNA1 AK1 PFKM TPRN SAE1 STAU1 SYAP1 NDUFA1 CCDC51 FDPS NSDHL MRPL54 SDHAF1 TBRG4 ZNF511 PAF1 SLC39A1 EDF1 NEDD8 BOLA1 PIPSL SERTAD1 TAX1BP3 BRI3 WDR34 COX8A DCI CD24 FKBP9
## G98 FOSB CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24
## G99 SMAGP PPAP2A HSD3B7 PLXNA1 AK1 PFKM TPRN SAE1 STAU1 SYAP1 NDUFA1 CCDC51 FDPS NSDHL MRPL54 SDHAF1 TBRG4 ZNF511 PAF1 SLC39A1 EDF1 NEDD8 BOLA1 PIPSL SERTAD1 TAX1BP3 BRI3 WDR34 COX8A DCI CD24 FKBP9 FOSB CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24
The plot can be also enriched with information about the most frequent
ontology domain for each node on the dendrogram.
findGO(exrtcga, topGO = 2, grouped = 'clustering', over.rep = TRUE)## DataFrame with 99 rows and 4 columns
## profile GOs
## <list> <list>
## G1 G1
## G2 G2
## G3 G3
## G4 G4 NA
## G5 G5
## ... ... ...
## G95 G95 c("GO:0000981", "GO:0006357")
## G96 G96
## G97 G97 GO:0003872
## G98 G98 GO:0000981
## G99 G99 GO:0003872
## p.values
## <list>
## G1 c(NA)
## G2 c(NA)
## G3 c(NA)
## G4 c(`NA` = NA)
## G5 c(NA)
## ... ...
## G95 c(`GO:0000981` = 0.011, `GO:0006357` = 0.038)
## G96 c(NA)
## G97 c(`GO:0003872` = 0.001)
## G98 c(`GO:0000981` = 0)
## G99 c(`GO:0003872` = 0.011)
## GENES
## <list>
## G1 SMAGP
## G2 NA NA NA
## G3 HSD3B7
## G4 PLXNA1
## G5 AK1
## ... ...
## G95 CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24
## G96 CD24 FKBP9
## G97 SMAGP PPAP2A HSD3B7 PLXNA1 AK1 PFKM TPRN SAE1 STAU1 SYAP1 NDUFA1 CCDC51 FDPS NSDHL MRPL54 SDHAF1 TBRG4 ZNF511 PAF1 SLC39A1 EDF1 NEDD8 BOLA1 PIPSL SERTAD1 TAX1BP3 BRI3 WDR34 COX8A DCI CD24 FKBP9
## G98 FOSB CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24
## G99 SMAGP PPAP2A HSD3B7 PLXNA1 AK1 PFKM TPRN SAE1 STAU1 SYAP1 NDUFA1 CCDC51 FDPS NSDHL MRPL54 SDHAF1 TBRG4 ZNF511 PAF1 SLC39A1 EDF1 NEDD8 BOLA1 PIPSL SERTAD1 TAX1BP3 BRI3 WDR34 COX8A DCI CD24 FKBP9 FOSB CCDC88A KIAA0226 LUC7L3 USP15 DENND4B ZCCHC11 ARID2 BPTF UBR2 ELF2 RCOR3 CLCN6 WASF2 ZNF592 FAM46A CEP110 RNF24