Retrieving Orthologs
Rohit Satyam, Alberto Maillo, David Gomez-Cabrero, Arnab Pain
King Abdullah University of Science & Technology, Saudi Arabia22 August, 2026
Source:vignettes/retrieving_orthologs.Rmd
retrieving_orthologs.RmdAbstract
This is a quick and dirty tutorial on how to fetch orthologs from OrthoMCL and InparanoiDB databases.
Introduction
Users might require Ortholog mapping to:
to assess phyletic profile of a set of genes and see if they have orthologs in other related organism or not.
integrate datasets obtained from two different species. Often such datasets are scRNASeq (eg: (Tebben et al. 2022))
borrow annotations such as localization for hypothesis generation.
perform co-expression analysis between two different species.
Whatever the use case maybe, we will see how we can exploit
getpairedOrthologs() to fetch orthologs of multiple
organisms of interest in batch manner and integrate them.
With release 71 of VEuPathDB, only users with subscription and API
key can access the database programmatically. Use
usethis::edit_r_environ(scope = "user") to add VEuPathDB
API using VEUPATHDB_API_KEY variable. By default, all the
functions look for VEUPATHDB_API_KEY variable for API key.
Should you store your key using any other variable, supply it using
api_key argument in VEuPathDB asscoiated functions.
# Load package and some other useful packages by using
suppressPackageStartupMessages(
suppressWarnings({
library(plasmoRUtils)
library(dplyr)
library(plyr)}))Here, we will try to get all Plasmodium falciparum 3D7 orthologs present in Toxoplasma gondii ME49, Plasmodium berghei ANKA and Plasmodium vivax Sal1.
#?getpairedOrthologs
ids <- listOrthomcl()
numids <- c(ids %>% dplyr::filter(Organism=="Plasmodium vivax Sal-1") %>% dplyr::pull(ID),
ids %>% dplyr::filter(Organism=="Toxoplasma gondii ME49") %>% dplyr::pull(ID),
ids %>% dplyr::filter(Organism=="Plasmodium berghei ANKA") %>% dplyr::pull(ID))
res <- lapply(numids, function(x){
getpairedOrthologs(from = ids %>% dplyr::filter(Organism=="Plasmodium falciparum 3D7") %>% dplyr::pull(ID),
to=x,
db="orthomcl",transform = FALSE)
}) %>% setNames(c("Pvivax","Tgme49","Pb"))
merged_df <- Reduce(function(x, y) merge(x, y, by = "Accession", all = TRUE), res)
merged_df %>% head()
#> Accession Target ID.x Target ID.y Target ID
#> 1 10181 No Ortholog Found No Ortholog Found No Ortholog Found
#> 2 10190 No Ortholog Found No Ortholog Found No Ortholog Found
#> 3 10214 No Ortholog Found No Ortholog Found No Ortholog Found
#> 4 1064799 No Ortholog Found No Ortholog Found No Ortholog Found
#> 5 1067437 No Ortholog Found No Ortholog Found No Ortholog Found
#> 6 1067464 No Ortholog Found No Ortholog Found No Ortholog Found
## Tidying up the dataframe by removing organisms prefixes
strip_prefix <- function(x) sub("^[^|]+\\|", "", x) # per-token
clean_list <- function(s) {
toks <- trimws(strsplit(s, ",", fixed = TRUE)[[1]])
paste(vapply(toks, strip_prefix, character(1)), collapse = ",")
}
merged_df[] <- lapply(merged_df, function(col) vapply(col, clean_list, character(1)))
## Changing column names
colnames(merged_df) <- c("P falciparum 3D7","P vivax Sal1","T gondii ME49","Plasmodium berghei ANKA")
merged_df %>% head()
#> P falciparum 3D7 P vivax Sal1 T gondii ME49 Plasmodium berghei ANKA
#> 1 10181 No Ortholog Found No Ortholog Found No Ortholog Found
#> 2 10190 No Ortholog Found No Ortholog Found No Ortholog Found
#> 3 10214 No Ortholog Found No Ortholog Found No Ortholog Found
#> 4 1064799 No Ortholog Found No Ortholog Found No Ortholog Found
#> 5 1067437 No Ortholog Found No Ortholog Found No Ortholog Found
#> 6 1067464 No Ortholog Found No Ortholog Found No Ortholog FoundSince OrthoMCL and other VEuPathDB database are not updated simultaneously, chances are that you might be using old OrthoMCL ID. If that’s the case above function will not be much helpful. The workaround is to fetch the old OrthoMCL IDs for each organisms and their respective databases and use it as an anchor to combine and collapse gene IDs from all the organisms. The code snippet below demonstrate the same.
## List all organisms
list <- listVeupathdb(customFields = c("primary_key","project_id","species",'species_ncbi_tax_id'))
## Subset organisms I am interested in
dbs <- list[grep("3D7|ME49$|Sal-1", list$Organism),]
dbs
#> # A tibble: 3 × 4
#> Organism project_id Species Species NCBI taxon I…¹
#> <chr> <chr> <chr> <dbl>
#> 1 Plasmodium falciparum 3D7 PlasmoDB Plasmodium falcip… 5833
#> 2 Plasmodium vivax Sal-1 PlasmoDB Plasmodium vivax 5855
#> 3 Toxoplasma gondii ME49 ToxoDB Toxoplasma gondii 5811
#> # ℹ abbreviated name: ¹`Species NCBI taxon ID`
df <- lapply(1:nrow(dbs), function(x){
plasmoRUtils::getTable(org=dbs[x,]$Organism, db=tolower(dbs[x,]$project_id),customFields = c("primary_key" ,"gene_orthomcl_name"))
})
## Retain protein coding genes
df2 <- lapply(df, function(x){
x[!(stringr::str_detect(pattern = "N/A",string = x$`Ortholog Group`)),]
})
## Combine all the tables
merged_df2 <- Reduce(function(x, y) merge(x, y, by = "Ortholog Group", all = TRUE), df2)
merged_df2 %>% head(n = 10)
#> Ortholog Group Gene ID.x Gene ID.y Gene ID
#> 1 OG8_0000001 PF3D7_1315100 PVX_122575 TGME49_243500
#> 2 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_237890
#> 3 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_255710
#> 4 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_262540
#> 5 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_305860
#> 6 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_225490
#> 7 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_235750
#> 8 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_239420
#> 9 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_316150
#> 10 OG8_0000003 PF3D7_1123100 PVX_095165 TGME49_301440
## Combining all the gene IDs in each column so that we have 1 row per orthogroup
merged_df2 <- merged_df2 %>%
dplyr::group_by( `Ortholog Group`) %>%
dplyr::summarise(dplyr::across(dplyr::everything(),
~paste(unique(.x), collapse = ",")), .groups = "drop")
## Changing column names
colnames(merged_df2) <- c("Orthogroup ID","P vivax Sal1","P falciparum 3D7","T gondii ME49")
merged_df2 %>% head(n = 10)
#> # A tibble: 10 × 4
#> `Orthogroup ID` `P vivax Sal1` `P falciparum 3D7` `T gondii ME49`
#> <chr> <chr> <chr> <chr>
#> 1 OG8_0000001 PF3D7_1315100 PVX_122575 TGME49_243500
#> 2 OG8_0000003 PF3D7_1123100,PF3D7_07045… PVX_095165,PVX_08… TGME49_237890,…
#> 3 OG8_0000008 PF3D7_1454300,PF3D7_14413… PVX_117765,PVX_11… TGME49_291050,…
#> 4 OG8_0000011 PF3D7_0517300,PF3D7_09330… PVX_091955,PVX_09… TGME49_208970,…
#> 5 OG8_0000016 NA NA TGME49_226300
#> 6 OG8_0000023 PF3D7_0922900 PVX_099555 TGME49_217740
#> 7 OG8_0000045 NA NA TGME49_306450
#> 8 OG8_0000051 NA NA TGME49_229140,…
#> 9 OG8_0000056 PF3D7_0422000 PVX_090140 TGME49_271888
#> 10 OG8_0000060 NA NA TGME49_234570,…Fetching orthologs from InParanoiDB
InParanoiDB uses NCBI taxon IDs as unique identifiers. For P. falciparum (36329), T. gondii (5811), P. vivax (126793) and P. berghei (5823), we will fetch pairwise orthologs by keeping P. falciparum as constant species. We can recycle the code above as follows.
#?getpairedOrthologs
# listipdb()
res <- lapply(c(126793,5811,5823), function(x){
getpairedOrthologs(from = 36329,
to=x,
db="ipdb",transform = TRUE)
}) %>% setNames(c("Pvivax","Tgme49","Pb"))
res[[1]] %>% head()
#> # A tibble: 6 × 11
#> `query_Group-id` query_Bitscore query_Species query_Inparalog-scor…¹
#> <dbl> <chr> <chr> <chr>
#> 1 1 8359 Plasmodium falciparum … 1
#> 2 2 8119 Plasmodium falciparum … 1
#> 3 3 7587 Plasmodium falciparum … 1
#> 4 4 7155 Plasmodium falciparum … 1
#> 5 5 6891 Plasmodium falciparum … 1
#> 6 6 6507 Plasmodium falciparum … 1
#> # ℹ abbreviated name: ¹`query_Inparalog-score`
#> # ℹ 7 more variables: `query_Protein-name` <chr>, `query_Seed-score` <chr>,
#> # target_Bitscore <chr>, target_Species <chr>,
#> # `target_Inparalog-score` <chr>, `target_Protein-name` <chr>,
#> # `target_Seed-score` <chr>Combining the results across InParanoiDB, is not as straightforward as for OrthoMCL because of lack of unique Orthogroup ID. However, user can use Query Uniprot ID as anchor, combine the results, convert Uniprot IDs to gene IDs and then collapse rows since two or more Uniprot IDs might correspond to same gene IDs.
Fetching Preconfigured tables from OrthoMCL-DB and map old OG IDs to new IDs
OrthoMCL gets frequently updated and with every update come new
Ortholog IDs. When using published dataset, its highly likely that users
would want to map the old orthogroups to new orthogroups or assess which
proteins changed the orthogroups and got clubbed, which orthogroup got
split etc. getPreconfiguredTableOrthomcl can help you do
this by accessing preconfigured table used by OrthoMCL to do that.
## Fetching all OG IDs that maps to IDs of interest
## New to old IDs. Old Ids are present in second and third column
ids = c("OG7_0008348", "OG7_0003896")
getPreconfiguredTableOrthomcl(ids,customField = "previousGroups")
#> # A tibble: 5 × 3
#> `Ortholog Group...1` `Ortholog Group...2` `Previous Ortholog Groups`
#> <chr> <chr> <chr>
#> 1 OG8_0001453 OG8_0001453 OG7_0008348
#> 2 OG8_0002935 OG8_0002935 OG7_0003896
#> 3 OG8_0010463 OG8_0010463 OG7_0008348
#> 4 OG8_0010910 OG8_0010910 OG7_0008348
#> 5 OG8_0033410 OG8_0033410 OG7_0008348
## Old Id. New IDs are present in first column
getPreconfiguredTableOrthomcl("OG3_10277")
#> # A tibble: 1 × 3
#> `Ortholog Group...1` `Ortholog Group...2` `Previous Ortholog Groups`
#> <chr> <chr> <chr>
#> 1 OG8_0002935 OG8_0002935 OG3_10277By changing customField arguments, you can also access
other preconfigured tables from the database. Eg. if you want keyword
frequencies and their associated orthogroups, this can be achieved
using
getPreconfiguredTableOrthomcl(customField = "DomainFrequency")
#> # A tibble: 229,403 × 3
#> `Ortholog Group` keyword frequency
#> <chr> <chr> <dbl>
#> 1 OG8_0000000 Protein kinase domain 0.909
#> 2 OG8_0000000 Protein tyrosine and serine/threonine kinase 0.0909
#> 3 OG8_0000001 Protein kinase domain 0.948
#> 4 OG8_0000001 Protein tyrosine and serine/threonine kinase 0.0174
#> 5 OG8_0000002 NAF domain 0.538
#> 6 OG8_0000002 Protein kinase domain 0.985
#> 7 OG8_0000002 Protein tyrosine and serine/threonine kinase 0.00769
#> 8 OG8_0000003 AAA-ATPase Vps4-associated protein 1 0.0000566
#> 9 OG8_0000003 AAA domain 0.000226
#> 10 OG8_0000003 Acetyltransferase (GNAT) family 0.000962
#> # ℹ 229,393 more rowsSession Info
utils::sessionInfo()
#> R version 4.4.1 (2024-06-14 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_India.utf8 LC_CTYPE=English_India.utf8
#> [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C
#> [5] LC_TIME=English_India.utf8
#>
#> time zone: Africa/Khartoum
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] plyr_1.8.9 dplyr_1.2.1 plasmoRUtils_1.1.2 rlang_1.3.0
#> [5] readr_2.2.0 janitor_2.2.1 BiocStyle_2.32.1
#>
#> loaded via a namespace (and not attached):
#> [1] segmented_2.2-1 fs_2.1.0
#> [3] ProtGenerics_1.36.0 matrixStats_1.5.0
#> [5] bitops_1.1-0 lubridate_1.9.5
#> [7] pRoloc_1.44.1 httr_1.4.8
#> [9] RColorBrewer_1.1-3 doParallel_1.0.17
#> [11] tools_4.4.1 MSnbase_2.30.1
#> [13] utf8_1.2.6 R6_2.6.1
#> [15] lazyeval_0.2.3 withr_3.0.3
#> [17] prettyunits_1.2.0 gridExtra_2.3.1
#> [19] preprocessCore_1.66.0 cli_3.6.6
#> [21] Biobase_2.64.0 textshaping_1.0.5
#> [23] gt_1.3.0 sass_0.4.10
#> [25] topGO_2.56.0 mvtnorm_1.4-2
#> [27] S7_0.2.2 randomForest_4.7-1.2
#> [29] proxy_0.4-29 pkgdown_2.2.1
#> [31] Rsamtools_2.20.0 systemfonts_1.3.2
#> [33] txdbmaker_1.0.1 AnnotationForge_1.46.0
#> [35] dichromat_2.0-1 parallelly_1.48.0
#> [37] limma_3.60.6 rstudioapi_0.19.0
#> [39] impute_1.78.0 RSQLite_3.53.3
#> [41] FNN_1.1.4.1 generics_0.1.4
#> [43] BiocIO_1.14.0 vroom_1.7.1
#> [45] gtools_3.9.5 dendextend_1.19.1
#> [47] GO.db_3.19.1 Matrix_1.7-5
#> [49] MALDIquant_1.22.3 drawProteins_1.24.0
#> [51] S4Vectors_0.42.1 abind_1.4-8
#> [53] lifecycle_1.0.5 yaml_2.3.12
#> [55] snakecase_0.11.1 SummarizedExperiment_1.34.0
#> [57] recipes_1.3.3 SparseArray_1.4.8
#> [59] BiocFileCache_2.12.0 grid_4.4.1
#> [61] blob_1.3.0 promises_1.5.0
#> [63] crayon_1.5.3 PSMatch_1.8.0
#> [65] lattice_0.22-9 beachmat_2.20.0
#> [67] GenomicFeatures_1.56.0 annotate_1.82.0
#> [69] chromote_0.5.1 mzR_2.38.0
#> [71] KEGGREST_1.44.1 pillar_1.11.1
#> [73] knitr_1.51 GenomicRanges_1.56.2
#> [75] rjson_0.2.23 lpSolve_5.6.23
#> [77] future.apply_1.20.2 codetools_0.2-20
#> [79] mgsub_2.0.0 glue_1.8.1
#> [81] pcaMethods_1.96.0 data.table_1.18.4
#> [83] MultiAssayExperiment_1.30.3 vctrs_0.7.3
#> [85] png_0.1-9 gtable_0.3.6
#> [87] kernlab_0.9-33 cachem_1.1.0
#> [89] gower_1.0.2 xfun_0.60
#> [91] prodlim_2026.03.11 S4Arrays_1.4.1
#> [93] coda_0.19-4.1 survival_3.8-6
#> [95] ncdf4_1.24 timeDate_4052.112
#> [97] SingleCellExperiment_1.26.0 iterators_1.0.14
#> [99] hardhat_1.4.3 lava_1.9.2
#> [101] statmod_1.5.2 MLInterfaces_1.84.0
#> [103] ipred_0.9-15 nlme_3.1-169
#> [105] bit64_4.8.2 progress_1.2.3
#> [107] filelock_1.0.3 LaplacesDemon_16.1.8
#> [109] GenomeInfoDb_1.40.1 bslib_0.12.0
#> [111] affyio_1.74.0 irlba_2.3.7
#> [113] rpart_4.1.27 otel_0.2.0
#> [115] colorspace_2.1-3 BiocGenerics_0.50.0
#> [117] DBI_1.3.0 nnet_7.3-21
#> [119] tidyselect_1.2.1 processx_3.9.0
#> [121] bit_4.6.0 compiler_4.4.1
#> [123] curl_7.1.0 rvest_1.0.5
#> [125] httr2_1.3.0 graph_1.82.0
#> [127] SparseM_1.84-2 xml2_1.6.0
#> [129] desc_1.4.3 DelayedArray_0.30.1
#> [131] plotly_4.12.1 bookdown_0.47
#> [133] rtracklayer_1.64.0 scales_1.4.0
#> [135] hexbin_1.28.6 affy_1.82.0
#> [137] rappdirs_0.3.4 stringr_1.6.0
#> [139] digest_0.6.39 mixtools_2.0.0.1
#> [141] rmarkdown_2.31 XVector_0.44.0
#> [143] htmltools_0.5.9 pkgconfig_2.0.3
#> [145] SingleR_2.6.0 sparseMatrixStats_1.16.0
#> [147] MatrixGenerics_1.16.0 dbplyr_2.6.0
#> [149] fastmap_1.2.0 htmlwidgets_1.6.4
#> [151] UCSC.utils_1.0.0 DelayedMatrixStats_1.26.0
#> [153] farver_2.1.2 jquerylib_0.1.4
#> [155] jsonlite_2.0.0 BiocParallel_1.38.0
#> [157] mclust_6.1.3 mzID_1.42.0
#> [159] ModelMetrics_1.2.2.2 BiocSingular_1.20.0
#> [161] RCurl_1.98-1.19 magrittr_2.0.5
#> [163] scuttle_1.14.0 GenomeInfoDbData_1.2.12
#> [165] Rcpp_1.1.2 viridis_0.6.5
#> [167] MsCoreUtils_1.16.1 vsn_3.72.0
#> [169] pROC_1.19.0.1 stringi_1.8.7
#> [171] zlibbioc_1.50.0 MASS_7.3-65
#> [173] listenv_1.0.0 parallel_4.4.1
#> [175] Biostrings_2.72.1 splines_4.4.1
#> [177] hms_1.1.4 igraph_2.3.3
#> [179] QFeatures_1.14.2 reshape2_1.4.5
#> [181] biomaRt_2.60.1 stats4_4.4.1
#> [183] ScaledMatrix_1.12.0 XML_3.99-0.23
#> [185] evaluate_1.0.5 BiocManager_1.30.27
#> [187] tzdb_0.5.0 foreach_1.5.2
#> [189] tidyr_1.3.2 purrr_1.2.2
#> [191] future_1.75.0 clue_0.3-68
#> [193] bio3d_2.4-5 ggplot2_4.0.3
#> [195] rsvd_1.0.5 xtable_1.8-8
#> [197] restfulr_0.0.17 AnnotationFilter_1.28.0
#> [199] easyPubMed_3.1.6 e1071_1.7-17
#> [201] later_1.4.8 viridisLite_0.4.3
#> [203] class_7.3-24 ragg_1.5.2
#> [205] tibble_3.3.1 websocket_1.4.4
#> [207] memoise_2.0.1 AnnotationDbi_1.66.0
#> [209] GenomicAlignments_1.40.0 IRanges_2.38.1
#> [211] cluster_2.1.8.2 globals_0.19.1
#> [213] timechange_0.4.0 caret_7.0-1
#> [215] sampling_2.11References
Tebben, Kieran, Aliou Dia, and David Serre. 2022. “Determination
of the Stage Composition of Plasmodium Infections
from Bulk Gene Expression Data.” mSystems 7 (4). https://doi.org/10.1128/msystems.00258-22.