R/OPNMF_correlation.R

In anbai106/NeuroStatisticR:

#' This is a pipeline to run GLM for group comparison based on the ROI measures for diffusion MRI
#' @param input_tsv the tsv file containing all the ROI measurement. The file should strictly contian these columns: "participant_id	session_id	group	age_bl	prevdemals_family_code	sex", then following each column for each ROI
#' @param output_dir the folder containing all the outputs
#' @param roi_start the first ROI name
#' @param roi_end the last ROI name
#' @param sep the separator for the tsv file, by default is '\t'
#' 
#' @return the table containing the results and also the png image for the meta-analysis of the data.
#' @export
#' @examples
#' stats_pipeline(input_tsv, output_dir, roi_start, roi_end)
#' 


stats_pipeline <- function(input_tsv, output_dir, roi_start, roi_end, sep='\t', groupPalette=c("blue","darkorange", "red"), groupPalette_CN_PS=c("blue", "darkorange")){
  
  ### Import des packages
  library(ggplot2)
  library(tidyr)
  library(car)
  library(MASS)
  library(reshape2)
  library(tableone)
  library(xtable)
  library(FactoMineR)
  library(parallel)
  library(xlsx)
  library(gridExtra)
  library(gtable)
  library(grid)
  library(lmtest)
  library(sandwich)
  library(dplyr)
  library(multcomp)
  library(lme4)
  library("factoextra")
  library(MuMIn)
  library(lsr)
  library(compute.es)
  library(factoextra)
  library("ggpubr")
  library("olsrr")
  
  ### Checkout the result folder, if not exist, creat it
  dir.create(output_dir, showWarnings = FALSE)
  
  # Th??me GGplot
  theme_set(theme_grey())
  
  options(xtable.comment = FALSE)
  
  ### Source the file in order to call functions from other r scripts
  source("/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/OPNMF_GLM_utils.R")
  
  ### Importation
  dataPrepres <- dataPrep(input_tsv, roi_start, roi_end)
  data_analyse <- dataPrepres[[1]] 
  nom_ROIs_volume <- dataPrepres[[2]] 
  
  variables <- c('disease_duration', 'cpz', 'age_disease_onset', 'panss_pos', 'panss_neg')
  
  ### loop over the clinical variables 
  for(vari in variables) {
    #######################
    ########### GLM for correlation 
    #######################
    #### only extract the data not empty for each variable
    data_used = data_analyse[!(is.na(data_analyse[,c(vari)]) | data_analyse[,c(vari)]=='NA'), ]
    
    model_pvalue <- general_linear_model_correlation(data = data_used, predictor =  c(vari), model_independent_var = c('age', 'sex', vari), ROIs = nom_ROIs_volume, model = 'fixe', correction = 'BH')
    write.table(model_pvalue, paste(output_dir, paste('GLM_', vari, '.csv', sep = ""), sep='/'), sep = ';', row.names = F)
    
    ### compute the pearson coefficient for each variable and plot the scatter plot
    for(ROI in nom_ROIs_volume) {
      print("ROI is:")
      print(ROI)
      ggscatter(data_used, x = ROI, y = vari, add = "reg.line", conf.int = TRUE, 
                add.params = list(color = "blue", fill = "lightgray"), cor.coef = TRUE, 
                cor.method = "pearson", xlab = ROI, ylab = vari)
      print("Stop here")
    }
  }
}

########################
########### Run the func
#########################
## CNSZ GM C17 
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CNSZ_GM_C17normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CNSZ_GM_C17"
roi_start <- 'component_1'
roi_end <- 'component_17' 
stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CNSZ WM C6 
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CNSZ_WM_C6normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CNSZ_WM_C6"
roi_start <- 'component_1'
roi_end <- 'component_6' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CNSZ CSF C6
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CNSZ_CSF_C6normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CNSZ_CSF_C6"
roi_start <- 'component_1'
roi_end <- 'component_6' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CN GM C28 
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CN_GM_C28normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CN_GM_C28"
roi_start <- 'component_1'
roi_end <- 'component_28' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CN WM C7 
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CN_WM_C7normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CN_WM_C7"
roi_start <- 'component_1'
roi_end <- 'component_7' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CN CSF C12
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CN_CSF_C12normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CN_CSF_C12"
roi_start <- 'component_1'
roi_end <- 'component_12' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)

## CNSZ GM C17 for mean signal not loading coefficient
input_tsv <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/R_PSCs_cognitive_scores_CNSZ_GM_C17normality_transformed.tsv"
output_dir <- "/Users/junhao.wen\ 1/test/SCZ_CLINICAL_OPNMF/clinical_correlation/CNSZ_GM_C17_mean_signal"
roi_start <- 'component_1'
roi_end <- 'component_17' 
#stats_pipeline(input_tsv, output_dir, roi_start, roi_end)