3.1 Gene Lists

The key input data that goSTAG requires is a set of gene lists. These gene lists can be derived from any type of gene-based analysis, such as differential expression analysis or ChIP-seq peak targets. They can be generated within R using another package, or using an external tool and then loaded into R.

Once you have generated your gene lists and loaded them into R, you must convert them into the format that goSTAG requires, which is a list of vectors. Each element of the list corresponds to a gene list, and the vectors contain the genes in the gene list.

We provide an example set of gene lists in the goSTAG package.

data( goSTAG_example_gene_lists )
head( lapply( goSTAG_example_gene_lists, head ) )

## $Day1
## [1] "CNN2"    "GHDC"    "GALNT10" "LIPT1"   "TPP1"    "TBPL2"  
## 
## $Day2
## [1] "GPR37L1"  "TNFSF13B" "OLAH"     "WWC1"     "NYNRIN"   "MRGPRX4" 
## 
## $Day3
## [1] "MAD2L1BP"      "MRPL47"        "SLIRP"         "MEF2BNB-MEF2B"
## [5] "TMEM245"       "UBA6"

For convienience, we have provided a function for loading external gene lists.

The first option is to put your gene lists into a single GMT file. This GMT file is a tab-delimited file containing one gene list per line. The first entry in the line is the name of the gene list, the second entry in the line is a description of the gene list (which goSTAG ignores), and the rest of the entries are the genes. For more information on the GMT format, see here.

Once you have created your GMT file, you can load it using the loadGeneLists function.

gene_lists <- loadGeneLists( "gene_lists.gmt" )

The second option is to have your gene lists in separate tab-delimited text files within a directory.

You can load the gene lists within this directory using the loadGeneLists function with the type argument set to "DIR".

gene_lists <- loadGeneLists( "/gene/lists/directory", type = "DIR" )

By default, the loadGeneLists function expects the files in this directory to have no header and to have the genes in the first column. You can modify this behavior using the header argument and the col argument.

In this example, the files contain a header that needs to be ignored, and the genes are in the 7th column.

gene_lists <- loadGeneLists( "/gene/lists/directory", type = "DIR", header = TRUE, col = 7 )

3.2 GO Terms

The second input data required by goSTAG is a set of GO terms and the genes associated with those GO terms. The GO terms should all belong to the same GO domain (either biological process, cellular component, or molecular function). The expected format is a list of vectors, where each list corresponds to a GO term, and the vectors contain the genes associated with the GO terms. The names of the lists should be the GO IDs of the GO terms. The list must also have an additional element named “ALL”, which is a vector that contains all annotated genes.

The goSTAG package provides the loadGOTerms function for loading the GO terms and gene associations from BioMart using the biomaRt package.

go_terms <- loadGOTerms()
head( lapply( go_terms, head ) )

## $`GO:0006810`
## [1] "SLC25A26" "GJA3"     "GJB2"     "HERC1"    "SLC35A2"  "BAG6"    
## 
## $`GO:0048011`
## [1] "BTC"   "FGF7"  "FURIN" "APH1B" "NDN"   "AGO2" 
## 
## $`GO:0007173`
## [1] "BTC"     "FGF7"    "AGO2"    "FGF9"    "MAPKAP1" "FGF2"   
## 
## $`GO:0008543`
## [1] "BTC"     "FGF7"    "AGO2"    "FGF9"    "MAPKAP1" "FGF2"   
## 
## $`GO:0048015`
## [1] "BTC"     "FGF7"    "UBE2C"   "AGO2"    "FGF9"    "MAPKAP1"
## 
## $`GO:0008284`
## [1] "BTC"   "FGF7"  "EDN3"  "RPS4X" "AVP"   "FGF9"

head( go_terms[["ALL"]] )

## [1] "SLC25A26" "BTC"      "GAST"     "HAP1"     "JUP"      "FKBP10"

Since loading the GO terms from BioMart can take several minutes, the default behavior of the loadGOTerms function is to use a previously generated version that has been archived in the goSTAG package, which is much faster.

You can load the most recent version of the GO terms and gene associations directly from BioMart by setting the use_archived argument to FALSE.

go_terms <- loadGOTerms( use_archived = FALSE )

By default, the loadGOTerms function will associate genes based on the annotation for human. It will also load GO terms for the biological process domain and remove any GO terms that has fewer than 5 genes associated with it. To modify this behavior, use the species, domain, and min_num_genes arguments. Currently, the supported options for the species argument are “human”, “mouse”, or “rat”. Acceptable options for the domain argument are “BP” (biological process), “CC” (cellular component), or “MF” (molecular function).

This example will load GO terms for mouse, using the molecular function domain, and remove any GO terms with fewer than 10 genes.

go_terms_mouse <- loadGOTerms( species = "mouse", domain = "MF", min_num_genes = 10 )

4.1 Generating the Enrichment Matrix

The first step is to create a matrix of GO enrichment scores. The goSTAG package calculates an enrichment score for each gene list and GO term using the hypergeometric test. It will remove any GO term that has no samples with significant enrichment. The enrichment score is calculated as the -log10 p-value of the hypergeometric test. For tests without significant enrichment, the enrichment score is set to 0.

To generate an enrichment matrix from your gene lists and GO terms, use the performGOEnrichment function.

enrichment_matrix <- performGOEnrichment( goSTAG_example_gene_lists, go_terms )
head(enrichment_matrix)

##                Day1     Day2     Day3
## GO:0007165 0.000000 0.000000 1.343007
## GO:0007369 0.000000 0.000000 1.389698
## GO:0042127 0.000000 0.000000 1.729184
## GO:0007249 1.862509 1.766861 1.654304
## GO:0008380 1.345807 0.000000 1.306134
## GO:0016049 0.000000 1.436532 0.000000

This function produces a matrix of enrichment scores, where columns are gene lists and rows are GO terms.

By default, this function will remove any GO term where none of the samples have a p-value of less than 0.05. You can specify whether to use a nominal p-value or an adjusted p-value using the filter_method argument, and you can change the significance threshold using the significance_threshold argument. This function uses the p.adjust method to adjust the p-value, and you can specify the method used to adjust the p-values with the p.adjust_method argument. See the method argument for the p.adjust function for available options.

This example generates the enrichment matrix and uses a Benjamini-Hochberg FDR of 0.3 to filter the GO terms.

enrichment_matrix_FDR <- performGOEnrichment( goSTAG_example_gene_lists, go_terms, filter_method = "p.adjust", significance_threshold = 0.3, p.adjust_method = "BH" )

4.2 Hierarchical Clustering

After generating the enrichment matrix, the next step is to cluster the GO terms. Before clustering the GO terms into groups, you must first perform hierarchical clustering.

You can perform the hierarchical clustering using the performHierarchicalClustering function.

hclust_results <- performHierarchicalClustering( enrichment_matrix )
hclust_results

## 
## Call:
## hclust(d = as.dist(1 - abs(cor(t(final_matrix)))), method = clustering_method)
## 
## Cluster method   : average 
## Number of objects: 155

This function uses the hclust function and returns an object of class hclust containing the hierarchical clustering tree.

By default, this function uses one minus the absolute value of the correlation to measure distance. You can change this behavior using the distance_method argument. Specifying anything other than correlation will cause this function to use the dist function to measure distance. The method specified in the distance_method argument is the method that the dist function will use. See the method argument for the dist function for available options.

This function uses the hclust function to perform the hierarchical clustering, with a default agglomeration method of average. You can change this by specifying the clustering_method argument.

This example performs hierarchical clustering on the GO terms using Euclidean distance and complete agglomeration.

hclust_results_euclidean <- performHierarchicalClustering( enrichment_matrix, distance_method = "euclidean", clustering_method = "complete" )

You can also perform hierarchical clustering on the samples by setting the feature argument to col.

sample_hclust_results <- performHierarchicalClustering( enrichment_matrix, feature = "col" )

4.3 Grouping the Clusters

After the hierarchical clustering is complete, you can group the GO terms into clusters using the groupClusters function.

clusters <- groupClusters( hclust_results_euclidean )
lapply( head( clusters ), head )

## $Cluster1
## [1] "GO:0007249"
## 
## $Cluster2
## [1] "GO:0043928"
## 
## $Cluster3
## [1] "GO:0043537" "GO:0000288"
## 
## $Cluster4
## [1] "GO:0006417"
## 
## $Cluster5
## [1] "GO:0045760" "GO:0090316"
## 
## $Cluster6
## [1] "GO:0002224"

This function produces a list of vectors, where each list corresponds to a cluster, and the vectors are the GO terms in the clusters.

By default, this function will group GO terms whose distance to each other is less than 0.2. This threshold can be adjusted by changing the distance_threshold argument. A larger distance threshold will produce fewer clusters with more GO terms, whereas a smaller one will produce more clusters with fewer GO terms.

The default distance threshold of 0.2 produces a moderate number of clusters.

length( clusters )

## [1] 35

Increasing the distance threshold to 0.5 produces fewer clusters, but with more GO terms on average.

clusters_larger_threshold <- groupClusters( hclust_results_euclidean, distance_threshold = 0.5 )
length( clusters_larger_threshold )

## [1] 21

Whereas decreasing the distance threshold to 0.05 produces more clusters, but with fewer GO terms on average.

clusters_smaller_threshold <- groupClusters( hclust_results_euclidean, distance_threshold = 0.05 )
length( clusters_smaller_threshold )

## [1] 60

4.4 Annotating the Clusters

The final step is to annotate each of the clusters using the GO term within each cluster’s subtree that has the most paths back to the root.

You can annotate the clusters using this method with the annotateClusters function.

cluster_labels <- annotateClusters( clusters )
head( cluster_labels )

##                                                                                            GO:0007249 
##                                                                 "I-kappaB kinase/NF-kappaB signaling" 
##                                                                                            GO:0043928 
## "exonucleolytic nuclear-transcribed mRNA catabolic process involved in deadenylation-dependent decay" 
##                                                                                            GO:0000288 
##                           "nuclear-transcribed mRNA catabolic process, deadenylation-dependent decay" 
##                                                                                            GO:0006417 
##                                                                           "regulation of translation" 
##                                                                                            GO:0090316 
##                                              "positive regulation of intracellular protein transport" 
##                                                                                            GO:0002224 
##                                                                "toll-like receptor signaling pathway"

This returns a vector containing the descriptions of the GO terms tagged to each cluster.