Building a community-driven bioinformatics platform to facilitate Cannabis sativa multi-omics research
Southern Cross University, Military Road, Lismore New South Wales 2480, Australia
International Rice Research Institute, Pili Drive, Los Baños Laguna 4031, Philippines
*ramil.mauleon@scu.edu.auAbstract
Global changes in Cannabis legislation after decades of stringent regulation, and heightened demand for its industrial and medicinal applications have spurred recent genetic and genomics research. An international research community emerged and identified the need for a web portal to host Cannabis-specific datasets that seamlessly integrates multiple data sources and serves omics-type analyses, fostering information sharing.
The Tripal platform was used to host public genome assemblies, gene annotations, QTL and genetic maps, gene and protein expression, metabolic profile and their sample attributes. SNPs were called using public resequencing datasets on three genomes. Additional applications, such as SNP-Seek and MapManJS, were embedded into Tripal. A multi-omics data integration web-service API, developed on top of existing Tripal modules, returns generic tables of sample, property, and values. Use-cases demonstrate the API’s utility for various -omics analyses, enabling researchers to perform multi- omics analyses efficiently.
Article notes
Competing Interest Statement
The authors have declared no competing interest.
Footnote Group
STATEMENT OF NEED
Background on Cannabis
Cannabis sativa (cannabis) is one of the earliest domesticated crops for food, fiber and medicine [1]. As hemp, Cannabis was a significant source of fiber in ropes traditionally used for shelter, clothing, traps, harnesses, sails and halyards. Its seeds (hempseed) are relatively rich in nutritious oils of high omega-3 and omega-6 fatty acid content [2], proteins, and contain a variety of micronutrients for human diet [3]. Its female flowers produce a mixture of secondary metabolites, some of which have psychotropic effects and have been used in religious rituals and recreational use. Because of these properties and uses, cannabis was also regarded as an ancient ‘wonder plant’ coined with various polarizing names including ‘devil’s weed’, ‘penicillin of Ayurvedic medicine’, ‘donator of joy’ [4]. The secondary metabolites found in high abundance in glandular trichomes on the surface of female flowers for pathogen resistance [5], insect deterrent [6], and UV protection [7] also have pharmacological and nutritional potential [8]. Cannabis glandular trichomes produce large amount of terpenoids, and cannabinoids. Calvi et al. [9] identified 554 compounds in cannabis including 113 cannabinoids and 120 terpenes. More than 200 metabolites from different genotypes are volatile, with 58 monoterpenes and 38 sesquiterpenes being reported. The most abundant monoterpenes are limonene, b-myrcene, a-pinene and linalool, while sesquiterpenes are E-caryophyllene, b-farnesene, and b-caryophyllene.
Cannabis is dioecious and diploid with 9 autosomes and a pair of sex chromosomes (female XX, male XY). The diploid genome size estimated by flow cytometry is 1636Mb for female (818Mb haploid size) and 1683Mb (haploid size 843Mb) for male plants [10]. The first draft sequence reported a haploid size of 534 Mb [11], while the accepted standard reference genome to date is 876Mb [12] Cannabis had been classified as narcotic for most of the 20th century (Mead, 2014), but recent revisions of legal status and recognition of defined medicinal benefits has supported it as an emerging economic crop with the World Health Organization scheduling recommendation [13]. Improving cannabis for the food, fiber and medicinal traits has thus recently been the focus of both public and private sector research (Small, 2017). The global market projection for medicinal cannabis is $176 billion by 2030 [14]. At 21% cumulative growth rate, the global cannabis cultivation market including hemp, could reach USD 1.8T by 2030 [15].
The need for Cannabis multi-omics database and analyses platform
The assembly of reference genomes [12,16–19] for hemp and medicinal types has facilitated an abundance of publications and associated publicly available genomic, transcriptomic, proteomic, and metabolomic data sets. Data are available from general web portals also include cannabis genomic and other omics data (e.g. NCBI, ENA, DDBJ, CoGe)[20–22], plant-specific databases (e.g. EnsemblPlant, Phytozome, PlantCyc)[23–26], and Cannabis-specific databases. Kannapedia [27] is a commercial site from Medicinal Genomics, which provides summaries on the genomic variation of cannabinoid biosynthesis genes for nearly two thousand cultivars. Bulk downloads of most of the data, however, are not available as a built-in feature. CannabisDraftMap.org [28] provides an inventory of raw proteomics and metabolomics data. The site has no analysis tools or results available, and only features short descriptions of various cannabis genomics projects. The Cannabis Genome Browser [29] hosts the results from the first cannabis draft sequence publication [11] but has not been updated since. Leafly [30] is a commercial site with cannabinoid and qualitative description of terpenes. SeedFinder [31] provides an interface for cannabis breeders, seedbanks, growers, and enthusiasts by collecting and standardizing information about cannabis. As of February 2024, the database has 31,823 different cultivars listed. While it is not an online shopping platform, it is supported by related advertisements. CannabisGDB [32] is an integrated cannabis functional genomics database with four modules: varieties, gene loci, metabolites, and proteins. ’Varieties’ provides description for 8 assembled genomes available from NCBI. ’Gene loci’ contains their in-house gene predictions and functional annotation on these 8 genomes, and these are displayed using JBrowse [33]. ‘Proteins’ includes the UniProt plant proteins matched to the in-house protein model. ’Metabolite’ displays various bar charts showing percent cannabinoid content, as collected from external publications.
Most of the Cannabis-specific websites and databases described above host proprietary data; hence the focus is mainly on specific cultivars and trait characterization. Some have omics-type data but are usually from unpublished or proprietary sources. These limitations provided the motivation for this project to build an authoritative ‘open-science’ reference bioinformatics portal to support cannabis multi-omics research with capabilities that we believe will be of great benefit to the research community, including: (a) tools to perform quick analyses using the loaded datasets, (b) tools allowing data/contents curation by the community (c) creation and promotion on the use of research data standards, and (d) a platform for collaboration and communication by the cannabis research community.
The features and capabilities were initially intended to fulfill the needs defined by the International Cannabis Genomics Research Consortium (ICGRC). Formed during the 2020 Plant and Animal Genome conference by an international community of cannabis researchers from the public and private/commercial sectors, ICGRC aims to provide an authoritative source of rigorous and peer- reviewed scientific information relating to Cannabis genetics, genomics and phenotyping. The scope of ICGRC activities encompass (a) Information and data sharing, (b) establishing standards and reference information, (c) facilitating the unambiguous description and availability of genetic resources, (d) contributions to consistent annotation of gene identity and function, and (e) describing common genetic resources and where necessary, virtual gene-banks [34].
IMPLEMENTATION
Tripal chosen as technology for the Cannabis omics Platform
There is a range of widely-adopted systems in the public domain suitable for genomic and genetic web databases/portals, such as InterMine [35], MetaCyc [36], EnsemblPlant [24], and Tripal [37]. Tripal is a content management system for genomics data based on the Drupal platform (drupal.org). The data architecture is based on Chado [38], a highly normalized schema that uses biological ontologies to define bioentities and their relationships. Tripal architecture is modular, allowing plug-in features such as ready-to-use specialized ‘omics modules for loading, visualization, and analysis, built-in support for user-driven content management, and is GMOD [39] standards-compliant. Tripal is used in many plant genomic research communities [37,40], and various genome hub projects [41] including several from the SouthGreen bioinformatics platform [42]. Its major advantage is flexibility for community-wide contribution and curation of content, thus Tripal was chosen as the software platform for the project. Tripal setup for the project is described in protocols.io, referred to as protocol from onwards (reserved DOI: dx.doi.org/10.17504/protocols.io.n2bvj3nz5lk5/v3)[43], of which the major steps are: Cannabis dataset download, reanalysis, and formatting for Tripal loading (protocol steps 1-8), Tripal docker container instantiation (protocol step 9) and loading of the reference genomes and annotation (protocol step 10). Modules have then been installed, enabled and loaded with the corresponding data: gene_search (protocol step 11), expression (protocol step 12), phenotype (protocol step 13), map (protocol step 14), and finally synteny (protocol step 15).
Cannabis –omics, genetic, and phenotypic datasets in the public domain
Published-omics datasets are sourced from databases such as NCBI Sequence Read Archive (SRA) (RRID:SCR_004891) (raw sequences)[44], NCBI RefSeq (RRID:SCR_003496) (genome assemblies and annotation)[45], Gene Expression Omnibus (GEO) (RRID:SCR_005012) (gene expression omnibus)[46], UniProt (RRID:SCR_002380) (protein sequences and annotations)[47], OBO Foundry biological ontologies and controlled vocabulary OBO (RRID:SCR_007083) [48], Crop Ontology Crop Ontology (RRID:SCR_010299) [49] crop-specific controlled vocabulary,. Datasets have also been manually curated from the published article tables, figures, and supplementary data and sites, which are common sources, especially for non-model or less studied species like cannabis.
Genome assemblies and annotations
The first peer-reviewed publication for cannabis genomes were for Purple Kush (PK), a medicinal type (THC-dominant) marijuana cultivar, and Finola (FN), a hemp-type cultivar developed for oil seed [11]. Both PK and Finola assemblies were updated and improved in Laverty et al. [16]. CBDRx (cs10) is a CBD-type cultivar derived from a cross from Skunk#1 (marijuana) with Carmen (hemp) and is estimated to be 89% marijuana and 11% hemp by admixture analysis. This was the first Cannabis assembly with NCBI RefSeq gene annotations (cs10, GCF_900626175.2) for three years and used by many publications for gene assignment. Pink Pepper (ASM2916894v1) is the second assembly annotated by RefSeq just recently (November 2023).
Three primary genome assemblies are currently loaded to the portal: cs10 (aka CBDRx, a CBD- dominant line) [12], PK and FN [16]. They were chosen to represent the major cannabis end-uses and have transcript sequences for gene annotation. Other Cannabis genome assemblies currently in the public domain, such as cultivars First Light [12], Jamaican Lions [19], Cannbio-2 [18], JL wild [17], and Pink Pepper were not included in this release of the portal.
The annotations from RefSeq for cs10, and our FINDER [50] annotation for PK and FN, were loaded and searchable using the ‘Gene Function Search’ module interface and displayed in the JBrowse browser (RRID:SCR_001004). The query constraints are cultivar/assembly, annotation analysis, genome region, gene or transcription accessions, and keyword in the functional annotations. Query by sequence is also possible using BLAST (RRID:SCR_004870) interface. A summary of annotations available in the portal is in Table 1.
Transcriptomes, comparative expression, and proteomics datasets
The public availability of transcript data has been the primary consideration for our choice of genomes to load in the portal (Table 2). For comparative analyses either absolute expression levels across multiple samples or relative fold-change of transcript levels between two samples is available. Currently three (3) Cannabis transcript assemblies are publicly available: canSat3, finola1 and cannbio2. All were used to annotate their respective genome, and in differential expression studies across tissues [11,18,51]. These transcript assemblies are viewable in the JBrowse browser. Functional annotation was assigned for each transcript according to homology with RefSeq mRNA sequences, and the genome location by splice-aware alignment with the different genomes. As with gene loci, the transcript assemblies can be queried in the ‘Gene Function Search’ page, and their locations and transcript levels are also available to be displayed in JBrowse.
Comparative expression studies for Cannabis are limited, mostly on changes related to secondary metabolism and fiber synthesis, comparing different cultivars, treatment conditions, tissues or developmental stages. A gene expression study for bast fiber development in the textile hemp variety Santhica 27 included tissues sampled from the bottom, mid and top stem region [57] carried out. In another study, Finola trichomes were separated from mid-stage flowers and categorized into bulbous, sessile, and mature stalks [52]. Terpene and cannabinoid profiles were measured and correlation with synthase gene transcript levels assessed [53,58]. Relevant expression datasets in NCBI GEO were also loaded. The list of absolute and differential expression datasets available within the ICGRC Tripal resource is shown in Table 3 and Table 4 respectively. Expression data are displayed as transcript heatmaps of genes by samples. It is important for this module to have traceable assignment of sample attributes from the source. Expression data are also displayed overlaid on the generic Biopathways module.
As of completion of this project, only one Cannabis proteome expression dataset has been available in GEO. Using quantitative time of flight mass spectrometry, 1240 proteins from glandular trichomes heads, 396 from trichome stalks, and 1682 from whole flower tissue isolates were identified by Conneely et al. [61].
Variants discovery using resequencing datasets
Details on the reanalysis of Cannabis Next Generation Sequencing NGS datasets are described in an accompanying paper (DOI: 10.25918/preprint.367) about the CannSeek SNPs database [62]. Data available in NCBI SRA as of December 2022 consisting of 26 trichome RNA-Seq, 383 genomic DNA, and the Phylos amplicons were reanalyzed using the GATK (RRID:SCR_001876) [63] and Parabricks [64] variant calling pipelines. The results are displayed using the CannSeek software, derived from Rice SNP-Seek [65], hosted in a Podman, or Docker, container [66,67] and embedded as a Tripal page.
Maps
The interactive Map Viewer, a part of the Tripal platform [68], can display the entire genome or chromosome, and features such as markers and QTLs in a single view. Maps positions may either be physical (use base position - bp), or genetic (use centi-Morgan - cM). Common markers are linked between the genetic and physical maps when displayed side by side. Public Cannabis genetic linkage data are scarce; only two recent publications of QTL maps are currently loaded; a study on gene duplication and divergence affecting drug content (Weiblen et al 2015) and QTL mapping for agronomic and biochemical traits (Woods et al. 2021).
Phenotypes
The Phenotype Viewer displays the distribution of measured traits including metabolites (cannabinoid and terpenoid content and their relative composition) across samples and experiments, visualized in a ‘violin plot’ [69]. Table 5 lists the phenotypic datasets loaded. This feature requires three genus- specific controlled vocabularies for trait, method, and unit. During setup, they can be defined from existing ontologies or community customized. New terms may be added as new datasets are loaded. Currently these ontologies use ChEBI [73] for chemical traits, PECO [74] for method and UO [75] for units. In the absence of a cannabis-specific vocabulary for methods and units in Crop Ontology, publication-specific terms were created when no suitable match in the existing ontologies were found.
Preparing Datasets for the Tripal modules
Gene model prediction and annotations
All resequencing data were trimmed using Trimmomatic (RRID:SCR_011848) [76] (protocol step 1). Cs10 genes and annotations were downloaded from NCBI RefSeq. (protocol step 2) The FINDER [50] pipeline was used to predict gene models for updated assemblies of Purple Kush and Finola using cultivar-specific RNA-Seq data. GffRead (RRID:SCR_018965) [77] was used to cluster overlapping gene loci predictions into a common locus. (protocol step 3) Proteins were annotated for putative function/structure/location by homology search against the UniProt [47] databases using mmseqs2 [78]. Gene Ontology (RRID:SCR_002811) and InterPro (RRID:SCR_006695) IPR protein domains annotation of these protein were done using InterProScan (RRID:SCR_005829) [79]
Gene expression
Expression data from GEO [46] were loaded using the loader for expression and biosample. To build a system that hosts and integrates multiple -omics datasets for the same sample sets, secondary analyses were performed on public data to complete these for all -omics data types. Metabolite measurements were loaded and RNA-Seq sequences were reanalyzed. Gene expression using cs10 mRNAs were quantified on 21 trichome-specific public RNA-Seq datasets [51–54] using Salmon (RRID:SCR_017036) (protocol step 5). Transcripts from the 21-trichomes RNA-Seq samples were assembled following (protocol step 4) using rnaSPAdes (RRID:SCR_016992) [80], then mapped to the three references using GMAP (RRID:SCR_008992) [81] and ORFs were identified using TransDecoder (RRID:SCR_017647) [82]. Quantification of transcript abundance was then performed using Salmon following (protocol step 5).
Synteny block discovery
This follows the Tripal synteny viewer instructions of generating synteny blocks using gene annotations (protocol step 8). BLASTN (RRID:SCR_001598) and MCScanX (RRID:SCR_022067) [83] were used on the cs10 gene annotation and our generated predictions for Purple Kush and Finola.
BLAST database
BLAST databases for the 3 genomes and the annotated genes are installed in the portal, allowing genes and protein sequences from the users to be searched against genome assembly, coding sequence and proteins for cs10, Purple Kush PK, and Finola FN. Cs10 gene model sequences are from NCBI Refseq, while PK and FN are from the FINDER gene prediction pipelines.
Data preparation for external, non-Tripal tools
Some data are hosted and visualized using modules that are not part of Tripal (e.g. genome browser, pathways visualization, genotype management). Although Tripal might have some equivalent module, or the module might be implemented differently in other Tripal sites, the following tools were used for user benefit.
JBrowse Genome Browser
The JBrowse Genome Browser [33] displays genome features along the reference chromosome or scaffolds. The data is loaded from interval files like GFF or BED. It is highly interactive and web-based making it popular and used by many genomic database sites. Genes and transcript alignments for cs10, PK and FN genomes are available.
Biopathways mapping
MapManJS, a web version of MapMan (RRID:SCR_003543) [84], is a gene expression and biological pathway visualization tool. It uses curated, organism-specific sets of genes for each pathway. (protocol step 7) Since MapMan does not have Cannabis-specific pathways, we lifted over (via best match with Cannabis cs10 genes using mmseq2 reciprocal best hit on proteins) genes and pathways from other plant maps, specifically A. thaliana, tomato and eucalyptus. The expression of the cs10 genes is visualized over the selected pathway mapping and study loaded in the Expression module.
SNPs discovery and genotype data management
Variant data of the 21 trichomes RNA seq were discovered using the GATK-RNASeq pipeline (protocol step 6). Genomic resequencing data for 398 samples were also used to discover SNP variants using the GATK germline SNP calling workflow [85] and Parabricks [64].
Genotype data results from variant discovery pipelines are presented as large genotype matrix of samples and base position. The Tripal Natural Diversity Genotypes (DOI:10.5281/zenodo.3731337) module used in KnowPulse [69] is able to display this matrix. However, it uses web technologies that limit its interactivity and has performance issues on larger genotype matrix queries. The Rice SNP-Seek database [65] can display thousands of SNP positions and samples while maintaining interactivity and allows more analysis on the result and is used in this project.
RESULTS
Setup and customizations of the Tripal platform, and re-engineering modules for multi-omics data integration
The detailed workflow to set up and customize the platform is enumerated in the protocol [86], and illustrated in Figure 1. Documentations about instantiating various Tripal modules and loading data are notoriously scattered across publications, and this workflow aims to collate and summarize the steps required to format files to simplify preparation of the various modules. The workflow generates the files necessary to initially populate and test the different Tripal modules. This can also guide installations of future Tripal sites, especially for plants with scarce published data, or where only raw sequences are available.
We added features that are not part of Tripal that were essential to enable multi-omics analysis, including pathway and expression visualization (MapMan JS), and genotype mining and visualization (SNP-Seek). NCBI BLAST was also included in the platform. (protocol step 18) provides details on the installation of these external applications (JBrowse, CannSeek, MapManJS, BLAST). In addition, we developed a novel module specifically for multi-omics integration. This is in the form of a web-service with a predefined Application Programming Interface (API) that queries and merges the different Tripal modules, facilitating rapid retrieval of different data types from multiple sources. The capabilities of this API are demonstrated with use-cases typically performed in -omics studies.
Challenges encountered in the existing Tripal and Chado design and solutions implemented
Although existing Tripal modules can host and display various biological -omics datasets, integrating data from the different modules is not straightforward. The underlying Chado schema is highly flexible and normalized. Unlike the more common Entity-Relation, Chado follows the Entity-Attribute-Value Model database model[87]. The data model relies on biological ontologies, which are loaded as data, to define datatypes and relationships. The different loaders for each module populate database tables, which may be scattered, duplicated and confusing in their utility and terminologies.
An additional complication is that Tripal introduces a layer of abstraction for biological entities on top of the storage tables in the Chado schema. Using ontology terms as type identifiers, a Chado table can store records for different biological entity types. Table 6 lists the default Tripal entities and their corresponding Chado tables and types. The Tripal BioEntity objects and Chado relational table mappings are illustrated in Figure 2. The feature, analysis and stock tables are used in common by multiple Tripal BioEntity types. Data management can be challenging when mapping the relationships between distinct bioentities because they are not equivalent to the entity-relationship defined by the Chado schema. Another challenge is the requirement to manually synchronize the Chado records and the Tripal entities when loaders are not used There are synchronization buttons in Tripal to publish Chado records to Tripal, unpublish deleted records, or delete unpublished bioentities. Using these features requires thorough understanding and experience on the datasets and platform.
Several solutions are implemented to overcome the challenges encountered. The integration process is built on top of existing Tripal modules, specifically for the ‘chado_gene_search’, ‘expression’, and ‘phenotype’ modules. These modules generate materialized views from complex queries used in the visualizations. The -omics integration API queries these materialized views by SQL. Merging from the different modules is implemented using table pivot. All queries return a triplet of property, biomaterial and value. A pivot function then transforms them into a table with the biomaterials as column labels, and properties as row labels. Figure 3 illustrates how the modules are linked for integration. Biomaterial labels can be from ‘dbxref’, ‘biomaterial’ or ‘stock’ tables. Properties are from ‘cvterm’ (sample metadata, traits) or ‘feature’ (genes, transcripts, proteins). The values are collected from any of the other tables following their references. Controlled terms are defined in the ‘cvterm’ table, but their assignment to a specific biomaterial and property are in their datatype-specific tables. The generic -omics query result table should be able to handle any of the dataset types, factors, and variable and values in Table 7.
The datasets are linked by sample IDs defined by their NCBI BioSample whenever available. Publications are required to define SAMN sample IDs when they submit sequences to NCBI, including NGS genotyping and RNA-Seq studies.
We also identified several ambiguous usage of terminology while integrating the Tripal modules. Table 8 summarizes the exact definition of tables mentioned as defined in Chado, or BioEntities as defined in Tripal.
For the ‘stock’ and ‘biomaterial’ tables, these appear to hold similar data. However, from the Chado definition, ‘stock’ represents the physical entity of an organism identified by unique name and stock type that is held within a living or preserved collection. In contrast, ‘biomaterial’ is a biological material (tissue, cells or serum) that may have been processed. Biomaterials should be traceable to source raw biomaterial, recognizing that unless vegetatively cloned, individual plants of a dioecious outcrossing species are likely to have different genetic makeup (haplotype). The Expression module uses the Biomaterial type, while the Phenotype module uses the Germplasm Accession type to represent the samples. Merging expression with phenotype data needs a link between stock and biomaterial which is not explicitly defined in Chado. The workaround implemented in the project is to use the same property and value (NCBI BioSample ID) in the stockprop and biomaterialprop tables. However, BioSample (SAMN) IDs may not be available since only sequences submitted to NCBI or other INDSC repositories have them. In this situation we assigned stocks/biomaterial ID using the simplest unique name deducible from the source publication.
Issues in the ‘organism’ and ‘genome assembly’ entities were also encountered. Although the portal handles only Cannabis sativa L. as the target organism, the ‘organism’ entity is used in Chado and Tripal to define a set of contigs and genes for a given genome. Thus, the different genome assemblies for cs10, PK and FN are stored as distinct records in the organism table. In Tripal, the Genome Assembly entity is defined as a computational analysis for a genome assembly in the analysis table. As implemented in the project, we loaded the different assemblies as separate organisms in the organism table, then modified all Organism labels to Assembly.
For the entities ‘Project’, ‘Study’, ‘Experiment’ and ‘Analysis’, by Chado definition, “Study” represents an experiment that produced microarray data; “Project” is a flexible property table for projects; “Analysis” is a single unit of computational analysis; “Assay” consists of a physical instance of an array, when combined with conditions create an Array.
The different types of computational analysis are reflected in Tripal by using the specific analysis entities. Currently defined are Genome Annotation, Genome Assembly, BLAST Results, InterPro Results. The phenotype module loads the experiment (Experiment is used in the label) to the Project table and is associated to a genus and the samples to Germplasm Accessions. A genus is defined in the phenotype configuration to consist of a controlled vocabulary of Trait, Method and Unit, or optionally by Crop Ontology. “Germplasm Accessions” is a type of stock in Tripal so phenotype samples are loaded to the stock table. Phenotype studies are added to the Project table. The Expression and Biomaterials loader in the expression module associates the data to an Analysis and Organism, which was used earlier for assemblies.
The interchanging use of terms made it confusing to merge samples between different modules, and our workaround was to use common terms for both the stocksprops and biomaterialprops table.
The Tripal genotype module, which uses the Natural Diversity tables [88], was tested and found to be performing sub optimally especially when loading large genotyping data matrices generated from NGS variant calling. To solve the speed issue, the SNP-Seek software system [65], developed to host large SNPs datasets, was instantiated in a podman container and embedded as a Tripal page. This feature, now named CannSeek database, is described in the accompanying paper. For the API calls to query SNPs, BCFtools (RRID:SCR_005227) [89] was used directly on the VCF files stored separately for each reference and dataset. The genotyped samples metadata are loaded to the ‘biomaterial’ table so they can be autoloaded from NCBI BioSample.
Some issues for the entities ‘gene’, ‘mRNA’, and ‘protein’ were also resolved. The central dogma of molecular biology relations is well-defined in the Tripal model and all are stored in the Chado feature table by different feature types. Their relationships are explicitly related as: mRNA part_of Gene, Exon part_of mRNA, Polypeptide derives_from mRNA. BLAST and InterPro Analysis Results were assigned to the specific molecule type used (i.e. DNA, mRNA or polypeptide). In the ‘chado_gene_search’ module, keyword search was modified to include the annotations assigned to all genes, mRNA and proteins by tracing their relationships. However, it became challenging to cross- reference genes between external databases and applications. The accessions or IDs will be different for the same gene but different molecule types. In addition, a gene locus can be associated with multiple mRNAs and proteins. To avoid the complicated query, a materialized view was created to map the corresponding genes, mRNA and polypeptides.
Replicates are challenging to handle when samples are measured, analyzed, and reported on multiple data types. Replicates may be biological, (i.e., multiple studies using the same plant source/genebank accession, multiple individuals in a study) or technical (i.e., multiple samples from a tissue with the same treatments). In NCBI, a biomaterial in BioSample (SAMN) can have multiple sequence replicates in SRA (SRR). Replicated NGS reads can be variant-called separately or jointly for a given SAMN. The measurement of metabolites for a SAMN can also have technical replicates. Some phenotype datasets loaded from publications report measurements per replicate [53], while others the mean and standard deviation [54]. The Tripal phenotype module documentation suggests loading the individual replicate data [88]. However, their proposed use-case is focused on geolocational and seasonal factors, which may be more applicable for a single project or institute, but not for a community portal that integrates from various data sources. These replicates create a Cartesian join between biomaterials and properties, that is, rows or columns will be duplicated except on one element. For this reason, we use average and standard deviation for technical replicates from published datasets. We pre-calculated these statistics before loading when only replicate values were available. It is also possible to have an array of values for a given biomaterial and property without a traceable replicate number. For publications that report in this format, the system can still store them and report in JSON format.
The -omics integration APIs: database query generation and web service implementation
The normalized structure of Chado makes database SQL queries complex and inefficient. To circumvent this, the Tripal modules prepare materialized views of the necessary data ready for user query and visualization. The same strategy was used to integrate the datasets from the different modules where they shared common BioEntities. Details of the data integration queries are described in (protocol step 16). In brief, five datatypes are queried separately, each returning four columns (i.e. datatype, property, sample, value). The datatypes are ID, PROP, PHEN, EXP and SNP. ID refers to sample identifications, PROP to sample attributes and metadata, PHEN are phenotype variables, EXP are transcripts or protein IDs, and SNPs are SNP IDs. The results are merged by a UNION operation, pivoted into a DATATYPE, PROPERTY by SAMPLE table, and returned to the client.
The API for omics data, implemented as a web service, usually returns JSON format text. They are returned for separate experiments, datasets, or omics type like BrAPI [90], or Tripal Web-Service [91]. Here we return an efficient and intuitive table format ready for further processing, which can be loaded easily into spreadsheets or any csv reader library. This API will work in complement with the Tripal Web Service module, which is also enabled in the ICGRC site [92]. The web-services module can supply the metadata of the biological entities, but not the values of the different -omics measurements. Specifically, the values this API queries are measurements in the materialized views taken from the Chado tables elementresult.signal (expression), phenotype.value (phenotype), stockprop.value (stock property values), biomaterialprops.value (biomaterial property values). The allele values for variant queries are read directly from the VCF files due to their large size if loaded to the Postgres database. The samples and metadata of the VCF are in the ‘biomaterial’ and ‘biomaterialprop’ tables.
To avoid Cartesian join duplicacy for replicates, complex element object representation like JSON was used for list of values. Being able to present values with replicates, their statistics, methods and units as a table element, the table still has the equivalent representation power of a JSON object typical in web-services APIs.
The Cannabis Tripal Web portal
The instantiated platform is available at https://icgrc.info [93]. Figure 4 shows the main features of the site. The functional menu includes Genomics, Phenome, Tools, and Infome. ’Genomics’ includes datasets pertaining to sequences, while ’Phenome’ includes phenotypes and gene expression data. ’Infome’ allows for management of Tripal Bioentities metadata, manuals, and collaborative information. The ’Genetic Resources’ section is reserved for data related to gene banks and other genetic resource collections. The ’Tools’ section is for performing computations using the datasets, currently BLAST is available. Details on the use are provided in the User’s Manual in the menu.
Drupal administrative modules allow site and user administration. Users are assigned roles of public, registered, data curator, site curator, or administrator. Registered users can view additional information including the metadata of all samples, biomaterials, analyses, or any bioentities made private to a group of users. They can also post comments on some pages and data. To register, only a name, affiliation and a business email address are required to be provided.
Most of the modules retain standard Tripal functionality, but some were revised to expand their capability. In ‘Gene function search’, keyword search now includes all annotation assigned to gene, mRNA, or protein. Using a Gene Ontology accession (GO:nnnnnn) includes all sequence features annotated with any of the transitive closure of the GO term. The expression module has an added feature to generate a file intended for the Biological Pathway viewer. The 21-trichomes transcript expression levels were also converted to bigwig format and visualized as JBrowse tracks at the coordinates identified with TransDecoder. Lastly, unique to this portal are the MapMan biological pathway expression viewer [84], and the CannSeek SNPs viewer. The gene search materialized view SQL definition and MapManJS HTML page revisions are in the protocols.io steps 11 and 18, respectively.
The protocol (DOI:10.17504/protocols.io.n2bvj3nz5lk5/v3) developed can also serve for future site update, and in setting-up new Tripal sites especially for crops where only sequence datasets are readily available.
The-omics API query
We designed and implemented an API to serve programmatic queries to the omics datasets. The available queries and parameters are described in the ICGRC Omics API documentation [94]. The general API query format is, https://icgrc.info/api/{user}/{datatype_list}/{constraints_list}?parameters=..&..&& where {user} could be an authenticated user with login and privileges, or just ‘user’ for public access. {datatype_list} is a comma separated set of datatypes to query with possible values of [gene, phenotype, expression, variant]. {constraints_list} is a comma separated list of genomic regions [chromosome:start-end] or keywords. Genomic region is used when Variant is queried. Keywords are queried on the gene functional annotations, as used by the ‘Gene Function Search’ page, to get genes for Expression or Variant queries. The available datasets and summaries are available at https://icgrc.info/api/user/{datatype}/dataset. The results of queries with multiple datatypes are combined by SAMPLE. For example, for samples with both expression and phenotype data, the sample column will have both genes and traits/metabolite rows.
DISCUSSION
Standardization for Cannabis vocabulary and nomenclatures
Standardization is mandatory to effectively use the platform, when hosting diverse datatypes and sources. Here we summarize terms that require community-wide standardization, their impact in storage and analyses, the current best practices from various research groups, possible immediate workarounds, and suggestions to prompt the establishment of accepted community standards.
Sequence and sequence features accession IDs
NCBI assigns unique sequence identifiers for all datasets it hosts and analyses. Similar ID assignments are implemented by other databases like Ensembl, UniProt, etc. External references for these identifiers are handled in the ‘dbxref’ table in Chado.
However, prior to publication, different research groups tend to generate their own IDs and nomenclatures. It would be highly beneficial if a standard in generating nomenclatures was adopted by the community. The standard alphanumeric coding system may consider the increments, versioning, feature types, relationships between feature types, and relationships with other relevant ontologies.
Stocks, cultivars, biomaterials, sample
This is a prevailing gene bank issue [95]; unfortunately there is no authoritative gene bank for Cannabis [96]. NCBI assigns sample SAMN IDs for submissions in BioSample. However, these IDs are more applicable for biomaterials, specifically laboratory materials. They can’t trace the pedigrees of genetic stocks.
In absence of a physical gene bank, the creation of a virtual gene bank was suggested [34]. A possible model for this is the Genesys platform[97]. Under this system, physical gene banks worldwide publish their passport and characterization data, which are then merged by Genesys into a common query interface. This or a similar entity would then have the authority to assign Stock or Germplasm accessions for Cannabis.
Traits, variables, methods, units for Cannabis
Detailed crop ontologies are available for well-studied crops, although these have their limitations [98]. These terms are often derived from legacy dictionaries used by gene banks in their passport transactions and characterization of collections. A standardized Cannabis crop ontology would benefit both the industry and scientific community to facilitate data curation and management. Moreover, using them during submission to public repositories simplifies their integration as mentioned above.
Multi-omics query using the API: use-cases
With the different omics data types loaded in the database, we are able to integrate these using API queries. We created Jupyter (RRID:SCR_018315) notebooks and Python modules to demonstrate the potential uses of the API for multi-omics analyses, all available in (protocol step 17). A static export of the notebook with outputs is in the ICGRC Omics Demo [99]. In the notebook, the python modules (module names are in courier font) perform specific tasks and use Pandas (RRID:SCR_018214) Dataframe for data handling and manipulation. Depending on the intended analysis, additional software or modules need to be installed in the client as described in the notebook.
Use case 1: Get phenotypes, normalize units, plot distributions and batch effects correction
https://icgrc.info/api/user/phenotype/list. This query returns a table of phenotype values from multiple datasets, in mixed units and format. When merging values from multiple sources, compatible variables need conversion of units. Units and standard deviation information for each value are included with flag with_stdunits=1. The values with different units can then be converted to a uniform unit (convert_units_to), for plotting and further processing. The main purpose of this feature is mainly for storage, integration, and retrieval. For the statistical validity of merging data from different sources and studies, we included a minimal feature (check_batcheffects) to detect and correct for batch effects and impute missing values using the samples metadata. This module uses ski-learn SimpleImputer and PCA [100], ppca [101] and pyComBat (RRID:SCR_010974)[102]. The metadata and sample attributes can also be used as covariates in various analyses.
Use case 2: Co-expression analysis using WGCNA
https://icgrc.info/api/user/expression,phenotype/list. This query returns a table with phenotype and gene expression values for each sample in the selected datasets. These data are used in Weighted correlation network analysis, WGCNA (RRID:SCR_003302) [103] in the (do_wgcna_prepare) and (do_wgcna_modulesgenes) modules. These modules return the top gene-trait pairs, with gene significance and p-value.
Use case 3: Matrix eQTL
https://icgrc.info/api/user/expression,variant/{querygenes}. This query returns a table of gene expressions and SNP allele data for genes and regions that meet the querygenes annotations constraints. The module (do_matrixeqtl) uses Matrix eQTL [104] for the top SNP-gene pairs by p-value. If provided with SNP and genes position information, Matrix eQTL can also analyze for local (cis) and distant (trans) pairs separately. The gene information is also queried using the API https://icgrc.info/api/user/gene, constrained by position for the SNPs, or by accession for the expression genes.
Use case 4: SNP-trait association study
https://icgrc.info/api/user/variant,phenotype/{querygenes}. This query returns a table of traits, mostly metabolite contents for cannabis, and SNP alleles for genes and regions in querygenes. The module (do_gwas) uses PLINK (RRID:SCR_001757)[105] –assoc function to get the top SNP- phenotype associations by p-value. When performed on all SNPs for a genome, and with enough samples this is basically a GWAS. To-date, there are only few public Cannabis samples that have both phenotype and genotype data publicly available, hence the limitation on datasets loaded in the portal. The -omics tools presented above used the API to retrieve multiple datasets for various analyses. Each analysis returns the top associations by p-value between genes, genotype variants, biomolecules, traits, or other factors. From the multiple evidence collected, it is logical next step to overlap the relationships into networks (union) or Venn diagrams (intersection) for further interpretation. The merge module (merge_matrixeqtl_wgcna_gwas) takes the results from WGCNA, eQTL, and plink -assoc, intersect the genes and traits, and merge the pairs to identify SNP gene-traits associations supported by two evidence, as illustrated in Figure 5.
The example use-cases illustrate the flexibility of the API for multi-omics, multi-source analyses. We showed two ways of merging multiple datasets. First is by merging datasets as input to an analysis tool. This requires shared samples, batch effect correction, covariates information, and most likely heavy imputations. Their important link is the shared genotype. In the second approach, the input for an analysis tool is from one dataset, but the results from different analyses are merged by shared biomolecules. Each approach has its own advantages and applications. Integration of sequence data is performed in the backend during variant calling, transcript assembly, or gene expression quantification.
CONCLUSION
We setup the Tripal platform to host available Cannabis genomics data available in the public domain. We also provided a protocol to generate various datasets from sequence data, to load, and test the different modules and visualizations. The protocol simplifies the process of regularly updating the site. One year since our last upload, a recent check on NCBI and PubMed returned new publications and datasets for potential reanalysis and loading to the site. This shows that Cannabis research is active and rapidly growing. The setup of this web-portal is timely, and ready to serve the research community. We also highlight the importance of using community standards especially when depositing datasets in public repositories. Standards minimize third-party curation and encourage reuse for meta-analysis. The multi-omics integration API and the example Jupyter notebooks demonstrate the convenience of doing analysis when the data are readily available in table format. Since they are built on top of Tripal modules, they can be readily adapted to any Tripal site for other research communities.
AVAILABILITY OF SOURCE CODE AND REQUIREMENTS
The source codes in this project can be grouped into: b. and c. are embedded or attached in the sections of the protocol described below. a. are available from the authors upon reasonable request.
- The ICGRC web application which is an instance of Tripal v3, the various Tripal/Drupal modules installed, site customizations/settings saved in the database, and the custom -omics API web server module developed in this project.
- the scripts and command lines part of the computational protocol for dataset generation from public datasets.
- Jupyter notebooks and python client library for -omics analyses which call the omics API
Project name: ICGRC Portal Tripal Data Generation and Setup
Project home page: https://www.protocols.io/private/34114B4120C511EF83720A58A9FEAC02
Operating system(s): Linux
Programming language: bash, PHP, Python
Other requirements: listed in the steps in the protocol
License: CC BY 4.0 International
DOI: dx.doi.org/10.17504/protocols.io.n2bvj3nz5lk5/v3 (reserved)
ABBREVIATIONS
API, Application Programming Interface; cs10, CBDRx cultivar assembly; FN, Finola cultivar assembly; GEO, Gene Expression Omnibus; ICGRC, International Cannabis Genomics Research Consortium; NCBI, National Center for Biotechnology Information; PK, Purple kush cultivar assembly; QTL, Quantitative Trait Loci; RefSeq, Reference Sequence Database; SAMN, BioSample accession; SNP, Single Nucleotide Polymorphism; SRA, Sequence Read Archive; SRR, SRA accession
DECLARATIONS
Ethical approval
The authors declare that ethical approval was not required for this type of research.
Competing Interests
None
Acknowledgements
The authors acknowledge the provision of computing and data resources provided by the Australian BioCommons Leadership Share (ABLeS) program. This program is co-funded by Bioplatforms Australia (enabled by NCRIS), the National Computational Infrastructure and Pawsey Supercomputing Centre.
Funding
This study was funded by the Australian Research Council (ARC) Linkage project LP210200606. In addition, first author Locedie Mansueto received a stipend from Southern Cross University (SCU). The ICGRC and CannSeek web servers are hosted and funded by SCU.