Fused inverse-normal method for integrated differential expression analysis of RNA-seq data

Background Use of next-generation sequencing technologies to transcriptomics (RNA-seq) for gene expression profiling has found widespread application in studying different biological conditions including cancers. However, RNA-seq experiments are still small sample size experiments due to the cost. Recently, an increased focus has been on meta-analysis methods for integrated differential expression analysis for exploration of potential biomarkers. In this study, we propose a p-value combination method for meta-analysis of multiple independent but related RNA-seq studies that accounts for sample size of a study and direction of expression of genes in individual studies. Results The proposed method generalizes the inverse-normal method without an increase in statistical or computational complexity and does not pre- or post-hoc filter genes that have conflicting direction of expression in different studies. Thus, the proposed method, as compared to the inverse-normal, has better potential for the discovery of differentially expressed genes (DEGs) with potentially conflicting differential signals from multiple studies related to disease. We demonstrated the use of the proposed method in detection of biologically relevant DEGs in glioblastoma (GBM), the most aggressive brain cancer. Our approach notably enabled the identification of over-expressed tumour suppressor gene RAD51 in GBM compared to healthy controls, which has recently been shown to be a target for inhibition to enhance radiosensitivity of GBM cells during treatment. Pathway analysis identified multiple aberrant GBM related pathways as well as novel regulators such as TCF7L2 and MAPT as important upstream regulators in GBM. Conclusions The proposed meta-analysis method generalizes the existing inverse-normal method by providing a way to establish differential expression status for genes with conflicting direction of expression in individual RNA-seq studies. Hence, leading to further exploration of them as potential biomarkers for the disease. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-022-04859-9.


Background
RNA sequencing (RNA-seq) technologies are now increasingly considered for whole transcriptome gene expression quantification studies as compared to traditional microarray technologies due to its high technical reproducibility and greater resolution [1].Over the last decade, it has found widespread application in studying different biological conditions including cancers.For instance, sequencing data archived on The Cancer Genome Atlas (https://portal.gdc.cancer.gov/)have been used in a number of studies to explore potential biomarkers and mechanisms in oncogenesis [2] [3].Despite its advantages and few large RNAseq datasets [4] [5], RNA-seq experiments are still small sample size experiments because of its high cost.This leads to a problem of reduced statistical power in studies such as differential expression analysis where thousands of genes are studied at a time but only have tens to hundreds of samples.Combination of data or results from multiple independent but related studies (referred to as meta-analysis) have been widely used to increase available sample size and consequently the statistical power to obtain a precise estimate of gene expression differentials [6] [7].In the context of differential expression analysis, several different metaanalysis approaches have been proposed for integrating microarray studies [8] [9] and some of them have later been adapted for RNA-seq data [10] [11].
For microarray gene expression studies, apart from vote-counting and direct merging of datasets, meta-analysis methods can mainly be classified into three types based on the combined statistic [7].First are methods based on effect sizes combination in which a combined effect (for instance, strength of differential expression between two conditions for a gene) is obtained based on the calculated effect sizes and its variance where two possible models namely, fixed and random effects model are used to obtain the combined effect [12].Second are approaches based on integration of p-values obtained from per-study analysis into a single combined p-value per gene [13].Lastly, are approaches based on rank combination which are non-parametric and allow for integration of studies based on a statistic that can be ordered, e.g., fold change of a gene [14].However, RNA-seq data are counts data, i.e., normalized number of sequenced reads within a certain gene or transcript, unlike the microarray data which are continuous, e.g., normalized signal intensity of image [15].Hence, the methods initially proposed for microarray data are not suited to be applied directly to RNAseq data in many cases [10].
In case of RNA-seq data, Poisson or Negative-Binomial distributions are typically used to model gene counts [16].Kulinskaya et al. (2008) [17] described an effect-size combination method using an Anscombe transformation of Poisson distributed data.However, as highlighted by Rau et al. (2014) [10], this effect-size combination approach is not appropriate for RNA-seq data due to over-dispersion among biological replicates and presence of zeroinflation.Rau et al. (2014) considered two p-value combination methods, namely Fisher and inverse normal (IN) or Stouffer's methods, previously proposed and used for meta-analysis of microarray studies [8][9] [13] and demonstrated how these can be adapted in RNA-seq data analysis.Their results illustrated that Fisher and IN methods were very similar to each other in terms of performance but were better than the global and per-study differential analysis [10].These two (Fisher and IN) p-value combination approaches have been implemented in several R packages, e.g., metaRNASeq [10], metaseqR [18] and metaSeq [19] and are most widely used methods due to its statistical simplicity and ease of direct application for meta-analysis of RNA-seq studies for differential expression.Among all the existing meta-analysis methods for RNA-seq data discussed above, only few of the p-value combination methods (e.g., IN and PANDORA p-value [18]) allow for incorporation of information regarding the number of replicates in different studies to be combined through specification of a set of weights.However, information related to the direction of expression (up-or down-regulated) of a gene across different studies is not accounted for or included in any of these meta-analysis methods for RNA-seq data.Underand over-expressed genes are analysed together and genes exhibiting conflicting direction of expression across studies are either removed prior to meta-analysis or are suggested to be identified and removed post-hoc [10] [11].Hence, no conclusion can be drawn with regards to differential expression for the genes that have conflicting direction of expression across different studies.Given that a significant proportion of genes may exhibit conflicting direction of expression across different gene expression studies [20], particularly when more and more RNA-seq data are publicly available and included into integration, emphasis is warranted on including this important prior information in a meta-analysis setting.
Recently, importance of inclusion of direction of expression information for genes in RNA-seq meta-analysis has been recognized leading to a generalization of some existing p-value combination methods such as Fisher method and Bayesian Hierarchical Model [21] [22] [23].However, these generalizations come at a cost of increased statistical and computational complexity which discourages their widespread application to transcriptomic studies.In this study, we aimed to develop a new approach for integrated differential meta-analysis of RNAseq data which accounts for both the sample size and direction of gene regulation in each study.The proposed approach leads to a generalization of the IN method without introducing additional cost and hence is simple and intuitive for real data application.First, we propose a modified inverse-normal (MIN) approach for p-value combination and assess its performance by comparing it with the IN method based on an extensive simulation study.Next, to overcome the conservative nature of MIN method, we further propose a fused inverse normal (FIN) method for p-value combination and assess its performance by comparing it to IN and MIN methods in a simulation study.Then an application to a set of real glioblastoma (GBM, the most aggressive type of brain cancer) studies has been conducted.Moreover, we assessed the relevance of the identified differentially expressed genes (DEGs) for GBM by using Ingenuity Pathway Analysis (IPA, www.qiagen.com/ingenuity)for pathway analysis and upstream regulator analysis.

Methods
Let   be the observed count for gene  ( = 1, 2, … , ) in condition  ( = 1, 2) of biological replicate  ( = 1, 2, … ,   ) in study  ( = 1, 2, … , ).For an integrated differential analysis of gene expression across multiple studies, we first conducted the differential expression analysis within a given study  using edgeR package (version 3.26.5) in R version 3.6.0[24] with likelihood ratio test as the test for differential expression.Let   be the raw p-value for per-gene and per-study obtained using the individual differential expression analysis within a given study  for gene .The null hypothesis tested in the individual differential analysis is that the gene is non-differentially expressed in the particular study.For notational convenience, the notations similar to the ones used in Rau et al. (2014) [10] were adopted in this study.

Modified inverse-normal method
Let   be a Bernoulli random variable which takes values 1 and -1 when a gene  is over and under expressed respectively in a study .A gene can be assessed as over-or under-expressed based on the fold change values (>1 or <1) of the gene in a study.Then, for a gene , we define a combined statistic (1) where   are a set of study specific weights described by Marot and Mayer [25] as follows: Here, ∑    is the total number of biological replicates in a study  for all condition c and ∑ ∑     indicates the total number of biological replicates in all studies.Moreover,   can be considered as a weighted z-score.An advantage of this weighting criteria is that larger weights are attributed to studies with larger sample sizes.Φ is the standard normal cumulative distribution function and   is the raw p-value obtained for gene  by differential analysis for study .
It is assumed that   are uniformly distributed under the null hypothesis and Φ −1 (1 −   ) is standard normal in the above formula (1), but this assumption of   is not automatically satisfied when dealing with RNA-seq data [10].However, filtering of very low expressed genes in each study results in p-values which are roughly uniformly distributed under the null hypothesis [10].Then, we have that   |Φ −1 (1 −   )| ~ (0, 1) (see Theorem 1).
Theorem 1: Let  and  be two independent random variables where ~(0, 1) and  is a Bernoulli random variable taking values 1 and −1.Then,  = || is standard normal distributed.Proof: Using the first principle, Now, if  < 0, the RHS of (3) becomes For  ≥ 0, the RHS of (3) becomes . Hence, by symmetry of the normal distribution, we have Thus,  ∼ (0, 1).■ Hence,   in equation ( 1) is a linear combination of independent standard normal variables.Thus, is also standard normal.A two-sided test can then be performed with  0 being that the gene  is not differentially expressed between two conditions (case vs control) and combined p-value is given by (  = ℙ(|| ≥   ), i.e.

Simulation study: MIN and IN comparison
To investigate the performance and compare the MIN method with state-of-the-art p-value combination method (IN with post-hoc filtering), we performed a simulation study.An extensive set of RNA-seq data was generated using the negative binomial distribution for the counts   and method described in Rau et al. (2014) [10] (see Supplementary 1, section: Simulation study model).Parameters for the simulation study were estimated from a real RNAseq dataset for Alzheimer's disease (AD) study downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/)[27] with accession number GSE125583 which contains data for 219 AD patients and 70 normal control subjects.The method used for estimation of mean and dispersion parameters from GSE125583 were as described in Rau et al. (2014) with BH p-value < 0.05 being used to classify a gene as a DEG.This dataset was chosen as it has considerable number of samples for both biological conditions, namely case and control.Simulation settings for inter-study variability parameter (), number of samples per condition and number of studies have been detailed in Table 1.The inter-study variability parameter represents the amount of variability between the studies considered for metaanalysis.In practice, the observed variability between human data studies is considerable (~0.5)[10].We chose two different values of  (0.15 and 0.5) to represent small and large amount of inter-study variability respectively.For each setting described in  For each simulation setting, individual p-values obtained from differential expression analysis using edgeR were combined using both MIN and IN methods.A gene was considered differentially expressed if the BH adjusted combined p-value (FDR) < 0.05.Based on area under the ROC curves (AUC) (Table 1, Figure 1), both meta-analysis methods performed well in terms of detection power in identifying DEGs under all simulation settings.For both low ( = 0.15) and high ( = 0.5) inter-study variability we observed that the MIN method was more conservative (AUC was smaller than IN) than the IN method.However, as the inter-study variability and the number of studies to be combined increased, both meta-analysis methods were found to have comparable performance (Figure 1c-1d).Although slightly conservative in its performance with respect to the IN approach, MIN method has the advantage of using direction of expression information leading to identification of DEGs among genes with conflicting direction of expression across studies.The conservative behaviour of the MIN method can be attributed to the fact that a two-sided hypothesis testing is performed as compared to a one-sided test on right-hand tail of the distribution in case of IN method.Hence, next we proposed the FIN method as a mixture of IN and MIN methods to circumvent the issue of conservativeness of MIN method.As expected, with increase in inter-study variability and number of studies to be combined, the number of genes with mismatched direction of expression was significantly higher (see Supplementary 1: Table 1).We also note that the FDR for all simulation settings was controlled well below 5% threshold (Figure 2a).In terms of uniquely identified DEGs by the MIN method as compared to IN method, the proportion of true positives (TP) was higher than 80% (Figure 2b) in all simulation settings.A large proportion of TPs among the unique DEGs identified by the MIN method indicates that the MIN approach can lead to DEGs that are biologically relevant to a disease in a real application.Moreover, as the inter-study variability, number of studies or both increased, there was an increase in the number of uniquely identified DEGs by the MIN method and proportion of TPs among them (Figure 2b).More importantly, a high percentage of these unique DEGs (>75% in all settings) were observed to have the true direction of expression (Figure 2c) suggesting that a significantly high percentage of uniquely identified DEGs by the MIN method in real data applications will have true direction of expression as their effective direction of expression.The effective observed direction of expression was determined by the sign of   as defined in equation (1).

Fused inverse-normal method
To address the conservative nature of MIN method, we propose a mixture method which is a mixture of IN and MIN method for integrated differential analysis.In contrast to formula (1) we define   as follows: , if  has same direction of expression across s Here,   , Φ and   have their usual meaning as described previously.As   follows a standard normal distribution given the assumption that   is uniformly distributed under the null hypothesis, a one-sided test on the right-hand tail of the distribution (as proposed in [10]) can be performed for genes with same direction of expression across studies.For the genes with conflicting direction of expression across studies, a two-sided test can be performed. 0 being the same in both the cases.Multiple testing correction to control the overall FDR can then be carried out using the BH method.A detailed interpretation of the FIN model in terms of differential expression of a gene and its direction of expression in individual studies can be found in Supplementary 1.

Simulation study: FIN, IN and MIN comparison
In addition to the simulation study for comparing MIN with IN method, we assess and compare the performance of FIN method to that of IN and MIN methods by using the same simulated data and settings described in Table 1.Based on AUC (Table 1, Figure 1), FIN performed similar or better than IN method and had better performance than MIN under all simulation settings.As with MIN, FIN method also has the advantage of using direction of expression information and hence identified DEGs among genes with conflicting direction of expression in contrast with IN method.More importantly, we observed that FIN significantly improved detection power for true DE genes with concordant differential expression patterns across studies as compared to MIN method and does not lead to increased number of false positive detections overall (Figure 3a).As compared to IN, the proportion of TPs among the uniquely identified DEGs by FIN method was higher than 90% (Figure 3b) indicating that FIN method can lead to DEGs that are biologically relevant to a disease in a real application.Similar to MIN, as the inter-study variability, number of studies or both increased, there was an increase in the number of uniquely identified DEGs by the FIN method as compared to IN method and proportion of TPs among them (Figure 3b).In addition, a high percentage of these unique DEGs (>80% in all settings) were observed to have the true direction of expression (Figure 3c) suggesting that a significantly high percentage of uniquely identified DEGs by the FIN method in real data applications will have true direction of expression as their effective direction of expression.
The effective observed direction of expression was determined by the sign of   for genes with conflicting direction of expression across studies.In case of same direction of expression of a gene across studies, the consistent direction of expression was kept as the effective direction of expression.

Application to brain cancer data
To demonstrate how the MIN and FIN method can be adapted in practice for differential metaanalysis of RNA-seq data and compare it with IN method, an application to real GBM studies has been conducted.

Data description
GBM RNA-seq datasets were searched in GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA databases (https://portal.gdc.cancer.gov/).Datasets were selected based on a selection criterion that at least 3 GBM patients and 3 normal tissue samples are available for analysis.Three different GBM RNA-seq datasets, two from GEO (with accession ID: GSE123892 and GSE151352) and one from TCGA (TCGA-GBM) matched our selection criteria and were considered for analysis (for details, see Table 2).Raw gene or transcript counts data (where available) was directly downloaded for TCGA-GBM and GSE123892 datasets.For GSE151352, raw FASTQ files were downloaded and processed using Galaxy web platform via the European UseGalaxy server (https://usegalaxy.eu/)[28] to obtain raw counts.The quality of the raw reads was assessed (using FastQC) and the specified adapter sequence ATCACCGACTGCCCATAGAGAGGCTGAGAC was removed with Cutadapt (version 1.16) [29].The parameters used for this step were the parameters provided by the submitter of the dataset on GEO.The adapter trimmed reads were aligned to the reference genome (GRCh37.p13)using sequence aligner RNA STAR (Galaxy version 2.7.5b) [30] where other parameters used were default settings.Following alignment, the generated BAM files were processed using the featureCounts tool (Galaxy version 1.6.4+galaxy2) to get raw counts for each RNA-seq data sample.More details of the processing pipeline used for GSE151352 can be found in Supplementary 1: Figure 1.  2. Information about GBM RNA-seq datasets used for integrated analysis using different p-value combination methods in our study.Up and down DEGs refer to the up and downregulated DEGs obtained in per-study differential analysis.

Per-study differential expression analysis
Each of these datasets were processed separately for quality control and differential expression analysis using edgeR package (version 3.26.5) in R. Raw counts data (transcript) were annotated by mapping Ensembl IDs to Entrez Gene IDs and gene symbols (org.Hs.eg.db package, version 3.8.2 in R [31]).Ensembl IDs with no Entrez ID mapping were filtered out.For those with multiple matchings, the one with highest aggregated count was selected.Counts per million (CPM) threshold was carefully selected to reduce the number of low expressed transcripts.We removed a transcript if number of samples equal to or more than the minimum number of samples in a condition had less than 0.85 CPM for that transcript.Although subjective, this choice of threshold seems to work well for the uniform distribution assumption for the p-values under  0 .Only genes left after low expression filtering were considered for individual differential expression analysis in order to satisfy the uniformity assumption on pvalues under the null hypothesis (Figure 4a).The remaining transcripts were then normalized using the TMM method [32].Common and tag-wise dispersion were estimated and a negative-binomial generalized log-linear model was fitted to the read counts using the glmFit function under the edgeR package.Raw p-values were then obtained from the differential analysis for case/control conditions.More information about the DEGs identified in per-study differential analysis based on the criteria |log 2 | > 1 and FDR p-value < 0.05 can be found in Table 2.We note that TCGA-GBM dataset has a much larger library size (~ 47 million reads, Illumina HiSeq 2000 v2 sequencer) as compared to GSE151352 (~ 4 million reads, Ion Torrent S5 sequencer) and GSE123892 (~ 35 million reads, Illumina HiSeq 2500 sequencer).Hence, we observed a differing number of genes left after filtering and consequently a much larger number of DEGs being observed for TCGA-GBM dataset as compared to the other two (Table 2) [33] in per-study differential expression analysis.As the sequencing output gets larger, the smaller count differences between samples are declared significant by models for differential expression in edgeR.A more detailed treatment of differential expression in RNA-seq data and how it is affected by sequencing depth and other factors can be found in Tarazona et al. (2011) [33].
Moreover, we also considered individual differential analysis for TCGA-GBM RNA-seq data by randomly selecting 20 cases together with available 5 normal samples in order to make all three datasets (GSE123892, GSE151352 and TCGA-GBM) comparable in terms of number of replicates for the meta-analysis (see Supplementary 1: Table 2).Hence, we considered two different meta-analysis scenarios.

Meta-analysis
Once the raw p-values were obtained from the individual differential expression analysis for each dataset, IN, MIN and FIN methods were then applied for p-value combination.Since the TCGA-GBM dataset (160 GBM vs 5 normal samples) is much larger in terms of number of samples as compared to GSE123892 (4 GBM vs 3 normal samples) and GSE151352 (12 GBM vs 12 normal samples), we considered two different combination scenarios.First, all TCGA-GBM samples were used for individual analysis to obtain the raw p-values.Second, 20 cases and 5 normal samples randomly selected from TCGA-GBM dataset were considered for individual analysis to get raw p-values and then considered for meta-analysis with the other two datasets (GSE123892 and GSE151352). 10 different random selections were made and individual differential expression analysis were conducted respectively.Second scenario ensured that the datasets included in meta-analysis had comparable sample sizes.
In scenario one, a total of 18325 unique gene pool was considered for meta-analysis which was the combination of genes identified in each RNA-seq data analysis after quality control and filtering (Table 2).13056 out of 18325 genes (~71%) were found to have the same direction of expression across the studies in which they were present whereas 5259 (~29%) of genes had conflicting or mismatched direction of expression.The direction of expression for a gene in an individual study was determined based on the sign of log 2  obtained for that gene in perstudy differential analysis.Hence, only 13056 genes were effectively considered for IN method as compared to 18325 genes for MIN and FIN methods for identifying DEGs because of posthoc removal of DEGs with conflicting direction of expression in the IN method.
For each of the combination methods, we assessed the number of DEGs based on average absolute log fold-change ∑ |log 2   |/  =1 > 1 and FDR p-value < 0.05 criteria.Here,  denotes the number of datasets in which a particular gene was present.In case a gene was absent in a dataset, the weights in the combination methods were estimated only using the number of replicates in datasets in which the gene was present.The three p-value combination methods were then compared based on number of DEGs identified and unique DEGs identified by each method.A total of 5918, 5892 and 6138 DEGs were identified by the IN, MIN and FIN methods respectively.Of the DEGs detected by all these meta-analysis methods, more than 90% of them were in common (Figure 4b) with FIN method having a higher detection power than the other two methods.Moreover, MIN and FIN method identified a total of 233 and 235 DEGs with mismatched direction of expression across studies by incorporating the direction of expression information.More importantly, in the subset of DEGs which were present in all three datasets, 5.26% of DEGs had conflicting direction of expression across studies.Although, small in proportion, this would be of importance in case a gene of interest for the disease being studied has conflicting direction of expression across different studies.Particularly when more datasets are included in meta-analysis, the number of genes considered in IN approach can be massively reduced.
Given that the FIN method has the highest power of DEG detection, we further explore the DEGs obtained using this meta-analysis procedure.Top 10 up and down-regulated DEGs identified by FIN method are presented in Table 3.For full list of DEGs identified by different meta-analysis methods, see Supplementary 2, 3 and 4. In terms of effective direction of expression of DEGs, 2914 DEGs with same direction of expression across studies and 180 DEGs with mismatched direction of expression were up-regulated.Similarly, 2989 (same direction) and 55 (mismatched direction) DEGs were down-regulated.Results for scenario 2 for random selection have been detailed in Supplementary 1: Table 3

Pathway analysis and biological significance
DEGs obtained by the FIN method were further explored to assess their biological relevance to GBM.QIAGEN's Ingenuity Pathway Analysis (IPA) (www.qiagen.com/ingenuity)tool was used to identify biological pathways in which DEGs were enriched and upstream regulator analysis (URA) identified upstream regulators for GBM.We performed pathway analysis and URA separately for DEGs that were up-regulated and down-regulated and were present in all three datasets.We also note that not all identified DEGs by the meta-analysis methods are present in all three studies considered.Number of DEGs present in one, two or all three datasets have been detailed in Table 4.In addition, the URA tool in IPA identified potential upstream regulators (transcription factors, genes or other small molecules) that has been experimentally observed to affect gene expression.It identifies these regulators by analysing linkage to DEGs through coordinated expression [36].Among the up-regulated DEGs, TGFB1 and TP53, which are also DEGs and important in GBM pathogenesis [37][38] are predicted to be the top two upstream regulators.293 up-regulated DEGs were identified as potential upstream regulators of gene upregulation out of a total of 2215 (BH corrected p-value < 0.01, see Supplementary 1: Table 7a) predicted upstream regulators.Out of 2215 predicted, 764 of these significant upstream regulators were activated and 112 were also observed as DEGs in our analysis.

Total
On the contrary, for the down-regulated DEGs, IPA identified 32 potential upstream regulators (BH corrected p-value < 0.01, see Supplementary 1: Table 7b) with TCF7L2 and MAPT as the top two.14 of the 32 upstream regulated were predicted to be inhibited and two among the inhibited are DEGs.TCF7L2 is a diabetes risk-associated gene which plays a key role in the Wnt-signaling pathway and is shown to be frequently mutated in colorectal cancer [39] and promote cell proliferation [40].However, exploration of its role in GBM pathogenesis warrant further studies.Interestingly, MAPT is also a DEG observed in our analysis and is one of the two hallmarks of AD [41].Gargini et al. (2019) [42] observed a strong correlation of Tau/MAPT expression and indicators of survival in glioma patients.Moreover, it has been found to be epigenetically controlled by balance between IDH1/2 wild-type and mutation in human gliomas [43].Thus, providing further evidence and reaffirming the involvement of MAPT in central nervous system disorders.
Of the DEGs with conflicting direction of expression across studies, 182 out of 235 DEGs are present in all three datasets.Among them CMTM6, RAD51, NOS1AP, MSANTD1, PGM2, PSD3, GPR82, SPTBN4, TSPAN6 and ARHGEF28 were identified as top 10 DEGs based on the absolute value of   (see Table 5).Interestingly, RAD51 and ARHGEF28 have previously been identified as a tumour suppressor and an oncogene respectively [44].More importantly, RAD51 was found to be effectively overexpressed in GBM in our study and have recently been shown as a target for inhibition to enhance radiosensitivity of GBM cells during treatment [45][46].On the other hand, ARHGEF28 was found to be effectively down-regulated in our study.It is an intracellular kinase that functions either as a Rho guanine exchange factor or a scaffolding protein to initiate FAK activation and cell contractibility [47].Furthermore, the RhoA-FAK pathway has been shown to be involved in colon cancer cell proliferation and migration [48].ARHGEF28 mRNA levels have also been found to be elevated in late-stage ovarian cancer and associated with decreased progression free and overall survival [49].However, its role in GBM growth and progression is yet to be elucidated and requires exploration in future studies.

Discussion
Although the implementation of MIN and FIN p-value combination methods are straight forward, they require some additional considerations.First, the used weighting criteria leads to a larger weight being given to a study with larger sample sizes.Intuitively, this is expected as a study with a larger sample size might be more robust than studies with lower sample sizes.However, importance must also be given to the quality of the RNA-seq data in each study.It must be assessed in case this information is available and other weights more appropriate as per the quality of the data may be specified.In our study, since TCGA-GBM has a much larger sample size as compared to the other two datasets, we compared the number of DEGs obtained by FIN method by considering all 165 samples and 10 random selections of 20 cases and 5 normal samples.On average, about 94% of the DEGs obtained when randomly selected subset was considered were also found in DEGs identified using the full TCGA-GBM dataset (see Supplementary 1: Table 8).Hence, suggesting that the identification of DEGs was stable across these two settings.
Next, the MIN and FIN are adaptive in a sense that they allow for consideration of genes that may not be present in all studies that are considered for integrated differential expression analysis.In case a gene is not present in some of the studies, the weights (  ) in the combination method can only be estimated using the number of replicates in the datasets in which the gene is present.However, for genes that are just present in one study, it would mean that the results from the meta-analysis for these genes would be the same as the per-study differential analysis.Hence, a careful consideration about the quality of the RNA-seq data and library size is required in case only the genes that are common among studies are considered.
For datasets of similar quality and library size, a large proportion of genes would not be excluded from meta-analysis if only common genes are used.However, a large number of genes might be excluded from meta-analysis in case of dissimilar library sizes and quality which could lead to potentially missing out on important genes for the disease.For instance, only 12345 out of 18315 unique genes are present in all 3 studies in our application where the library sizes are not similar.Thus, a balanced approach is suggested.
Finally, we used edgeR for per-study differential analysis in our study but other popular packages such as DESeq2 [50] and NOIseq [33] can be applied.Moreover, the FIN model can be extended to multi-group comparisons apart from a two-group comparison discussed in this study.The proposed meta-analysis method relies on the fact that the same test statistics are used for per-study differential expression analysis to obtain individual p-values and all studies under consideration have the same experimental considerations.For instance, in case DESeq2 is used for multi-group differential expression analysis in each study, a likelihood ratio test is used rather than Wald statistics being used for two group differential expression analysis.

Conclusions
In this study, we proposed MIN and consequently FIN method for meta-analysis of RNA-seq data.The developed methods account for both the sample size of study and direction of expression of a gene in each study allowing for detection of potentially robust biologically significant DEGs even when they have conflicting direction of expression across studies.In contrast with the existing IN method, the proposed methods have the advantage of identifying DEGs among genes with conflicting direction of expression across studies.For the genes with concordant differential expression patterns across studies the MIN method exhibited a similar DEG detection power and performance as compared to IN method particularly when there was high inter-study variability and increased number of studies were considered.FIN method exhibited a similar or improved DEG detection power as compared to IN method and was significantly better in performance as compared to MIN method.More importantly, in a real data application, we demonstrated the use of FIN method in detection of biologically relevant DEGs to GBM.Hence, this meta-analysis method provides a way to establish differential expression status for genes with conflicting direction of expression in individual RNA-seq studies and further exploration of them as potential biomarkers for the disease.With lowering costs and increase in the number of RNA-seq studies being archived on public databases, this method might provide a way to integrate a greater number of studies without losing much prior information and consequently considering all the genes in the analysis irrespective of their direction of expression.Tables Table 1.Simulation settings for inter-study variability parameter (), number of studies and number of replicates per study.Area under the ROC curves (AUC) for IN, MIN and FIN methods computed using 100 trials for each simulation setting.Table 2. Information about GBM RNA-seq datasets used for integrated analysis using different p-value combination methods in our study.Up and down DEGs refer to the up and downregulated DEGs obtained in per-study differential analysis.Table 3. Top 10 up-and down-regulated DEGs identified by FIN method.The DEGs have been sorted based on the value of the statistic   and the mean of absolute value of the log 2  have been reported.Effect signifies the direction of expression of DEGs in the per-study differential analysis.Table 4. Number of DEGs found in one, two or all three datasets.Same and mismatched represents if the direction of expression of a DEG was consistent across a study or not respectively.Table 5. Top 10 DEGs with mismatched direction of expression across datasets identified by FIN method.The DEGs have been sorted based on the absolute value of the statistic   and the mean of absolute value of the log 2  have been reported.Effect signifies the direction of expression of DEGs in the per-study differential analysis for GSE123892, GSE151352 and TCGA-GBM respectively.Supplementary Table 1.Total number of genes (# genes), number of genes with conflicting direction of expression (# conf.genes), number of differentially expressed genes (# DEGs), number of genes with conflicting direction of expression (# conf.DEGs) and number of true DEGs (# true DEGs) present for all simulation settings averaged over 100 trials.Supplementary Table 6.IPA canonical pathways for a. up and b. down-regulated DEGs present in all three GBM datasets and identified by the FIN method.Ratio denotes the number of DEGs enriched in a pathway to the total number of genes in that pathway.

Figure 1 .
Figure 1.Performance comparison of modified inverse-normal, inverse-normal and fused inverse-normal methods.Plots of receiver operating characteristics (ROC) curves averaged over 100 trials for each simulation setting.Simulation settings are represented by rows (from top to bottom): corresponding to low ( = 0.15) and high ( = 0.5) inter-study variability and columns (from left to right): corresponding to 3 (S=3) and 5 studies (S=5) combined.The red, turquoise and magenta ROC curves represent the modified inverse-normal, inverse-normal and fused inverse-normal methods respectively.

Figure 2 .
Figure 2. Characteristics of modified inverse-normal method.a).False discovery rates for modified inverse-normal method for all simulation settings.b).Proportion of true-positives (TPs) among unique differentially expressed genes (DEGs) identified by modified inversenormal (MIN) method as compared to inverse-normal (IN) method.c).Proportion of unique DEGs (MIN) with the observed effective direction of expression as the true direction of expression.

Figure 3 .
Figure 3. Characteristics of fused inverse-normal method.a).False discovery rates for fused inverse-normal method for all simulation settings.b).Proportion of true-positives (TPs) among unique differentially expressed genes (DEGs) identified by fused inverse-normal (FIN) method as compared to inverse-normal (IN) method.c).Proportion of unique DEGs (FIN) with the observed effective direction of expression as the true direction of expression.

Figure 4 .
Figure 4. Comparison of results from meta-analysis methods.a).Histograms of raw pvalues obtained from per-study differential analysis of GSE123892 and GSE151352 GBM datasets from gene expression omnibus used in real data application.b).Venn diagram of the differentially expressed genes identified using inverse-normal (IN DEGs), modified inversenormal (MIN DEGs) and fused inverse-normal (FIN DEGs) methods.

Figure 5 .
Figure 5. Significant pathways identified by IPA.The top ten significant pathways based on BH p-value among the canonical pathways identified by IPA for the up-regulated DEGs (orange bar) and down-regulated DEGs (green bar).The numbers on the bar plot show the ratio between the numbers of DEGs enriched and total number of genes in each of these pathways.

Figures Figure 3 .
FiguresFigure 3. Performance comparison of modified inverse-normal, inverse-normal and fused inverse-normal methods.Plots of receiver operating characteristics (ROC) curves averaged over 100 trials for each simulation setting.Simulation settings are represented by rows (from top to bottom): corresponding to low ( = 0.15) and high ( = 0.5) inter-study variability and columns (from left to right): corresponding to 3 (S=3) and 5 studies (S=5) combined.The red, turquoise and magenta ROC curves represent the modified inverse-normal, inverse-normal and fused inverse-normal methods respectively.

Figure 4 .
Figure 4. Characteristics of modified inverse-normal method.a).False discovery rates for modified inverse-normal method for all simulation settings.b).Proportion of true-positives (TPs) among unique differentially expressed genes (DEGs) identified by modified inversenormal (MIN) method as compared to inverse-normal (IN) method.c).Proportion of unique DEGs (MIN) with the observed effective direction of expression as the true direction of expression.

Figure 3 .
Figure 3. Characteristics of fused inverse-normal method.a).False discovery rates for fused inverse-normal method for all simulation settings.b).Proportion of true-positives (TPs) among unique differentially expressed genes (DEGs) identified by fused inverse-normal (FIN) method as compared to inverse-normal (IN) method.c).Proportion of unique DEGs (FIN) with the observed effective direction of expression as the true direction of expression.

Figure 4 .
Figure 4. Comparison of results from meta-analysis methods.a).Histograms of raw pvalues obtained from per-study differential analysis of GSE123892 and GSE151352 GBM

Figure 5 .
Figure 5. Significant pathways identified by IPA.The top ten significant pathways based on BH p-value among the canonical pathways identified by IPA for the up-regulated DEGs (orange bar) and down-regulated DEGs (green bar).The numbers on the bar plot show the ratio between the numbers of DEGs enriched and total number of genes in each of these pathways.

Figure 1 .
Figure 1.Flow of different steps used in processing of the raw RNA-seq fastq files for GSE151352 using Galaxy to get raw counts.

Supplementary Table 5 .
Number of DEGs identified by the fused inverse-normal (FIN) method in the GBM data application for both scenarios and number of these DEGs present in one, two or all three datasets considered.a.For DEGs with the same direction of expression across all three studies.b.For DEGs with conflicting direction of expression across studies a.For DEGs with the same direction of expression across all 3

Table 1 .
1, 100 independent trials were considered.Simulation settings for inter-study variability parameter (), number of studies and number of replicates per study.Area under the ROC curves (AUC) for IN, MIN and FIN methods computed using 100 trials for each simulation setting.

Table 3 .
Top 10 up-and down-regulated DEGs identified by FIN method.The DEGs have been sorted based on the value of the statistic   and the mean of absolute value of the log 2  have been reported.Effect signifies the direction of expression of DEGs in the per-study differential analysis.

Table 5 .
Top 10DEGs with mismatched direction of expression across datasets identified by FIN method.The DEGs have been sorted based on the absolute value of the statistic   and the mean of absolute value of the log 2  have been reported.Effect signifies the direction of expression of DEGs in the per-study differential analysis for GSE123892, GSE151352 and TCGA-GBM respectively.

Table 2 .
Total number of up and down-regulated DEGs identified in perstudy differential analysis in 10 random selections of 20 GBM cases and 5 controls from TCGA-GBM dataset.

Table 3 .
Number of DEGs identified by the inverse-normal (IN) method in the GBM data application for both scenarios and number of these DEGs present in one, two or all three datasets considered.

Table 4 .
Number of DEGs identified by the modified inverse-normal (MIN) method in the GBM data application for both scenarios and number of these DEGs present in one, two or all three datasets considered.a.For DEGs with the same direction of expression across all three studies.b.For DEGs with conflicting direction of expression across studies. a.

Table 7 .
IPA upstream regulator analysis results for a. up and b. downregulated DEGs present in all three GBM datasets and identified by the FIN method.

Table 8 .
Number of common DEGs identified by the FIN method in both scenarios considered for TCGA-GBM data.On average about 94% of the DEGs obtained when randomly selected subset was considered were also found in DEGs identified using the full TCGA-GBM dataset.