Whole transcriptome analysis of Arabidopsis thaliana

Overview
Creative Commons License: CC-BY Questions:
  • Which miRNAs are upregulated in response to brassinosteroids?

  • Which genes are potential target of brassinosteroid-induced miRNAs?

Objectives:
  • Perform miRNA differential expression analysis

  • Understand the quasi-mapping-based Salmon method for quantifying the expression of transcripts using RNA-Seq data

  • Idenfity potential miRNAs involved in brassinosteroid-mediated regulation networks

Requirements:
Time estimation: 2 hours
Supporting Materials:
Published: Apr 8, 2021
Last modification: Nov 3, 2023
License: Tutorial Content is licensed under Creative Commons Attribution 4.0 International License. The GTN Framework is licensed under MIT
purl PURL: https://gxy.io/GTN:T00292
rating Rating: 4.0 (0 recent ratings, 4 all time)
version Revision: 90

As sessile organisms, the survival of plants under adverse environmental conditions depends, to a large extent, on their ability to perceive stress stimuli and respond appropriately to counteract the potentially damaging effects. Coordination of phytohormones and reactive oxygen species are considered a key element for enhancing stress resistance, allowing fine-tuning of gene expression in response to environmental changes (Planas-Riverola et al. 2019, Ivashuta et al. 2011). These molecules constitute complex signalling networks, endowing with the ability to respond to a variable natural environment.

Brassinosteroids (BRs) are a group of plant steroid hormones essential for plant growth and development, as well as for controlling abiotic and biotic stress. Structurally, BRs are polyhydroxylated sterol derivatives with close similarity to animal hormones (Figure 1). This group of phytohormones is comprised of around 60 different compounds, of which brassinolide (BL), 24-epibrassinolide (EBR), and 28-homobrassinolide (HBR) are considered the most bioactive.

fig1:Esteroid hormones. Open image in new tab

Figure 1: Structure of various plant and animal steroid hormones.

Several recent studies suggest that the BR-mediated gene regulatory networks have the potential to reshape the future of agriculture, not only by alleviating the antagonistic effect diverse abiotic stress conditios, such as drought, but also by enhancing plant growth and yield. For instance, in tomato (Solanum lycopersicum), EBR treatment enhances drought tolerance, improving photosynthetic capacity, leaf water status, and antioxidant defense (Wang et al. 2018). In pepper (Capsicum annuum), BL treatment increased the efficiency of light utilization under drought (Hu et al. 2013). Gram (Cicer arietinum) plants exposed to drought stress and treated with BL showed significant increases in weight (Anwar et al. 2018). However, the mechanisms of BRs action in enhancing plant tolerance to abiotic stresses still remain largely unknown.

MicroRNAs (miRNAs), mainly 20–22 nucleotide small RNAs (sRNAs), are characterized for regulating gene expression at the post-transcriptional level. miRNAs are distinguished from other sRNAs by being generated from precursor harboring an imperfect stem‐loop structure. Unlike in animals, the pre-processing of plant miRNA occurs in the nucleus (Figure 2). The pre-miRNAs are then exported to the cytoplasm after methylation and incorporated into the Argonaute 1 protein to form RISC (RNA-induced silencing complex). The miRNA itself does not have the ability to cleave mRNAs or interfere with translation, but it plays a role in scanning the appropriate target.

fig2:miRNA biosynthesis. Open image in new tab

Figure 2: Plant miRNA biosynthesis, homeostasis and mechanisms of action (Wang et al. 2018).

miRNAs have been found to be important regulators of many physiological processes, such as stress and hormonal responses. Four factors justify the miRNAs to be considered as master regulators of the plant response to the surrounding environment:

  • Multiple miRNA genes are regulated under given environmental conditions
  • Computational predictions estimate that each miRNA regulates hundreds of genes
  • The majority of plant miRNAs regulate genes encoding for transcription factors (TFs)
  • Targets include not only mRNAs but also long noncoding RNAs (lncRNAs)

In plants, miRNAs can silence targets through RNA degradation as well as translational repression pathways, and unlike animals, a large proportion of miRNA and their targets have less than four mismatches. This feature has been exploited for developing miRNAs target prediction tools, providing an efficient approach to elucidate the miRNA-mediated regulatory networks, which can contribute to biotechnological solutions to improve crops productivity.

In this tutorial, inspired by Park et al. 2020, we aim to explore the interplay between brassinosteroids and the miRNA-gene silencing pathway, considered one of the most versatile regulatory mechanisms in response to stressful situations in plants.

Agenda

In this tutorial, we will cover:

  1. Experimental design
  2. Background on data
  3. Get data
  4. miRNA data analysis
    1. Quality assessment of miRNA reads
    2. miRNA quantification: MiRDeep2
    3. Differential expression analysis of miRNAs: DESeq2
    4. Filter significantly differentially expressed miRNAs
  5. mRNA data analysis
    1. Quality assessment of mRNA reads
    2. Quantification of gene expression: Salmon
    3. Differential expression analysis of mRNAs: DESeq2
    4. Filter significantly differentially expressed mRNAs
  6. Identification of miRNA targets
    1. miRNA target prediction using TargetFinder
  7. Optional exercise
  8. Conclusion

Experimental design

The main objective of this training is to identify potential targets of miRNAs whose expression is induced by brassinosteroids. Our starting hypothesis is that there must be brassinosteroid-induced miRNAs that have high sequence complementarity with mRNAs whose expression is inhibited in the presence of these hormones (Figure 3).

fig3:Experimental design. Open image in new tab

Figure 3: Experimental design

Background on data

The datasets to be used in this training can be classified into three groups: miRNA reads, mRNA reads and additional datasets.

miRNA reads

The miRNA datasets consist of six FASTQ files, obtained by using the Illumina GAxII sequencing platform. The plant samples were obtained from wild-type Ws-2 seedlings treated with mock or 1 μM EBR for 90 min before harvest. The original datasets are available in the NCBI SRA database, with the accession number SRP258575. As in the previous case, for this tutorial, we will use a reduced version of the data.

mRNA reads

The mRNA datasets consist of four FASTQ files, generated through the Illumina HiSeq 2000 sequencing system. The samples were obtained from wild-type Columbia (Col-0) seedlings treated with mock or 100 nM BL for 4 hours. The original datasets are available in the NCBI SRA database, with the accession number SRP032274. For this tutorial, subsets from the original data were generated in order to reduce the analysis run time.

Additional datasets

In addition to the RNA-Seq reads obtained from the NCBI database, we will use datasets from two sources:

  • AtRTD2 is a high-quality transcript reference dataset developed to exploit the accuracy of transcript quantification tools such as Salmon and Kallisto in analyzing Arabidopsis RNA-Seq data.
  • PmiREN is a comprehensive functional plant miRNA database that includes more than 20,000 annotated miRNAs diverse plant species.

Get data

The first step of our analysis consists of retrieving the miRNA-Seq datasets from Zenodo and organizing them into collections.

Hands-on: Retrieve miRNA-Seq and mRNA-Seq datasets
  1. Create a new history for this tutorial
  2. Import the files from Zenodo:

    • Open the file galaxy-upload upload menu
    • Click on Rule-based tab
    • “Upload data as”: Collection(s)
    • Copy the following tabular data, paste it into the textbox and press Build

      SRR11611349	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611349_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
      SRR11611350	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611350_MIRNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
      SRR11611351	Control miRNA	https://zenodo.org/record/4710649/files/SRR11611351.MIRNASEQ_CTLfastqsanger.gz	fastqsanger.gz
      SRR11611352	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611352_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
      SRR11611353	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611353_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
      SRR11611354	BR treated miRNA	https://zenodo.org/record/4710649/files/SRR11611354_MIRNASEQ_BL.fastqsanger.gz	fastqsanger.gz
      SRR1019436	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019436_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
      SRR1019437	Control mRNA	https://zenodo.org/record/4710649/files/SRR1019437_RNASEQ_CTL.fastqsanger.gz	fastqsanger.gz
      SRR1019438	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019438_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz
      SRR1019439	BR treated mRNA	https://zenodo.org/record/4710649/files/SRR1019439_RNASEQ_BL.fastqsanger.gz	fastqsanger.gz
      
    • From Rules menu select Add / Modify Column Definitions
      • Click Add Definition button and select List Identifier(s): column A

        Then you’ve chosen to upload as a ‘dataset’ and not a ‘collection’. Close the upload menu, and restart the process, making sure you check Upload data as: Collection(s)

      • Click Add Definition button and select Collection Name: column B
      • Click Add Definition button and select URL: column C
      • Click Add Definition button and select Type: column D
    • Click Apply and press Upload
  1. Click on the dataset to expand it
  2. Click on Add Tags galaxy-tags
  3. Add a tag starting with #
    • Tags starting with # will be automatically propagated to the outputs of tools using this dataset.
  4. Press Enter
  5. Check that the tag appears below the dataset name

Next we will retrieve the remaining datasets.

Hands-on: Retrieve the additional datasets
  1. Import the files from Zenodo:

    • Open the file galaxy-upload upload menu
    • “Upload data as”: Datasets
    • Once again, copy the tabular data, paste it into the textbox and press Build

      annotation_AtRTD2.gtf	https://zenodo.org/record/4710649/files/annotation_AtRTD2_19April2016.gtf.gz
      transcriptome.fasta	https://zenodo.org/record/4710649/files/transcriptome_AtRTD2_12April2016.fasta.gz
      star_miRNA_seq.fasta	https://zenodo.org/record/4710649/files/star_miRNA_seq.fasta
      mature_miRNA_AT.fasta	https://zenodo.org/record/4710649/files/mature_miRNA_AT.fasta
      miRNA_stem-loop_seq.fasta	https://zenodo.org/record/4710649/files/miRNA_stem-loop_seq.fasta
      
    • From Rules menu select Add / Modify Column Definitions
      • Click Add Definition button and select Name: column A
      • Click Add Definition button and select URL: column B
    • Click Apply and press Upload

As indicated above, for this tutorial the depth of the samples was reduced in order to speed up the time needed to carry out the analysis. This was done as follows:

Hands-on: Dataset subsampling
  1. Sub-sample sequences files ( Galaxy version 0.2.5) with the following parameters:
    • param-files “Multiple datasets”: Each of the fastq files
    • “Subsampling approach”: Take every N-th sequence (or pair e.g. every fifth sequence)
    • “N”: 100

In this way, we will only take 1% of reads at a random sampling rate.

miRNA data analysis

Once we have imported the data, we can begin to study how brassinosteroid exposure alters gene expression patterns.

Quality assessment of miRNA reads

Due to technical limitations, sequencing is considered an error-prone process. In Illumina sequencing platforms, substitution type miscalls are the dominant source of errors, which can lead to inconsistent results. Another factor that can interfere with our analyses is the presence of adapter contaminants, which can result in an increased number of unaligned reads, since the adapter sequences are synthetic and do not occur in the genomic sequence.

Sequence quality control is therefore an essential first step in your analysis. We will use two popular tools for evaluating the quality of our raw reads: FastQC and MultiQC.

Comment

In order to visualize the data from both collections together in the MultiQC tool, it will be necessary to combine the results generated by FastQC. For more information on the topic of quality control, please see our dedicated training materials.

Hands-on: Initial quality check
  1. FastQC ( Galaxy version 0.72+galaxy1) with the following parameters:
    • param-collection “Dataset collection”: Control miRNA
  2. Rename the outputs as FastQC unprocessed control miRNA: RawData and FastQC unprocessed control miRNA: Webpage
  3. Repeat the previous steps with the BR treated miRNA collection.
  4. Merge Collections with the following parameters:
    • In “Input collections”:
      • “1.Input Collections”: FastQC unprocessed control miRNA: RawData
      • “2.Input Collections”: FastQC unprocessed BR treated mRiNA: RawData
  5. MultiQC ( Galaxy version 1.8+galaxy1) with the following parameters:
    • In “Results”:
      • “Which tool was used generate logs?”: FastQC
      • param-collection “Dataset collection”: select the output generated in the previous step.
    • In “Report title”: miRNA initial quality check
  6. Click on the galaxy-eye (eye) icon and inspect the generated HTML file
Question

Based on the information provided by MultiQC, is it necessary to trimming/filtering the reads?

The report generated by MultiQC indicates that three quality parameters show values outside the recommended limits: per sequence G/C content, overrepresented sequences and adapter content.

fig4:GC content. Open image in new tab

Figure 4: Per sequence GC content of the miRNA samples.
fig5:Overrepresented sequences. Open image in new tab

Figure 5: Overrepresented sequences in miRNA samples.
fig6:Adaptor content. Open image in new tab

Figure 6: Adaptor content of miRNA samples.

Particularly remarkable is the content of adapters, which are up to 80% of the reads in some samples. Given the abundance of adapters, one would expect that if we eliminate this source of contamination, the rest of the parameters would show a noticeable improvement.

To remove the adapters contamination, we will employ the Trim Galore tool, a wrapper script around Cutadapt and FastQC which automates quality and adapter trimming.

Hands-on: Trimming of adapter sequences
  1. Trim Galore! ( Galaxy version 0.4.3.1) with the following parameters:
    • “Is this library paired- or single-end?”: Single-end
      • param-collection “Reads in FASTQ format”: Control miRNA
      • “Adapter sequence to be trimmed”: Illumina small RNA adapters
  2. Rename the output collection as Control miRNA trimmed
  3. Repeat the previous steps with the BR treated miRNA collection.

Next, we will reassess the quality of the sequences to check if the adapters have been removed.

Hands-on: Post-processing quality check
  1. FastQC ( Galaxy version 0.72+galaxy1) with the following parameters:
    • param-collection “Dataset collection”: Control miRNA trimmed
  2. Rename the outputs as FastQC processed control miRNA: RawData and FastQC processed control miRNA: Webpage
  3. Repeat the previous steps with the BR treated miRNA trimmed collection.
  4. Merge Collections with the following parameters:
    • In “Input collections”:
      • “1.Input Collections”: FastQC processed control miRNA: RawData
      • “2.Input Collections”: FastQC processed BR treated miRNA: RawData
  5. MultiQC ( Galaxy version 1.8+galaxy1) with the following parameters:
    • In “Results”:
      • “Which tool was used generate logs?”: FastQC
      • param-collection “Dataset collection”: select the output generated in the previous step.
    • In “Report title”: miRNA trimmed quality check
  6. Click on the galaxy-eye (eye) icon and inspect the generated HTML file

The evaluation of the report generated by MultiQC after having processed the samples by Trim Galore indicates that the G/C content has been successfully corrected, and that the adapter contamination has been eliminated. However, the samples still show a high degree of over-represented sequences (Figure 7).

fig7:Overexpressed sequences. Open image in new tab

Figure 7: Report of overexpressed sequences in the miRNA reads
Question

What do you think could be the cause of the high number of over-represented sequences?

Two of the factors that may determine the abundance of overrepresented sequences are the existence of highly expressed miRNAs (Seco-Cervera et al. 2018), the existence of highly conserved sequence motifs within the miRNA (Glazov et al. 2008), and contamination with foreign sequences.

miRNA quantification: MiRDeep2

Quantification of miRNAs requires to use two different tools:

  • The MiRDeep2 Mapper tool for preprocessing the reads.
  • The MiRDeep2 Quantifier tool for mapping the deep sequencing reads to predefined miRNA precursors and determining the expression of the corresponding miRNAs. It is carried out in two steps: firstly, the predefined mature miRNA sequences are mapped to the predefined precursors (optionally, predefined star sequences can be mapped to the precursors too). And second, the deep sequencing reads are mapped to the precursors.
Hands-on: Quantification of miRNAs
  1. MiRDeep2 Mapper ( Galaxy version 2.0.0) with the following parameters:
    • param-collection “Deep sequencing reads”: Control miRNA trimmed
    • “Remove reads with non-standard nucleotides”: Yes
    • “Clip 3’ Adapter Sequence”: Don't Clip
    • “Collapse reads and/or Map”: Collapse
  2. Rename the output as Collapsed control miRNA
  3. Repeat the previous stages by providing BR treated miRNA trimmedas input, and rename it as Collapsed BR treated miRNA.

  4. MiRDeep2 Quantifier ( Galaxy version 2.0.0) with the following parameters:
    • param-collection “Collapsed deep sequencing reads”: Collapsed control miRNA
    • param-file “Precursor sequences”: miRNA_stem-loop_seq.fasta
    • param-file “Mature miRNA sequences”: mature_miRNA_AT.fasta
    • param-file “Star sequences”: star_miRNA_seq.fasta
  5. Rename the outputs as MiRDeep2 control miRNA and MiRDeep2 control miRNA (html).
  6. Repeat the fourth step by providing Collapsed BR treated miRNA as input, and rename the outputs as MiRDeep2 BR treated miRNA and MiRDeep2 BR treated miRNA (html)

To use the outputs generated by MiRDeep2 Quantifier in the differential expression analysis, it is necessary to modify the datasets.

Hands-on: Edition of MiRDeep2 Quantifier outputs
  1. Cut columns from a table with the following parameters:
    • “Cut columns”: c1,c2
    • “Delimited by”: Tab
    • param-collection “From”: MiRDeep2 control miRNA
  2. Rename the output as control miRNA counts
  3. Cut columns from a table with the following parameters:
    • “Cut columns”: c1,c2
    • “Delimited by”: Tab
    • param-collection “From”: MiRDeep2 BR treated miRNA
  4. Rename the output as BR treated miRNA counts

Differential expression analysis of miRNAs: DESeq2

DESeq2 is a tool for differential gene expression analysis based on a negative binomial generalized linear model. DESeq2 internally corrects the differences in library size, due to which no preliminary normalization of input datasets is required.

Comment

It is desirable to use at least three replicates of each experimental condition to ensure sufficient statistical power.

Hands-on: Differential expression analysis using DESeq2
  1. DESeq2 ( Galaxy version 2.11.40.6+galaxy1) with the following parameters:
    • “How”: Select datasets per level
      • In “Factor”:
        • param-repeat “Insert Factor”
          • “Specify a factor name, e.g. effects_drug_x or cancer_markers”: effects_brassinolide
          • In “Factor level”:
            • param-repeat “Insert Factor level”
              • “Specify a factor level”: brassinolide
              • param-collection “Counts file(s)”: BR treated miRNA counts
            • param-repeat “Insert Factor level”
              • “Specify a factor level”: control
              • param-collection “Counts file(s)”: control miRNA counts
    • “Files have header?”: No
    • “Choice of Input data”: Count data (e.g. from HTSeq-count, featureCounts or StringTie)
  2. Rename the outputs as DESeq2 results miRNA and DESeq2 plots miRNA
  3. Click on the galaxy-eye (eye) icon and inspect the DESeq2 plots miRNA file

DESeq2 generated 2 outputs: a table with the normalized counts and a graphical summary of the results. To evaluate the similarity of our samples, we are going to inspect the Principal Component Analysis (PCA) plot. PCA allows evaluating the dominant directions of the highest variability in the data. Thus, the samples subjected to the same conditions should cluster together.

fig8:PCA miRNA. Open image in new tab

Figure 8: PCA plot of expression data from control and BR treated miRNA samples.

As can be seen, the main axes account for only 47% and 19% of the total variation. This suggests that the effect of brassinosteroids on miRNA regulation is limited (Figure 8).

Filter significantly differentially expressed miRNAs

Next, we will extract those genes whose expression is statistically significantly differentially expressed (DE) due to BR treatment by selecting those whose adjusted p-value is less than or 0.05. A cut-off value of 0.05 indicates that the probability of a false positive result is less than 5%.

The p-value is a measure of the probability that an observed difference could have occurred just by random chance. A small p-value indicates that there is a small chance of getting this data if no real difference existed. A p-value threshold of 0.05 indicates that there is a 5% chance that the result is a false positive. The p-adj (also known as q-value) is an adjusted p-value which taking in to account the false discovery rate (FDR). Applying a FDR becomes necessary when we’re measuring thousands of variables from a small sample set.

Hands-on: Filter differentially expressed miRNAs
  1. Filter data on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: DESeq2 results miRNA
    • “With following condition”: c7<0.05
Question

How many miRNAs show statistically significant differential expression?

Unfortunately, we have not detected any differentially expressed miRNAs. In this case, this is caused by the fact that the subsampled datasets don’t have sufficient read depth to test for differential expression.

To get any sensible results, it is worth analyzing the full dataset. You can analyze the full datasets following the above tutorial. Otherwise, you can import the DESeq2 analysis results that we generated from full dataset into your history.

Hands-on: Retrieve the DESeq2 analysis results on full miRNA dataset
  1. Import the files from Zenodo:

    • Open the file galaxy-upload upload menu
    • Click of the Paste/Fetch button
    • Copy the Zenodo links and press Start

      https://zenodo.org/record/4663389/files/miRNA_DESeq2_results_complete_dataset.tabular
      
  2. Rename each dataset according to the sample id (e.g. miRNA_DESeq2_results_complete_dataset)
  3. Add all miRNA data analysis related tags: #miRNA #BR #control

Repeat the filtering step on the imported DESeq2 result.

Hands-on: Filter and sort differentially expressed miRNAs from the full dataset
  1. Filter data on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: miRNA_DESeq2_results_complete_dataset
    • “With following condition”: c7<0.05
  2. Rename the output as Differentially expressed miRNAs
  3. Filter data on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: Differentially expressed miRNAs
    • “With following condition”: c3>0
  4. Rename the output as Upregulated miRNAs
  5. Sort data in ascending or descending order (Galaxy Version 1.1.0) with the following parameters:
    • param-file “Sort Dataset”: Upregulated miRNAs
    • “on column”: Column: 3
    • “everything in”: Descending order
Question
  1. How many miRNAs are differentially expressed?
  2. How many miRNAs show statistically significant upregulation and what is the most upregulated miRNA?
  1. Out of 442 miRNAs, 39 show significant differential expression.
  2. There are 16 upregulated miRNAs and Ath-miR156g is the most upregulated one.

mRNA data analysis

Once the differential expression analysis of miRNAs has been carried out, the next stage is to analyze how mRNA expression is altered in response to brassinosteroids.

Quality assessment of mRNA reads

As in the previous section, we shall begin by assessing the quality of our sequences.

Hands-on: Initial quality check
  1. FastQC ( Galaxy version 0.72+galaxy1) with the following parameters:
    • param-collection “Dataset collection”: Control mRNA
  2. Rename the outputs as FastQC unprocessed control mRNA: RawData and FastQC unprocessed control mRNA: Webpage
  3. Repeat the previous step with the BR treated mRNA collection.
  4. Merge Collections with the following parameters:
    • In “Input collections”:
      • “1.Input Collections”: FastQC unprocessed control mRNA: RawData
      • “2.Input Collections”: FastQC unprocessed BR treated mRNA: RawData
  5. MultiQC ( Galaxy version 1.8+galaxy1) with the following parameters:
    • In “Results”:
      • “Which tool was used generate logs?”: FastQC
      • param-collection “Dataset collection”: select the output generated in the previous step.
    • “Report title”: mRNA initial quality check
  6. Click on the galaxy-eye (eye) icon and inspect the generated HTML file
Question

Are there any stats that indicate the need to process the samples in order to improve their quality?

All the stats are within acceptable limits. However, the adapter content shows the presence of Illumina universal adapters in our reads, which can be removed to avoid possible interferences at later stages (Figure 10).

fig10:Adapter content of mRNA samples. Open image in new tab

Figure 9: Quality assessment of mRNA samples
Comment

Although our samples have adapters, we are not going to trim them in this case. We will explain the reason in the next section.

Quantification of gene expression: Salmon

After performing the quality assessment of the reads, we can move on to quantifying the gene expression. The aim of this step is to identify which transcript each read comes from and the total number of reads associated with each transcript. In this tutorial, we will use the Salmon tool (Patro et al. 2017) for the quantification of mRNA transcripts.

fig11:Salmon logo.

One of the characteristics of Salmon is that it doesn’t require performing a base-to-base alignment, which is the time-consuming step of tools such as STAR and HISAT2. Salmon relies on the quasi-mapping concept, a new mapping technique that allows the rapid and accurate mapping of RNA-Seq reads to a target transcriptome. Rather than a standard alignment, quasi-mapping seeks to find the best mappings for each read, and does so by finding minimal collections of dynamically sized, right-maximal, matching contexts between target and query positions (Srivastava et al. 2016)

The quasi-mapping approach by Salmon requires a reference index to determine the position and orientation information for accurate mapping prior to quantification. It allows providing the transcriptome in a format that optimizes its use in transcript identification and quantification.

After determining the best mapping for each read, Salmon generates the final transcript abundance estimation after modeling sample-specific parameters and biases. Reads that map equally well to more than one transcript will have the count divided between all of the mappings, thus avoiding the loss of information on the different gene isoforms.

The quasi-mapping algorithm makes use of two main data structures, the generalized suffix array (SA) of the transcriptome T, and a hash table (h) that maps each k-mer in T to its SA interval (by default k = 31). During the quasi-mapping procedure, a read is scanned from left to right until a k-mer (ki, starting at position i on the read) is encountered that appears in h. The ki is looked up in the hash table and the SA intervals are retrieved, giving all suffixes containing the k-mer ki. Then, the Maximal Mappable Prefix (MMPi) is computed by finding the longest substring of the read that matches the reference suffixes. Owing to sequencing errors, the MMPs may not span the complete read. In this case, the next informative position (NIP) is determined based on the longest common prefix of the SA intervals of MMPi. Suffix array search continues from k-bases before NIP until the final NIP is less than k-bases before the read end. Finally, for each read, the algorithm reports the transcripts it mapped to, location and strand information (Srivastava et al. 2016).

fig12:Quasi-mapping algorithm. Open image in new tab

Figure 10: Illustration of the quasi-mapping of a read using k=3. Hash table lookup of ki returns the suffix array interval [b, e). The base 'G' at position 6 in the read is a result of sequencing error. Hence, the MMPi is 'ATTGA' and the SA interval of MMPi is [b', e'). The next k-mer starts at k-bases prior to the NIP which is base after the longest common prefix of the interval [b',e'). In the end, the read in the above example most likely map to the suffix array at [e'-1, e').
Hands-on: Quantify gene expression with Salmon
  1. Salmon quant ( Galaxy version 0.14.1.2+galaxy1) with the following parameters:
    • “Select salmon quantification mode:”: Reads
      • “Select a reference transcriptome from your history or use a built-in index?”: Use one from the history
        • In “Salmon index”:
          • param-file “Transcripts fasta file”: transcriptome.fasta
      • In “Data input”:
        • “Is this library mate-paired?”: Single-end
          • param-collection “FASTQ/FASTA file”: Control mRNA
      • “Validate mappings”: Yes
    • param-file “File containing a mapping of transcripts to genes”: annotation_AtRTD2.gtf
  2. Rename the outputs as Salmon control mRNA (Quantification) and Salmon control mRNA (Gene Quantification)
  3. Repeat the previous procedure by using the BR treated mRNA dataset
Comment: Quasi-mapping sequence requirements

Trimming the reads is not required when using this method, since if there are k-mers in the reads that are not in the hash table, they are not counted. Quantification of the reads is only as good as the quality of the reference transcriptome.

Salmon generates two outputs for each input sample:

  • Quantification: summarizes the quantifications by transcript
  • Gene quantification: summarizes the quantification by gene

Each output consists of a tabular dataset with five columns:

Field Description
Name The name of the target transcript provided in the input transcriptome.
Length The length of the target transcript.
EffectiveLength The computed effective length.
TPM The relative abundance of this transcript in units of Transcripts Per Million.
NumReads The number of reads mapping to each transcript that was quantified.
Question

Could you explain why we did not trim the reads before?

The reason is that, since the k-mers which contain the adaptor sequences are not present in the transcriptome (from which the hash table is generated), they are omitted.

Differential expression analysis of mRNAs: DESeq2

Now, let’s analyze which genes show statistically significant differential expression in response to brassinosteoids.

Hands-on: Differential expression analysis using DEseq2
  1. DESeq2 ( Galaxy version 2.11.40.6+galaxy1) with the following parameters:
    • “How”: Select datasets per level
      • In “Factor”:
        • param-repeat “Insert Factor”
          • “Specify a factor name, e.g. effects_drug_x or cancer_markers”: effects_brassinolide
          • In “Factor level”:
            • param-repeat “Insert Factor level”
              • “Specify a factor level”: brassinolide
              • param-collection “Counts file(s)”: Salmon BR treated mRNA (Gene Quantification)
            • param-repeat “Insert Factor level”
              • “Specify a factor level”: control
              • param-collection “Counts file(s)”: Salmon control mRNA (Gene Quantification)
    • “Choice of Input data”: TPM values (e.g. from kallisto, sailfish or salmon)
      • “Program used to generate TPMs”: Salmon
      • “Gene mapping format”: GTF/GFF3
        • param-file “GTF/GFF3 annotation file”: annotation_AtRTD2.gtf
  2. Rename the outputs as DESeq2 results mRNA and DESeq2 plots mRNA
fig13:PCA plot of mRNA samples. Open image in new tab

Figure 11: PCA plot of differential expression data from control and BR treated samples.
Question

What conclusions about the similarity of the samples can be derived from the PCA plot?

From the information provided by the plot, it is possible to state that there is a high similarity between the samples belonging to the same experimental conditions, since the first dimension (x-axis) allows to explain 81% of the variability, and the samples are located at opposite ends of the x-axis.

Filter significantly differentially expressed mRNAs

To conclude the analysis of the differential expression of mRNAs, we will extract those transcripts that show a significant induction of expression in response to brassinosteroids. Before continuing with further analysis, similar to miRNA data analysis, import the DESeq2 results generated from full mRNA datasets.

Hands-on: Retrieve the DESeq2 analysis results on full mRNA dataset
  1. Import the files from Zenodo:

    • Open the file galaxy-upload upload menu
    • Click of the Paste/Fetch button
    • Copy the Zenodo links and press Start

      https://zenodo.org/record/4663389/files/mRNA_DESeq2_results_complete_dataset.tabular
      
  2. Rename each dataset according to the sample id (e.g. mRNA_DESeq2_results_complete_dataset)
  3. Add all mRNA data analysis related tags: #mRNA #BR #control

Now we continue with the DE genes analysis.

Hands-on: Extract significantly DE genes
  1. Filter data on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: mRNA_DESeq2_results_complete_dataset
    • “With following condition”: c7<0.05
  2. Rename the output as Differentially expressed mRNAs

  3. Filterdata on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: Differentially expressed mRNAs
    • “With following condition”: c3>1
  4. Rename it as Upregulated mRNAs

  5. Filter data on any column using simple expressions (Galaxy Version 1.1.1) with the following parameters:
    • param-file “Filter”: Differentially expressed mRNAs
    • “With following condition”: c3<-1
  6. Rename it as Downregulated mRNAs
Question
  1. How many genes show statistically significant differential expression?
  2. How many of them are significantly upregulated (at least two-fold change)? And downregulated?
  3. What is the most significantly DE downregulated gene and what is it biological function?
  1. The total number of genes whose expression differential between the two experimental conditions is 4176.
  2. Of them, 328 are significantly downregulated by the BR treatment and 778 are upregulated.
  3. The most significantly DE gene is AT3G30180 (BR60X2). BR60X2 encodes a cytochrome p450 enzyme that catalyzes the last reaction in the production of brassinolide. It is capable of converting 6-deoxocastasterone into castasterone, a C-6 oxidation, as well as the further conversion of castasterone into brassinolide (source: TAIR database).

Identification of miRNA targets

To predict which miRNAs target which mRNAs, first we need their transcriptomic sequences in FASTA format. Now we will obtain the sequences of miRNAs whose expression is induced by brassinosteroids.

Hands-on: Obtaining the sequences of upregulated miRNAs
Comment

The following tools can be found in Text Manipulation and Filter and Sort sections.

  1. Cut columns from a table with the following parameters:
    • “Cut columns”: c1
    • “Delimited by”: Tab
    • param-collection “From”: Upregulated miRNAs
  2. Rename the output as Upregulated miRNA ids
  3. Filter FASTA ( Galaxy version 2.1) with the following parameters:
    • param-file “FASTA sequences”: mature_miRNA_AT.fasta
    • “Criteria for filtering on the headers”: List of IDs
      • param-file “List of IDs to extract sequences for”: Upregulated miRNA ids
      • “Match IDs by”: Default: ID is expected at the beginning: >ID
  4. Filter FASTA ( Galaxy version 2.1) with the following parameters:
    • param-file “FASTA sequences”: star_miRNA_seq.fasta
    • “Criteria for filtering on the headers”: List of IDs
      • param-file “List of IDs to extract sequences for”: Upregulated miRNA ids
      • “Match IDs by”: Default: ID is expected at the beginning: >ID
  5. Concatenate datasets tail-to-head (Galaxy Version 1.0.0) with the following parameters:
    • param-file “Concatenate Dataset”: output of Filter FASTA tool on mature_miRNA_AT.fasta
      • Insert Dataset
        • param-file “Select”: output of Filter FASTA tool on star_miRNA_seq.fasta
  6. Rename it as Upregulated miRNA sequences
  7. Click on the galaxy-eye (eye) icon and inspect the Upregulated miRNA sequences file

To identify putative targets of upregulated miRNAs, it is necessary to obtain the sequences of all downregulated mRNAs in FASTA format.

Hands-on: Obtaining the gene sequences of downregulated mRNAs
  1. Cut columns from a table with the following parameters:
    • “Cut columns”: c1
    • “Delimited by”: Tab
    • param-collection “From”: Downregulated mRNAs
  2. Rename the output as Downregulated mRNA ids
  3. Filter FASTA ( Galaxy version 2.1) with the following parameters:
    • param-file “FASTA sequences”: transcriptome.fasta (Input dataset)
    • “Criteria for filtering on the headers”: List of IDs
      • param-file *“List of IDs to extract sequences for: Downregulated mRNA ids
      • “Match IDs by”: Custom regex pattern
        • “Regex search pattern for ID”: GENE=(AT.{7})
  4. Rename it as Downregulated mRNA sequences
  5. Click on the galaxy-eye (eye) icon and inspect the Downregulated mRNA sequences file

miRNA target prediction using TargetFinder

We are now ready to launch the search for miRNA target genes. For this we will use the TargetFinder tool that is implemented and used for miRNA target prediction in plants (Srivastava et al. 2014, Ma et al. 2017).

Hands-on: identification of potential miRNA targets
  1. TargetFinder ( Galaxy version 1.7.0+galaxy1) with the following parameters:
    • param-file “Input small RNA sequences file”: Upreguled miRNA sequences
    • param-file “Target sequence database file”: Downregulated mRNA sequences
    • “Prediction score cutoff value”: 5.0
    • “Output format”: Tab-delimited format
  2. Click on the galaxy-eye (eye) icon and inspect the output of TargetFinder.

Congratulations! We have identified the following 5 potential genes involved in the brassinosteroid-miRNA regulatory network.

fig15:TargetFinder results. Open image in new tab

Figure 12: TargetFinder results.

Finally, we can access all the information available on the genes identified in the TAIR database:

Both AT4G14365 and AT1G26890 are not well characterized genes. In the case of AT5G50570, experimental research have demonstrated that this gene is involved in flooding tolerance in Medicago sativa, through a signaling path mediated by miR156 (Feyissa et al. 2021). On the other hand, according to Gao et al. 2017, SPL13 is likely a negative regulator of plant vegetative growth. The interaction between miR156 and SPL transcription factors has been suggested for diverse Poaceae family members (Yue et al. 2021).

Comment: Hypothesis testing

One of the hypotheses that we could propose from our results is: the inhibition of the AT2G46850 gene can result in plants with improved resistance to drought conditions. Is it possible to validate it? Yes! We propose this approach: to acquire AT2G46850 mutant seeds and wild type seeds, grow them under two controlled conditions: watered and drought stress, and analyze plant weight after 33 days (Figure 13).

fig16:Plant growth. Open image in new tab

Figure 13: Arabidopsis growth conditions protocol (de Ollas et al. 2019).

Optional exercise

As additional activity, you can try to repeat the workflow by using the sequences stored in the NCBI GEO database with the accession number GSE119382. In that case, we will compare gene expression patterns of mutants overexpressing the brassinosteroid receptor BRL3 under two experimental conditions: control and drought-stress. The required datasets are available in the data library:

Hands-on: Import data from the Data Libraries
  1. Go into Shared data (top panel) and click on Data Libraries
  2. In the search box enter the following identifier: 4710649
  3. Select the following files:
    https://zenodo.org/record/4710649/files/SRR7779222_BRL3_mRNA_drought.fastqsanger.gz
    https://zenodo.org/record/4710649/files/SRR7779223_BRL3_mRNA_drought.fastqsanger.gz
    https://zenodo.org/record/4710649/files/SRR7779224_BRL3_mRNA_drought.fastqsanger.gz
    
  4. Click on Export to History and as a Collection
  5. In the pop-up window press Continue
  6. Provide it the name BRL3 mRNA drought and push Create list
  7. Repeat the previous procedure with the remaining files:
    https://zenodo.org/record/4710649/files/SRR7779228_BRL3_mRNA_watered.fastqsanger.gz
    https://zenodo.org/record/4710649/files/SRR7779229_BRL3_mRNA_watered.fastqsanger.gz
    
  8. Finally provide it the name BRL3 mRNA control and push Create list

We will use the upregulated miRNAs obtained in the previous analysis in order to identify potential targets.

Upregulated miRNA   https://zenodo.org/record/4710649/files/upregulated_miRNA_BR_complete_dataset.fasta

Conclusion

In this tutorial, we have analyzed RNA sequencing data to extract information about potential genes regulated by brassinosteroids. For this purpose, the approach chosen was the identification of genes complementary to miRNAs upregulated in response by brassinosteroids. The final result has been the identification of five potential miRNA targets.