Describes how to use the BioCor package to answer several biological questions and how to use functional similarities with other measures.
BioCor 1.2.1
In this vignette we assume that the reader is already familiarized with the other vignette. But want 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.
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.0000000 0.72472885 0.09903930
## 675 0.7247289 1.00000000 0.05156284
## 10 0.0990393 0.05156284 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.0001213771
## 10 0.0000000 0.0001213771 0.2100000000
similarities(sim, prod)
## 672 675 10
## 672 1.0000000 0.0824295011 0.0000000000
## 675 0.0824295 1.0000000000 0.0005779864
## 10 0.0000000 0.0005779864 1.0000000000
similarities(sim, mean)
## 672 675 10
## 672 1.00000000 0.54121475 0.07074236
## 675 0.54121475 1.00000000 0.03918538
## 10 0.07074236 0.03918538 1.00000000
This functions are similar to weighted.mean
, except that first the multiplication by the weights is done and then the NA
s 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
In this example we will use the functional similarities in a classical differential study.
We start using data from the RNAseq workflow and following the analysis:
suppressPackageStartupMessages(library("airway"))
data("airway")
library("DESeq2")
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")
As we can see here there are around 4000 genes differentially expressed genes, some of them differentially expressed above 2^2.5.
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"))
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, somre 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 73 genes from the original 127.
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?
Here we will use a subset of 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, ]
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")
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"]
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.1156056
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.4575425
hist(b$t, main = "Distribution of scores", xlab = "Similarity score")
abline(v = scoreDEG, col = "red")
Comparing the genes differentially expressed and those who aren’t doesn’t show that they are non-randomly selected (p-value 0.4575425). The genes with a p-value below the threshold are not more closely functionally related than all the other genes1 From 400 genes there are 111 with pathway information and only 19 where significantly differentially expressed in this subset..
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"))
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")
The maximal functional similarity between genes up-regulated and down-regulated are at 0.9 log2 fold change.
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 ilustration 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")
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 combineSources
3 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")
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")
Now we can appreciate that most of the functional similarity is brought by Reactome database
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).4 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-100 hsa-miR-451 hsa-miR-184 hsa-miR-615-3p hsa-miR-191
## 92 36 54 28 118
## hsa-miR-99b hsa-miR-296-3p hsa-miR-487b hsa-miR-99a hsa-miR-551a
## 92 125 23 92 17
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")
hsa-miR-100 | hsa-miR-451 | hsa-miR-184 | hsa-miR-615-3p | hsa-miR-191 | hsa-miR-99b | hsa-miR-296-3p | hsa-miR-487b | hsa-miR-99a | hsa-miR-551a | |
---|---|---|---|---|---|---|---|---|---|---|
hsa-miR-100 | 1.00 | 0.79 | 0.78 | 0.58 | 0.75 | 1.00 | 0.82 | 0.72 | 1.00 | 0.50 |
hsa-miR-451 | 0.79 | 1.00 | 0.75 | 0.66 | 0.76 | 0.79 | 0.80 | 0.69 | 0.79 | 0.53 |
hsa-miR-184 | 0.78 | 0.75 | 1.00 | 0.60 | 0.69 | 0.78 | 0.74 | 0.61 | 0.78 | 0.47 |
hsa-miR-615-3p | 0.58 | 0.66 | 0.60 | 1.00 | 0.66 | 0.58 | 0.60 | 0.72 | 0.58 | 0.73 |
hsa-miR-191 | 0.75 | 0.76 | 0.69 | 0.66 | 1.00 | 0.75 | 0.77 | 0.66 | 0.75 | 0.51 |
hsa-miR-99b | 1.00 | 0.79 | 0.78 | 0.58 | 0.75 | 1.00 | 0.82 | 0.72 | 1.00 | 0.50 |
hsa-miR-296-3p | 0.82 | 0.80 | 0.74 | 0.60 | 0.77 | 0.82 | 1.00 | 0.72 | 0.82 | 0.50 |
hsa-miR-487b | 0.72 | 0.69 | 0.61 | 0.72 | 0.66 | 0.72 | 0.72 | 1.00 | 0.72 | 0.71 |
hsa-miR-99a | 1.00 | 0.79 | 0.78 | 0.58 | 0.75 | 1.00 | 0.82 | 0.72 | 1.00 | 0.50 |
hsa-miR-551a | 0.50 | 0.53 | 0.47 | 0.73 | 0.51 | 0.50 | 0.50 | 0.71 | 0.50 | 1.00 |
So for instance hsa-miR-99b, hsa-miR-99a, hsa-miR-100 are functionally related despite being from different families.
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.4.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")
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 compartimentalization of said genes.
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] GOSemSim_2.4.0 targetscan.Hs.eg.db_0.6.1
## [3] Hmisc_4.0-3 ggplot2_2.2.1
## [5] Formula_1.2-2 survival_2.41-3
## [7] lattice_0.20-35 boot_1.3-20
## [9] DESeq2_1.18.0 airway_0.112.0
## [11] SummarizedExperiment_1.8.0 DelayedArray_0.4.0
## [13] matrixStats_0.52.2 GenomicRanges_1.30.0
## [15] GenomeInfoDb_1.14.0 reactome.db_1.62.0
## [17] org.Hs.eg.db_3.4.2 AnnotationDbi_1.40.0
## [19] IRanges_2.12.0 S4Vectors_0.16.0
## [21] Biobase_2.38.0 BiocGenerics_0.24.0
## [23] BioCor_1.2.1 BiocStyle_2.6.0
##
## loaded via a namespace (and not attached):
## [1] bit64_0.9-7 splines_3.4.2 highr_0.6
## [4] latticeExtra_0.6-28 blob_1.1.0 GenomeInfoDbData_0.99.1
## [7] yaml_2.1.14 RSQLite_2.0 backports_1.1.1
## [10] digest_0.6.12 RColorBrewer_1.1-2 XVector_0.18.0
## [13] checkmate_1.8.5 colorspace_1.3-2 htmltools_0.3.6
## [16] Matrix_1.2-11 plyr_1.8.4 XML_3.98-1.9
## [19] pkgconfig_2.0.1 genefilter_1.60.0 bookdown_0.5
## [22] zlibbioc_1.24.0 xtable_1.8-2 GO.db_3.4.2
## [25] scales_0.5.0 BiocParallel_1.12.0 annotate_1.56.0
## [28] tibble_1.3.4 htmlTable_1.9 nnet_7.3-12
## [31] lazyeval_0.2.1 magrittr_1.5 memoise_1.1.0
## [34] evaluate_0.10.1 foreign_0.8-69 tools_3.4.2
## [37] data.table_1.10.4-3 stringr_1.2.0 locfit_1.5-9.1
## [40] munsell_0.4.3 cluster_2.0.6 compiler_3.4.2
## [43] rlang_0.1.2 grid_3.4.2 RCurl_1.95-4.8
## [46] htmlwidgets_0.9 bitops_1.0-6 base64enc_0.1-3
## [49] rmarkdown_1.6 gtable_0.2.0 codetools_0.2-15
## [52] DBI_0.7 gridExtra_2.3 knitr_1.17
## [55] bit_1.1-12 rprojroot_1.2 stringi_1.1.5
## [58] Rcpp_0.12.13 geneplotter_1.56.0 rpart_4.1-11
## [61] acepack_1.4.1