library(data.table) rlog::log_info("Connecting to Ensembl API") # Object to access the Ensembl API. ensembl <- biomaRt::useEnsembl("ensembl", version = 105) # Retrieve species information. rlog::log_info("Retrieving species information") ensembl_datasets <- data.table(biomaRt::listDatasets(ensembl)) # Filter out species ID and name from the result. species <- ensembl_datasets[, .( id = stringr::str_match(dataset, "(.*)_gene_ensembl")[, 2], name = stringr::str_match(description, "(.*) genes \\(.*\\)")[, 2] )] # List of assemblies that the Ensembl Rest API advertises as chromosomes. # Mitochondrial DNA has been manually removed. Unfortunately, species IDs from # the Ensembl REST API don't map to dataset names in the BioMart interface. # Because of that, we can't programatically filter chromosome names. # # See get_all_chromosomes() valid_chromosome_names <- c( "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "Z", "1A", "4A", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "X", "25LG1", "25LG2", "LGE22", "Y", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "LG34", "LG35", "2A", "2B", "LG1", "LG2", "LG3", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10", "LG11", "LG12", "LG13", "LG14", "LG15", "LG16", "LG17", "LG18", "LG19", "LG20", "LG21", "LG22", "LG23", "W", "LG24", "LG25", "LG26", "LG27", "LG28", "LG29", "LG30", "LG01", "LG02", "LG03", "LG04", "LG05", "LG06", "LG07", "LG08", "LG09", "A1", "A2", "A3", "B1", "B2", "B3", "B4", "C1", "C2", "D1", "D2", "D3", "D4", "E1", "E2", "E3", "F1", "F2", "LGE64", "LG7_11", "a", "b", "c", "d", "f", "g", "h", "LG28B", "LG30F", "LG36F", "LG37M", "LG42F", "LG44F", "LG45M", "LG48F", "LG49B", "ssa01", "ssa02", "ssa03", "ssa04", "ssa05", "ssa06", "ssa07", "ssa08", "ssa09", "ssa10", "ssa11", "ssa12", "ssa13", "ssa14", "ssa15", "ssa16", "ssa17", "ssa18", "ssa19", "ssa20", "ssa21", "ssa22", "ssa23", "ssa24", "ssa25", "ssa26", "ssa27", "ssa28", "ssa29", "2a", "2b", "7a", "7b", "I", "II", "III", "IV", "V", "VI", "VII", "VIII", "IX", "XI", "XII", "XIII", "XIV", "XV", "XVI", "LGE22C19W28_E50C23", "1a", "22a", "sgr01", "sgr02", "sgr03", "sgr04", "sgr05", "sgr06", "sgr07", "sgr08", "sgr09", "sgr10", "sgr11", "sgr12", "sgr13", "sgr14", "sgr15", "sgr16", "sgr17", "sgr18", "sgr19", "XVII", "XVIII", "XIX", "XX", "XXI", "XXII", "XXIII", "XXIV", "groupI", "groupII", "groupIII", "groupIV", "groupV", "groupVI", "groupVII", "groupVIII", "groupIX", "groupX", "groupXI", "groupXII", "groupXIII", "groupXIV", "groupXV", "groupXVI", "groupXVII", "groupXVIII", "groupXIX", "groupXX", "groupXXI", "2L", "2R", "3L", "3R", "MIC_1", "MIC_10", "MIC_11", "MIC_2", "MIC_3", "MIC_4", "MIC_5", "MIC_6", "MIC_7", "MIC_8", "MIC_9", "X1", "X2", "X3", "X4", "X5" ) #' Get all chromosome names for an Ensembl dataset. #' #' The function tries to filter out valid chromosome names from the available #' assemblies in the dataset. get_chromosome_names <- function(dataset) { chromosome_names <- biomaRt::listFilterOptions(dataset, "chromosome_name") chromosome_names[chromosome_names %chin% valid_chromosome_names] } # Retrieve information on human genes. This will only include genes on # assembled chromosomes. Chromosomes are filtered using get_chromosome_names(). rlog::log_info("Retrieving information on human genes") dataset <- biomaRt::useDataset("hsapiens_gene_ensembl", mart = ensembl) human_data <- data.table(biomaRt::getBM( attributes = c( "ensembl_gene_id", "hgnc_symbol", "chromosome_name", "start_position", "end_position" ), filters = "chromosome_name", values = get_chromosome_names(dataset), mart = dataset )) # Remove duplicated gene IDs (at the time of writing, there are a handful). human_data <- unique(human_data, by = "ensembl_gene_id") # Only keep relevant information on genes. genes <- human_data[, .( id = ensembl_gene_id, name = hgnc_symbol, chromosome = chromosome_name )] # Retrieve gene distance data across species. rlog::log_info("Retrieving distance data") distances <- data.table() #' Handle data for one species. handle_species <- function(species_id, species_data) { chromosomes <- species_data[, .(chromosome_length = max(end_position)), by = chromosome_name ] # Store the number of chromosomes in the species table. species[id == species_id, n_chromosomes := nrow(chromosomes)] # Store the median chromosome length in the species table. species[ id == species_id, median_chromosome_length := chromosomes[, median(chromosome_length)] ] # Precompute the genes' distance to the nearest telomere. species_distances <- species_data[ chromosomes, .( species = species_id, gene = ensembl_gene_id, distance = pmin( start_position, chromosome_length - end_position ) ), on = "chromosome_name" ] # Add species distances to the distances table. distances <<- rbindlist(list(distances, species_distances)) } # Handle the human first, as we already retrieved the data and don't need to # filter based on orthologies. handle_species("hsapiens", human_data) # Iterate through all other species and retrieve their distance data. for (species_id in species[86:nrow(species), id]) { rlog::log_info(sprintf("Loading species \"%s\"", species_id)) dataset <- biomaRt::useDataset( sprintf("%s_gene_ensembl", species_id), mart = ensembl ) # Besides the attributes that are always present, we need to check for # human orthologs. Some species don't have that information and will be # skipped. if (!"hsapiens_homolog_ensembl_gene" %chin% biomaRt::listAttributes(dataset, what = "name")) { rlog::log_info("No data on human orthologs") species <- species[id != species_id] next } chromosome_names <- get_chromosome_names(dataset) # Skip the species, if there are no assembled chromosomes. if (length(chromosome_names) <= 0) { rlog::log_info("No matching chromosome assemblies") species <- species[id != species_id] next } # Retrieve information on all genes of the current species, that have # human orthologs. This is called "homolog" in the Ensembl schema. species_distances <- data.table(biomaRt::getBM( attributes = c( "hsapiens_homolog_ensembl_gene", "chromosome_name", "start_position", "end_position" ), filters = c("with_hsapiens_homolog", "chromosome_name"), values = list(TRUE, chromosome_names), mart = dataset )) # Only include human genes that we have information on. species_distances <- species_distances[ hsapiens_homolog_ensembl_gene %chin% genes$id ] # Only include one ortholog per human gene. species_distances <- unique( species_distances, by = "hsapiens_homolog_ensembl_gene" ) # Rename gene ID column to match the human data. setnames( species_distances, "hsapiens_homolog_ensembl_gene", "ensembl_gene_id" ) handle_species(species_id, species_distances) } # Save data in the appropriate place. usethis::use_data(species, overwrite = TRUE) usethis::use_data(genes, overwrite = TRUE) usethis::use_data(distances, overwrite = TRUE)