Methodology

We used six different methods to find interesting genes in our study: eggNOG-mapper, Orthofinder, TapScan, Gene Ontology, best blast hit against Araport11 and InterProScan. Each of them approach a similar question with a different method and we believe the collective look can help us to target interesting genes better than individual approaches. Please read the corresponding paper for each tool if you need a comprehensive explanation of the tools.

In the context of the analyses described in this section, it is vital to emphasize that only protein sequences from the indicated species were used. Consequently, this collection does not include non-coding RNA sequences. Alternative methodology or customized criteria-based studies may be necessary for investigating non-coding RNA or identifying specific genes of interest, potentially necessitating independent study or the use of distinct analytical approaches.

OrthoFinder is a software for identifying and categorizing orthologous genes across species or genomes. It employs an algorithm that considers both sequence similarity and gene evolutionary relationships. By creating a similarity graph based on pairwise comparisons of protein sequences, OrthoFinder clusters genes into orthogroups, which represent genes descended from a common ancestral gene. This approach ensures accurate grouping, even for genes with complex evolutionary histories.To obtain the best results from Orthofinder, we must include a phylogenetically diverse sample of species for orthogroup inference. As a result, this section contains more species.

EggNOG-mapper is a tool for assigning functional annotations to genes by mapping them to orthologous groups and functional categories in the EggNOG database. Users submit gene sequences, and EggNOG-mapper detects orthologous groups and provides functional insights based on known functions in the database. This technology simplifies gene annotation and helps us understand gene functions and evolutionary links across species.

With the increased availability of genomic data, TapScan is a tool designed specifically for plant genome analysis to identify transcription factors with high precision.

Gene ontology (GO) is a database and standardized approach for classifying and categorizing genes and gene products based on their molecular functions, biological processes, and cellular components. It establishes a common vocabulary and framework for identifying gene features and behaviors across various animals, making it easier to evaluate Genomic data and better understand gene function in biological systems. GO organizes genes into hierarchical categories and links, allowing researchers to perform comparative analysis, identify functional similarities, and gain biological insights across species. We integrated the InterProScan and eggNOG-mapper results and filtered the “obsolute” GO IDs to generate comprehensive GO-Gene-Term sets.

InterProScan helps also with the functional annotation of protein sequences. By entering protein sequences into InterProScan, we may detect conserved domains, motifs, and functional signatures. This information contributes to our understanding of these proteins’ potential roles and activities in biological processes. Furthermore, InterProScan can predict protein families and connect sequences to known biological pathways, which is critical for understanding the molecular mechanisms that underpin many biological occurrences.

Unfortunately, the multiple sequence alignment viewer is not compatible with the Firefox browser. It’s a bug in the package we use for MSA visualization, and we can’t do anything about it. We checked, and MSA viewer works fine with Chrome, Brave, and Safari. Other browsers have not been checked.

The parameters that were used for each program are listed below.

Orthofinder

We used Orthofinder with two different settings, and used the results of the second run since we believe it is more accurate. The species tree that we used for the second run is shown below.

# 1st run
orthofinder.py -S diamond -M msa -A mafft -T fasttree -t 50 -a 6 -y -n run_1
# 2nd run
orthofinder.py -t 50 -a 6 -y -n run_2 -ft Results_run_1 -s SpeciesTree_input.txt

eggNOG-mapper

emapper.py -m diamond --itype proteins --data_dir eggnog-mapper-data/ --dmnd_iterate yes --dbmem --cpu 0 --evalue 1e-10 --sensmode ultra-sensitive --tax_scope 33090 --dmnd_db eggnog-mapper-data/eggnog_proteins_default_viridiplantae.dmnd

InterProScan

interproscan.sh -cpu 150 -pa -goterms  &> iprsc.log 

Introduction

We utilized principal component analysis (PCA) to extract the core of the data and find the most significant sources of variation between samples. Multidimensional Scaling (MDS) goes one step further, allowing us to see the subtle interactions between samples, discovering underlying structures that would otherwise go undiscovered. Finally, our Hierarchical Clustering reveals natural data groupings, allowing us to identify physiologically significant groups.

Methods

We employed a range of R packages for our exploratory data analysis, and you can find the comprehensive list in the Methods section. To give an overview, our process involved importing Kallisto quantification files using the txImport tool and applying a lengthScaledTPM transformation. We then filtered the data, retaining reads with a minimum Count-Per-Million (CPM) of at least 10 across a minimum of three samples. Ensuring robustness, we conducted Smooth Quantile Normalization, accounting for varying experimental conditions, and subsequently transformed the data using the voom function from the limma package. To perform hierarchical clustering, we computed distances using the Euclidean method and opted for the “ward.D”” method for agglomeration.

A note on abbreviations

Each treatment has a name composed of two components separated by an underscore. The first portion is the treatment, and the second part is the time (in hours) since the start of the experiment. We had two highlight settings: stress (“s”) and recovery (“r”).

Visualization of a single gene expression profile

Here, we provide the chance to see how particular genes are expressed. We think it’s critical to analyze gene expression using a variety of visualization techniques, such as the boxplot, dot plot, and heatmap. By providing a distinct viewpoint, each approach makes it possible to decipher the intricate workings of gene regulation. Because of its simplicity, the boxplot shows central tendencies and outliers in a clear description of the distribution of gene expression. Conversely, the dot plot provides a finer perspective by displaying individual data points together with their distribution. In the meantime, the heatmap presents gene expression across samples in a clear and understandable way, painting a complete picture. In addition, we applied the Z-score transformation, a statistical method that normalizes gene expression levels to facilitate comparisons between various genes and datasets.

A note on abbreviations

Each treatment has a name composed of two components separated by an underscore. The first portion is the treatment, and the second part is the time (in hours) since the start of the experiment. We had two highlight settings: stress (“s”) and recovery (“r”).

Introduction

This section covers differential gene expression analysis. Our primary goal in this part was to compare individual treatments to appropriate control samples (at the same time point if possible) in order to better understand how the treatments affected gene expression patterns.

Methodology

To begin our differential gene expression analysis, we used pre-processed count data (normalized with qsmooth) transformed with the voom function from the limma package. Our analytical technique starts with the formulation of a contrast matrix. When we compared the treatment group to a properly selected control sample, we estimated the fold change and adjusted p-value values using gene expression shifts. Our research was supported by a number of essential functions from the limma package, such as lmFit, contrast.fit, eBayes, and decideTests. These methods helped us discover genes with differential expression (DEGs). We next visualized the dynamic gene expression profiles using ggplot2 and a hierarchical clustering technique. For each comparison, genes were divided into two unique groups: “up-regulated” and “down-regulated” using the cutree algorithm. For the majority of our experimental samples, only one treatment and control sample combination was available. However, in cases where this matching was impossible due to experimental restrictions, we used a cautious approach. We used the nearest accessible time point as a surrogate reference for comparison. After that, we looked at over-representation analysis, which is a critical step in understanding the functional implications of DEGs. Using the clusterProfiler package, we attempted to identify Gene Ontology (GO) concepts that were strongly enriched inside each genetic cluster. Our background set included all expressed genes in our samples.

Parameter selection

The default values are suitable, but you can adjust them. Please be patient while a new plot or table is loading. It takes time to do the calculation. If there are no enhanced terms, you will see an empty plot/table. Then perhaps you might choose a different set of parameters or a different comparison.

Introduction

Co-expression network analysis with RNA-Seq data provides insights into the complex regulatory mechanisms that control gene expression. We can find functional relationships, identify important regulatory genes, and discover biological pathways by analyzing the correlation patterns between gene expression profiles across diverse situations. This method enables the discovery of co-regulated genes that may be implicated in shared biological processes or pathways, providing a comprehensive understanding of gene expression dynamics. Furthermore, co-expression networks provide a systematic framework for selecting candidate genes for additional experimental validation, making it possible to uncover novel biomarkers in a variety of biological situations.

Dirichlet process and Gaussian process (DPGP)

Dirichlet and Gaussian process models have emerged as promising tools for capturing temporal dependencies and dynamic interactions among genes in co-expression network research using time series data. Dirichlet process models provide a versatile non-parametric framework for clustering genes into co-expression modules while supporting different cluster sizes and architectures. Gaussian process models, on the other hand, offer a robust probabilistic framework for modeling time-varying gene expression trajectories, incorporating both deterministic and stochastic variations. Integrating these Bayesian modeling tools with RNA-Seq time series data allows for the identification of dynamic gene regulatory networks, which sheds insight on the temporal dynamics of gene expression and regulatory interactions in biological systems. If you want to find out more about this tool, please read this paper. Please remember that if no terms get enriched in an analysis there will be no table and plot for that specific setting. You may get a table or plot by tuning the parameters below.

Please be patient. Loading the results in this section take a few seconds.

Introduction

Co-expression network analysis with RNA-Seq data provides insights into the complex regulatory mechanisms that control gene expression. We can find functional relationships, identify important regulatory genes, and discover biological pathways by analyzing the correlation patterns between gene expression profiles across diverse situations. This method enables the discovery of co-regulated genes that may be implicated in shared biological processes or pathways, providing a comprehensive understanding of gene expression dynamics. Furthermore, co-expression networks provide a systematic framework for selecting candidate genes for additional experimental validation, making it possible to uncover novel biomarkers in a variety of biological situations.

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA is a widely used method in bioinformatics for creating and analyzing co-expression networks. WGCNA identifies modules of highly connected genes, allowing to discover biologically significant gene clusters associated with specific phenotypes or conditions. By weighting gene-gene correlations based on their significance, WGCNA provides a strong framework for finding co-expression patterns and choosing genes with high relevance to the biological context under study. For more information about this tool, see this paper. For network construction, we used the following settings: merge threshold = 0.20, network type = signed, TOM type = signed, Min. module size = 30, Max. P. outliers = 0.05. We used 20 as soft threshold for M. endlicherianum and P. patens and 13 for Z. circumcarinatum based on our screening for scale-free network properties.

Please use the following acronym table for measuring metabolites:

Please be patient. Loading the results in this section take a few seconds.

Introduction

Gene regulatory networks are the intricate webs of interactions between genes and regulatory elements that control biological activities. Understanding these networks is essential for determining the underlying mechanisms that control biological systems. Time-series RNA-Seq data provide a unique chance to dynamically infer gene regulatory networks by capturing the temporal dynamics of gene expression changes in response to different stimuli or perturbations. By studying gene expression profiles across time, we can discover regulatory relationships, identify important regulatory genes, and figure out the dynamic interactions that shape gene expression programs. Time-series RNA-Seq data are thus an effective tool for deciphering the temporal characteristics of gene regulation networks and understanding how they influence cellular behaviours and responses.

Sliding Window Inference for Network Generation (SWING)

Despite the availability of time-resolved, high-throughput data, many algorithms ignore the temporal delays inherent in regulatory systems, resulting in unreliable network inferences. SWING is used to address this issue by only taking temporal information into account when identifying time-delayed edges. SWING’s tolerance to user-defined parameters allows for the successful identification of regulatory mechanisms from time-series gene expression data. SWING uses multivariate Granger causality to capture the regulatory relationships between genes throughout time. SWING, unlike traditional Granger approaches, uses a sliding window approach to evaluate numerous upstream regulators at the same time over a range of time delays. If you would like to investigate more about this tool, please see this paper.

The algorithm is O(2^N), indicating that the number of input genes has a significant impact on the computing time. It was not possible to use all expressed genes as input for this tool. We filtered for transcription factors as well as stress response gene homologs of A. thaliana (identified using sequence similarity search). For approximately 3,000 genes and 5 to 9 time points, the calculation took 4 to 16 days to complete. As a result, these databases do not include all genes or their interactions. We did two analyses. First, we include all time points where we measured RNA-Seq and metabolite levels for each species under treatment. Second, we only considered RNA-Seq data for each species subjected to a treatment. For M. endlicherianum in both instances and transcript-only data sets of P. patens and Z. circumcarinatum-1b, the following settings were used as input: k_min = 0, k_max = 1, w = 4, method = RandomForest, trees = 500, and lag_method='mean_mean'. For P. patens and Z. circumcarinatum-1b, we utilized the following settings for the data set of RNA-Seq and metabolites: k_min = 0, k_max = 1, w = 2, method = RandomForest, trees = 500, lag_method="mean_mean".

Please use the following acronym table for measuring metabolites: