Genome assemblies and annotations

The first peer-reviewed publication for cannabis genomes was on Purple Kush (PK), a medicinal type tetrahydrocannabinol (THC)-dominant marijuana cultivar, and Finola (FN), a hemp-type cultivar developed for oil seeds [11]. Both PK and FN assemblies were updated and improved in Laverty et al. [18]. CBDRx (cs10) is a cannabidiol (CBD)-type cultivar derived from a cross of Skunk#1 (marijuana) and 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 was used by many publications for gene assignment. Pink Pepper (ASM2916894v1) is the second assembly annotated by RefSeq recently (November 2023).

Three primary genome assemblies are currently loaded to the portal: cs10 (also known as CBDRx, a CBD-dominant line) [12], PK and FN [18]. 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 the cultivars First Light [12], Jamaican Lions [21], Cannbio-2 [20], JL wild [19], and Pink Pepper, were not included in this release of the portal.

The annotations from RefSeq for cs10 and our FINDER [53] annotation for PK and FN were loaded and are searchable using the ‘Gene Function Search’ module interface, and are displayed in the JBrowse browser. The query constraints in the functional annotations are cultivar/assembly, annotation analysis, genome region, gene or transcription accessions, and keyword. Query by sequence is also possible using the BLAST interface (RRID:SCR_004870). 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 are 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, 20, 54]. 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 different genomes. As with gene loci, the transcript assemblies can be queried on 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 and mostly explore changes relative 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 regions [58]. In another study, FN trichomes were separated from mid-stage flowers and categorized into bulbous, sessile and mature stalks [55]. Terpene and cannabinoid profiles were measured, and the correlation with synthase gene transcript levels was assessed [56, 57]. 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 Tables 3 and 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.

At the completion of this project, only one cannabis proteome expression dataset had been made available in GEO. Using quantitative time of flight mass spectrometry, 1,240 proteins from glandular trichome heads, 396 from trichome stalks, and 1,682 from whole flower tissue isolates were identified by Conneely et al. [59].

Variants discovery using resequencing datasets

Details on the reanalysis of Next Generation Sequencing (NGS) datasets of cannabis are described in an accompanying paper [60] about the CannSeek single nucleotide polymorphisms (SNPs) database (RRID:SCR_025579) [61]. Data available in NCBI SRA as of December 2022 consisting of 26 trichome RNA sequencing (RNA-Seq), 383 genomic DNA, and the Phylos amplicons were reanalyzed using GATK (RRID:SCR_001876) [62] and Parabricks [63] variant calling pipelines. The results are displayed using the CannSeek software, derived from Rice SNP-Seek [68], hosted in a Podman or Docker container [69, 70] and embedded as a Tripal page.

Maps

The interactive Map Viewer, a part of the Tripal platform [71], can display an entire genome or chromosome, and features such as markers and quantitative trait loci (QTLs) in a single view. Maps positions may either be physical (using base position, or bp) or genetic (using centiMorgan, or 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 [72] and QTL mapping for agronomic and biochemical traits [73].

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’ [74]. 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 [75] for chemical traits, PECO [76] for method and UO [77] 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 matches in the existing ontologies were found.

Gene model prediction and annotations

All resequencing data were trimmed using Trimmomatic (RRID:SCR_011848) [78] (protocol step 1). Cs10 genes and annotations were downloaded from NCBI RefSeq (protocol step 2). The FINDER [53] pipeline was used to predict gene models for updated assemblies of PK and FN using cultivar-specific RNA-Seq data. GffRead (RRID:SCR_018965) [79] 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 [50] databases using MMseqs2 (RRID:SCR_022962) [80]. Gene Ontology (RRID:SCR_002811) and InterPro (IPR) (RRID:SCR_006695) protein domain annotation of these proteins were done using InterProScan (RRID:SCR_005829) [81].

Gene expression

Expression data from GEO [49] 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 was quantified on 21 trichome-specific public RNA-Seq datasets [54, 55, 57, 85] 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) [86], then mapped to the three references using GMAP (RRID:SCR_008992) [87], and ORFs were identified using TransDecoder (RRID:SCR_017647) [88]. Quantification of transcript abundance was then performed using Salmon (RRID:SCR_017036) following (protocol step 5).

Synteny block discovery

This follows the Tripal synteny viewer instructions for generating synteny blocks using gene annotations (protocol step 8). BLASTN (RRID:SCR_001598) and MCScanX (RRID:SCR_022067) [89] were used on the cs10 gene annotation and our generated predictions for PK and FN.

BLAST database

BLAST databases for the three 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, PK, and FN. Cs10 gene model sequences are from NCBI Refseq, while PK and FN are from the FINDER gene prediction pipelines.

JBrowse Genome Browser

The JBrowse Genome Browser [35] displays genome features along the reference chromosome or scaffolds. The data is loaded from interval files like General Feature Format (GFF) or Browser Extensible Data (BED) format. 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) [90], 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 the study is loaded into the Expression module.

SNPs discovery and genotype data management

Variant data of the 21 trichomes RNA seq were discovered using the GATK-RNASeq (RRID:SCR_001876) pipeline (protocol step 6). Genomic resequencing data for 398 samples were also used to discover SNP variants using the GATK germline SNP calling workflow [91] and Parabricks [63].

Genotype data results from variant discovery pipelines are presented as a large genotype matrix of samples and base position. The Tripal Natural Diversity Genotypes [92] module used in KnowPulse [74] 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 [68] can display thousands of SNP positions and samples while maintaining interactivity, allows more analyses of the result and was thus used in this project.

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 [93]. 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 a thorough understanding and experience with 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 using 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, variables and values in Table 7.

The datasets are linked by sample IDs defined by their NCBI BioSample whenever available. Publications are required to define BioSample accessions (SAMNs) 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.

The ‘stock’ and ‘biomaterial’ tables appear to hold similar data. However, according to Chado’s definition, ‘stock’ represents the physical entity of an organism identified by a 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, 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 IDs using the simplest and most unique name deduced 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’s 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 which, when combined with specific conditions, creates 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 and InterPro Results. The phenotype module loads the experiment (Experiment is used in the label) to the Project table and is associated with 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 into 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 [94], 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 [68], developed to host large SNP 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) [95] 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, the keyword search was modified to include the annotations assigned to all genes, mRNAs 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 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 [57], while others report mean and standard deviation [85]. The Tripal phenotype module documentation suggests loading the individual replicate data [94]. 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 various data sources. These replicates create a Cartesian join between biomaterials and properties; that is, rows or columns will be duplicated except for 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 and value). The data types 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, omics types like BrAPI [96] or Tripal Web-Service [97]. 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 [98]. The web-services module can supply the metadata of the biological entities but not the values of the different omics measurements. Specifically, the values queried by this API 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 duplication for replicates, complex element object representation, like JSON, was used to manage the 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.

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 and UniProt. 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 for generating nomenclatures were 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, samples

This is a prevailing gene bank issue [101]; unfortunately, there is no authoritative gene bank for cannabis [102]. NCBI assigns sample SAMN IDs for submissions in BioSample. However, these IDs are more applicable for biomaterials, specifically laboratory materials. They cannot trace the pedigrees of genetic stocks.

In the absence of a physical gene bank, the creation of a virtual gene bank was suggested [36]. A possible model for this is the Genesys platform [103]. 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.

NCBI BioSample metadata submissions

Related to biomaterials, the NCBI BioSample SAMN metadata can be loaded automatically into Tripal. Thus, terms used in NCBI submissions are reflected in the Tripal ’Biomaterial’ pages. The sample metadata are also in the ID, and PROP rows are returned to API queries. Therefore, it is beneficial to use standardized terminologies during submissions to public repositories, like NCBI, to use the information efficiently and accurately in subsequent analyses.

Querying PubMed for publications relating to cannabis plant biology or agronomy is complicated by a wealth of papers pertaining to clinical or pharmacological studies. We therefore recommend using standard keywords or MESH terms specific to the study of the plant itself, not its clinical/pharmacological use.

Traits, variables, methods and units for cannabis

Detailed crop ontologies are available for well-studied crops, although these have their limitations [104]. 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 in facilitating data curation and management. Moreover, using them during submission to public repositories simplifies their integration, as mentioned above.

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 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 [106], ppca [107] and pyComBat (RRID:SCR_010974) [108]. 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, or WGCNA (RRID:SCR_003302) [109], 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 [110] for the top SNP-gene pairs by p-value. If provided with SNP and gene 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 gene expression set.

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 PLINK1 (RRID:SCR_001757) [111] –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 genome-wide association study (GWAS). To date, there are only a 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, the logical next step is 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, intersects the genes and traits, and merges the pairs to identify SNP gene-traits associations supported by two pieces of evidence, as illustrated in Figure 5.

The example use cases illustrate the flexibility of the API for multi-omics and multi-source analyses. We showed two ways of merging multiple datasets. First, datasets are merged as input for an analysis tool. This requires shared samples, batch effect correction, covariates information, and most likely heavy imputations. Their essential 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 advantages and applications. Integration of sequence data is performed in the backend during variant calling, transcript assembly, or gene expression quantification.