1 Introduction

In this vignette we assume that the reader is already familiarized with the introduction vignette, but wants to know how can it help to answer other questions in other situations

We follow the same convention about names of the objects used genesReact and genesKegg are the list with information about the pathways the genes are involved in.

2 Merging similarities

If one calculates similarities with KEGG data and Reactome or other input for the same genes or clusters BioCor provides a couple of functions to merge them.

We can set a weight to each similarity input with weighted.sum, multiply them also using a weight for each similarity (with weighted.prod), doing the mean or just adding them up. Similarities allow us to apply a function to combine the matrices of a list. Here we use some of the genes used in the first vignette:

kegg <- mgeneSim(c("672", "675", "10"), genesKegg)
react <- mgeneSim(c("672", "675", "10"), genesReact)
## We can sum it adding a weight to each origin
weighted.sum(c(kegg["672", "675"], react["672","675"]), w = c(0.3, 0.7))
## [1] 0.7247289

## Or if we want to perform for all the matrix
## A list of matrices to merge
sim <- list("kegg" = kegg, "react" = react)
similarities(sim, weighted.sum, w = c(0.3, 0.7))
##            672       675         10
## 672 1.00000000 0.7247289 0.09926055
## 675 0.72472885 1.0000000 0.05211140
## 10  0.09926055 0.0521114 1.00000000
similarities(sim, weighted.prod, w = c(0.3, 0.7))
##           672          675           10
## 672 0.2100000 0.0173101952 0.0000000000
## 675 0.0173102 0.2100000000 0.0001227335
## 10  0.0000000 0.0001227335 0.2100000000
similarities(sim, prod)
##           672          675           10
## 672 1.0000000 0.0824295011 0.0000000000
## 675 0.0824295 1.0000000000 0.0005844451
## 10  0.0000000 0.0005844451 1.0000000000
similarities(sim, mean)
##            672        675         10
## 672 1.00000000 0.54121475 0.07090039
## 675 0.54121475 1.00000000 0.03957722
## 10  0.07090039 0.03957722 1.00000000

This functions are similar to weighted.mean, except that first the multiplication by the weights is done and then the NAs are removed:

weighted.mean(c(1, NA), w = c(0.5, 0.5), na.rm = TRUE)
## [1] 1
weighted.mean(c(1, 0.5, NA), w = c(0.5, 0.25, 0.25), na.rm = TRUE)
## [1] 0.8333333
weighted.sum(c(1, NA), w = c(0.5, 0.5))
## [1] 0.5
weighted.sum(c(1, 0.5, NA), w = c(0.5, 0.25, 0.25))
## [1] 0.625
weighted.prod(c(1, NA), w = c(0.5, 0.5))
## [1] 0.5
weighted.prod(c(1, 0.5, NA), w = c(0.5, 0.25, 0.25))
## [1] 0.0625

3 Assessing a differential study

In this example we will use the functional similarities in a classical differential study.

3.1 Obtaining data

We start using data from the RNAseq workflow and following the analysis comparing treated and untreated:

suppressPackageStartupMessages(library("airway"))
data("airway")
library("DESeq2")
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, "untrt")
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res <- results(dds, alpha = 0.05)
summary(res)
## 
## out of 33469 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2211, 6.6%
## LFC < 0 (down)     : 1817, 5.4%
## outliers [1]       : 0, 0%
## low counts [2]     : 16676, 50%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plot(res$log2FoldChange, -log10(res$padj), pch = 16, xlab = "log2FC", 
     ylab = "-log10(p.ajd)", main = "Untreated vs treated")
logFC <- 2.5
abline(v=c(-logFC, logFC), h = -log10(0.05), col = "red")

Figure 1: Volcano plot
The airway data

As we can see here there are around 4000 genes differentially expressed genes, some of them differentially expressed above 2^2.5.

3.2 Selecting differentially expressed genes

Usually in such a study one selects genes above certain logFC or fold change threshold, here we use the absolute value of 2.5:

fc <- res[abs(res$log2FoldChange) >= logFC & !is.na(res$padj), ]
fc <- fc[fc$padj < 0.05, ]
# Convert Ids (used later)
genes <- select(org.Hs.eg.db, keys = rownames(res), keytype = "ENSEMBL", 
                column = c("ENTREZID", "SYMBOL"))
## 'select()' returned 1:many mapping between keys and columns
genesFC <- genes[genes$ENSEMBL %in% rownames(fc), ]
genesFC <- genesFC[!is.na(genesFC$ENTREZID), ]
genesSim <- genesFC$ENTREZID
names(genesSim) <- genesFC$SYMBOL
genesSim <- genesSim[!duplicated(genesSim)]
# Calculate the functional similarity
sim <- mgeneSim(genes = genesSim, info = genesReact, method = "BMA")
## Warning in mgeneSim(genes = genesSim, info = genesReact, method = "BMA"):
## Some genes are not in the list provided.

Once the similarities for the selected genes are calculated we can now visualize the effect of each method:

nas <- apply(sim, 1, function(x){all(is.na(x))})
sim <- sim[!nas, !nas]

MDSs <- cmdscale(1-sim)
plot(MDSs, type = "n", main = "BMA similarity", xlab = "MDS1", ylab = "MDS2")
up <- genes[genes$ENSEMBL %in% rownames(fc)[fc$log2FoldChange >= logFC], "SYMBOL"]
text(MDSs, labels =  rownames(MDSs), col = ifelse(rownames(MDSs) %in% up, "black", "red"))
abline(h = 0, v = 0)
legend("top", legend = c("Up-regulated", "Down-regulated"), fill = c("black", "red"))

Figure 2: Functional similarity of genes with logFC above 2,5
Similar genes cluster together.

This plot illustrate that some differentially expressed genes are quite similar according to their pathways. Suggesting that there might be a relationship between them. Furthermore, some up-regulated genes seem functionally related to down-regulated genes indicating a precise regulation of the pathways where they are involved.

Note that here we are only using 72 genes from the original 127.

3.3 Are differentially expressed genes selected by their functionality?

In the previous section we have seen that some differentially expressed genes are functionally related and that they have a high logFC value. Are genes differentially expressed more functional related than those which aren’t differential expressed?

For simplicity we will use a subset of 400 genes represented again in a volcano plot and we will look for the functional similarities between those genes:

set.seed(250)
subRes <- res[!is.na(res$log2FoldChange), ]
subs <- sample.int(nrow(subRes), size = 400)
subRes <- subRes[subs, ]
g <- genes[genes$ENSEMBL %in% rownames(subRes), "ENTREZID"]
gS <- mgeneSim(g[g %in% names(genesReact)], genesReact, "BMA")
deg <- rownames(subRes[subRes$padj < 0.05 & !is.na(subRes$padj), ])
keep <- rownames(gS) %in% genes[genes$ENSEMBL %in% deg, "ENTREZID"]
plot(subRes$log2FoldChange, -log10(subRes$padj), pch = 16, xlab = "log2FC", 
     ylab = "-log10(p.ajd)", main = "Untreated vs treated")
abline(v=c(-logFC, logFC), h = -log10(0.05), col = "red")

Figure 3: Volcano plot of the subset of 400 genes
This subset will be used in the following sections

We can answer this by testing it empirically:

library("boot")
# The mean of genes differentially expressed
(scoreDEG <- mean(gS[!keep, keep], na.rm = TRUE)) 
## [1] 0.1188639
b <- boot(data = gS, R = 1000, statistic= function(x, i){
    g <- !rownames(x) %in% rownames(x)[i]
    mean(x[g, i], na.rm = TRUE)
})
(p.val <- (1 + sum(b$t>scoreDEG))/1001)
## [1] 0.9220779
hist(b$t, main = "Distribution of scores", xlab = "Similarity score")
abline(v = scoreDEG, col = "red")

Figure 4: Distribution of scores between differentially expressed genes and those who aren’t
The line indicates the mean score of the similarity between differentially expressed genes and those which aren’t differentially expressed.

Comparing the genes differentially expressed and those who aren’t doesn’t show that they are non-randomly selected (p-value 0.9220779). The genes with a p-value below the threshold are not more closely functionally related than all the other genes11 From 400 genes there are 123 with pathway information and only 23 where significantly differentially expressed in this subset..

3.4 Are functionally related the selected differentially expressed genes?

We have seen that the genes differentially expressed aren’t selected by their functionality. However they could be more functionally related than the other genes. Are the differentially expressed genes more functionally similar than it would be expected ?

(scoreW <- combineScores(gS[keep, keep], "avg"))
## [1] 0.1562942
b <- boot(data = gS, R = 1000, statistic= function(x, i){
    mean(x[i, i], na.rm = TRUE)
})
(p.val <- (1 + sum(b$t>scoreW))/1001) # P-value
## [1] 0.01098901
hist(b$t, main = "Distribution of scores", xlab = "Similarity score")
abline(v = scoreW, col = "red")

Figure 5: Distribution of the similarity within differentially expressed genes
The line indicates the mean funtional similarity whitin them.

If we selected randomly the genes from our pool we would expect a score around 0.1562942 with a probability of 0.010989. That means that the differentially expressed genes is highly different compared to the other genes if we use a significance threshold of 0.05 22 Remember that this is a small subset..

3.5 Influence of the fold change in the functionally similarity of the genes

We have seen that the genes differentially expressed are not selected by functional similarity but they are functionally related. Now we would like to know if selecting a fold change threshold affects the functional similarity between them.

To know the relationship between the fold change and the similarity between genes we have several methods:

s <- seq(0, max(abs(subRes$log2FoldChange))+0.05, by = 0.05)
l <- sapply(s, function(x){
    deg <- rownames(subRes[abs(subRes$log2FoldChange) >= x, ])
    keep <- rownames(gS) %in% genes[genes$ENSEMBL %in% deg, "ENTREZID"]
    BetweenAbove <- mean(gS[keep, keep], na.rm = TRUE)
    AboveBelow <- mean(gS[keep, !keep], na.rm = TRUE)
    BetweenBelow <- mean(gS[!keep, !keep], na.rm = TRUE)
    c("BetweenAbove" = BetweenAbove, "AboveBelow" = AboveBelow, 
      "BetweenBelow" = BetweenBelow)
})
L <- as.data.frame(cbind(logfc = s, t(l)))
plot(L$logfc, L$BetweenAbove, type = "l", xlab = "abs(log2) fold change", 
     ylab = "Similarity score", 
     main = "Similarity scores along logFC", col = "darkred")
lines(L$logfc, L$AboveBelow, col = "darkgreen")
lines(L$logfc, L$BetweenBelow, col = "black")
legend("topleft", 
       legend = c("Between genes above and below threshold", 
                  "Whitin genes above threshold", 
                  "Whitin genes below threshold"), 
       fill = c("darkgreen", "darkred", "black"))

Figure 6: Similarity of genes along abs(logFC)
Assessing the similarity of genes according to their absolute log2 fold change.

The functional similarity of the genes above the threshold increases with a more restrictive threshold, indicating that a logFC threshold selects genes by their functionality, or in other words that genes differentially expressed tend to be of related pathways. The similarity between those genes above the threshold and below remains constant as well as within genes below the threshold.

l <- sapply(s, function(x){
    # Names of genes up and down regulated
    degUp <- rownames(subRes[subRes$log2FoldChange >= x, ])
    degDown <- rownames(subRes[subRes$log2FoldChange <= -x, ])
    
    # Translate to ids in gS
    keepUp <- rownames(gS) %in% genes[genes$ENSEMBL %in% degUp, "ENTREZID"]
    keepDown <- rownames(gS) %in% genes[genes$ENSEMBL %in% degDown, "ENTREZID"]
    
    # Calculate the mean similarity between each subgrup
    between <- mean(gS[keepUp, keepDown], na.rm = TRUE)
    
    c("UpVsDown" = between)
})
L <- as.data.frame(cbind("logfc" = s, "UpVsDown" = l))
plot(L$logfc, L$UpVsDown, type = "l", xlab = "abs(log2) fold change threshold", ylab = "Similarity score", 
     main = "Similarity scores along logFC")
legend("topright", legend = "Up vs down regulated genes",
       fill = "black")

Figure 7: Functional similarity between the up-regulated and down-regulated genes

The maximal functional similarity between genes up-regulated and down-regulated are at 1.15 log2 fold change.

4 Assessing a new pathway

Sometimes the top differentially expressed genes or some other key genes are selected as a signature or a potential new group of related genes. In those cases we can test how does the network of genes change if we add them. Here we create a new pathway named deg, and we see the effect on the functional similarity score for all the genes:

# Adding a new pathway "deg" to those genes
genesReact2 <- genesReact
diffGenes <- genes[genes$ENSEMBL %in% deg, "ENTREZID"]
genesReact2[diffGenes] <- sapply(genesReact[diffGenes], function(x){c(x, "deg")})
plot(ecdf(mgeneSim(names(genesReact), genesReact)))
curve(ecdf(mgeneSim(names(genesReact2), genesReact2)), color = "red")

This would take lot of time, for a illustration purpose we reduce to some genes:

library("Hmisc")
genesReact2 <- genesReact
diffGenes <- genes[genes$ENSEMBL %in% deg, "ENTREZID"]
# Create the new pathway called deg
genesReact2[diffGenes] <- sapply(genesReact[diffGenes], function(x){c(x, "deg")})
ids <- unique(genes[genes$ENSEMBL %in% rownames(subRes), "ENTREZID"])
Ecdf(c(mgeneSim(ids, genesReact, method = "BMA"),  
       mgeneSim(ids, genesReact2, method = "BMA")), 
     group = c(rep("Reactome", length(ids)^2), rep("Modified", length(ids)^2)), 
     col = c("black", "red"), xlab = "Functional similarities", main = "Empirical cumulative distribution")

Figure 8: The effect of adding a new pathway to a functional similarity
In red the same network as in black but with the added pathway.

5 Merging sources of information

Sometimes we have several origin of information, either several databases, or information from other programs… We can merge this in the single object required by the function in BioCor using the function combineSources33 See the help page of combineSources.

This functions helps to evaluate what happens when we add more pathway information. For instance here we add the information in Kegg and the information in Reactome and we visualize it using the same procedure as previously:

genesKegg <- as.list(org.Hs.egPATH)
gSK <- mgeneSim(rownames(gS), genesKegg)
mix <- combineSources(genesKegg, genesReact)
## Warning in combineSources(genesKegg, genesReact): More than 50% of genes identifiers of a source are unique
## Check the identifiers of the genes
gSMix <- mgeneSim(rownames(gS), mix)
Ecdf(c(gS, gSK, gSMix), 
     group = c(rep("Reactome", length(gS)), rep("Kegg", length(gSK)), 
               rep("Mix", length(gSMix))), 
     col = c("black", "red", "blue"), xlab = "Functional similarities", main = "ecdf")
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

Figure 9: Comparison of functional similarity in different databases

When mixed, there is a huge increase in the genes that share a pathway using the max method. Observe in next figure (10) how does the method affects to the results:

gSK2 <- mgeneSim(rownames(gS), genesKegg, method = "BMA")
gS2 <- mgeneSim(rownames(gS), genesReact, method = "BMA")
gSMix2 <- mgeneSim(rownames(gS), mix, method = "BMA")
Ecdf(c(gS2, gSK2, gSMix2), 
     group = c(rep("Reactome", length(gS)), rep("Kegg", length(gSK)), 
               rep("Mix", length(gSMix))), 
     col = c("black", "red", "blue"), xlab = "Functional similarities (BMA)", main = "ecdf")
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

Figure 10: Comparison of functional similarity in different gene sets

Now we can appreciate that most of the functional similarity is brought by Reactome database

6 miRNA analysis

miRNAs are RNAs that interact with many genes and transcripts and is subject to change with location and time, thus defining the effect of an miRNA is difficult. In this section we try to answer how functionally similar are miRNA between them to provide a help for potentially closely related miRNA.

First we look for human miRNAs and prepare them for using them as the input for the cluster functions (Restricting to 10 miRNA with less than 150 genes).44 This preparation has been adapted from a previous discussion.

library("targetscan.Hs.eg.db")
## select human mirna 
humanMirnaIdx <- grep("hsa", mappedkeys(targetscan.Hs.egMIRNA)) 
## select seed-based families for human mirna
humanMirna <- mappedkeys(targetscan.Hs.egMIRNA)[humanMirnaIdx]
## select targets of families 
humanMirnaFamilies <- unlist(mget(humanMirna, targetscan.Hs.egMIRBASE2FAMILY)) 
humanMirnaTargets <- mget(humanMirnaFamilies, revmap(targetscan.Hs.egTARGETS))
names(humanMirnaTargets) <- humanMirna
# Restrict to miRNA with more than one target and less than 150
miRNAs <- sample(humanMirnaTargets[lengths(humanMirnaTargets) > 1 &
                                       lengths(humanMirnaTargets) < 150], 10)
lengths(miRNAs)
##    hsa-miR-187 hsa-miR-127-3p   hsa-miR-551a    hsa-miR-191 hsa-miR-615-3p 
##             13             17             17            118             28 
##    hsa-miR-99b   hsa-miR-551b   hsa-miR-487b    hsa-miR-184    hsa-miR-451 
##             92             17             23             54             36

Now we calculate the functional similarity of those miRNAs using a cluster approach.

cluster1 <- mclusterSim(miRNAs, genesReact, method = "BMA")
## Warning in mclusterSim(miRNAs, genesReact, method = "BMA"): Some genes are
## not in the list provided.
knitr::kable(round(cluster1, 2), caption = "The similarity between miRNA", format = "html")

So for instance hsa-miR-551b, hsa-miR-551a are functionally related despite being from different families.

7 Comparing with GO similarities

As suggested in the main vignette functional similarities can be compared to semantic similarities such as those based on GO. Here a comparison using the biological process from the gene ontologies is done:

library("GOSemSim")
## GOSemSim v2.10.0  For help: https://guangchuangyu.github.io/GOSemSim
## 
## If you use GOSemSim in published research, please cite:
## Guangchuang Yu, Fei Li, Yide Qin, Xiaochen Bo, Yibo Wu, Shengqi Wang. GOSemSim: an R package for measuring semantic similarity among GO terms and gene products Bioinformatics 2010, 26(7):976-978
## 
## Attaching package: 'GOSemSim'
## The following objects are masked from 'package:BioCor':
## 
##     clusterSim, combineScores, geneSim, mclusterSim, mgeneSim
BP <- godata('org.Hs.eg.db', ont="BP", computeIC=TRUE)
## preparing gene to GO mapping data...
## preparing IC data...
gsGO <- GOSemSim::mgeneSim(rownames(gS), semData = BP, measure = "Resnik", verbose = FALSE)
keep <- rownames(gS) %in% rownames(gsGO)
hist(as.dist(gS[keep, keep]-gsGO), 
     main = "Difference between functional similarity and biological process", 
     xlab = "Functional similarity - biological process similarity")

Figure 11: Comparison of similarities
Functional similarities compared to biological process semantic similarity.

On this graphic we can observe that some genes have a large functional similarity and few biological similarity. They are present together in several pathways while they share few biological process, indicating that they might be key elements of the pathways they are in. On the other hand some other pairs of genes show higher biological process similarity than functional similarity, indicating specialization or compartmentalization of said genes.