Add partial correlations for batch correction if needed

Author: mhurtado13

R/cell_deconvolution.R

@@ -579,12 +579,28 @@ compute.cell.types = function(data, cells_extra = NULL){
 #' - Highly correlated features found
 #' - Cell type name
 #'
-removeCorrelatedFeatures <- function(data, threshold, name, n_seed, corr_method = "spearman") {
+removeCorrelatedFeatures <- function(data, threshold, name, n_seed, corr_method = "spearman", batch = NULL) {
 
   features_high_corr = c()
   cell_name = c()
+
   # Compute correlation matrix
-  corr_matrix <- stats::cor(data, method = corr_method)
+  if(is.null(batch)){
+    corr_matrix <- stats::cor(data, method = corr_method)
+  } else {
+    if(is.factor(batch) || is.character(batch)) batch <- as.numeric(as.factor(batch))
+    # Pairwise partial correlation
+    corr_matrix <- matrix(NA, ncol=ncol(data), nrow=ncol(data))
+    colnames(corr_matrix) <- colnames(data)
+    rownames(corr_matrix) <- colnames(data)
+    for(i in 1:(ncol(data)-1)){
+      for(j in (i+1):ncol(data)){
+        pc <- ppcor::pcor.test(data[[i]], data[[j]], batch)
+        corr_matrix[i,j] <- pc$estimate
+        corr_matrix[j,i] <- pc$estimate
+      }
+    }
+  }
   # Find highly correlated features
   contador = 1
   while(nrow(corr_matrix)>0){
@@ -667,7 +683,7 @@ remove_subgroups = function(groups){
 #' - Cell subgroups obtained by proportionality correlation
 #' - Discard cell features either because of low variance or high zero number
 #'
-compute_subgroups = function(deconvolution, thres_corr, corr_type, file_name){
+compute_subgroups = function(deconvolution, thres_corr, corr_type, file_name, batch = NULL){
   data = data.frame(deconvolution)
   cell_subgroups = list()
   #cell_groups_similarity = list()
@@ -757,7 +773,7 @@ compute_subgroups = function(deconvolution, thres_corr, corr_type, file_name){
     terminate = FALSE
     iteration = 1
     while (terminate == FALSE) {
-      corr_df <- correlation(data.matrix(data), corr_type = corr_type)
+      corr_df <- correlation(data, corr_type = corr_type, batch = batch)
       vec = colnames(data)
       indice = 1
       subgroup = list()
@@ -876,28 +892,46 @@ compute_subgroups = function(deconvolution, thres_corr, corr_type, file_name){
 #'
 #' @return Dataframe containing all significant correlations (pvalue < 0.05)
 #'
-correlation <- function(data, corr_type = "spearman") {
-
-  M <- Hmisc::rcorr(as.matrix(data), type = corr_type)
-  Mdf <- purrr::map(M[c("r", "P", "n")], ~data.frame(.x))
-
-  corr_df = Mdf %>%
-    purrr::map(~tibble::rownames_to_column(.x, var="measure1")) %>%
-    purrr::map(~tidyr::pivot_longer(.x, -measure1, names_to = "measure2")) %>%
-    dplyr::bind_rows(.id = "id") %>%
-    tidyr::pivot_wider(names_from = id, values_from = value) %>%
-    dplyr::rename(p = P) %>%
-    dplyr::mutate(sig_p = ifelse(p < .05, T, F),
-                  p_if_sig = ifelse(sig_p, p, NA),
-                  r_if_sig = ifelse(sig_p, r, NA))
+correlation <- function(data, corr_type = "spearman", batch = NULL) {
+  if (!is.null(batch)) {
+    # Convert batch to numeric if factor or character
+    if(is.factor(batch) || is.character(batch)){
+      batch <- as.numeric(as.factor(batch))
+    }
 
-  corr_df = stats::na.omit(corr_df)  #remove the ones that are the same Tfeatures (pval = NA)
-  corr_df <- corr_df[which(corr_df$sig_p==T),]  #remove not significant
-  corr_df <- corr_df[order(corr_df$r, decreasing = T),]
-  corr_df$AbsR =  abs(corr_df$r)
+    # Compute all pairwise partial correlations controlling for batch
+    vec <- colnames(data)
+    corr_df <- data.frame(measure1 = character(0), measure2 = character(0),
+                          r = numeric(0), p = numeric(0))
+    for(i in 1:(length(vec)-1)){
+      for(j in (i+1):length(vec)){
+        pc <- ppcor::pcor.test(data[[vec[i]]], data[[vec[j]]], batch, method = corr_type)
+        corr_df <- rbind(corr_df, data.frame(measure1 = vec[i],
+                                             measure2 = vec[j],
+                                             r = pc$estimate,
+                                             p = pc$p.value))
+      }
+    }
+  } else {
+    # Original correlation using Hmisc::rcorr
+    M <- Hmisc::rcorr(data.matrix(data), type = corr_type)
+    Mdf <- purrr::map(M[c("r", "P", "n")], ~data.frame(.x))
+    corr_df <- Mdf %>%
+      purrr::map(~tibble::rownames_to_column(.x, var = "measure1")) %>%
+      purrr::map(~tidyr::pivot_longer(.x, -measure1, names_to = "measure2")) %>%
+      dplyr::bind_rows(.id = "id") %>%
+      tidyr::pivot_wider(names_from = id, values_from = value) %>%
+      dplyr::rename(p = P) %>%
+      dplyr::mutate(sig_p = ifelse(p < 0.05, TRUE, FALSE),
+                    p_if_sig = ifelse(sig_p, p, NA),
+                    r_if_sig = ifelse(sig_p, r, NA))
+    corr_df <- stats::na.omit(corr_df)
+    corr_df <- corr_df[which(corr_df$sig_p == TRUE), ]
+    corr_df <- corr_df[order(corr_df$r, decreasing = TRUE), ]
+  }
 
+  corr_df$AbsR <- abs(corr_df$r)
   return(corr_df)
-
 }
 
 
@@ -970,7 +1004,7 @@ remove_low_variance <- function(data, plot = FALSE) {
 #'
 #' processed_deconvolution = compute.deconvolution.analysis(deconvolution, cells_extra = "mesenchymal")
 #'
-compute.deconvolution.analysis <- function(deconvolution, corr = 0.7, corr_type = "pearson", seed = NULL, cells_extra = NULL, file_name = NULL, return = FALSE, verbose = FALSE){
+compute.deconvolution.analysis <- function(deconvolution, corr = 0.7, corr_type = "spearman", seed = NULL, batch = NULL, cells_extra = NULL, file_name = NULL, return = FALSE, verbose = FALSE){
   deconvolution.mat = deconvolution
 
   #####Unsupervised filtering
@@ -1020,7 +1054,7 @@ compute.deconvolution.analysis <- function(deconvolution, corr = 0.7, corr_type
     if(is.null(ncol(data))==T){
       cells[[i]] = data
     }else if(ncol(data)>1){
-      data = removeCorrelatedFeatures(data, 0.9, names(cells)[i], seed, corr_method = corr_type)
+      data = removeCorrelatedFeatures(data, 0.9, names(cells)[i], seed, corr_method = corr_type, batch = batch)
       cells[[i]] = data[[1]]
       if(length(data[[2]])>0 && is.null(data[[3]])==F){
         features_high_corr[[j]] = data[[2]]
@@ -1036,7 +1070,7 @@ compute.deconvolution.analysis <- function(deconvolution, corr = 0.7, corr_type
   #groups_similarity = list()
   groups_discard = list()
   for (i in 1:length(cells)) {
-    x = compute_subgroups(cells[[i]], file_name = names(cells)[i], thres_corr = corr, corr_type = corr_type)
+    x = compute_subgroups(cells[[i]], file_name = names(cells)[i], thres_corr = corr, corr_type = corr_type, batch = batch)
     res = c(res, x[1])
     groups = c(groups, x[2])
     #groups_similarity = c(groups_similarity, x[3])
@@ -2698,7 +2732,7 @@ prepare_multideconv_folds <- function(
 #' @importFrom factoextra fviz_nbclust hcut
 #'
 #' @export
-deconvolution_dictionary = function(deconv_subgroups, pathway_matrix){
+deconvolution_dictionary = function(deconv_subgroups, pathway_matrix, batch_id = NULL){
   pathway_matrix = pathway_matrix[,!colnames(pathway_matrix)%in%c("Androgen", "Estrogen")]
   cell_subgroups = deconv_subgroups[["Deconvolution subgroups per cell types"]]
   cell_clusters = list()
@@ -2709,6 +2743,7 @@ deconvolution_dictionary = function(deconv_subgroups, pathway_matrix){
     deconv_subgroups[["Deconvolution matrix"]],
     pathway_matrix,
     return = TRUE,
+    batch = batch_id,
     plot = FALSE
   )
 
@@ -2749,7 +2784,7 @@ deconvolution_dictionary = function(deconv_subgroups, pathway_matrix){
       rownames(cell_subgroups[[cell]]) <- rownames(deconv_subgroups[["Deconvolution matrix"]])
 
       #Compute module correlation between cell deconvolution features and PROGENy pathways
-      x <- CellTFusion::compute.modules.relationship(cell_subgroups[[cell]], pathway_matrix, return = TRUE, plot = FALSE)
+      x <- CellTFusion::compute.modules.relationship(cell_subgroups[[cell]], pathway_matrix, return = TRUE, batch = batch_id, plot = FALSE)
       corr_matrix <- data.frame(x[[1]])
 
       for (k in seq_along(clusters_global)) {

vignettes/Results/Distribution_variances_deconvolution.pdf

@@ -0,0 +1,31 @@
+"","Androgen","EGFR","Estrogen","Hypoxia","JAK-STAT","MAPK","NFkB","PI3K","TGFb","TNFa","Trail","VEGF","WNT","p53"
+"SAM7f0d9cc7f001",0.501202327084813,0.0481403051897762,0.175521670439651,-0.932496396633032,-0.092175193044574,0.006307947520786,1.72858261966861,-1.00961757184729,0.356033124198143,0.335880865957153,0.700174564846386,-0.362059214320684,-0.780070024826196,0.300823867976749
+"SAM4305ab968b90",-0.615956795831225,-0.764123264923606,-0.35461161550043,0.34357586506766,-1.75885764817319,0.626192147223636,-1.91475574555241,-0.0655723168766853,-2.42759525527634,-0.911500034266245,-0.853026736109109,1.71082124358451,-0.228807172772938,-0.304542276410043
+"SAMcf018fee2acd",0.461030774150299,0.318016962563996,-1.13732028467355,-0.747583165637903,0.273875533102506,-0.330634433399792,-0.0648672645091573,-2.4702705512752,-0.560645591002033,-0.164819955734901,2.03439801789168,-0.105596970398538,0.0298227548846749,1.67823428774881
+"SAMcc4675f394a1",0.708340651336115,-0.229449109739751,-0.879667167910258,-0.45101467950415,-0.478251917650965,-1.96120331038764,-1.12802415213619,-0.67496164305272,-0.0431076957260577,0.967063015289295,0.972153464158632,-0.7958874309622,0.705154625716835,0.908103786422819
+"SAM49f9b2e57aa5",-0.0525077324829043,0.0238999094748491,0.0906155535213645,-0.168563350893111,-1.33216402626442,0.785099310064374,0.135988617767126,0.144188778629872,-0.550286778322233,-1.1834230564947,-0.505437690984083,-0.745678639127147,0.114206152274328,0.101036453453709
+"SAM2e7aa8fa0ab3",-1.73112856733671,-1.32886292288447,-0.893879815820183,0.860734530704065,0.0393097128572464,-1.48801041505734,0.140247458574208,-0.00671281061819084,-0.131494729920062,-0.645402621171688,-0.118102017942277,0.932505964397447,1.44098140356897,0.147146159738315
+"SAMdf3e42c8672a",0.43014591487667,1.50829424022276,1.26294632729694,1.86471783642638,0.608844028964011,0.852366602652316,0.614438532853088,1.23604028204655,-0.383480853524494,0.773536742469119,-1.01187649153445,0.132268903803536,0.389648404207173,-0.973536700945024
+"SAMd027124354ce",-0.486884003961032,-0.428957153119357,1.98277679446853,1.86939270667137,-0.119402626043375,0.356331595115521,-0.400718670477007,1.86640721399838,-0.216136216442976,-0.418526520901296,-1.4644917015642,-1.16393174803946,-0.708565319114249,-1.42808741111167
+"SAMe7bf6c015192",-0.157585425042266,0.480212391069881,-0.340705028111347,0.173828516431692,0.704014141141414,-0.59936752823094,1.00719498334749,0.11245193624228,0.171285571365987,0.535999674146456,1.56538205084502,-0.542378131457993,-0.175758880239431,0.107128265822987
+"SAM6dd7ad1d797d",-1.47427184981032,1.79440079468457,2.26574084088425,-0.860587951578688,2.00813973568238,0.43275335173471,1.55749987974988,2.28781574290141,1.51940072767977,1.34949069693871,0.879019807989398,2.12232258624472,0.730157149943444,-0.927286654275007
+"SAM18039827e1b9",0.321322086834494,-0.665161830430496,-1.03826786328342,0.510819735031411,-1.1764099036649,0.822900564848599,-0.392770838285755,-0.688390180130047,-0.896394686890645,-0.477443884862368,-0.135611621772851,0.782411736750316,0.904368113894293,-0.166774620761731
+"SAMc692536a795a",0.0938145999731334,0.497648256011942,-0.273220653307457,-1.20979070501855,0.630422552924659,-1.81272168255174,0.128854795562907,-0.913699831159416,0.992384773113748,0.0333012258186793,0.863487105111853,-0.887947352337348,-0.861004082819969,-0.159079180815172
+"SAM9a2cf3c06fb3",0.916005056553408,0.829181700419048,1.33957559794147,-1.49402011184786,0.885634987865809,1.696838919168,1.10004137523273,0.54994450402526,0.698479439659968,0.334052880894688,-0.0439206634017195,-0.343133715879574,-1.66380139271114,-0.604120510327014
+"SAM557dde1b9f3e",-0.160642804981361,-0.157766757664521,-0.165607966915137,-0.352357811819535,0.265411967985976,-0.48378020407324,0.285948879001779,-1.12436359674278,1.51699274078716,0.299899958830839,0.582649165016261,-1.54764236304185,-0.866855528741966,1.16533724555277
+"SAM23aa15d4a0b0",-0.355653487639992,-0.056341805112318,-0.329126892949966,0.160711725408663,-0.351525838152739,-0.12170952922031,0.408782075163163,-0.130363527073975,0.0523969167961646,0.0056954098903415,-0.366007594142631,-1.25623072490406,-0.620772803207441,-0.119459082316105
+"SAM468a9e1dc821",-0.477300873528647,1.19256307568451,0.492701904522933,0.372158931345949,-0.145893584338381,-0.221867335635066,-0.762147531089719,0.825747803076648,0.423882063021743,0.357050271192651,0.0402716195561987,0.622956208079934,1.43236705033637,1.99158354860332
+"SAM81b71522417a",0.0544883497043489,-1.32971491799097,-0.660868520973388,0.433347330904573,-1.33077010226212,1.1985086079805,-0.974459502569838,0.783474821470852,-1.12507863068551,-1.72456219774716,-1.08080386227248,0.158063209548548,0.693489088696513,-1.25411917677895
+"SAMb963dda93cfd",-0.643459810509838,0.148026566895117,1.72451999998291,-0.545526705756573,0.534839152108086,1.4608125792785,0.484559061423702,-0.514787356653263,0.947527110646201,0.648747169980428,-1.54369405459326,0.623893973291663,-0.639697939599069,-1.95547214529753
+"SAMbcbc7957c264",-0.174732882417787,-0.0393808998029008,-0.592292887909991,0.662737604947458,-0.774706391594431,1.11787584346812,-0.273487590674765,0.545804077223036,-0.787421416928816,-0.616154052650002,-1.46762946062021,0.272048917698994,0.959932218111358,0.320611056374883
+"SAM7fb6987514a4",0.470985796465443,0.0740636083523108,1.62448861333991,-1.11908515027197,1.61318184480753,0.521374016790038,2.19197584629197,0.179872020878216,0.90971533400136,0.162396410691025,0.0933355293793949,-0.561342652872902,-0.485934243668331,-0.198903211565923
+"SAM63405b04ab2d",-0.58788115549129,-0.211171176694518,-0.0206968055120434,-0.396123870853721,1.26815604540648,-0.724156856398065,1.10306397758686,0.458062088062592,1.53531414405831,-0.504371403360444,0.578768300042518,0.00159980721684351,1.28023086405201,0.0746110206062119
+"SAM18bc1078bc15",-0.184293895967879,-1.68277377828351,-1.7606494669777,-0.874778951455845,0.293935720282401,-1.08310248274646,-0.887367898031778,-0.823903811307423,-0.462081468625133,0.211684930518473,0.282733858437624,0.532609352963711,-0.310513298830893,-1.49426337095219
+"SAMd1bd63734394",-0.565298257096645,1.08661471126978,-0.485200246285454,1.29058643702759,1.78837280623672,0.603692225938965,-0.893361424408252,0.126834330635557,-2.01005941596333,1.38157753603666,-1.6347505604278,1.33205758021056,-1.94380452839539,-0.592036666569154
+"SAMe9ae8beb82fa",0.120273152887152,-0.170109669487559,-0.640471631381582,-0.209304293360986,1.33825119190994,-1.94844381470598,-1.48240370804962,0.64615698317165,0.77436211618826,1.21575185737726,1.5494109545972,-0.612171559111438,0.557691253855057,1.26651208718389
+"SAMba7176afe070",-0.399882917080415,-1.10739898371845,-1.00562804494344,-1.40230728624507,-0.591437046382024,0.210802373831932,-0.47035401224863,-1.37396633930082,-1.74970426601584,-2.20401546261127,-0.0985335132426545,1.51774085097194,1.39003193828094,1.71680918304397
+"SAMbe83eae4026e",0.578483268945518,-0.51037762634421,0.0419348423737584,1.60876046964542,-0.398319422588754,0.514561849621346,0.967407358696127,1.15318983769793,0.0323499103318224,0.282973289497091,0.0440981058279232,0.226002540644956,-0.751982746552687,0.151865490508125
+"SAMe5bc41772bc9",0.0312858969835789,2.94422026557683,0.348277890635478,0.557415407184958,-0.714397168007047,-0.424581421278836,-0.994142628711901,0.438397783232794,0.648194169770616,2.30475514956641,1.41350147966279,-0.188605773576151,0.313759026621819,1.28721028770409
+"SAM23095936e611",-0.463692354845003,-0.688803022704654,-0.687603441501235,1.48861670102848,-1.34732741133346,1.0937423964987,-0.829563467867355,-0.805166847228467,0.216243984112538,-0.183531501438348,-1.22708857338757,0.616826952682477,0.500818713225474,-0.283946590514282
+"SAM7114d99032ec",-0.446966488478995,-0.746776369921522,-0.276314010811397,0.17056631746865,-0.808704243218696,-0.409747462819467,-0.852614247228216,0.212263010155344,0.798931008507836,-0.15555031884445,0.17299982103872,-0.0609469795818617,0.927038531139419,0.256001758251038
+"SAMdb3f50c9129c",4.29076142670735,-0.828113498592561,0.193032309360767,-1.60442968441734,-0.832046898556087,-0.69083385523118,0.466453220920956,-0.964874830182083,-0.250006128916178,-2.01055607501239,-0.22140930240631,-2.41057657247895,-2.33212932732897,-1.0113869003519