---
title: "End-to-End De Novo Transcriptomics Workflow"
canonical: "https://help.biobam.com/space/OED0324/4824793089/End-to-End+De+Novo+Transcriptomics+Workflow"
format: markdown
---
- [Introduction](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
- [Bioinformatic Analysis](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
- [1. Data Preprocessing](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
- [2. De Novo Assembly & Refinement](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
- [3. Differential Expression Analysis](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
- [4. Enrichment Analysis](https://biobam.atlassian.net/wiki/spaces/~712020bf3a1778a7ec4a9d97c54cc083d01c86/pages/4767088641)
## **Introduction**
In non-model transcriptomics, the ultimate goal is to translate raw sequencing data into clear biological insights. While establishing a high-quality *de novo* transcriptome is a fundamental requirement, the true value of the experiment lies in accurately quantifying transcript abundance and identifying the specific biological processes that respond to external stimuli. This tutorial demonstrates a complete, high-resolution workflow in OmicsBox, moving from raw data preprocessing and assembly to robust differential expression and gene set enrichment analysis. By following this pipeline, we will identify the molecular mechanisms behind the tolerance of *Rorippa indica* to the mustard aphid, identifying not just which genes are present, but which ones drive the plant's immune response.
### **Dataset Description**
The dataset used in this tutorial comes from an mRNA-Seq study of the wild mustard *Rorippa indica*, a plant known for its significant tolerance to the mustard aphid (*Lipaphis erysimi*). The data consists of four paired-end libraries representing two experimental conditions: uninfested control plants and aphid-infested plants harvested at 24 hours post-infestation.
These data were originally generated and published by [Bandopadhyay et al. (2024)](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA792073) as part of a study to perform a *de novo* transcriptome assembly and global analysis of differential gene expression in response to biotic stress.
- **Organism:** *Rorippa indica*.
- **Experimental Design:** Two conditions (control and infested), each with two libraries prepared by pooling three biological replicates.
- **Sequencing:** 150 bp paired-end reads generated on the Illumina HiSeq 2000 platform.
- **Data Access:** [PRJNA792073](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA792073).
## **Bioinformatic Analysis**
### **1. Data Preprocessing**
In a *de novo* assembly project, the absence of a reference genome makes data purity the most critical factor for success. Non-target sequences or low-quality bases can lead to chimeric contigs or an artificially inflated transcriptome. For the *Rorippa indica* dataset, we implemented a rigorous preprocessing pipeline to ensure only high-quality plant mRNA is processed downstream.
### **1.1. FASTQ Quality Check**
This tool performs an initial assessment of the raw sequencing reads to identify potential issues such as per-base quality distribution drops or adapter contamination. Looking at these results is essential to deciding the parameters for further read preprocessing.
#### **Application**
The FASTQ Quality Check tool can be found under General Tools → FASTQ Tools → Short Read Quality Assessment with FASTQC. For further information, refer to the OmicsBox User Guide section[ FastQC Quality Check](https://docs.omicsbox.biobam.com/latest/FastQC-Quality-Check/).
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/1R2k_BszQ3SgwGb3RUGCjWRmx4iOr7w7l): Raw sequencing data from 4 libraries (2 control, 2 infested) in FASTQ format.
**Parameters**
Provide the FASTQ files and leave all parameters by default.
**Output**
The analysis generates two primary outputs that must be reviewed to guide downstream preprocessing:
- **FASTQC Summary Report:** This report evaluates the library using nine quality modules, assigning a status of "PASS" (normal), "WARNING" (slightly abnormal), or "FAIL" (unusual).
[Summary Report](https://drive.google.com/file/d/1E3RjccwTv6VCjH56mVulzt1a7Iz7Bcbx/view).
- **FASTQC Results Table:** This table provides a high-level overview of the library composition statistics. To generate the specific charts shown in this tutorial, you can use two methods:
- Directly from the Results Table: Right-click on a sample row and navigate to the Show Statistics submenu to select a specific module (e.g., Per Base Sequence Quality).
- From the Individual Report: Open the report of a sample (Right-click → Show report). Inside this new tab, you will find a summary table where the Details column contains buttons to open each specific chart.
For this specific dataset, the Per Base Sequence Quality profiles (Figure 1) showed median and mean quality scores consistently above 30 across all positions. As the values remain within the "green zone" of the plot, the base calling is highly reliable throughout the entire sequencing run. However, the Adapter Content module (Figure 2) revealed a cumulative increase of Illumina Universal Adapter sequences towards the 3' ends of the reads. The detection of these sequences is a critical indicator that an adapter trimming step is mandatory.
[Results Table](https://drive.google.com/file/d/1QFVRRdyb8mHWHtORGeqx9S1J0kHChf8s/view).
**Figure 1.** Per-base sequence quality profile for *R. indica* reads (sample SRR17346648).
**Figure 2. **Adapter Content for R. indica reads (sample SRR17346648).
### **1.2. Trimming**
After identifying unusual metrics in the data, it is highly recommended to remove them to prevent interference with the assembly process. This step uses Trimmomatic for targeted adapter removal and safety pruning.
#### **Application**
The FASTQ Quality Check tool can be found under General Tools → FASTQ Tools → Short Read Preprocessing with Trimmomatic. For further information, refer to the OmicsBox User Guide section [Fastq Preprocessing with Trimmomatic](https://docs.omicsbox.biobam.com/latest/Fastq-Preprocessing-with-Trimmomatic/).
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/1R2k_BszQ3SgwGb3RUGCjWRmx4iOr7w7l): Same input as the previous assessment step.
**Parameters**
- Paired-End: Illumina sequencing reads.
- Remove Adapters: True.
- Quality Trimming: True.
Leave the remaining parameters by default.
**Output**
The preprocessing step generates a set of cleaned files and a report summarizing the filtering efficiency.
- **Processed FASTQ files:** For each sample, four files (.fastq.gz) are produced: two containing paired reads (where both members of the pair passed the filters) and two containing unpaired reads (where only one member survived). These cleaned paired reads are the primary input for the Trinity assembly.
[Paired reads](https://drive.google.com/drive/folders/1MhkMi80pCivG5aqkmeRlXUXd9aI3_8Ld).
- **Trimmomatic Summary Report:** This report displays the number of "surviving" reads. In this dataset, more than 97% of the raw reads were retained after clipping, indicating that the library preparation was robust and most sequences were of high biological value.
[Summary Report](https://drive.google.com/file/d/1tMqwgqVdfN1Fit2_8JUQVIjHH9KgqWe6/view).
- **Validation via Post-Trimming FastQC:** To confirm the effectiveness of the process, a second quality assessment was performed on the output files using the same configuration as the initial FASTQ Quality Check step (1.1). More importantly, the Adapter Content plot (Figure 3) shows a flat 0% baseline across all positions, confirming that the residual Illumina adapters identified in the raw data were successfully clipped.
[Summary Report](https://drive.google.com/file/d/10PeD19Mzs0cUq6H8I3WvYIGxZ7Ag1v1D/view).
[Results Table](https://drive.google.com/file/d/1wg9cs6ByflyLqyteHkjWGOvYvfpOfbkB/view).
**Figure 3. **Adapter Content for *R. indica* reads after trimming (sample SRR17346648).
### **1.3. Contaminant Removal**
In this study, working with an infested plant means the library is a mixture of *Rorippa indica* and the mustard aphid (*Lipaphis erysimi*). Furthermore, human-derived sequences (*Homo sapiens*) can be amplified during PCR, potentially contaminating the final library. This tool ensures the assembly pipeline only processes plant-specific mRNA by removing these "off-target" sequences. Since a reference genome for the mustard aphid was unavailable, a proxy-species strategy was employed using the genome of the cabbage aphid (*Brevicoryne brassicae*).
#### **Application**
The Contaminant Removal tool can be found under Metagenomics → Taxonomic Classification → Contaminant Removal. For further information, refer to the OmicsBox User Guide section [Contaminant Removal](https://docs.omicsbox.biobam.com/latest/Contaminant-Removal/).
#### **1.3.1. Contaminant Removal with *****Homo sapiens***** genome**
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/1MhkMi80pCivG5aqkmeRlXUXd9aI3_8Ld): Files obtained from the previous preprocessing step (trimmed).
**Parameters**
- Sequencing Data: Paired-End Reads.
- Reads: Illumina sequencing reads.
- Database Index: Homo sapiens (grch38).
Leave the remaining parameters by default.
**Output**
This analysis generates four primary outputs:
- **Cleaned FASTQ files**: These files contain the sequences that have not mapped against the Human genome ("non-contaminant" reads).
[FASTQ files](https://drive.google.com/drive/folders/18H9QW7o8g81xEnEdHTAjggPZAzj4qkNz).
- **Summary Report**: A table summarizing the number of reads processed and the filtering efficiency.
[Summary Report](https://drive.google.com/file/d/1fWLJ6eUEoVZpqwIQyLRMwVxeL3zZ776N/view).
- **Alignments per Category**: A bar chart showing absolute counts of mapped (contaminant) and unmapped reads.
[ Alignments per Category](https://drive.google.com/file/d/10dngh6E7TRG8EtwkKXuDid18NKKsJafa/view).
- **Relative Alignments per Category**: A percentage-based visualization (Figure 4). For this dataset, the Relative Alignments per Category plot shows that nearly 100% of the reads across all samples are categorized as "Unmapped" (represented in red). This indicates that human DNA contamination is negligible in these libraries, confirming high sample purity and allowing the unmapped reads to proceed to the next filtering stage.
[Relative Alignments per Category](https://drive.google.com/file/d/1b7A6b7W-IWzVjerpBDd_tLJ5Pbuca0Vt/view).
**Figure 4. **Read alignment distribution against *H. sapiens* genome.
#### **1.3.2. Contaminant Removal with *****Brevicoryne brassicae***** genome**
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/18H9QW7o8g81xEnEdHTAjggPZAzj4qkNz): Filtered FASTQ files from the contaminant removal step against the *Homo sapiens* genome.
**Parameters**
- Sequencing Data: Paired-End Reads.
- Reads: Illumina sequencing reads.
- Database Index: Create database from genome.
- Target genome:* **[Brevicoryne brassicae](https://drive.google.com/file/d/1rPh4BFzENHffwOBuRcxjq5gBNZ4UIy9V/view)** *genome assembly.
Leave the remaining parameters by default.
**Output**
Similar to the previous step, this analysis produces four outputs:
- **Filtered FASTQ files:** Plant-specific reads, free from aphid contamination (the sequences that have not mapped against the aphid’s genome).
[FASTQ Files](https://drive.google.com/drive/folders/1NelfF_KQD1CA3GDWVCQ_AgTNDC1Iqmaa).
- **Summary Report:** Detailed mapping statistics for each library against the aphid genome.
[Summary Report](https://drive.google.com/file/d/1Ypov-9k3aZAXab2cDKgQbVaDH-ZVKNVf/view).
- **Alignments per Category:** Absolute distribution of aphid-aligned (green) vs. plant-specific reads (red).
[Alignments per Category](https://drive.google.com/file/d/1jJTF4J0zCRkkHxtyxnR5vhAGr9Lg7hFt/view).
- **Relative Alignments per Category:** Percentage-based distribution plot (Figure 5). The Relative Alignments per Category plot reveals that the majority of the reads remain "Unmapped" (red category), meaning they do not belong to the aphid genome. These filtered sequences represent the high-quality, host-specific mRNA required for a reliable Trinity assembly. By removing the small fraction of aphid-aligned reads (green category), we prevent the formation of chimeric transcripts during the assembly process.
[Relative Alignments per Category](https://drive.google.com/file/d/1jUHsntECe5kj5BfHti7e2-lJsIVkeIgJ/view).
**Figure 5.** Read alignment distribution against *B. brassicae *genome.
### **1.4. rRNA Removal**
Ribosomal RNA (rRNA) often comprises up to 90% of the total RNA pool and can contaminate libraries despite experimental enrichment kits. SortMeRNA is used to deplete these fragments, saving computational resources and reducing noise.
#### **Application**
The Ribosomal RNA Removal tool can be found under Metagenomics → Metatranscriptomics Tools → rRNA Removal. For further information, refer to the OmicsBox User Guide section [rRNA Removal with SortMeRNA](https://docs.omicsbox.biobam.com/latest/rRNA-Removal-with-SortMeRNA/).
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/1NelfF_KQD1CA3GDWVCQ_AgTNDC1Iqmaa): Filtered FASTQ files from the contaminant removal step against *Brevicoryne brassicae* genome.
**Parameters**
- Sequencing Data: Paired-End Reads.
- Reads: Illumina sequencing reads.
Leave the remaining parameters by default.
**Output**
The SortMeRNA analysis generates two primary outputs:
- **Cleaned (mRNA-enriched) FASTQ files:** Files containing only the reads that did not match any ribosomal database. These are the final sequences used for the assembly.
[FASTQ Files](https://drive.google.com/drive/folders/1voV4mM83myRBL3P-ofPwAqcFzZ8-KR2H).
- **rRNA Removal Summary Report:** A statistical table summarizing the depletion efficiency for each sample. As shown in Table 1, the analysis identified varying levels of residual rRNA across the libraries. While the control samples showed relatively low contamination (3.12% and 5.31%), the infested samples exhibited significantly higher levels, reaching up to 15.14% in library SRR17327158. Removing these millions of non-coding reads significantly reduced potential noise for the subsequent *de novo* assembly.
[Summary Report](https://drive.google.com/file/d/1X5rRd3FDAMiVmJXnAeYNFEBVG4VQYPiK/view).
**Table 1. **Alignment statistics of residual ribosomal RNA (rRNA) using SortMeRNA.
### **2. De Novo Assembly & Refinement**
Reconstructing a transcriptome without a reference genome requires grouping raw reads into clusters and resolving them into isoforms. For the *Rorippa indica* project, we followed a pipeline of assembly, redundancy reduction, and biological validation.
### **2.1. RNA-Seq Assembly**
Trinity is considered the gold standard for *de novo* transcriptomics. It reconstructs transcripts using de Bruijn graphs to resolve the underlying isoforms.
#### **Application**
The Trinity tool can be found under Transcriptomics → RNA-Seq Assembly → RNA-Seq Assembly (Trinity). For further information, refer to the OmicsBox User Guide section [RNA-Seq de novo Assembly](https://docs.omicsbox.biobam.com/latest/RNA-Seq-de-novo-Assembly/).
**Input **
- [FASTQ files](https://drive.google.com/drive/folders/1voV4mM83myRBL3P-ofPwAqcFzZ8-KR2H): Preprocessed (cleaned and filtered) FASTQ files.
**Parameters**
- Sequencing Data: Paired-End Reads.
- Reads: Illumina sequencing reads.
Leave the remaining parameters by default.
**Output**
The assembly process produces two primary outputs:
- **Assembled Transcripts (FASTA):** A FASTA file containing the nucleotide sequences of all reconstructed transcripts. The accession names (e.g., TRINITY_DN…) encode the cluster, gene, and isoform information.
[Transcripts File](https://drive.google.com/file/d/1UkKDIgNCZZ9lO2_fY6ljbfcvOKulgQM2/view).
- **Transcript-to-Gene Mapping File:** A tab-delimited file that links each individual isoform to its corresponding gene identifier.
[Transcript-to-Gene Mapping File](https://drive.google.com/file/d/1GWItfxGIF7GP8RUrUTIhIFZxR2QHKy4o/view).
- **Summary Report:** A comprehensive statistical overview that informs about the total number of genes and isoforms and contiguity metrics (Nx statistics). The assembly yielded a total of 56,650 transcripts grouped into 21,713 genes. To evaluate the success of the assembly, we analyzed the length distribution and contiguity statistics (Table 2). The N50 value is a critical metric for assembly quality, representing the length at which at least 50% of the total assembled bases are contained in contigs of that size or longer. For this plant transcriptome, an N50 of 2,003 bp and an average length of 1,831.4 bp are excellent results. These values indicate a highly contiguous assembly, where more than half of the sequence information is contained in transcripts longer than 2 kb.
[ Summary report](https://drive.google.com/file/d/1ivYtn9Uu6E9V4HL2vkq2jKX4KMGDcCKu/view).
**Table 2.** Length distribution and contiguity statistics for the *de novo* transcriptome assembly.
### **2.2. Sequence Clustering (CD-HIT)**
*De novo* assemblers often produce redundant transcripts or fragments. CD-HIT is used to "collapse" these sequences into a representative set of high-quality isoforms.
#### **Application**
The CD-HIT tool can be found under Transcriptomics → RNA-Seq Assembly → Sequence Clustering with CD-HIT. For further information, refer to the OmicsBox User Guide section [Clustering](https://docs.omicsbox.biobam.com/latest/Clustering/).
**Input **
- [FASTA file](https://drive.google.com/file/d/1UkKDIgNCZZ9lO2_fY6ljbfcvOKulgQM2/view): The assembled FASTA file from the Trinity step.
**Parameters**
Provide the assembled FASTA file and leave all parameters by default.
### **2.3. Completeness Assessment (BUSCO)**
To evaluate the biological quality of the assembly, we search for evolutionarily conserved genes expected to be present in any complete transcriptome.
#### **Application**
The BUSCO tool can be found under Transcriptomics → RNA-Seq Assembly → ** **Completeness Assessment with BUSCO.** **For further information, refer to the OmicsBox User Guide section [Completeness Assessment](https://docs.omicsbox.biobam.com/latest/Completeness-Assessment/).
**Input **
- [FASTA file](https://drive.google.com/file/d/1To2KNYSSMbg6f0aknG0BWXzEzhubbPlD/view): The refined transcriptome in FASTA format obtained in the previous step.
**Parameters**
- Input Sequences: The refined transcriptome.
- Lineage: Viridiplantae (kingdom).
Leave the remaining parameters by default.
**Output**
The analysis generates three primary outputs to evaluate the gene space coverage:
- **Main Results Viewer:** A detailed table listing each individual BUSCO ID and its matching sequence in the transcriptome (Figure 7).
[BUSCO Project](https://drive.google.com/file/d/1GGxpORIju0i_6JZvftBSnnpqj5I1iGhH/view?usp=drive_open).
- **Completeness Assessment Report:** A summary providing a quick evaluation and ID lists for each category:
- Complete: Transcripts that match the expected score and length of the BUSCO profile. These can be Single-copy or Duplicated.
- Fragmented: Sequences that meet the score threshold but are significantly shorter than the expected length, indicating incomplete transcripts.
- Missing: Genes that could not be identified, either because they are absent from the dataset or too fragmented to meet the minimum criteria.
[Summary Report](https://drive.google.com/file/d/1ZVUIaesnv3H93T6RbzoCNBKCNq_TjibT/view).
- **Completeness Assessment Summary Chart:** A visual representation (pie chart) of the distribution of orthologs. As shown in Figure 8, the *R. indica* transcriptome achieved an outstanding completeness score of 96.47% (combining single-copy and duplicated orthologs). Only a small fraction of genes was reported as fragmented (1.41%) or missing (2.12%). It is important to note that a significant portion of the "Complete" genes were classified as duplicated. While high duplication levels can be a concern in genome assemblies, they are expected and normal in *de novo* transcriptomes. This occurs because Trinity reconstructs multiple isoforms (splicing variants) for the same genomic locus; BUSCO correctly identifies these multiple versions of the same conserved gene as duplicates. This high overall score confirms that the assembled isoforms capture nearly the entire essential gene space of the plant, ensuring a high-quality dataset for functional annotation.
[Summary Chart](https://drive.google.com/file/d/1-ml03k_giognqQtvlmEEVw-rnLjKGaG1/view).
**Figure 7.** BUSCO results table for the *R. indica* transcriptome.
**Figure 8. **BUSCO completeness assessment showing 96.47% coverage of essential orthologs against the Viridiplantae database.
### **3. Differential Expression Analysis**
After the *de novo* transcriptome is assembled and refined, the next objective is to estimate the expression levels of each transcript across samples and identify which sequences show significant changes in response to aphid infestation.
### **3.1. Transcript-level Quantification (RSEM)**
To estimate the abundance of each transcript across our samples, we employed RSEM (RNA-Seq by Expectation-Maximization). Unlike simple counting methods, RSEM is specifically designed to handle the complexity of *de novo* assemblies. It uses an expectation-maximization approach to accurately allocate multi-mapping reads among isoforms, providing a more reliable quantification of relative abundances.
#### **Application**
The RSEM tool can be found under Transcriptomics → RNA-Seq Read Quantification → Transcript-level Quantification.** **For further information, refer to the OmicsBox User Guide section [Transcript-level Quantification](https://docs.omicsbox.biobam.com/latest/Transcript-level-Quantification/).
**Input **
- [Sequencing Data](https://drive.google.com/drive/folders/1voV4mM83myRBL3P-ofPwAqcFzZ8-KR2H): The original preprocessed FASTQ files generated by the rRNA removal step (1.4).
- [Reference Transcriptome](https://drive.google.com/file/d/1To2KNYSSMbg6f0aknG0BWXzEzhubbPlD/view): The refined FASTA file generated by the clustering step (2.2).
- [Transcript to Gene Map File](https://drive.google.com/file/d/1GWItfxGIF7GP8RUrUTIhIFZxR2QHKy4o/view): The mapping file generated during the Trinity assembly stage (2.1) to enable gene-level estimations.
**Parameters**
- Input Reads: Sequencing Data.
- Input FASTA: Reference Transcriptome.
- Gene-level Estimations: Enabled.
- Transcript to Gene Mapping: Transcript to Gene Map File.
Leave the remaining parameters by default.
**Output**
The analysis generates several outputs that validate the assembly and provide the foundation for statistical testing:
- **Count Tables:** As the "Gene-level Estimations" were enabled, the tool produced two tables: one at the isoform-level and another at the gene-level (where gene counts are the sum of their respective transcripts).
[Count Table at gene-level](https://drive.google.com/file/d/1BeWamLJjb95tj40uyJ7Necp-VqOVndcN/view).
[Count Table at isoform-level](https://drive.google.com/file/d/1QlmxWCahfRi0cqhN-OgoxMaHukXTiQ4F/view?usp=drive_open).
- **Summary Report:** The alignment results (Table 3) show a highly consistent overall alignment rate across all samples, ranging from 70.28% to 72.58%. Approximately 39% of the reads mapped uniquely ("Aligned 1 Time"), while roughly 32% mapped to multiple isoforms ("Aligned Multiple Times"). This significant proportion of multi-mapping is expected in *de novo* transcriptomics, as it reflects the alternative splicing and transcript diversity captured by Trinity.
[Summary Report](https://drive.google.com/file/d/1PrqcTzgm0O_Z1uVKD1vmqS7qgCuoN7BI/view?usp=drive_open).
- **PCA Visualization.** To assess the global structure of the data, a Principal Component Analysis (PCA) was performed (Figure 10). This visualization was generated directly from the isoform-level count table** **via Side Panel → Charts → PCA Plot by providing the [experimental design file](https://drive.google.com/file/d/1VVE9g4oyJOJBcvJ40JwbPDEw4eRBwbln/view). The first principal component (PC1) accounts for 50% of the total variance and perfectly separates the infested samples from the controls. However, the second principal component (PC2) accounts for 28% of the variance and appears to separate the samples by replicate rather than by condition. This suggests that Replicate 1 of both groups might have been processed, sequenced, or handled differently than Replicate 2 (batch effect). Although a thorough review of the original study's metadata provided no explicit information regarding library preparation, it is worth noting that OmicsBox allows for the correction of such technical noise. If this experimental information were available, it could be corrected during the differential expression analysis by setting the "batch" factor as a “blocking factor” in the statistical model.
[PCA Plot](https://drive.google.com/file/d/1VJHEbtOdwYwKVgL8HsM4fDf4pqZJHihT/view).
**Figure 9. **Isoform-level estimation count table.
**Table 3. **Read distribution for the R. indica transcriptome quantification using RSEM.
**Figure 10. **PCA plot based on transcript-level quantification.
### **3.2. Pairwise Differential Expression (edgeR)**
To identify the specific genes driving the resistance in *Rorippa indica*, we performed a statistical analysis using the edgeR package, which is optimized for experiments with a limited number of replicates.
#### **Application**
With the Count Table project open, the edgeR tool can be easily accessed via Side Panel → Actions → Diff. Expression Analysis → Differential Expression Analysis with edgeR. Alternatively, the tool can be found under Transcriptomics → Differential Expression Analysis → Count Table Project/Count Table File → Differential Expression Analysis with edgeR.** **For further information, refer to the OmicsBox User Guide section [Pairwise Differential Expression Analysis with EdgeR](https://docs.omicsbox.biobam.com/latest/Pairwise-Differential-Expression-Analysis/).
**Input **
- [Count Table Project](https://drive.google.com/file/d/1QlmxWCahfRi0cqhN-OgoxMaHukXTiQ4F/view): The isoform-level count table project generated in the previous RSEM step.
- [Experimental Design File](https://drive.google.com/file/d/1VVE9g4oyJOJBcvJ40JwbPDEw4eRBwbln/view): A tab-delimited file associating each sample with its condition (Control or Infested).
**Parameters**
- Count Table Project: The isoform-level count table project.
- Experimental Design: Experimental Design File.
- CPM Filter: 1
- Samples reaching CMP Filter: 2
- Primary Contrast Condition: Infested.
- Primary Reference Condition: Control.
- Select a Statistical Test: Exact Test.
Leave the remaining parameters by default.
**Output**
Several outputs were generated by the analysis:
- **Summary Report:** Provides a comprehensive overview of gene filtering, total library sizes pre-filtering and post-filtering, and the total count of upregulated and downregulated features.
[Summary Report](https://drive.google.com/file/d/1iqME3KUikhbcZp3QlJgUOsKrDtl9pNQZ/view).
- **edgeR Project: **The results table that contains the differential expression statistics (Figure 11). Genes that did not pass the initial filtering step are excluded from this view. Each row corresponds to a genomic feature with the following statistics:
- **logFC:** A measure of the magnitude of expression change between conditions (log2-fold-changes).
- **logCPM:** The average log2-counts-per-million across all samples.
- **LR / F:** Depending on the test, these represent the Likelihood Ratio or the Quasi-likelihood F-statistic.
- **FDR:** False Discovery Rate calculated by the Benjamini-Hochberg method to correct for multiple testing.
- **Tags:** Indicators that label a gene as Upregulated (FDR **≤** 0.05, log2FC ≥ 0) or Downregulated (FDR ≤ 0.05, log2FC ≤ 0) based on the active thresholds.
While OmicsBox establishes a default threshold of False Discovery Rate (FDR) **≤** 0.05 and a |log2FC| ≥ 1, we increased the stringency to |log2FC| ≥ 2 via Side Panel → Actions → Set Up/Down Tags. Under these criteria, edgeR identified 697 Differentially Expressed Genes (DEGs), with 477 upregulated and 220 downregulated in response to aphid infestation.
[Results table](https://drive.google.com/file/d/1CSkOd87suzsd62YS36ll0cZj8JJ6VkNA/view).
- **Volcano Plot (Figure 12):** This plot visualizes the high-confidence DEGs. The wide horizontal spread (log2FC) and high vertical significance (-log10FDR) demonstrate a massive transcriptomic shift triggered by the aphid. This visualization was generated directly from the Results Table via Side Panel → Charts → Volcano Plot.
[Volcano Plot](https://drive.google.com/file/d/1j1nPN1Vcax-Kdtz_XRbBnrEv-ue8HxR0/view).
- **Hierarchical Clustering Heatmap (Figure 13):** This plot standardizes expression values using Z-scores to highlight relative changes. The top dendrogram shows perfect sample clustering by condition, and the distinct "red vs. green" switch confirms the biological impact of the treatment. This visualization was generated directly from the Results Table via Side Panel → Charts → Heatmap.
[Heatmap](https://drive.google.com/file/d/1MNDwdhocMMdWp2kwvImOhniJd09moMI1/view).
**Figure 11. **edgeR pairwise differential expression results table.
**Figure 12. **Volcano plot of the differentially expressed genes obtained by edgeR, highlighting upregulated (red) and downregulated (green).
**Figure 13. **Hierarchical clustering heatmap of differentially expressed genes obtained by edgeR.
### **3.3. Pairwise Differential Expression (DESeq2)**
To complement and validate our findings, we performed a secondary analysis using DESeq2, a statistical package known for its conservative approach to identifying DEGs.
#### **Application**
With the Count Table project open, the edgeR tool can be easily accessed via Side Panel → Actions → Diff. Expression Analysis → Differential Expression Analysis with DESeq2. Alternatively, the tool can be found under Transcriptomics → Differential Expression Analysis → Count Table Project/Count Table File → Differential Expression Analysis with DESeq2. For further information, refer to the OmicsBox User Guide section [Pairwise Differential Expression Analysis with DESeq2](https://docs.omicsbox.biobam.com/latest/Pairwise-Diffexp-with-DESeq2/).
**Input **
- [Count Table Project](https://drive.google.com/file/d/1QlmxWCahfRi0cqhN-OgoxMaHukXTiQ4F/view): The isoform-level count table project generated in the previous RSEM step (3.1).
- [Experimental Design File](https://drive.google.com/file/d/1v_RomriQjkY25K7_QUlcuTKe82LqF9VN/view): A tab-delimited file associating each sample with its condition (Control or Infested).
**Parameters**
- Count Table Project: The isoform-level count table project.
- Experimental Design: Experimental Design File.
- CPM Filter: 1
- Samples reaching CMP Filter: 2
- Primary Contrast Condition: Infested.
- Primary Reference Condition: Control.
- Select a Statistical Test: Wald Significance Test.
Leave the remaining parameters by default.
**Output**
Upon completion, DESeq2 generates a Summary Report and a Results Table containing specific statistical metrics:
- **Summary Report:** A result page that provides a global overview of the analysis.
[Summary Report](https://drive.google.com/file/d/11fwP0VL2K7AgdYc-7Xm8HlVwaTRYwlI6/view).
- **DESeq2 project: **The results table that contains the differential expression statistics. Using the same stringent thresholds as previously (FDR ≤ 0.05 and |log2FC| ≥ 2), DESeq2 identified 558 Differentially Expressed Genes (DEGs), consisting of 402 upregulated and 156 downregulated transcripts. Although slightly more conservative than edgeR (697 DEGs), the core biological signals remain highly consistent.
** ** [Results Table](https://drive.google.com/file/d/1CbLPgrhkKl_kkbOmLjiB1fPimwwRcz46/view).
- **Volcano Plot (Figure 14):** Visualizes the distribution of DEGs, highlighting those with the highest significance and magnitude of change.
[Volcano Plot](https://drive.google.com/file/d/1VkDS0gwuaTG_kMuBojkLXJf_eSBgks-T/view).
- **Hierarchical Clustering Heatmap (Figure 15):** The resulting dendrograms reproduce the distinct separation between infested and control plants, confirming that the biological treatment is the primary driver of the transcriptomic profile.
[Heatmap](https://drive.google.com/file/d/1CI1WBaUBYJMTloplyt9qZ7gziZOTYnlx/view).
**Figure 14. **Volcano plot of the differentially expressed genes obtained by DESeq2, highlighting upregulated (red) and downregulated (green).
**Figure 15. **Hierarchical clustering heatmap of differentially expressed genes obtained by DESeq2.
### **3.4. Venn Diagram**
To ensure the robustness of our results across different statistical models, we used the Venn Diagram tool to visualize the overlap between the DEGs identified by edgeR and DESeq2.
#### **Application**
The Venn Diagram tool can be found under General Tool → Venn Diagram.** **For further information, refer to the OmicsBox User Guide section [Venn Diagram](https://docs.omicsbox.biobam.com/latest/Venn-Diagram/).
**Input **
- [Lists of identifiers](https://drive.google.com/drive/folders/1gahuYUyxImJUXwCYKE_CnfoSKEyS8bCZ): The Venn tool requires a list of identifiers (IDs) in plain text, BOX, or B2G format. To obtain these, we first filtered the results tables from both edgeR and DESeq2 by their** **Tag (Up or Down). Then, using Side Panel → Export → Export Table, we generated .txt files containing the sequence IDs for each category. To ensure the files contained only the IDs, the export was configured as follows:
- Export Table Configuration: In the first page of the wizard, uncheck the Include Header box and check the Select Specific Columns box.
- Select Columns for Export: On the following page, move only the Feature column (which refers to the gene/isoform ID) to the "Selected Columns for Export" list.
**Parameters**
Provide the ID lists of up or down-regulated genes from both differential expression analyses (two separate analyses were performed).
**Output**
- **Venn Diagram. **The comparison revealed a high agreement between both packages: 392 out of the 402 upregulated genes identified by DESeq2 (Figure 16) were also identified by edgeR. For downregulated genes (Figure 17), the overlap was near-perfect (154 out of 156). Using the Proportional check box in the Side Panel, the size of the circles was adjusted to reflect the number of elements in each list.
[Venn Diagram](https://drive.google.com/drive/folders/18Ce0itDhxZPpwHfd8v9-0DrIDrHmou6z).
**Figure 16.** Venn Diagram showing the overlap between upregulated genes identified by edgeR and DESeq2.
**Figure 17. **Venn Diagram showing the overlap between downregulated genes identified by edgeR and DESeq2.
### **4. Enrichment Analysis**
This final phase identifies coordinated shifts in functional categories to interpret the molecular defense strategy of the plant. To gain a multi-angled view, we performed two complementary analyses using the outputs from both edgeR and DESeq2, utilizing annotated sequences from [Functional Analysis](https://help.biobam.com/space/OED0324/3525083363/De+Novo+Transcriptome+Characterization#5--Functional-Annotation).
### **4.1. Fisher’s Exact Test**
Fisher’s Exact Test (FET) is used to identify GO terms that are significantly over-represented or under-represented within discrete lists of differentially expressed genes (DEGs) compared to a reference group.
#### **Application**
The Fisher's Exact Test tool can be accessed directly from the Pairwise Differential Expression results table (either edgeR or DESeq2) by clicking the Fisher's Exact Test button located in the Side Panel → Actions. For further information, refer to the OmicsBox User Guide section [Fisher Exact Test](https://docs.omicsbox.biobam.com/latest/Fisher-Exact-Test/).
**Input **
- Test-set Files: Pairwise Differential Expression results project obtained in the previous [edgeR](https://drive.google.com/file/d/1CSkOd87suzsd62YS36ll0cZj8JJ6VkNA/view) or [DESeq2](https://drive.google.com/file/d/1CbLPgrhkKl_kkbOmLjiB1fPimwwRcz46/view) steps.
- [Reference-set](https://drive.google.com/file/d/1ppG2gOGYi9Cl3yytjxLkmLlIYDX0_-wA/view): The complete set of functionally annotated sequences from the transcriptome project, providing the necessary biological context.
**Parameters**
- Test-set Genes: Subset of significant genes to be considered as the test set. It allows to pick between up-regulated or down-regulated genes. We analyzed upregulated and downregulated genes separately to identify direction-specific biological responses. While ID lists can be imported as text files, the tool also allows selecting "Up-regulated" or "Down-regulated" subsets directly from the DE results table.
- Reference annotation: The pre-annotated transcriptome file.
Leave the remaining parameters by default.
**Output**
- **Enrichment Results Table**: The analysis generates a new table containing the following key metrics:
- Tag: Indicates OVER (more frequent in DEGs) or UNDER (less frequent in DEGs) representation.
- Adj. P-Value: The statistical significance after correction.
- Nr Test vs. Nr Reference: Compares how many genes have the term in our list versus the background.
Following the generation of this table, we applied a GO Redundancy Reduction step using the Side Panel → Actions → Reduce to Most Specific option. This refined our results by collapsing redundant parental terms into their most informative child nodes.
[FET Results](https://drive.google.com/drive/folders/1Z526dmjUGX9407KqDv7MZVnqxls2LHAG).
- **Bubble Plots**: Visualizations were generated via Side Panel → Charts → Bubble Plot. The results from both the edgeR and DESeq2 perspectives show a highly robust defensive signature:
- Upregulated Genes: Both analyses identified a significant over-representation of chitin-based extracellular matrix and chitin binding. The enrichment of these terms suggests the plant is active in recognizing the aphid’s structural components.
- Downregulated Genes: Conversely, suppressed genes were almost exclusively associated with the cytoskeleton and cell cycle progression. Terms such as microtubule-based movement, kinesin complex, and mitotic spindle assembly dominated the results.
[Bubble Plots](https://drive.google.com/drive/folders/1f71JSjSkZW6cR4D-OfX7mVl4YUCqO16X).
**Figure 18. **Fisher's Exact Test results for upregulated genes (DESeq2).
**Figure 19. **Bubble plot of significantly over-represented GO terms among upregulated genes (DESeq2).
**Figure 20. **Bubble plot of significantly over-represented GO terms among downregulated genes (DESeq2).
### **4.2. Gene Set Enrichment Analysis (GSEA)**
Unlike Fisher’s test, GSEA considers the entire expression profile, detecting subtle but coordinated shifts in functional groups that might be missed when looking only at a subset of DEGs. GSEA ranks all isoforms from "most upregulated" to "most downregulated." It then moves down this list, increasing a "running-sum" statistic whenever it encounters a gene annotated to a certain GO term. If the genes of a GO term accumulate at the very top (positive enrichment) or the very bottom (negative enrichment) of the list, the GO term is considered significantly enriched.
**Application**
The GSEA tool is accessed through the Pairwise Differential Expression results table by clicking the Gene Set Enrichment Analysis button in the Side Panel → Actions.** **For further information, refer to the OmicsBox User Guide section [Gene Set Enrichment Analysis (GSEA)](https://docs.omicsbox.biobam.com/latest/Gene-Set-Enrichment-Analysis/).
**Input **
- Pairwise Differential Expression results: Pairwise Differential Expression results project obtained in the previous [edgeR](https://drive.google.com/file/d/1CSkOd87suzsd62YS36ll0cZj8JJ6VkNA/view) or [DESeq2](https://drive.google.com/file/d/1CbLPgrhkKl_kkbOmLjiB1fPimwwRcz46/view) steps.
- [Reference Annotation](https://drive.google.com/file/d/1ppG2gOGYi9Cl3yytjxLkmLlIYDX0_-wA/view): The complete set of functionally annotated sequences from the transcriptome project, providing the necessary biological context.
**Parameters**
- Reference Annotation: The pre-annotated transcriptome file.
- Filter Value: 0.05
Leave the remaining parameters by default.
**Output**
- **Enrichment Results Table:** The analysis generates a new table containing the following key metrics:
- Tags: Indicates if a GO term is enriched at the TOP (upregulated) or BOTTOM (downregulated) of the ranked list.
- NES (Normalized Enrichment Score): The enrichment score after it has been normalized across all analyzed gene sets.
- FDR q-val: The estimated probability that the NES represents a false positive finding.
Following the generation of this table, we applied a GO Redundancy Reduction step using the Side Panel → Actions → Reduce to Most Specific option. This refined our results by collapsing redundant parental terms into their most informative child nodes.
[GSEA Results](https://drive.google.com/drive/folders/17pwLM6ElsmVKFH6-jNRnAMyIkAZx2MWm).
- **Bubble Plots:** Visualizations were generated via Side Panel → Charts → Bubble Plot. The results from both the edgeR and DESeq2 perspectives provided a high-resolution map of the "metabolic sacrifice" occurring in *Rorippa indica* during aphid stress:
- Positive Signals (TOP): The most significant positive signal was the chitin binding, validating that the plant's Pattern-Triggered Immunity (PTI) is highly active. Additionally, the enrichment of indole glucosinolate metabolic processes confirmed a robust chemical defense, as these tryptophan-derived compounds are essential barriers against sucking herbivores.
- Negative Signals (BOTTOM): Significant negative Normalized Enrichment Scores (NES) for microtubule binding, microtubule-based movement, and kinesin complex provided statistical confirmation of a systematic "Growth-Defense Tradeoff." This massive, coordinated downregulation of cellular infrastructure confirms that *R. indica* prioritizes survival over development during biotic stress to optimize fitness.
[Bubble Plots](https://drive.google.com/drive/folders/1UB2efKMKEvsF8FebC_zvHCOpsgkfhjbl).
**Figure 21. **Gene Set Enrichment Analysis results for upregulated genes (DESeq2).
**Figure 22. **GSEA bubble plot for Gene Ontology terms, filtered to show the 15 most enriched categories (top NES) based on DESeq2 differential expression.
**Figure 23. **GSEA bubble plot for Gene Ontology terms, filtered to show the 15 most depleted categories (bottom NES) based on DESeq2 differential expression.
## **CloudUnits Summary**
OmicsBox utilizes [OmicsCloud](https://www.biobam.com/omicscloud-cloud-platform/), an AWS-based high-performance computing (HPC) infrastructure, to execute resource-intensive bioinformatics tasks. These costs are measured in [Cloud Units](https://www.biobam.com/cloud-computation/), which directly correspond to the compute infrastructure used, such as CPU seconds and data transfer volume. Most features within OmicsBox are included for free with an active subscription; only features related to Sequence Alignment and Assembly consume Cloud Units. Different tools need more or less computational resources, depending on the algorithm, input data, and used parameters.
The following table summarizes the resources consumed during this study. Please note that these are approximate values intended to serve as a reference:
| | | | |
| --- | --- | --- | --- |
| **Tool** | **Input Size** | **Execution Time** | **Cloud Units Charged** |
| Trinity *de novo* Assembly | ~8.5 GB (4 paired-end samples) | ~4 hours | ~250,000 |