Introduction

Because of the plummeting in the cost of sequencing technologies during the last decade, RNA massive sequencing (RNA-seq) has become mainstream to study the transcriptome [1]. Currently, short-read sequencing technologies, typically producing outputs of 30 million reads per sample, are the most extensively used methodologies for gene expression profiling [2]. This pace of data generation requires computational processing to produce interpretable results. Thus, the use of pipelines to perform the different steps of transcriptomic data processing have become a widespread practice. The core of these is usually composed by spliced aligners as STAR [3], HISAT2 [4] or Rail-RNA [5], which map reads against a reference genome, or by pseudo-alignment tools as Salmon [6] or Kallisto [7], that directly obtain a quantification for the regions of interest using probabilistic models. Additionally, there are pipelines which are intended to be run by the user in local computers or high-performance environments, as QuickRNASeq [8], or interactively in cloud-based platforms, after uploading raw data to an external service, as BioJupies [9] or RaNA-Seq [10]. Typically, the interpretation of the experiment involves differential expression analysis, carried out using count based or linear models, with packages as edgeR [11], DESeq2 [12] or limma [13], followed by methods, such as over representation analysis [14] or the gene set enrichment analysis [15], to extract functional information from the obtained results.

Despite different pipelines to perform the aforementioned tasks are available (Tables 1 and 2), most of them present two major drawbacks. First, the genomic information contained in the RNA-Seq reads usually remains unused. However, genomic variants, which may contain crucial information about the functionality and potential activity of the resulting proteins in the different processes where they participate, can be retrieved from such sequences. In this sense, it is well known that RNA-Seq has some limitations for DNA variant calling. There are two main points to consider: (i) lowly expressed genes include lower depth, so variant calling is harder in those regions and (ii) the detection of heterozygous variants can be limited due to allele-specific gene expression [16]. Despite these limitations, it has been demonstrated that variants can be called even for low expressed genes in deeper RNA-Seq sequence samples. Moreover, some studies have shown that RNA-Seq variant calling is able to provide a good sensitivity of 99.7%-99.8% in both heterozygous and homozygous variants whereas precision still reaches 97.6% in homozygous but 90% in heterozygous [17]. The second major drawback is that conventional functional analysis strategies are mainly descriptive, and very limited in providing biological insights of the underlying molecular mechanisms that produce the observed phenotypic responses. Recently, a new generation of methods, known as mechanistic pathway analyses, are outperforming traditional approaches in both biological explanatory power and interpretability [18]. Here we present MIGNON (Mechanistic InteGrative aNalysis Of rNa-seq), a complete and versatile workflow able to exploit all the information contained in RNA-Seq data and producing not only the conventional normalized gene expression matrix, but also an annotated VCF file per sample with the corresponding mutational profile. Moreover, MIGNON can combine both files to model signaling pathway activities through an integrative functional analysis using the mechanistic modeling algorithm Hipathia [19]. Signaling circuit outputs can further be easily linked to phenotypic features (e.g. disease outcome, drug response, etc.) [19–21]. Mechanistic modeling has been successfully applied to understand disease mechanisms in rare diseases [22,23], complex diseases [21], and, especially in cancer [19,24–26]. Specifically, the hiPathia algorithm has demonstrated to have a superior sensitivity and specificity than other similar algorithms available [27].

Table 1

Features of the workflows for RNA-seq data analysis.

Workflow	URL	Google scholar citations	Year	Implementation	Read pre-processing	(Pseudo) Alignment	Variant calling and annotation	Differential gene expression	Functional analysis	Omic integration
QuickRNASeq	https://sourceforge.net/projects/quickrnaseq/	26	2016	Shell, Perl and R scripts	-	STAR	VarScan2	-	-	No
SePIA	http://anduril.org/sepia	25	2016	Anduril workflow	FastX-Toolkit Trimmomatic TrimGalore	STAR TopHat Bowtie Bowtie2	Bambino ANNOVAR	Cuffdiff DESeq DESeq2 EdgeR	SPIA	No
Recount2	https://jhubiostatistics.shinyapps.io/recount/	154	2017	Shiny app and R package	-	Rail-RNA	-	-	-	No
RNACocktail	https://bioinform.github.io/rnacocktail/	71	2017	Python scripts	-	HISAT2	GATK	DESeq2	-	No
BioJupies (ARCHS4)	https://amp.pharm.mssm.edu/biojupies/	143	2018	Web service	-	Kallisto	-	Limma Characteristic direction	Enricher	No
GREIN	https://shiny.ilincs.org/grein	9	2019	Shiny app and R package	Trimmomatic	Salmon	-	EdgeR	-	No
VaP	https://modupeore.github.io/VAP/	1	2019	Perl scripts	AfterQC Trimmomatic	TopHat2 HISAT2 STAR	GATK ANNOVAR VEP	-	-	No
DEWE	http://www.sing-group.org/dewe/	3	2019	Java app	Trimmomatic	Bowtie2 HISAT2	-	Ballgown EdgeR	PathfindR	No
RaNa-Seq	https://ranaseq.eu/	0	2019	Web service	Fastp	Salmon	-	DESeq2	GOseq fgsea	No
MIGNON	https://github.com/babelomics/MIGNON	-	2020	WDL workflow	Fastp	Salmon STAR HISAT2	GATK VEP	EdgeR	hiPathia	Yes

Table 2

Analysis outputs of the workflows for RNA-seq data analysis.

Workflow	Normalized gene expression	Differential gene expression	Transcriptomic-based functional results	Genomic variants	Annotated genomic variants	Integrated Transcriptomic + genomic functional results
QuickRNASeq	Yes	-	-	Yes	-	-
SePIA	Yes	Yes	Yes	Yes	Yes	-
Recount2	Yes	-	-	-	-	-
RNACocktail	Yes	Yes	-	Yes	-	-
BioJupies (ARCHS4)	Yes	Yes	Yes	-	-	-
GREIN	Yes	Yes	-	-	-	-
VaP	-	-	-	Yes	Yes	-
DEWE	Yes	Yes	Yes	-	-	-
RaNa-Seq	Yes	Yes	Yes	-	-	-
MIGNON	Yes	Yes	Yes	Yes	Yes	Yes

Design and implementation

Quality control and alignment

First, using raw sequences as the input for the workflow, fastp (v0.20.0) [28] is applied to perform the quality trimming and filtering of reads using the default values for windows size and required mean quality and length. Then, FastQC (v0.11.5) can be used to create a quality report for each pre-processed read file. After the quality control step, five modes for the execution of the workflow can be selected (see Table 3). Each execution mode uses a different combination of “core” tools to perform the alignment or pseudo-alignment of pre-processed reads, as explained in the tool documentation (see also Fig 1). In brief, all of them make use of a combination of STAR (v2.7.2b), HISAT2 (v2.1.0), Salmon (v0.13.0) and FeatureCounts (v1.6.4) [29] to align (or pseudo-align) reads against a reference genome (or transcriptome) and subsequently obtain the counts per gene matrix. The hisat2 and star modes use a conventional counting strategy, employing FeatureCounts to summarize the number of sequences overlapping the genomic regions of interest (genes), as specified by a genome annotation file. On the other hand, the core component of salmon-hisat2, salmon-star and salmon consist of the pseudo-aligner Salmon, which directly obtains transcript level quantification using a probabilistic model. Note that in the salmon-hisat2 and salmon-star modes, the execution of STAR or HISAT2 is still necessary to generate the alignment files that feed the variant calling sub-workflow.

Fig 1

MIGNON Workflow.

Directed graph summarizing the tools employed by the workflow (blue boxes) and the strategy used by MIGNON to integrate genomic and transcriptomic information into signaling circuits. Gene expression and LoF variants are obtained from reads and integrated by doing an in-silico knockdown of genes that present a LoF variant. Then, this combined matrix is used as the input for hiPathia, that estimates the signaling circuits activation status by using expression values as proxies for protein signaling activities.

Table 3

MIGNON execution modes.

Execution mode	Alignment	Quantification	Variant calling	Computational profile
salmon-hisat2	HISAT2	Salmon	Yes	Low memory consumption. Slower than STAR.
salmon-star	STAR	Salmon	Yes	High memory consumption. Faster than HISAT2.
hisat2	HISAT2	featureCounts	Yes	Low memory consumption. Slower than STAR.
star	STAR	featureCounts	Yes	High memory consumption. Faster than HISAT2.
salmon	-	Salmon	No	Low memory consumption and fast.

Results

MIGNON integrative approach for the mechanistic interpretation of multi-omic information into a pathway framework

MIGNON is the first pipeline able to extract genomic and transcriptomic information from RNA-seq data and integrate them within a mechanistic framework. The ultimate protein activity is assessed from the transcriptional activity conditioned to the integrity of the gene. No matter its level of expression, a gene that harbors a deleterious mutation is in-silico knocked-down by the model to simulate the loss of function (Fig 1). To evaluate how the proposed strategy affects the predicted signaling circuit activities, two different runs of MIGNON were carried out over 462 unrelated human lymphoblastoid cell line samples from the 1000 Genomes sample collection, corresponding to the CEU, FIN, GBR, TSI and YRI populations [37]. In the reference run, only transcriptomic information (raw) was used, while in the case example run the knock-down strategy was applied. Fig 2A and 2B clearly depicts how the knock-down due to LoF mutations interrupts the transduction of the signal in three circuit/sample pairs. Moreover, Fig 2C shows that the overall predicted signaling circuit activities are significantly lower (paired Wilcoxon signed-rank test P value < 2.2x10^-16) when the genomic information is integrated in the model. This example clearly shows how the use of transcriptomics data alone produced an incomplete picture of the real signaling activity and proves the usefulness of multi-omic data integration.

Fig 2

In silico knock-downs effect on predicted signaling circuit activities.

A) Network representation of three signaling circuits that contain genes with loss of function variants for three subjects from the 1000 genomes cohort. The node color indicates whether a gene contained in it has a loss of function variant (yellow) or not (black). Red and blue arrows indicate stimulations and inhibitions, respectively. B) Predicted signaling activity for three circuit/sample pairs on the sub-figure A. Color represents signaling circuit activity with and without considering the genomic information. C) Violin plots showing all the predicted signaling circuit activities with and without the genomic information for the 1000 genomes cohort (paired Wilcoxon signed-rank test P value < 2.2x10^-16).

Workflow performance evaluation

To assess MIGNON performance and resource consumption, the workflow was executed over 6 different human datasets (S2 Table), comprising a total of 42 samples. It was tested with cromwell (v47) and singularity (v3.5), using 6 different CPU configurations on tasks allowing multi-threading. This analysis revealed that the slower components of the workflow are the aligners (HISAT2 and STAR) and the MarkDuplicates and HaplotypeCaller steps of the GATK sub-workflow. Fig 3 summarizes the time and memory consumption of the tools which allow multi-threading using 6 different CPU configurations. While HISAT2 is slower than STAR, the second one makes a more intensive use of available memory. Therefore, both aligners are available in MIGNON since this tradeoff should be considered if planning to deploy the workflow in cloud-computing based environments or, contrarily, in limited memory computing environments. Additionally, Fig 4 shows the time and memory consumption of the different steps that compose the GATK sub-workflow. Here, MarkDuplicates displays the highest memory consumption and HaplotypeCaller shows the longest runtime. Overall, the different tasks carried out by the workflow show a maximum memory usage under the 32 gigabytes, which makes the pipeline deployable under most computational environments. Finally, and due to the WDL, cromwell and docker combination, the workflow is something fast and easy to deploy and setup.

Fig 3

MIGNON performance results.

Multi-thread tasks. A) Memory consumption by task. Each boxplot represents the maximum memory consumption in Gigabytes (y axis) for each CPU configuration (X axis) and each multi-thread task (facets). Dashed lines indicate the following memory configurations: 4, 8, 16 and 32 gigabytes. B) Elapsed time by task. Each boxplot represents the elapsed time (Y axis) for each CPU configuration (X axis) and each task (facets). Dashed lines indicate time points: 30, 60, 120 and 240 minutes.

Fig 4

GATK sub-workflow performance results.

A) Memory consumption by task. Each boxplot represents the maximum memory consumption in Gigabytes (Y axis) for each task (X axis). Dashed lines indicate the following memory configurations: 4, 8, 16 and 32 gigabytes. B) Elapsed time by task. Each boxplot represents the elapsed time (Y axis) for each task (X axis). Dashed lines indicate the following time points: 30, 60, 120 and 240 minutes.

Conclusions

In summary, MIGNON represents an innovative concept of RNA-Seq data analysis that automates the sequence of steps that leads from the uninformative raw reads to the ultimate sophisticated functional interpretation of the experiment, providing, for the first time, a user-friendly framework for integration of genomic and transcriptomic data.

MIGNON makes use of several popular methods to perform the initial processing of reads and utilize the HiPathia mathematical model to provide a mechanistic interpretation of the experiment in the context of human signaling. MIGNON has an enormous application potential in personalized medicine, especially in the analysis of cancer transcriptomes, given its ability to interpret putative driver mutations along with gene expression in the context of signaling activity, a process highly relevant in tumorigenesis.

MIGNON can be easily deployed in different computer environments making an optimal use of the resources. Additionally, the modularity with which the workflow has been designed makes its upgrade and maintenance a straightforward task.