From 1d4e0dd75e621ba1fc9ef105e053ee63981e5471 Mon Sep 17 00:00:00 2001 From: "soel.micheletti" Date: Wed, 25 Oct 2023 02:09:30 +0200 Subject: [PATCH 1/4] Feature: included pcor decomposition natively in cobra --- R/COBRA.R | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/R/COBRA.R b/R/COBRA.R index 14276b72..ec9b357c 100644 --- a/R/COBRA.R +++ b/R/COBRA.R @@ -10,7 +10,7 @@ #' Inputs: #' @param X : design matrix of size (n, q), n = number of samples, q = number of covariates #' @param expressionData : gene expression as a matrix of size (g, n), g = number of genes -#' @param standardize : boolean flag to standardize the gene expression as a pre-processing step +#' @param method : if pearson, the decomposition of the co-expression matrix is compouted. If pcor, cobra decompose the (regularized) partial co-expression #' #' Outputs: #' @return psi : impact of each covariate on the eigenvalues as a matrix of size (q, n) @@ -31,21 +31,26 @@ #' #' @export -cobra <- function(X, expressionData, standardize=T){ +cobra <- function(X, expressionData, method = "pearson"){ + + if(!(method %in% c("pearson", "pcor"))){ + stop("Only Pearson and pcor methods are supported. Please make sure to provide a valid method argument.") + } + numSamples <- ncol(expressionData) N <- min(ncol(expressionData),nrow(expressionData)) - - if (standardize){ + C <- 0 + if(method == "pearson"){ G_star <- expressionData-rowMeans(expressionData) G <- (G_star/sqrt(rowSums(G_star^2))) G <- as.matrix(G) - } else { - G <- expressionData - G <- (G/sqrt(rowSums(G^2))) - G <- as.matrix(G) + C <- tcrossprod(G) + } + if(method == "pcor"){ + C <- pcor.shrink(t(expressionData)) } - eigenG <- rARPACK::eigs_sym(tcrossprod(G),N) + eigenG <- rARPACK::eigs_sym(C,N) Q <- eigenG$vectors D <- diag(eigenG$values) From 02313d8d3e34e7602956225e15918ed236c9a0b6 Mon Sep 17 00:00:00 2001 From: "soel.micheletti" Date: Fri, 27 Oct 2023 17:27:03 +0200 Subject: [PATCH 2/4] Feature: included pcor decomposition natively in cobra --- DESCRIPTION | 3 ++- R/COBRA.R | 26 +++++++++++++++++--------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d10b7a98..50323bd3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -73,7 +73,8 @@ Imports: cmdstanr, GeneNet, loo, - rARPACK + rARPACK, + corpcor License: GPL-3 Encoding: UTF-8 LazyData: false diff --git a/R/COBRA.R b/R/COBRA.R index ec9b357c..c5b9bced 100644 --- a/R/COBRA.R +++ b/R/COBRA.R @@ -31,9 +31,9 @@ #' #' @export -cobra <- function(X, expressionData, method = "pearson"){ +cobra <- function(X, expressionData, method = "pearson", layer2 = NULL){ - if(!(method %in% c("pearson", "pcor"))){ + if(!(method %in% c("pearson", "pcorsh", "dragon"))){ stop("Only Pearson and pcor methods are supported. Please make sure to provide a valid method argument.") } @@ -42,12 +42,20 @@ cobra <- function(X, expressionData, method = "pearson"){ C <- 0 if(method == "pearson"){ G_star <- expressionData-rowMeans(expressionData) - G <- (G_star/sqrt(rowSums(G_star^2))) - G <- as.matrix(G) - C <- tcrossprod(G) + expressionData <- (G_star/sqrt(rowSums(G_star^2))) + expressionData <- as.matrix(expressionData) + C <- tcrossprod(expressionData) } - if(method == "pcor"){ - C <- pcor.shrink(t(expressionData)) + if(method == "dragon"){ + if(is.null(layer2)){ + method <- "pcorsh" + } + else{ + C <- dragon(t(expressionData), t(layer2), pval=FALSE)$cov + } + } + if(method == "pcorsh"){ + C <- matrix(as.numeric(pcor.shrink(t(expressionData))), dim(expressionData)[1], dim(expressionData)[1]) } eigenG <- rARPACK::eigs_sym(C,N) @@ -57,11 +65,11 @@ cobra <- function(X, expressionData, method = "pearson"){ hatmat <- ginv(crossprod(X))%*%t(X) Qinv <- ginv(Q) - QinvG <- Qinv%*%(G) + QinvG <- Qinv%*%(expressionData) est <- t(sapply(seq_len(nrow(hatmat)), function(hatmatRow){ diag(QinvG%*%(numSamples*diag(hatmat[hatmatRow,]))%*%t(QinvG)) })) - list(psi=est, Q=Q, D=eigenG$values, G=G) + list(psi=est, Q=Q, D=eigenG$values, G=expressionData) } \ No newline at end of file From 82482193cb6554d96efcf3dcf67492082273bf6d Mon Sep 17 00:00:00 2001 From: "soel.micheletti" Date: Mon, 30 Oct 2023 23:48:41 +0100 Subject: [PATCH 3/4] Feature: included pcor decomposition natively in cobra --- R/COBRA.R | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/R/COBRA.R b/R/COBRA.R index c5b9bced..3bbba4c2 100644 --- a/R/COBRA.R +++ b/R/COBRA.R @@ -31,14 +31,16 @@ #' #' @export -cobra <- function(X, expressionData, method = "pearson", layer2 = NULL){ +cobra <- function(X, expressionData, method = "pearson"){ if(!(method %in% c("pearson", "pcorsh", "dragon"))){ stop("Only Pearson and pcor methods are supported. Please make sure to provide a valid method argument.") } numSamples <- ncol(expressionData) - N <- min(ncol(expressionData),nrow(expressionData)) + if(method != "dragon"){ + N <- min(ncol(expressionData),nrow(expressionData)) + } C <- 0 if(method == "pearson"){ G_star <- expressionData-rowMeans(expressionData) @@ -47,12 +49,17 @@ cobra <- function(X, expressionData, method = "pearson", layer2 = NULL){ C <- tcrossprod(expressionData) } if(method == "dragon"){ - if(is.null(layer2)){ - method <- "pcorsh" + if(length(expressionData) != 2){ + stop("Dragon needs two layers, please provide a list of two matrices in expressionData or consider using the pcorsh argument") } - else{ - C <- dragon(t(expressionData), t(layer2), pval=FALSE)$cov + if(dim(expressionData[[1]])[2] != dim(expressionData[[2]])[2]){ + stop("Dragon layers in expressionData list must have the same number of samples") } + C <- dragon(t(expressionData[[1]]), layer2 = t(expressionData[[2]]), pval=FALSE)$cov + N <- dim(expressionData[[1]])[2] + numSamples <- N + expressionData <- rbind(expressionData[[1]], expressionData[[2]]) + print(dim(expressionData)) } if(method == "pcorsh"){ C <- matrix(as.numeric(pcor.shrink(t(expressionData))), dim(expressionData)[1], dim(expressionData)[1]) From 894ba78610e801925fdd0c206b832c8143102284 Mon Sep 17 00:00:00 2001 From: Marouen Date: Mon, 30 Oct 2023 20:47:47 -0400 Subject: [PATCH 4/4] Update COBRA.R --- R/COBRA.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/COBRA.R b/R/COBRA.R index 3bbba4c2..d07eec8f 100644 --- a/R/COBRA.R +++ b/R/COBRA.R @@ -10,7 +10,7 @@ #' Inputs: #' @param X : design matrix of size (n, q), n = number of samples, q = number of covariates #' @param expressionData : gene expression as a matrix of size (g, n), g = number of genes -#' @param method : if pearson, the decomposition of the co-expression matrix is compouted. If pcor, cobra decompose the (regularized) partial co-expression +#' @param method : if pearson, the decomposition of the co-expression matrix is compouted. If pcorsh, COBRA decomposes the regularized partial co-expression #' #' Outputs: #' @return psi : impact of each covariate on the eigenvalues as a matrix of size (q, n) @@ -79,4 +79,4 @@ cobra <- function(X, expressionData, method = "pearson"){ })) list(psi=est, Q=Q, D=eigenG$values, G=expressionData) -} \ No newline at end of file +}