avidaR: an R library to perform complex queries on an ontology-based database of digital organisms

avidaDB

We have built a comprehensive database that encompasses the genomes, phenotypes, and transcriptomes of over one million digital organisms resulting from evolution experiments carried out using Avida. A significant portion of the data, more than half, was obtained by randomly sampling the sequence space of genomes that were 100 instructions long (Fortuna et al., 2017). These genomes encoded the same phenotype within a single environment. Furthermore, in a subsequent study (Fortuna, 2022), the phenotype encoded by the genomes of these 512,000 organisms, as well as their transcriptomes, were computed in a set of 1000 distinct environments. The database includes the following key data from these two experiments:

• genome instruction sequence: a linear string of letters representing the instruction codes that make up the genome of a digital organism.

• phenotype: a unique combination of logic operations that a digital organism performs on 32-bit one- and two-binary input numbers.

• transcriptome instruction sequence: a linear string of letters representing the instruction codes that make up the transcriptome executed by a digital organism.

In a third study (Rico and Fortuna, unpublished), we extended the analysis by computing the phenotype and transcriptome of all single-point mutants for 200 digital organisms from the previous dataset. These single-point mutants were generated by introducing a single mutation at different positions in the genome. The computations were performed in the same 1,000 environments as in the previous study. As a result, we obtained an additional dataset comprising 500,000 digital organisms (2500 single-point mutants for each of the 200 organisms). In this expanded dataset, we have stored the same type of data as in the previous studies, including genomes, phenotypes, and transcriptomes. Moreover, we have also included information about the specific genomic locations of all single-point mutations:

• mutant of: a relation linking two digital organisms if their genomes differ in a single instruction code.

• position at genome: the position of an instruction code along the digital genome of a digital organism.

In addition to the genomic, phenotypic, and transcriptomic data, we have also computed and stored information on the viability and generation time of the 1,012,000 digital organisms in the same 1,000 environments:

• viable: the ability of a digital organism to produce an offspring able to replicate by executing its genome.

• generation time: number of instruction codes executed by a digital organism to produce a viable offspring (i.e., the inverse of fitness in the absence of any selective pressure other than reducing the number of instructions that must be executed for successful replication).

Tracing the provenance and reproducibility of data

Clarifying the provenance of data is crucial for traceability and reproducibility in scientific research. We have added functions into avidaR to trace the context and provenance of the data obtained from Avida experiments:

• get_experiment_id_from_organism_id() # establishing a connection between a digital organism and the experiment it was derived from enables researchers to easily retrieve the relevant experiment information associated with that specific organism.

• get_doi_from_experiment_id() # establishing a connection between an experiment and the identifier of the scientific publication it was reported in enables researchers to easily access the publication, which provides additional details on the experimental design.

• get_docker_image_from_experiment_id() # establishing a connection between and experiment and a docker image that contains the configuration files and options used in the experiment enables researchers to reproduce the experimental conditions and obtain repeatable results.

Providing semantics to avidaDB

avidaDB has been implemented as an RDF store (also known as graph database or triple store) to express relationships between resources using subject–predicate–object triples. The most relevant advantage of RDF stores is that they represent, store and query data as a graph. The second characteristic is that they are semantic, which means that they can store not only data but also explicit descriptions of the meaning of that data. These explicit descriptions are known by the RDF and linked data community as ontologies. An ontology is a machine-readable description of a domain that typically includes a vocabulary of terms and some specification of how these terms relate to each other, imposing a structure on the data for such domain. This is also known as a schema.

We have developed the Ontology for Avida (OntoAvida) to provide semantics to avidaDB. The semantic relationships between the terms commonly used in Avida are expressed in the W3C standard ontology language OWL-DL. OntoAvida is part of the Open Biological and Biomedical Ontologies (https://obofoundry.org/ontology/ontoavida.html) and is already available to the scientific community. We have also customized pyLODE to obtain a readable version of OntoAvida, which is also available at https://owl.fortunalab.org/ontoavida.

Connecting to avidaDB

Semantic databases are databases that store RDF data and allow the querying of RDF data via the SPARQL query language (i.e., a programming language used to retrieve and manipulate data stored as subject–predicate–object triples). The library avidaR can connect to triple-stores that support the RDF4J server REST API such as GraphDB. The SPARQL endpoint of avidaDB (https://graphdb.fortunalab.org/) allows both basic connection (requiring no password or requiring basic HTTP user-password authentication) or connection secured with an API access token. The way avidaR connects avidaDB requires user-password authentication (i.e.,  public_avida as both username and password).

avidaR retrieves data from avidaDB without specific knowledge of SPARQL or RDF. The functions in avidaR can be grouped by the type of data requested: (i) genomes as sequences of letters representing code instructions (e.g., get genomes encoding specific phenotypes in particular environments), (ii) phenotypes as a combination of the logic operations computed (e.g., get the phenotype of organisms that computed specific logic operations), (iii) transcriptomes executed during the replication cycle as both sequences of letters representing code instructions and chord diagrams (e.g., get the transcriptomes of organisms whose genomes encoded specific phenotypes in specific environments), and (iv) single-point mutants at specific positions on the genome (e.g., get the genome of the single-point mutant that substituted the instruction coded by letter i by j at a particular position on the genome of a wild-type organism). Users can build programming workflows by piping multiple functions together so that the output of one becomes the input to another (e.g., get a genome encoding different phenotypes in distinct environments (i.e., plastic organism), and retrieve the generation time as a surrogate of fitness in the absence of any selective pressure other than reducing the number of instructions executed).

avidaDB access (required to establish a connection to the database)

• avidaDB <- triplestore_access$new() # create triplestore object.

• avidaDB$set_access_options(url = “https://graphdb.fortunalab.org”, user = “public_avida”, password = “public_avida”, repository = “avidaDB”) # set access options.

avidaDB summary

• get_db_summary(triplestore = avidaDB) # get a summary of the content of the database.

Genetic language used by the organisms in Avida

• instruction_set(inst_set = “heads”) # get the default Avida genetic language.

Logic operations computed by digital organisms

• logic_operation() # get the list of logic operations that a digital organism can compute (i.e., the default logic-9 environment).

• get_logic_operation_from_phenotype_id() # get the list of logic operations computed by a digital organism whose genome encodes a specific phenotype out of the 512 distinct phenotypes comprising the entire phenotype space (i.e., ranging from phenotype 0 to phenotype 511).

Genomes harbored by digital organisms

• get_genome_id_from_logic_operation() # get genome identifiers of organisms that computed specific logic operations.

• get_genome_id_from_phenotype_id() # get genome identifiers encoding specific phenotypes in particular environments.

• get_genome_id_from_transcriptome_id() # get genome identifiers of organisms that executed specific transcriptomes in particular environments.

• get_genome_id_from_genome_seq() # get genome identifiers of genome sequences.

• get_genome_seq_from_genome_id() # get genome sequence from genome identifiers.

Phenotypes encoded by genomes of digital organisms

• get_phenotype_id_from_genome_id() # get the phenotype identifier encoded by a genome identifier.

• get_phenotype_id_from_genome_seq() # get the phenotype identifier encoded by a genome sequence.

• get_phenotype_id_from_logic_operation() # get the phenotype identifier of organisms that computed specific logic operations.

• get_phenotype_id_from_transcriptome_id() # get the phenotype identifier of organisms that executed specific transcriptomes in specific environments.

Transcriptomes executed by digital organisms

• get_transcriptome_id_from_genome_id() # get the transcriptome identifier encoded by a genome identifier.

• get_transcriptome_id_from_genome_seq() # get the transcriptome identifier encoded by a genome sequence.

• get_transcriptome_id_from_logic_operation() # get the transcriptome identifier of organisms that computed specific logic operations.

• get_transcriptome_id_from_phenotype_id() # get the transcriptome identifier of organisms whose genomes encoded specific phenotypes in specific environments.

• get_transcriptome_seq_from_transcriptome_id() # get the transcriptome sequence from the transcriptome identifier.

• plot_transcriptome(inst_set = “heads”) = # plot a transcriptome as a chord diagram using the default instruction set.

Single-point mutants

• get_genome_id_of_wild_type_organisms() # get the genome identifier of the list of organisms for which genomic, transcriptomic, and phenotypic information of all their single-point mutants is available (i.e., wild-type organisms).

• get_mutant_at_pos() # get the genome sequence of specific single-point mutants at a particular position on the genome of a wild-type organism.

The fitness cost of phenotypic plasticity

In the first case study, we investigate the impact of phenotypic plasticity on the fitness of digital organisms. Phenotypic plasticity is not only a pervasive feature of biological organisms but also of digital ones (Fortuna, 2022). Specifically, the environment can modify the sequence of instructions executed from a digital organism’s genome (i.e., its transcriptome), which results in changes in its phenotype (i.e., the ability of the digital organism to perform Boolean logic operations). This epigenetic pathway for plasticity comes at a fitness cost to an organism’s viability and replication time: the longer the replication time (higher fitness cost), the more chances for the environment to modify the genetic execution flow control, and the higher the likelihood for the genome to encode novel phenotypes.

We obtain a set of genomes that exhibit phenotypic plasticity, where organisms can switch between different phenotypes in response to changes in their environment. We also collect a set of genomes that do not exhibit plasticity and maintain a fixed phenotype across all environments. By comparing the replication times of organisms with and without plasticity in different environments, we quantify the fitness cost associated with phenotypic plasticity. Here is the code to retrieve the data and perform the analysis.

 
 
# load required libraries: 
library(avidaR) 
library(tidyverse) 
# data reproducibility: 
set.seed(1000) 
# create triplestore object: 
avidaDB <- triplestore_access$new() 
# set access options: 
avidaDB$set_access_options( 
url = "https://graphdb.fortunalab.org", 
user = "public_avida", 
password = "public_avida", 
repository = "avidaDB" 
) 
# define function to get the data: 
data <- function(boolean_func, seeds, plastic) { 
    get_genome_id_from_logic_operation( 
        # logic operations that the target organisms must perform: 
        logic_operation = boolean_func, 
        # number of distinct environments where that phenotype will be computed: 
        seed_id = 1:seeds, 
        triplestore = avidaDB) %>% 
    # show data for each genome that computes those logic operations: 
    group_by(genome_id) %>% 
    # number of environments where the phenotype was encoded: 
    summarize(n_seeds = n()) %>% 
    # classify the genome as plastic or no plastic: 
    filter(if (plastic == TRUE) n_seeds == seeds else n_seeds == 1) %>% 
    # randomly select 10 organisms from those that meet the requirements: 
    sample_n(10) %>% 
    select(genome_id) %>% 
    # get the transcriptome from the genome of each selected organisms: 
    get_transcriptome_id_from_genome_id( 
        genome_id = as.integer(gsub("genome_", "", .$genome_id)), 
        transcriptome_seq = TRUE, 
        seed_id = 1:seeds, 
        triplestore = avidaDB) %>% 
    select(genome_id, transcriptome_id, transcriptome_seq) %>% 
    distinct %>% 
    # get the number of instructions that the organism executed to replicate as a surrogate for fitness: 
    mutate(generation_time = nchar(transcriptome_seq), 
        phen_group = boolean_func, plastic_group = plastic) %>% 
  
 
 
                select(-transcriptome_seq) 
            } 
            # get the data and join them in a single dataframe: 
            # select 10 non-plastic organisms performing the logic function NOT: 
            data_NOT_non_plastic <- 
                data(boolean_func = "not", seeds = 10, plastic = FALSE) 
            # select 10 plastic organisms performing the logic function NOT: 
            data_NOT_plastic <- 
                    data(boolean_func = "not", seeds = 10, plastic = TRUE) 
            # select 10 non-plastic organisms performing the logic function EQUALS: 
            data_EQU_non_plastic <- 
                    data(boolean_func = "equals", seeds = 10, plastic = FALSE) 
            # select 10 plastic organisms performing the logic function EQUALS: 
            data_EQU_plastic <- 
                    data(boolean_func = "equals", seeds = 10, plastic = TRUE) 
            # provide data as a single data frame: 
            df <- rbind( 
                data_NOT_non_plastic, data_NOT_plastic, 
                data_EQU_non_plastic, data_EQU_plastic) 
            # plot the results: 
            ggplot(df, aes(x = phen_group, y = generation_time, fill = plastic_group)) + 
            geom_boxplot() 
    

This preliminary finding suggests that plasticity may incur a fitness cost for simple phenotypes, but not for complex phenotypes. This is evidenced by an increased time for digital organisms to replicate, indicating a decrease in fitness (Fig. 1).

The fitness effects of single-point mutations

In the second case study, we retrieve data from avidaDB to identify the single-point mutation that has the largest positive effect on replication rate. We consider the replication rate as a surrogate for fitness, where a higher replication rate corresponds to higher fitness. Replication in digital organisms involves the copying of the genome instruction by instruction into a new memory region, a process that can occasionally result in errors known as mutations. Mutations occur when an instruction is inaccurately copied, leading to its replacement in the offspring genome with a randomly selected instruction. Mutations can have different effects on fitness, including being lethal (resulting in non-viable offspring), neutral (no change in phenotype or fitness), deleterious (decreasing fitness), or beneficial (increasing fitness).

The process of identifying the most significant single-point mutation on replication rate involves multiple steps:

• randomly sampling the genome sequence of a single organism (wild-type organism).

• retrieving the transcriptome of the wild-type organism and its single-point mutants in a specific environment.

• selecting the mutant with the highest fitness based on generation time as a surrogate measure.

• comparing the genome sequences of the wild-type organism and the fittest mutant to identify the specific mutated instruction and its position in the genome.

• finally, plotting the transcriptomes of the wild-type organism and the fittest mutant for comparison (Fig. 2).

Here is the code to retrieve the data and perform the analysis.

 
 
# load required libraries: 
library(avidaR) 
library(tidyverse) 
# data reproducibility: 
set.seed(10) 
# create triplestore object: 
avidaDB <- triplestore_access$new() 
# set access options: 
avidaDB$set_access_options( 
url = "https://graphdb.fortunalab.org", 
user = "public_avida", 
password = "public_avida", 
repository = "avidaDB" 
) 
# get data: 
data_wild_type <- 
    get_genome_id_of_wild_type_organisms( 
        triplestore = avidaDB) %>% 
    # randomly select the genome of an organism for which all of its single-point mutants are known: 
    sample_n(1) %>% 
    # get the transcriptome from the genome of the selected organism: 
    get_transcriptome_id_from_genome_id( 
        genome_id = as.integer(gsub("genome_", "", .$genome_id_wild_type)), 
        transcriptome_seq = TRUE, 
        seed_id = 1, 
        triplestore = avidaDB) %>% 
    # get the number of instructions that the organism executed to replicate as a surrogate for fitness: 
    mutate(generation_time = nchar(transcriptome_seq)) %>% 
    select(genome_id, transcriptome_id, generation_time) %>% 
    arrange(generation_time) 
# get all single-point mutants from the selected organism: 
data_mutant <- 
    # specificy the position of the mutation along the genome of the organism: 
    get_mutant_at_pos( 
        genome_id = as.integer( 
            gsub("genome_", "", data_wild_type$genome_id)), 
        triplestore = avidaDB) %>% 
    select(genome_id_mutant) %>% 
    # get the transcriptome from the genome of the mutant organisms: 
    get_transcriptome_id_from_genome_id( 
        genome_id = as.integer(gsub("genome_", "", .$genome_id_mutant)), 
        transcriptome_seq = TRUE, 
  
 
 
                    seed_id = 1, 
                    triplestore = avidaDB) %>% 
                select(genome_id, transcriptome_id, transcriptome_seq) %>% 
                distinct %>% 
                # get the number of instructions that the mutant executed to replicate as a surrogate for fitness: 
                mutate(generation_time = nchar(transcriptome_seq)) %>% 
                select(-transcriptome_seq) %>% 
                # sort the mutants by the number of instructions that the mutant executed to replicate: 
                arrange(generation_time) %>% 
                # select a single mutant organism, i.e., the one with the highest fitness: 
                head(1) 
            # join data for both organisms in a single dataframe: 
            df <- inner_join( 
                rbind(data_wild_type, data_mutant), 
                get_genome_seq_from_genome_id( 
                    genome_id = c(as.integer( 
                        gsub("genome_", "", data_wild_type$genome_id)), 
                                  as.integer( 
                        gsub("genome_", "", data_mutant$genome_id))), 
                    triplestore = avidaDB), 
                by = "genome_id") 
            # get the relative fitness of the mutant compared to that of the wildtype: 
            df$generation_time[1] / df$generation_time[2] 
            # get the position of the mutation in the genome of the mutant: 
            pos <- mapply(function(x, y) which(x != y)[1], 
                strsplit(df$genome_seq[1], ""), strsplit(df$genome_seq[2], "")) 
            # get the mutated instruction: 
            substr(df$genome_seq[1], pos, pos) 
            substr(df$genome_seq[2], pos, pos) 
            # plot the transcriptomes of both organisms: 
            plot_transcriptome( 
                inst_set = "heads", 
                transcriptome_id = as.integer( 
                    gsub("transcriptome_", "", df$transcriptome_id[1])), 
                seed_id = 1, 
                save = TRUE, save_path = getwd(), format = "pdf", 
                triplestore = avidaDB) 
    

A single-point mutation has the potential to significantly impact fitness by altering the execution flow of the instructions within a digital organism’s genome. In our analysis, we observed a notable fitness improvement resulting from a specific mutation. Specifically, the substitution of the math instruction n (dec) at position 54 in the genome with the control-flow instruction h (jump-head) led to a 29% increase in the replication rate of the digital organism (Fig. 2). This finding highlights the important influence that individual mutations can have on the performance and fitness of digital organisms.