--- a
+++ b/4-Multi-Omic Integration/scripts/5.LOSO.delta.r
@@ -0,0 +1,277 @@
+# files required:
+# 1) "4_MetaG.MetaB.modules.linked.txt":  resulting data generated from script "7.MetaG.MetaB.link.r"
+# 2) "1_metaG-combined.gct" MetaG module abundance table without excluding any species
+# 3) "1_metaB.module_eigengene.txt" MetaB abundance table generated by script "1.WGCNA_metabolomics.r"
+# 4) the directory "LOSO_metaG_DR" contains all ssGSEA output files named as "speciesX-combined.gct", one for each species removed
+
+# output:
+# "LOSO_delta.spearman.r.txt" file recording delta.spearman.r for each MetaG module - MetaB module - species combination
+
+ 
+
+
+
+
+## ##########################################################################################
+## 
+##  libraries
+## 
+## ##########################################################################################
+
+
+# if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
+# if (!requireNamespace("cmapR", quietly = TRUE)) BiocManager::install("cmapR")
+
+# library(cmapR)
+library(doParallel)
+library(dplyr)
+library(data.table)
+
+
+## ##########################################################################################
+## 
+## import data 
+## 
+## ##########################################################################################
+log.file = "LOSO.deltaR.log"
+error.file = "LOSO.deltaR.error.txt"
+
+cat(paste(as.character(Sys.time()), '\n'),  file=log.file, append=T)
+cat('Importing data: \n',  file=log.file, append=T)
+
+module.pairs = fread("4_MetaG.MetaB.modules.NEU.linked.txt", select = 1, data.table = F) %>% 
+  mutate(MetaG.module = sapply(strsplit(MetaG.MetaB_modulePair,"_", fixed = T),"[[", 1)) %>%
+  mutate(MetaB.module = sapply(strsplit(MetaG.MetaB_modulePair,"_", fixed = T),"[[", 2))
+
+
+# MetaG original module data
+input.ds <- "metaG-combined.gct"
+gct.unique <- NULL
+dataset <- try(parse.gctx(input.ds), silent = T)
+if(class(dataset) != 'try-error' ){
+  
+  m <- dataset@mat
+  gene.names <- dataset@rid
+  gene.descs <- dataset@rdesc
+  sample.names <- dataset@cid
+  sample.descs <- dataset@cdesc
+  
+} else {
+  
+  ## - cmapR functions stop if ids are not unique
+  ## - import gct using readLines and make ids unique
+  if(length(grep('rid must be unique', dataset) ) > 0) {
+    gct.tmp <- readLines(input.ds)
+    #first column
+    rid <- gct.tmp %>% sub('\t.*','', .)
+    #gct version
+    ver <- rid[1]
+    #data and meta data columns
+    meta <- strsplit(gct.tmp[2], '\t') %>% unlist() %>% as.numeric()
+    if(ver=='#1.3')
+      rid.idx <- (meta[4]+3) : length(rid)
+    else
+      rid.idx <- 4:length(rid)
+    
+    #check whether ids are unique
+    if(length(rid[rid.idx]) > length(unique(rid[rid.idx]))){
+      warning('rids not unique! Making ids unique and exporting new GCT file...\n\n')
+      #make unique
+      rid[rid.idx] <- make.unique(rid[rid.idx], sep='_')
+      #other columns
+      rest <- gct.tmp %>% sub('.*?\t','', .)
+      rest[1] <- ''
+      gct.tmp2 <- paste(rid, rest, sep='\t') 
+      gct.tmp2[1] <-  sub('\t.*','',gct.tmp2[1])
+      
+      #export
+      gct.unique <- sub('\\.gct', '_unique.gct', input.ds)
+      writeLines(gct.tmp2, con=gct.unique)
+      
+      #import using cmapR functions
+      dataset <- parse.gctx(fname = gct.unique)
+      
+      ## extract data 
+      m <- dataset@mat
+      gene.names <- sub('_[0-9]{1,5}$', '', dataset@rid)
+      gene.descs <- dataset@rdesc
+      sample.names <- dataset@cid
+      sample.descs <- dataset@cdesc
+    }
+    
+  } else { #end if 'rid not unique'
+    
+    ########################################################
+    ## display a more detailed error message if the import 
+    ## failed due to other reasons than redundant 'rid'
+    stop("\n\nError importing GCT file using 'cmapR::parse.gctx()'. Possible reasons:\n\n1) Please check whether you have the latest version of the 'cmapR' installed. Due to submission to Bioconductor the cmap team changed some naming conventions, e.g 'parse.gctx()' has been renamed to 'parse.gctx()'.\n2) The GCT file doesn't seem to be in the correct format! Please see take a look at https://clue.io/connectopedia/gct_format for details about GCT format.
+\nError message returned by 'cmapR::parse.gctx()':\n\n", dataset, '\n\n')
+  } 
+} #end if try-error
+MetaG.mod.bf <- m %>% t() %>% data.frame()
+#head(MetaG.mod.bf[,1:6])
+
+
+
+# MetaB module data
+m1 <- fread("metaB.module_eigengene.txt",data.table = F) %>%
+  dplyr::filter(!grepl("#",`#NAME`,fixed=T)) 
+m1[-1] <- sapply(m1[-1], as.numeric)
+m1 <- m1 %>% tibble::column_to_rownames("#NAME")
+feature.abb_df1 <- cbind.data.frame(feature = rownames(m1),
+                                    abb = paste("feature",seq(1,nrow(m1),1),sep = ""),
+                                    stringsAsFactors = F)
+rownames(m1) <- sapply(rownames(m1), function(x) feature.abb_df1$abb[which(feature.abb_df1$feature == x)])
+m1 <- t(m1) %>% as.data.frame(stringsAsFactors=F)
+colnames(m1) <- sapply(colnames(m1), function(x) feature.abb_df1$feature[which(feature.abb_df1$abb == x)])
+
+MetaB.mod <- m1
+#head(MetaB.mod[,1:6])
+
+#all(rownames(MetaB.mod) %in% rownames(MetaG.mod.bf))
+
+MetaG.B.bf <- merge(MetaG.mod.bf, MetaB.mod, by=0) 
+MetaG.B.bf <- MetaG.B.bf[complete.cases(MetaG.B.bf),] # remove NA because otherwise spearman r might be NA
+
+## ##########################################################################################
+## 
+## LOSO analysis 
+## 
+## ##########################################################################################
+cat('Starting LOSO analysis: \n',  file=log.file, append=T)
+
+MetaG.mod.after.dir = "LOSO_metaG_DR"
+AllSpecies = sub("\\-combined\\.gct","",  basename(list.files(MetaG.mod.after.dir,full.names = T, pattern = "combined.gct"))) 
+if(length(AllSpecies) == 0) stop("\n\nError importing GCT file, GCT file names must contain '-combined.gct' as exported by ssGSEA2 function\n\n")
+
+
+# combination of the module pairs and species 
+allCombs <- expand.grid(paste(module.pairs$MetaG.module, module.pairs$MetaB.module, sep = "_"),
+                        AllSpecies)
+
+
+
+moduelPair.res <- NULL
+for( i in c(1:nrow(allCombs))){
+  mdp = strsplit(as.character(allCombs$Var1[i]), "_", fixed = T)[[1]]
+  
+  metaG.md = sub("MetaG\\.", "", mdp[1]  )
+  metaB.md = sub("MetaB\\.", "",  mdp[2])
+  
+  r.bf = cor(MetaG.B.bf[,metaG.md], MetaG.B.bf[,metaB.md], method = "spearman")
+  
+  
+  spc = as.character(allCombs$Var2[i])
+  
+  spc.mod.files = list.files(MetaG.mod.after.dir, full.names = T, pattern = spc)
+  spc.mod.file = spc.mod.files[grepl("-combined.gct", spc.mod.files, fixed = T)]
+  
+  # cat(paste('module pairs: ', paste(mdp,collapse = " | "), " \nSpecies: ",spc,"\n", sep=""), file=log.file, append=T)
+  if(i %% 100 ==  0) cat(paste("-----------progress: ",i," out of ", nrow(allCombs), " combinations.", sep = ""), file=log.file, append=T)
+  
+  m <- NA
+  readGctx_success <- F
+  if(T){
+    input.ds <- spc.mod.file
+    gct.unique <- NULL
+    dataset <- try(parse.gctx(input.ds), silent = T)
+    if(class(dataset) != 'try-error' ){
+      
+      m <- dataset@mat
+      gene.names <- dataset@rid
+      gene.descs <- dataset@rdesc
+      sample.names <- dataset@cid
+      sample.descs <- dataset@cdesc
+      readGctx_success <- T
+    } else {
+      
+      ## - cmapR functions stop if ids are not unique
+      ## - import gct using readLines and make ids unique
+      if(length(grep('rid must be unique', dataset) ) > 0) {
+        gct.tmp <- readLines(input.ds)
+        #first column
+        #rid <- gct.tmp %>% sub('\t.*','', .)
+        rid <-  sub('\t.*','', gct.tmp)
+        #gct version
+        ver <- rid[1]
+        #data and meta data columns
+        meta <- as.numeric( unlist( strsplit(gct.tmp[2], '\t'))) 
+        if(ver=='#1.3'){rid.idx <- (meta[4]+3) : length(rid)} else {rid.idx <- 4:length(rid)}
+        
+        #check whether ids are unique
+        if(length(rid[rid.idx]) > length(unique(rid[rid.idx]))){
+          warning('rids not unique! Making ids unique and exporting new GCT file...\n\n')
+          #make unique
+          rid[rid.idx] <- make.unique(rid[rid.idx], sep='_')
+          #other columns
+          rest <- sub('.*?\t','', gct.tmp)
+          rest[1] <- ''
+          gct.tmp2 <- paste(rid, rest, sep='\t') 
+          gct.tmp2[1] <-  sub('\t.*','',gct.tmp2[1])
+          
+          #export
+          gct.unique <- sub('\\.gct', '_unique.gct', input.ds)
+          writeLines(gct.tmp2, con=gct.unique)
+          
+          #import using cmapR functions
+          dataset <- parse.gctx(fname = gct.unique)
+          
+          ## extract data 
+          m <- dataset@mat
+          gene.names <- sub('_[0-9]{1,5}$', '', dataset@rid)
+          gene.descs <- dataset@rdesc
+          sample.names <- dataset@cid
+          sample.descs <- dataset@cdesc
+          
+          readGctx_success <- T
+          
+          
+        }
+        
+      } else { #end if 'rid not unique'
+        
+        ########################################################
+        ## display a more detailed error message if the import 
+        ## failed due to other reasons than redundant 'rid'
+        cat(paste("\n\nError importing GCT file using 'cmapR::parse.gctx()' for species: ", spc,
+                  "\nPossible reasons:\n\n1) Please check whether you have the latest version of the 'cmapR' installed. Due to submission to Bioconductor the cmap team changed some naming conventions, e.g 'parse.gctx()' has been renamed to 'parse.gctx()'.\n2) The GCT file doesn't seem to be in the correct format! Please see take a look at https://clue.io/connectopedia/gct_format for details about GCT format.
+\nError message returned by 'cmapR::parse.gctx()':\n\n", sep=""),  file = error.file, append=T)
+        
+        # next
+        
+        
+      } 
+    } 
+    
+    
+  }# wrap the script of importing gct file  
+  
+  if(!readGctx_success){
+    delta.r_vec = c(metaG.md, metaB.md, spc, NA)
+    names(delta.r_vec) <- c("MetaG.module","MetaB.module","Species","delta.spearman.r")
+  }else{
+    
+    spc.metaG.mod.aft <- data.frame(t(m))
+    MetaG.B.aft <- merge(spc.metaG.mod.aft, MetaB.mod, by=0) 
+    MetaG.B.aft <- MetaG.B.aft[complete.cases(MetaG.B.aft),] # remove NA because otherwise spearman r might be NA
+    r.aft = cor(MetaG.B.aft[,metaG.md], MetaG.B.aft[,metaB.md], method = "spearman")
+    
+    delta.r = r.aft - r.bf
+    
+    delta.r_vec = c(metaG.md, metaB.md, spc, delta.r)
+    names(delta.r_vec) <- c("MetaG.module","MetaB.module","Species","delta.spearman.r")
+  }
+  
+  #delta.r_vec
+  
+  moduelPair.res <- bind_rows(moduelPair.res, delta.r_vec)
+}
+
+
+final_res <- moduelPair.res %>% mutate(delta.spearman.r = as.numeric(delta.spearman.r)) %>%
+  group_by(MetaG.module, MetaB.module) %>% mutate(zscore = (delta.spearman.r - mean(delta.spearman.r, na.rm=T))/sd(delta.spearman.r, na.rm=T))
+
+
+cat('Saving results as file: \n',  file=log.file, append=T)
+write.table(final_res, file = "5_LOSO_NEU.delta.spearman.r.txt", sep = "\t", quote = F, row.names = F)
+