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 and Sharma, et al.

For automatic Bioconductor package installation, start R, and run:

`BiocManager::install("projectR")`

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, Wu, and Fyfe (2009) . 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 generic projectR function is executed as follows:

```
library(projectR)
projectR(data, loadings, dataNames=NULL, loadingsNames=NULL, NP = NULL, full = false)
```

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`

`dataNames`

`rownames(data)`

. Order of names in vector must match order of rows in data.`loadingsNames`

`rownames(loadings)`

. Order of names in vector must match order of rows in loadings.`NP`

`full`

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).

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 *limma*.

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`

.```
# 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=""))
```

```
# 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"]])
```

```
## [1] "93 row names matched between data and loadings"
## [1] "Updated dimension of data: 93 9"
```

```
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=""))
```

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 Wang, Kossenkov, and Ochs (2006), the matrices \({\bf{A}}\) and \({\bf{P}}\) in CoGAPS are assumed to have the atomic prior described in Sibisi and Skilling (1997). 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 *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)
```

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`

`dataNames`

`loadingsNames`

`NP`

`full`

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 *limma*.

```
# get data
library(projectR)
AP <- get(data("AP.RNAseq6l3c3t")) #CoGAPS run data
AP <- AP$Amean
# heatmap of gene weights for CoGAPs patterns
library(gplots)
```

```
##
## Attaching package: 'gplots'
```

```
## The following object is masked from 'package:stats':
##
## lowess
```

```
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")
)
```

```
# data to project into PCs from RNAseq6l3c3t expression data
library(projectR)
data('p.ESepiGen4c1l4')
```

`## Warning in data("p.ESepiGen4c1l4"): data set 'p.ESepiGen4c1l4' not found`

```
data('p.RNAseq6l3c3t')
NMF2ESepi <- projectR(p.ESepiGen4c1l$mRNA.Seq,loadings=AP,full=TRUE,
dataNames=map.ESepiGen4c1l[["GeneSymbols"]])
```

```
## [1] "93 row names matched between data and loadings"
## [1] "Updated dimension of data: 93 9"
```

```
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=""))
```

```
## $x
## [1] "Projected PC1 (18.32% of varience)"
##
## $y
## [1] "Projected PC2 (17.12% of varience)"
##
## $title
## [1] "Encode RNAseq in target PC1 & PC2"
##
## attr(,"class")
## [1] "labels"
```

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`

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)
```

The inputs that must be set each time are the clusters and data.

The arguments are as follows:

** clusters** a clustering object

`NP`

`data`

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`

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)
```

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`

`pval`

`full`

`k`

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 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.

```
library(projectR)
correlateR(genes = NA, dat = NA, threshtype = "R", threshold = 0.7, absR = FALSE, ...)
```

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`

`threshtype`

`threshold`

`absR`

`...`

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.

`correlateR`

objects.```
# 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))
```

`## Using indx as id variables`

```
library(gplots)
library(ggplot2)
library(viridis)
```

`## Loading required package: viridisLite`

```
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()
)
```

```
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```
# 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)
```

```
## [1] "9 row names matched between data and loadings"
## [1] "Updated dimension of data: 9 9"
```

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)
```

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`

`loadings`

`loadingsNames`

`pattern_name`

`pvalue`

`display`

`normalize_pattern`

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.

```
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
```

```
## [1] "2996 row names matched between datasets"
## [1] "2996"
## [1] "Updated dimension of data: 2996"
```

```
conf_intervals <- drivers$mean_ci[drivers$significant_genes,]
str(conf_intervals)
```

```
## 'data.frame': 253 obs. of 2 variables:
## $ low : num 1.86 0.158 -0.562 -0.756 0.155 ...
## $ high: num 2.03943 0.26729 -0.00197 -0.18521 0.23239 ...
```

The arguments for plotConfidenceIntervals are:

** confidence_intervals** A dataframe of features x estimates

`interval_name`

`pattern_name`

`sort`

`genes`

`weights`

`weights_clip`

`weights_vis_norm`

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.

```
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)
```

```
## Warning: Removed 249 rows containing missing values
## (`geom_label_repel()`).
```