--- title: "projectR Vignette" author: - "Gaurav Sharma" - "Jared N. Slosberg" - "Loyal A. Goff" - "Genevieve L. Stein-O'Brien" date: "`r BiocStyle::doc_date()`" output: BiocStyle::html_document bibliography: projectR.bib description: | Functions for the Projection of Weights from PCA, CoGAPS, NMF, Correlation, and Clustering vignette: > %\VignetteIndexEntry{projectR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include=FALSE} knitr::opts_chunk$set(echo = TRUE) options(scipen = 1, digits = 2) set.seed(1234) ``` # Introduction Technological advances continue to spur the exponential growth of biological data as illustrated by the rise of the omics—genomics, transcriptomics, epigenomics, proteomics, etc.—each with there own high throughput technologies. In order to leverage the full power of these resources, methods to integrate multiple data sets and data types must be developed. The reciprocal nature of the genomic, transcriptomic, epigenomic, and proteomic biology requires that the data provides a complementary view of cellular function and regulatory organization; however, the technical heterogeneity and massive size of high-throughput data even within a particular omic makes integrated analysis challenging. To address these challenges, we developed projectR, an R package for integrated analysis of high dimensional omic data. projectR uses the relationships defined within a given high dimensional data set, to interrogate related biological phenomena in an entirely new data set. By relying on relative comparisons within data type, projectR is able to circumvent many issues arising from technological variation. For a more extensive example of how the tools in the projectR package can be used for *in silico* experiments, or additional information on the algorithm, see [Stein-O'Brien, et al](https://www.sciencedirect.com/science/article/pii/S2405471219301462) and [Sharma, et al](https://academic.oup.com/bioinformatics/article/36/11/3592/5804979). # Getting started with projectR ## Installation Instructions For automatic Bioconductor package installation, start R, and run: ``` BiocManager::install("projectR") ``` ## Methods Projection can roughly be defined as a mapping or transformation of points from one space to another often lower dimensional space. Mathematically, this can described as a function $\varphi(x)=y : \Re^{D} \mapsto \Re^{d}$ s.t. $d \leq D$ for $x \in \Re^{D}, y \in \Re^{d}$ @Barbakh:2009bw . The projectR package uses projection functions defined in a training dataset to interrogate related biological phenomena in an entirely new data set. These functions can be the product of any one of several methods common to "omic" analyses including regression, PCA, NMF, clustering. Individual sections focusing on one specific method are included in the vignette. However, the general design of the projectR function is the same regardless. ## The base projectR function The generic projectR function is executed as follows: ``` library(projectR) projectR(data, loadings, dataNames=NULL, loadingsNames=NULL, NP = NULL, full = false) ``` ### Input Arguments The inputs that must be set each time are only the data and loadings, with all other inputs having default values. However, incongruities in the feature mapping between the data and loadings, i.e. a different format for the rownames of each object, will throw errors or result in an empty mapping and should be checked before running. To overcoming mismatched feature names in the objects themselves, the ``dataNames`` and `loadingNames` arguments can be manually supplied by the user. The arguments are as follows: **`data`** a dataset to be projected into the pattern space **`loadings`** a matrix of continous values with unique rownames to be projected **`dataNames`** a vector containing unique name, i.e. gene names, for the rows of the target dataset to be used to match features with the loadings, if not provided by `rownames(data)`. Order of names in vector must match order of rows in data. **`loadingsNames`** a vector containing unique names, i.e. gene names, for the rows of loadings to be used to match features with the data, if not provided by `rownames(loadings)`. Order of names in vector must match order of rows in loadings. **`NP`** vector of integers indicating which columns of loadings object to use. The default of NP = NA will use entire matrix. **`full`** logical indicating whether to return the full model solution. By default only the new pattern object is returned. The `loadings` argument in the generic projectR function is suitable for use with any genernal feature space, or set of feature spaces, whose rows annotation links them to the data to be projected. Ex: the coeffients associated with individual genes as the result of regression analysis or the amplituded values of individual genes as the result of non-negative matrix factorization (NMF). ### Output The basic output of the base projectR function, i.e. `full=FALSE`, returns `projectionPatterns` representing relative weights for the samples from the new data in this previously defined feature space, or set of feature spaces. The full output of the base projectR function, i.e. `full=TRUE`, returns `projectionFit`, a list containing `projectionPatterns` and `Projection`. The `Projection` object contains additional information from the proceedure used to obtain the `projectionPatterns`. For the the the base projectR function, `Projection` is the full `lmFit` model from the package `r BiocStyle::Biocpkg("limma")`. # PCA projection Projection of principal components is achieved by matrix multiplication of a new data set by previously generated eigenvectors, or gene loadings. If the original data were standardized such that each gene is centered to zero average expression level, the principal components are normalized eigenvectors of the covariance matrix of the genes. Each PC is ordered according to how much of the variation present in the data they contain. Projection of the original samples into each PC will maximize the variance of the samples in the direction of that component and uncorrelated to previous components. Projection of new data places the new samples into the PCs defined by the original data. Because the components define an orthonormal basis set, they provide an isomorphism between a vector space, $V$, and $\Re^n$ which preserves inner products. If $V$ is an inner product space over $\Re$ with orthonormal basis $B = v_1,...,v_n$ and $v \epsilon V s.t [v]_B = (r_1,...,r_n)$, then finding the coordinate of $v_i$ in $v$ is precisely the inner product of $v$ with $v_i$, i.e. $r_i = \langle v,v_i \rangle$. This formulation is implemented for only those genes belonging to both the new data and the PC space. The **`projectR`** function has S4 method for class `prcomp`. ## Obtaining PCs to project. ```{r prcomp, warning=FALSE} # data to define PCs library(projectR) data(p.RNAseq6l3c3t) # do PCA on RNAseq6l3c3t expression data pc.RNAseq6l3c3t<-prcomp(t(p.RNAseq6l3c3t)) pcVAR <- round(((pc.RNAseq6l3c3t$sdev)^2/sum(pc.RNAseq6l3c3t$sdev^2))*100,2) dPCA <- data.frame(cbind(pc.RNAseq6l3c3t$x,pd.RNAseq6l3c3t)) #plot pca library(ggplot2) setCOL <- scale_colour_manual(values = c("blue","black","red"), name="Condition:") setFILL <- scale_fill_manual(values = c("blue","black","red"),guide = FALSE) setPCH <- scale_shape_manual(values=c(23,22,25,25,21,24),name="Cell Line:") pPCA <- ggplot(dPCA, aes(x=PC1, y=PC2, colour=ID.cond, shape=ID.line, fill=ID.cond)) + geom_point(aes(size=days),alpha=.6)+ setCOL + setPCH + setFILL + scale_size_area(breaks = c(2,4,6), name="Day") + theme(legend.position=c(0,0), legend.justification=c(0,0), legend.direction = "horizontal", panel.background = element_rect(fill = "white",colour=NA), legend.background = element_rect(fill = "transparent",colour=NA), plot.title = element_text(vjust = 0,hjust=0,face="bold")) + labs(title = "PCA of hPSC PolyA RNAseq", x=paste("PC1 (",pcVAR[1],"% of varience)",sep=""), y=paste("PC2 (",pcVAR[2],"% of varience)",sep="")) ``` ## Projecting prcomp objects ```{r projectR.prcomp, warning=FALSE} # data to project into PCs from RNAseq6l3c3t expression data data(p.ESepiGen4c1l) library(projectR) PCA2ESepi <- projectR(data = p.ESepiGen4c1l$mRNA.Seq,loadings=pc.RNAseq6l3c3t, full=TRUE, dataNames=map.ESepiGen4c1l[["GeneSymbols"]]) pd.ESepiGen4c1l<-data.frame(Condition=sapply(colnames(p.ESepiGen4c1l$mRNA.Seq), function(x) unlist(strsplit(x,'_'))[1]),stringsAsFactors=FALSE) pd.ESepiGen4c1l$color<-c(rep("red",2),rep("green",3),rep("blue",2),rep("black",2)) names(pd.ESepiGen4c1l$color)<-pd.ESepiGen4c1l$Cond dPCA2ESepi<- data.frame(cbind(t(PCA2ESepi[[1]]),pd.ESepiGen4c1l)) #plot pca library(ggplot2) setEpiCOL <- scale_colour_manual(values = c("red","green","blue","black"), guide = guide_legend(title="Lineage")) pPC2ESepiGen4c1l <- ggplot(dPCA2ESepi, aes(x=PC1, y=PC2, colour=Condition)) + geom_point(size=5) + setEpiCOL + theme(legend.position=c(0,0), legend.justification=c(0,0), panel.background = element_rect(fill = "white"), legend.direction = "horizontal", plot.title = element_text(vjust = 0,hjust=0,face="bold")) + labs(title = "Encode RNAseq in target PC1 & PC2", x=paste("Projected PC1 (",round(PCA2ESepi[[2]][1],2),"% of varience)",sep=""), y=paste("Projected PC2 (",round(PCA2ESepi[[2]][2],2),"% of varience)",sep="")) ``` ```{r, fig.show='hold', fig.width=10, fig.height=5, echo=FALSE, message= FALSE} library(gridExtra) #grid.arrange(pPCA,pPC2ESepiGen4c1l,nrow=1) ``` # NMF projection NMF decomposes a data matrix of $D$ with $N$ genes as rows and $M$ samples as columns, into two matrices, as $D ~ AP$. The pattern matrix P has rows associated with BPs in samples and the amplitude matrix A has columns indicating the relative association of a given gene, where the total number of BPs (k) is an input parameter. CoGAPS and GWCoGAPS seek a pattern matrix (${\bf{P}}$) and the corresponding distribution matrix of weights (${\bf{A}}$) whose product forms a mock data matrix (${\bf{M}}$) that represents the gene-wise data ${\bf{D}}$ within noise limits ($\boldsymbol{\varepsilon}$). That is, \begin{equation} {\bf{D}} = {\bf{M}} + \boldsymbol{\varepsilon} = {\bf{A}}{\bf{P}} + \boldsymbol{\varepsilon} ..............(1) \label{eq:matrixDecomp} \end{equation} The number of rows in ${\bf{P}}$ (columns in ${\bf{A}}$) defines the number of biological patterns (k) that CoGAPS/GWCoGAPS will infer from the number of nonorthogonal basis vectors required to span the data space. As in the Bayesian Decomposition algorithm @Ochs2006, the matrices ${\bf{A}}$ and ${\bf{P}}$ in CoGAPS are assumed to have the atomic prior described in @Sibisi1997. In the CoGAPS/GWCoGAPS implementation, $\alpha_{A}$ and $\alpha_{P}$ are corresponding parameters for the expected number of atoms which map to each matrix element in ${\bf{A}}$ and ${\bf{P}}$, respectively. The corresponding matrices ${\bf{A}}$ and ${\bf{P}}$ are found by MCMC sampling. Projection of CoGAPS/GWCoGAPS patterns is implemented by solving the factorization in (1) for the new data matrix where ${\bf{A}}$ is the fixed nonorthogonal basis vectors comprising the average of the posterior mean for the CoGAPS/GWCoGAPS simulations performed on the original data. The patterns ${\bf{P}}$ in the new data associated with this amplitude matrix is estimated using the least-squares fit to the new data implemented with the `lmFit` function in the `r BiocStyle::Biocpkg("limma")` package. The `projectR` function has S4 method for class `Linear Embedding Matrix, LME`. ``` library(projectR) projectR(data, loadings,dataNames = NULL, loadingsNames = NULL, NP = NA, full = FALSE) ``` ### Input Arguments The inputs that must be set each time are only the data and patterns, with all other inputs having default values. However, inconguities between gene names--rownames of the loadings object and either rownames of the data object will throw errors and, subsequently, should be checked before running. The arguments are as follows: **`data`** a target dataset to be projected into the pattern space **`loadings`** a CogapsResult object **`dataNames`** rownames (eg. gene names) of the target dataset, if different from existing rownames of data **`loadingsNames`** loadingsNames rownames (eg. gene names) of the loadings to be matched with dataNames **`NP`** vector of integers indicating which columns of loadings object to use. The default of NP = NA will use entire matrix. **`full`** logical indicating whether to return the full model solution. By default only the new pattern object is returned. ### Output The basic output of the base projectR function, i.e. `full=FALSE`, returns `projectionPatterns` representing relative weights for the samples from the new data in this previously defined feature space, or set of feature spaces. The full output of the base projectR function, i.e. `full=TRUE`, returns `projectionFit`, a list containing `projectionPatterns` and `Projection`. The `Projection` object contains additional information from the procedure used to obtain the `projectionPatterns`. For the the the base projectR function, `Projection` is the full `lmFit` model from the package `r BiocStyle::Biocpkg('limma')`. ## Obtaining CoGAPS patterns to project. ```{r} # get data library(projectR) AP <- get(data("AP.RNAseq6l3c3t")) #CoGAPS run data AP <- AP$Amean # heatmap of gene weights for CoGAPs patterns library(gplots) par(mar=c(1,1,1,1)) pNMF<-heatmap.2(as.matrix(AP),col=bluered, trace='none', distfun=function(c) as.dist(1-cor(t(c))) , cexCol=1,cexRow=.5,scale = "row", hclustfun=function(x) hclust(x, method="average") ) ``` ## Projecting CoGAPS objects ```{r} # data to project into PCs from RNAseq6l3c3t expression data library(projectR) data('p.ESepiGen4c1l4') data('p.RNAseq6l3c3t') NMF2ESepi <- projectR(p.ESepiGen4c1l$mRNA.Seq,loadings=AP,full=TRUE, dataNames=map.ESepiGen4c1l[["GeneSymbols"]]) dNMF2ESepi<- data.frame(cbind(t(NMF2ESepi),pd.ESepiGen4c1l)) #plot pca library(ggplot2) setEpiCOL <- scale_colour_manual(values = c("red","green","blue","black"), guide = guide_legend(title="Lineage")) pNMF2ESepiGen4c1l <- ggplot(dNMF2ESepi, aes(x=X1, y=X2, colour=Condition)) + geom_point(size=5) + setEpiCOL + theme(legend.position=c(0,0), legend.justification=c(0,0), panel.background = element_rect(fill = "white"), legend.direction = "horizontal", plot.title = element_text(vjust = 0,hjust=0,face="bold")) labs(title = "Encode RNAseq in target PC1 & PC2", x=paste("Projected PC1 (",round(PCA2ESepi[[2]][1],2),"% of varience)",sep=""), y=paste("Projected PC2 (",round(PCA2ESepi[[2]][2],2),"% of varience)",sep="")) ``` # Clustering projection As canonical projection is not defined for clustering objects, the projectR package offers two transfer learning inspired methods to achieve the "projection" of clustering objects. These methods are defined by the function used to quantify and transfer the relationships which define each cluster in the original data set to the new dataset. Briefly, `cluster2pattern` uses the corelation of each genes expression to the mean of each cluster to define continuous weights. These weights are output as a `pclust` object which can serve as input to `projectR`. Alternatively, the `intersectoR` function can be used to test for significant overlap between two clustering objects. Both `cluster2pattern` and `intersectoR` methods are coded for a generic list structure with additional S4 class methods for kmeans and hclust objects. Further details and examples are provided in the followin respecitive sections. ## cluster2pattern `cluster2pattern` uses the corelation of each genes expression to the mean of each cluster to define continuous weights. ``` library(projectR) data(p.RNAseq6l3c3t) nP<-5 kClust<-kmeans(t(p.RNAseq6l3c3t),centers=nP) kpattern<-cluster2pattern(clusters = kClust, NP = nP, data = p.RNAseq6l3c3t) kpattern cluster2pattern(clusters = NA, NP = NA, data = NA) ``` ### Input Arguments The inputs that must be set each time are the clusters and data. The arguments are as follows: **`clusters`** a clustering object **`NP`** either the number of clusters desired or the subset of clusters to use **`data`** data used to make clusters object ### Output The output of the `cluster2pattern` function is a `pclust` class object; specifically, a matrix of genes (rows) by clusters (columns). A gene's value outside of its assigned cluster is zero. For the cluster containing a given gene, the gene's value is the correlation of the gene's expression to the mean of that cluster. ## intersectoR `intersectoR` function can be used to test for significant overlap between two clustering objects. The base function finds and tests the intersecting values of two sets of lists, presumably the genes associated with patterns in two different datasets. S4 class methods for `hclust` and `kmeans` objects are also available. ``` library(projectR) intersectoR(pSet1 = NA, pSet2 = NA, pval = 0.05, full = FALSE, k = NULL) ``` ### Input Arguments The inputs that must be set each time are the clusters and data. The arguments are as follows: **`pSet1`** a list for a set of patterns where each entry is a set of genes associated with a single pattern **`pSet2`** a list for a second set of patterns where each entry is a set of genes associated with a single pattern **`pval`** the maximum p-value considered significant **`full`** logical indicating whether to return full data frame of signigicantly overlapping sets. Default is false will return summary matrix. **`k`** numeric giving cut height for hclust objects, if vector arguments will be applied to pSet1 and pSet2 in that order ### Output The output of the `intersectoR` function is a summary matrix showing the sets with statistically significant overlap under the specified $p$-value threshold based on a hypergeometric test. If `full==TRUE` the full data frame of significantly overlapping sets will also be returned. # Correlation based projection Correlation based projection requires a matrix of gene-wise correlation values to serve as the Pattern input to the `projectR` function. This matrix can be user-generated or the result of the `correlateR` function included in the projectR package. User-generated matrixes with each row corresponding to an individual gene can be input to the generic `projectR` function. The `correlateR` function allows users to create a weight matrix for projection with values quantifying the within dataset correlation of each genes expression to the expression pattern of a particular gene or set of genes as follows. ## correlateR ``` library(projectR) correlateR(genes = NA, dat = NA, threshtype = "R", threshold = 0.7, absR = FALSE, ...) ``` ### Input Arguments The inputs that must be set each time are only the genes and data, with all other inputs having default values. The arguments are as follows: **`genes`** gene or character vector of genes for reference expression pattern dat **`data`** matrix or data frame with genes to be used for to calculate correlation **`threshtype`** Default "R" indicates thresholding by R value or equivalent. Alternatively, "N" indicates a numerical cut off. **`threshold`** numeric indicating value at which to make threshold **`absR`** logical indicating where to include both positive and negatively correlated genes **`...`** addtion imputes to the cor function ### Output The output of the `correlateR` function is a `correlateR` class object. Specifically, a matrix of correlation values for those genes whose expression pattern pattern in the dataset is correlated (and anti-correlated if absR=TRUE) above the value given in as the threshold arguement. As this information may be useful in its own right, it is recommended that users inspect the `correlateR` object before using it as input to the `projectR` function. ## Obtaining and visualizing `correlateR` objects. ```{r correlateR-exp} # data to library(projectR) data("p.RNAseq6l3c3t") # get genes correlated to T cor2T<-correlateR(genes="T", dat=p.RNAseq6l3c3t, threshtype="N", threshold=10, absR=TRUE) cor2T <- cor2T@corM ### heatmap of genes more correlated to T indx<-unlist(sapply(cor2T,rownames)) indx <- as.vector(indx) colnames(p.RNAseq6l3c3t)<-pd.RNAseq6l3c3t$sampleX library(reshape2) pm.RNAseq6l3c3t<-melt(cbind(p.RNAseq6l3c3t[indx,],indx)) library(gplots) library(ggplot2) library(viridis) pCorT<-ggplot(pm.RNAseq6l3c3t, aes(variable, indx, fill = value)) + geom_tile(colour="gray20", size=1.5, stat="identity") + scale_fill_viridis(option="B") + xlab("") + ylab("") + scale_y_discrete(limits=indx) + ggtitle("Ten genes most highly pos & neg correlated with T") + theme( panel.background = element_rect(fill="gray20"), panel.border = element_rect(fill=NA,color="gray20", size=0.5, linetype="solid"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_blank(), axis.ticks = element_blank(), axis.text = element_text(size=rel(1),hjust=1), axis.text.x = element_text(angle = 90,vjust=.5), legend.text = element_text(color="white", size=rel(1)), legend.background = element_rect(fill="gray20"), legend.position = "bottom", legend.title=element_blank() ) ``` ```{r, fig.show='hold', fig.width=10, fig.height=5, echo=FALSE} pCorT ``` ## Projecting correlateR objects. ```{r} # data to project into from RNAseq6l3c3t expression data data(p.ESepiGen4c1l) library(projectR) cor2ESepi <- projectR(p.ESepiGen4c1l$mRNA.Seq,loadings=cor2T[[1]],full=FALSE, dataNames=map.ESepiGen4c1l$GeneSymbols) ``` # Differential features identification. ## projectionDriveR Given loadings that define the weight of features (genes) in a given latent space (e.g. PCA, NMF), and the use of these patterns in samples, it is of interest to look at differential usage of these features between conditions. These conditions may be defined by user-defined annotations of cell type or by differential usage of a (projected) pattern. By examining differences in gene expression, weighted by the loadings that define their importance in a specific latent space, a unique understanding of differential expression in that context can be gained. This approach was originally proposed and developed in (Baraban et al, 2021), which demonstrates its utility in cross-celltype and cross-species interpretation of pattern usages. ``` library(projectR) projectionDriveR(cellgroup1, cellgroup2, loadings, loadingsNames = NULL, pvalue, pattern_name, display = T, normalize_pattern = T) ``` ### Input Arguments The required inputs are two feature by sample (e.g. gene by cell) matrices to be compared, the loadings that define the feature weights, and the name of the pattern (column of feature loadings). If applicable, the expression matrices should already be corrected for variables such as sequencing depth. The arguments for projectionDriveR are: **`cellgroup1`** Matrix 1 with features as rows, samples as columns. **`cellgroup2`** Matrix 2 with features as rows, samples as columns. **`loadings`** Matrix or dataframe with features as rows, columns as patterns. Values define feature weights in that space **`loadingsNames`** Vector of names corresponding to rows of loadings. By default the rownames of loadings will be used **`pattern_name`** the column name of the loadings by which the features will be weighted **`pvalue`** Determines the significance of the confidence interval to be calculated between the difference of means **`display`** Boolean. Whether or not to plot the estimates of significant features. Default = T **`normalize_pattern`** Boolean. Whether or not to normalize the average feature weight. Default = T ### Output The output of `projectionDriveR` is a list of length five `mean_ci` holds the confidence intervals for the difference in means for all features, `weighted_ci` holds the confidence intervals for the weighted difference in means for all features, and normalized_weights are the weights themselves. In addition, `significant_genes` is a vector of gene names that are significantly different at the threshold provided. `plotted_ci` returns the ggplot figure of the confidence intervals, see `plotConfidenceIntervals` for documentation. ### Identifying differential features associated with learned patterns ```{r projectionDriver, message = F, out.width="100%"} options(width = 60) library(projectR) library(dplyr, warn.conflicts = F) #gene weights x pattern data("retinal_patterns") #size-normed, log expression data("microglial_counts") #size-normed, log expression data("glial_counts") #the features by which to weight the difference in expression pattern_to_weight <- "Pattern.24" drivers <- projectionDriveR(microglial_counts, #expression matrix glial_counts, #expression matrix loadings = retinal_patterns, #feature x pattern dataframe loadingsNames = NULL, pattern_name = pattern_to_weight, #column name pvalue = 1e-5, #pvalue before bonferroni correction display = T, normalize_pattern = T) #normalize feature weights conf_intervals <- drivers$mean_ci[drivers$significant_genes,] str(conf_intervals) ``` ## plotConfidenceIntervals ### Input The arguments for plotConfidenceIntervals are: **`confidence_intervals`** A dataframe of features x estimates **`interval_name`** names of columns that contain the low and high estimates, respectively. (default: c("low","high")) **`pattern_name`** string to use as the title for the plots **`sort`** Boolean. Whether or not to sort genes by their estimates (default = T) **`genes`** a vector with names of genes to include in plot. If sort=F, estimates will be plotted in this order (default = NULL will include all genes.) **`weights`** weights of features to include as annotation (default = NULL will not include heatmap) **`weights_clip`** quantile of data to clip color scale for improved visualization (default: 0.99) **`weights_vis_norm`** Which processed version of weights to visualize as a heatmap. One of c("none", "quantile"). default = "none" ### Output A list of the length three that includes confidence interval plots and relevant info. `ci_estimates_plot` is the point-range plot for the provided estimates. If called from within `projectionDriveR`, the unweighted estimates are used. `feature_order` is the vector of gene names in the order shown in the figure. `weights_heatmap` is a heatmap annotation of the gene loadings, in the same order as above. ### Customize plotting of confidence intervals ```{r} suppressWarnings(library(cowplot)) #order in ascending order of estimates conf_intervals <- conf_intervals %>% mutate(mid = (high+low)/2) %>% arrange(mid) gene_order <- rownames(conf_intervals) #add text labels for top and bottom n genes conf_intervals$label_name <- NA_character_ n <- 2 idx <- c(1:n, (dim(conf_intervals)[1]-(n-1)):dim(conf_intervals)[1]) gene_ids <- gene_order[idx] conf_intervals$label_name[idx] <- gene_ids #the labels above can now be used as ggplot aesthetics plots_list <- plotConfidenceIntervals(conf_intervals, #mean difference in expression confidence intervals sort = F, #should genes be sorted by estimates weights = drivers$normalized_weights[rownames(conf_intervals)], pattern_name = pattern_to_weight, weights_clip = 0.99, weights_vis_norm = "none") pl1 <- plots_list[["ci_estimates_plot"]] + ggrepel::geom_label_repel(aes(label = label_name), max.overlaps = 20, force = 50) pl2 <- plots_list[["weights_heatmap"]] #now plot the weighted differences weighted_conf_intervals <- drivers$weighted_mean_ci[gene_order,] plots_list_weighted <- plotConfidenceIntervals(weighted_conf_intervals, sort = F, pattern_name = pattern_to_weight) pl3 <- plots_list_weighted[["ci_estimates_plot"]] + xlab("Difference in weighted group means") + theme(axis.title.y = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank()) cowplot::plot_grid(pl1, pl2, pl3, align = "h", rel_widths = c(1,.4, 1), ncol = 3) ``` # References