A Brief Introduction to TReNA
Overview of TReNA
TReNA is a package for inferring relationships between target genes and their candidate regulators. By using TReNA, we can generate hypotheses around regulation of gene expression and enable the construction of gene network models. In the most basic sense, TReNA provides two operations using the data. First, it filters putitive transcriptional regulators based on user-specified criteria, such as presence of motifs or footprints. Second, given a matrix of gene expression data (that must include the expression levels of transcription factors), a target gene, and list of candidate regulators as predictors, it uses a user-specified solver to perform feature selection. Thus, the standard TReNA workflow is the following:
- Select an expression matrix (RNA-seq or microarray) of interest.
- Use one of the built-in filtering methods to determine your set of candidate regulators.
- Use one of the built-in solver methods to perform feature selection, identifying putative transcriptional regulators for the target gene.
The filtering step, while technically optional, can constrain the solution space and reduce computational time by reducing the number of predictors supplied to the solver. Rather than supplying hundreds of transcription factors that can obscure meaningful relationships between genes, filters allow the user to screen out less likely candidates, leaving only the most likely affectors.
Starting off, we’ll load the TReNA library.
Selecting an Expression Matrix
Running the solvers in TReNA requires a matrix of gene expression data; there are several example data sets included with the TReNA package. In this case, we will load data from a study on Alzheimer’s disease (Nat Sci Data, doi:10.1038/sdata.2016.89), in which 278 samples were obtained for 154 transcription factors.
This gives us a matrix, mtx.sub, that contains our gene expression data in units of TMM. Before proceeding further , it is wise to inspect the distribution of the matrix and consider using some sort of data transformation to normalize the mean-variance relationship. We can see, for instance, that our dataset is fairly skewed, with most of the expression values falling close to 0 and only a few values on the higher end of the spectrum.
The skewed nature of the matrix may profoundly affect our results. To overcome this, we will transform our matrix using the voom function in the limma package.
If we now plot our transformed matrix, we can see that it’s been scaled down substantially.
In general, we recommend a transformation of this nature to adjust skewed data. The voom transformation or something similar (e.g. asinh()) adjusts the mean-variance relationship and can greatly change the results we see. Moving forward, we will use this transformed data set in the rest of our analyses.
Basic Interface: The Trena Class
The most straightforward way to filter candidates and perform feature selection is by constructing a Trena object.
Filtering Candidate Regulators
TReNA currently supports 4 different candidate transcription factor inputs, all of which are subclasses of the CandidateFilter class:
VarianceFilter: returns transcription factors with expression variance within a certain range of the target gene’s expression varianceGeneOntologyFilter: given an organism database and a Gene Ontology term (GOterm), returns only those genes that match the GOtermHumanDHSFilter: given a genome database, DNAse hypersensitivity (DHS) database, and region specifications, return only those genes that fall within the given regionFootprintFilter: using a PostgreSQL or SQLite database of footprint data, returns transcription factors with footprints for a given target gene for a specified genomic range
In order to get a list of transcription factor candidates, you must construct the appropriate object by supplying the necessary parameters for that CandidateFilter, then call the getCandidates() method on the candidate inputs. In each case, the call to getCandidates() is identical and operates on only the CandidateFilter object.
The VarianceFilter Subclass
The VarianceFilter requires the specification of the expression matrix and a target gene, plus it allows you to optionally specify a variance size range. It finds the variance of all genes in the assay matrix, then returns all transcription factors with variance within the variance size range of the target gene. For instance, we can create a VarianceFilter and use it to find all transcription factors with variance within 50% of the target gene’s variance. This will return a named list with both the names of the transcription factors and their variances.
## List of 2
## $ tfs : chr [1:33] "BCL6B" "CREB5" "CUX2" "DRGX" ...
## $ tf.vars: Named num [1:33] 0.799 0.57 0.996 0.578 0.504 ...
## ..- attr(*, "names")= chr [1:33] "BCL6B" "CREB5" "CUX2" "DRGX" ...
The GeneOntologyFilter Subclass
The
Feature Selection: the Solver Object
Once you know your candidate regulators and you’ve got a target gene in mind, you’re ready to construct a Solver object. A Solver object is a construct that contains the expression matrix, candidate regulators, target gene, and a solver method of choice. . TReNA currently supports 9 different solvers. Their associated solvers strings for use in TReNA are (in no particular order):
- LassoSolver - an application of
glmnet using alpha = 0.9 as the default elastic net parameter
- RidgeSolver - an application of
glmnet using alpha = 0 as the default elastic net parameter
- RandomForestSolver - an application of `randomForest*, the standard Random Forest function
- BayesSpikeSolver - an application of
vbsr, the Bayes Spike function from the vbsr package
- SqrtLassoSolver - an application of the
slim function from the flare package, with q = 2
- LassoPVSolver - an application of
lassopv, the P-Value LASSO function
- PearsonSolver - an application of the
cor function using the default parameters
- SpearmanSolver - an application of the
cor function using method = "spearman"
- EnsembleSolver - a combination of 2 or more of the other solvers; by default, it uses all but
BayesSpikeSolver and SqrtLassoSolver
When it comes to selecting a solver, we recommend doing some critical thinking about which solver most suits your purposes before choosing. Each of solvers has its own strengths and weaknesses that we will not delve into here, but that are generally detailed within their own packages. To begin with, we will use the LassoSolver, which uses LASSO regularlized regression to select features. Once a solver has been chosen, the particular Solver object is simply specified:
We have now created an object of the LassoSolver class, which in this case uses “MEF2C” as the target gene and the candidate regulators we found using the FootprintFilter above. Regardless of the specific subclass, any Solver object can be executed using the run command:
## 'data.frame': 7 obs. of 2 variables:
## $ beta : num 1.625 -0.335 -0.229 -0.2 -0.186 ...
## $ intercept: num 0.739 0.739 0.739 0.739 0.739 ...
The output of each run method can vary with solver type, but in general it will return the coefficients for each feature, in this case contained in the beta field. Notably, the LassoSolver, SqrtLassoSolver, and EnsembleSolver often return subsets of the supplied candidate regulators, whereas the remaining methods return coefficients for all regulators. This behavior can be tuned as well; for instance, we can alter the LassoSolver object by adjusting lambda, a tuning parameter that corresponds to the penalty imposed on regression coefficients.
First, we will alter a different parameter of the LassoSolver object, the keep.metrics parameter. This optional Boolean argument will return several parameters of the fit if specified as TRUE (default = FALSE), including the lambda value determined via permutation testing. We can see this behavior by creating and running this new LassoSolver object:
## List of 3
## $ mtx.beta:'data.frame': 7 obs. of 2 variables:
## ..$ beta : num [1:7] 1.625 -0.335 -0.229 -0.2 -0.186 ...
## ..$ intercept: num [1:7] 0.739 0.739 0.739 0.739 0.739 ...
## $ lambda : num 0.0862
## $ r2 : Named num 0.649
## ..- attr(*, "names")= chr "1"
As shown, we now get a list where the first element is our data frame from the previous example, the second is the value determined for lambda, and the final element is the value determined for \(r^{2}\). Notably, the keep.metrics argument is only available for the LassoSolver and RidgeSolver objects at this time.
Now we can change the value of lambda to change how many features are returned by the solver. In this case, we’ll make lambda much smaller; this will relax the penalty imposed on regression coefficients, thus we expect to get more features returned to us:
## 'data.frame': 7 obs. of 2 variables:
## $ beta : num 1.7 -0.357 -0.248 -0.247 -0.245 ...
## $ intercept: num -0.0443 -0.0443 -0.0443 -0.0443 -0.0443 ...
As expected, we get nearly twice as many features as in our original example. Different Solver subclasses have different optional arguments corresponding to their specific solver methods. For more information on Solver object, check the documentation for your object of interest; for more information on a particular solver method, check the original package referenced in the help documentation for that Solver object.
Each solver will necessarily give a different result, neccessitating some amount of thoughtfulness on the part of the user. For instance, if the goal is to return a sparse list of only the most influential genes, then LASSO may be the best choice. By contrast, using the Random Forest method will return scores for all genes, and will thus result in inferring more relationships.
Using an Ensemble Solver Approach
In addition to the individual solvers, TReNA includes an ensemble approach as a solver option. The EnsembleSolver class allows the user to specify a vector of multiple TReNA solvers as part of the argument list (e.g. solverNames = c("lasso","ridge")). Alternatively, the user can leave this argument unspecified to run all solvers but Bayes Spike and square root LASSO (the most sensitive to outliers and the computationally intensive, respectively). Other solver-specific parameters may also be specified, such as the lambda value for LASSO; these additional arguments can be found in the help documentation for the EnsembleSolver
The EnsembleSolver also has a geneCutoff argument, whereby the run method cuts off a certain fraction of the total regulators supplied. For example, if 200 candidates are supplied and the geneCutoff is left at its default value (default = 0.1), only roughly 20 genes will be returned in the final data frame. The geneCutoff argument is designed to help the user hone in on only the most important features in a model, but it can be changed to 100% of the features if so desired. In cases of smaller sets of regulators, the cutoff will cease being imposed once the number of features drops below 10.
Besides computing and returning the individual scores of each solver, the ensemble solver generates an “extreme” score (pcaMax) and a “concordance” score (concordance) based on all solvers. Both of these new scores are computed using a principle component analysis on all individual solver scores and can perhaps best be understood as follows:
pcaMax essentially quantifies the extent to which dissimilar solvers agree on selecting a given feature. For instance, if Random Forest and LASSO both deem a particular feature important, then it is likely to have a high pcaMax score. By contrast, if Pearson and Spearman both deem a feature important, it would be less likely to have a high pcaMax score because the methods are quite similar to one another. The basic idea here is that if lots of dissimilar methods agree on the importance of a feature, it is more likely to actually be important.
concordance essentially weighs each solver equally and sums up the contributions from each one, regardless of method similarity. Thus, if Random Forst and LASSO both deem a particular feature important, it is equally likely to have a high concordance score as if Pearson and Spearman both deem a feature important. This score functions as more of a basic “voting” system of which solves consider a feature significant and ranges from 0-1.
We can see this in action by constructing and running an EnsembleSolver object with the default arguments for solverNames and geneCutoff:
## gene betaLasso lassoPValue pearsonCoeff rfScore betaRidge
## 3 FOXO3 1.63008787 2.696918e-38 0.67481919 72.555013 1.05142781
## 6 SOX1 0.06836318 4.597246e-04 0.25367217 16.405686 0.15069654
## 4 LHX2 0.02236146 9.008229e-03 0.03267022 28.070971 0.04759860
## 7 SOX8 -0.23022048 1.837246e-16 -0.54312583 37.091426 -0.24392258
## 2 FOXL1 0.00000000 4.122076e-01 -0.22176415 8.143075 -0.06977827
## 8 VSX1 -0.33642027 4.101894e-17 -0.44575194 22.936443 -0.26481503
## 5 POU4F1 -0.19048425 1.231688e-06 -0.39202043 20.115456 -0.32708785
## 1 FOXD2 -0.20292360 4.948668e-10 -0.36088504 25.887043 -0.19184874
## spearmanCoeff concordance pcaMax
## 3 0.7767056 0.40348652 2.4045903
## 6 0.1073104 0.32259095 0.6803177
## 4 -0.3232110 0.19837126 0.3385371
## 7 -0.5830969 0.19387110 0.3305646
## 2 -0.1999422 0.16615970 0.2739545
## 8 -0.4358512 0.14929858 0.2453073
## 5 -0.3767792 0.11664126 0.1861632
## 1 -0.4688832 0.08779695 0.1396837
Note that in the final data frame, regulators are returned in descending order of pcaMax score. If we so choose, we can change the solverNames and geneCutoff as well; for this example, we’ll change them to just 4 solver methods and 100% of the genes, respectively:
## 'data.frame': 8 obs. of 7 variables:
## $ gene : chr "FOXO3" "SOX1" "LHX2" "FOXL1" ...
## $ betaLasso : num 1.6345 0.0695 0.0257 0 -0.2312 ...
## $ pearsonCoeff: num 0.6748 0.2537 0.0327 -0.2218 -0.5431 ...
## $ rfScore : num 75.27 15.63 26.93 7.74 35.82 ...
## $ betaRidge : num 1.0289 0.1497 0.0455 -0.0719 -0.2406 ...
## $ concordance : num 0.446 0.321 0.231 0.167 0.157 ...
## $ pcaMax : num 2.648 0.633 0.42 0.275 0.257 ...
As demonstrated, the full set of regulators here is 64; when using a 10% cutoff, we only get 9 features in the final model as the cutoff is designed to go into effect only for 10 or more features.
The ensemble approach is an attractive option for a couple of reasons:
Solver choice can be a non-trivial task for a TReNA user, as it requires some working knowledge of the pros and cons of the different methods. The ensemble method alleviates the need to choose one solver.
Different solvers will give different predictions for important transcription factors; the ensemble method provides a way to create a composite score that accounts for all included solvers and can thus be seen as something of an overall metric.
The relationship between transcription factors and target genes is a young area of active investigation. How this relates to the various solvers is an open question. TReNA was born out of a desire to explore these relationships.
Thus, if you are a first-time TReNA user, we recommend that you start with the EnsembleSolver as your Solver object of choice.