1 Introduction

Methods to find similarities have been developed for several purposes, being Jaccard and Dice similarities the most known. In bioinformatics much of the research on the topic is centered around Gene Ontologies because they provide controlled vocabularies, as part of their mission:

The mission of the GO Consortium is to develop an up-to-date, comprehensive, computational model of biological systems, from the molecular level to larger pathways, cellular and organism-level systems.

However, there is another resource of similarities between genes: metabolic pathways. Metabolic pathways describe the relationship between genes, proteins, lipids and other elements of the cells. A pathway describes, to some extent, the function in which it is involved in the cell. There exists several databases about which gene belong to which pathway. Together with pathways, gene sets related to a function or to a phenotype are a source of information of the genes function. With this package we provide the methods to calculate functional similarities based on this information.

Here we provides functions to calculate functional similarities for pathways, gene sets, genes and clusters of genes.

As it development started aiming to improve clustering of genes by functionality in co-expression networks using WGCNA it also has some functions to combine similarities.

2 Citation

The main article describing the software and its usefulness is currently under writing.

3 Installation

The BioCor package is available at Bioconductor and can be downloaded and installed via biocLite:

source("http://bioconductor.org/biocLite.R")
biocLite("BioCor")

You can install the latest version of BioCor from Github with:

library("devtools")
install_github("llrs/BioCor")

4 Using BioCor

4.1 Preparation

We can load the package and prepare the data for which we want to calculate the similarities:

library("BioCor")
## 
## If you use BioCor in published research, please cite:
## Load libraries with the data of the pathways
library("org.Hs.eg.db")
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
##     as.data.frame, cbind, colMeans, colSums, colnames, do.call,
##     duplicated, eval, evalq, get, grep, grepl, intersect,
##     is.unsorted, lapply, lengths, mapply, match, mget, order, paste,
##     pmax, pmax.int, pmin, pmin.int, rank, rbind, rowMeans, rowSums,
##     rownames, sapply, setdiff, sort, table, tapply, union, unique,
##     unsplit, which, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
library("reactome.db")
genesKegg <- as.list(org.Hs.egPATH)
genesReact <- as.list(reactomeEXTID2PATHID)

To avoid having biased data it is important to have all the data about the pathways and genes associated to all pathways for organism under study. Here we assume that we are interested in human pathways. We use this two databases KEGG and Reactome as they are easy to obtain the data. However KEGG database is no longer free for large retrievals therefore it is not longer updated in the Bioconductor annotation packages.

However, one can use any list where the names of the list are the genes and the elements of the list the pathways or groups where the gene belong. One could also read from a GMT file or use GeneSetCollections in addition or instead of those associations from a pathway database and convert it to list using:

library("GSEABase")
paths2Genes <- geneIds(getGmt("/path/to/c3.all.v2.5.symbols.gmt",
                 collectionType=BroadCollection(category="c3"),
                 geneIdType=SymbolIdentifier()))

genes <- unlist(paths2Genes, use.names = FALSE)
pathways <- rep(names(paths2Genes), lengths(paths2Genes))
genes2paths <- split(pathways, genes)

gmt_sim <- mgeneSim(names(genes2paths), genes2paths)

4.2 Pathway similarities

We can compute similarities (Dice similarity, see question 1 of FAQ) between two pathways or between several pathways and combine them, or not:

(paths <- sample(unique(unlist(genesReact)), 2))
## [1] "R-HSA-1834949" "R-CEL-1834949"
pathSim(paths[1], paths[2], genesReact)
## [1] 0

(pathways <- sample(unique(unlist(genesReact)), 10))
##  [1] "R-SSC-2465910" "R-HSA-5658034" "R-SSC-349425"  "R-CEL-196299" 
##  [5] "R-HSA-69613"   "R-XTR-76005"   "R-RNO-425410"  "R-HSA-5689896"
##  [9] "R-DME-390471"  "R-XTR-2173789"
mpathSim(pathways, genesReact)
##               R-SSC-2465910 R-HSA-5658034 R-SSC-349425 R-CEL-196299
## R-SSC-2465910             1             0            0            0
## R-HSA-5658034             0             1            0            0
## R-SSC-349425              0             0            1            0
## R-CEL-196299              0             0            0            1
## R-HSA-69613               0             0            0            0
## R-XTR-76005               0             0            0            0
## R-RNO-425410              0             0            0            0
## R-HSA-5689896             0             0            0            0
## R-DME-390471              0             0            0            0
## R-XTR-2173789             0             0            0            0
##               R-HSA-69613 R-XTR-76005 R-RNO-425410 R-HSA-5689896
## R-SSC-2465910  0.00000000           0            0    0.00000000
## R-HSA-5658034  0.00000000           0            0    0.00000000
## R-SSC-349425   0.00000000           0            0    0.00000000
## R-CEL-196299   0.00000000           0            0    0.00000000
## R-HSA-69613    0.98245614           0            0    0.08510638
## R-XTR-76005    0.00000000           1            0    0.00000000
## R-RNO-425410   0.00000000           0            1    0.00000000
## R-HSA-5689896  0.08510638           0            0    1.00000000
## R-DME-390471   0.00000000           0            0    0.00000000
## R-XTR-2173789  0.00000000           0            0    0.00000000
##               R-DME-390471 R-XTR-2173789
## R-SSC-2465910            0             0
## R-HSA-5658034            0             0
## R-SSC-349425             0             0
## R-CEL-196299             0             0
## R-HSA-69613              0             0
## R-XTR-76005              0             0
## R-RNO-425410             0             0
## R-HSA-5689896            0             0
## R-DME-390471             1             0
## R-XTR-2173789            0             1

When the method to combine the similarities is set to NULL mpathSim returns a matrix of pathway similarities, otherwise it combines the values. In the next section we can see the methods to combine pathway similarities.

4.2.1 Combining values

To combine values we provide a function with several methods:

sim <- mpathSim(pathways, genesReact)
methodsCombineScores <- c("avg", "max", "rcmax", "rcmax.avg", "BMA")
sapply(methodsCombineScores, combineScores, scores = sim)
##       avg       max     rcmax rcmax.avg       BMA 
## 0.1015267 1.0000000 0.1085106 0.9982456 0.9982456

We can also specify the method to combine the similarities in mpathSim, geneSim, mgeneSim, clusterSim, mclusterSim, clusterGeneSim and mclusterGeneSim, argument method. By default the method is set to “max” to combine pathways (except in mpathSim where the default is to show all the pathway similarities) and “BMA” to combine similarities of genes or for cluster analysis. This function is adapted from GOSemSim package.

4.3 Gene similarities

To compare the function of two genes there is the geneSim function and mgeneSim function for several comparisons. In this example we compare the genes BRCA1 and BRCA2 and NAT2, which are the genes 672, 675 and 10 respectively in ENTREZID:

geneSim("672", "675", genesKegg)
## [1] 0.0824295
geneSim("672", "675", genesReact)
## [1] 1

mgeneSim(c("BRCA1" = "672", "BRCA2" = "675", "NAT2" = "10"), genesKegg)
##           BRCA1       BRCA2        NAT2
## BRCA1 1.0000000 0.082429501 0.000000000
## BRCA2 0.0824295 1.000000000 0.008241758
## NAT2  0.0000000 0.008241758 1.000000000
mgeneSim(c("BRCA1" = "672", "BRCA2" = "675", "NAT2" = "10"), genesReact)
##           BRCA1      BRCA2       NAT2
## BRCA1 1.0000000 1.00000000 0.14579350
## BRCA2 1.0000000 1.00000000 0.07011901
## NAT2  0.1457935 0.07011901 1.00000000

Note that for the same genes each database or list provided has different annotations, which result on different similarity scores. In this example BRCA1 has 3 and 13 pathways in KEGG and Reactome respectively and BRCA2 has 1 and 36 pathways in KEGG and Reactome respectively which results on different scores.

4.4 Gene cluster similarities

There are two methods:

  • Combining all the pathways for each cluster and compare between them.
  • Calculate the similarity between genes of a cluster and the other cluster.

4.4.1 By pathways similarities

As explained, in this method all the pathways of a cluster are compared with all the pathways of the other cluster. If a method to combine pathways similarities is not provided, all pathway similarities are returned:

clusterSim(c("672", "675"), c("100", "10", "1"), genesKegg)
## [1] 0.04210526
clusterSim(c("672", "675"), c("100", "10", "1"), genesKegg, NULL)
##            00230       01100       05340 00232 00983
## 04120 0.00000000 0.000000000 0.011764706     0     0
## 03440 0.04210526 0.006908463 0.000000000     0     0
## 05200 0.00000000 0.008241758 0.005540166     0     0
## 05212 0.00000000 0.001666667