mgpCalculationPerDatasetFinal.Rmd

---
title: "MGP Calculation (rCTP derivation) --> unique cell lists per dataset -- FINAL"
output:
  html_document:
    df_print: paged
---

```{r echo=FALSE, cache=FALSE, comment=FALSE, warning=FALSE}
#load necessary libraries
library(magrittr)
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(limma)
library(edgeR)
library(rms)
library(useful)
library(stringr)
library(markerGeneProfile)
library(rlist)
library(kableExtra)
```

Now that we've loaded in the necessary libraries we can set up the marker gene lists to run the MGP algorithm. These marker gene lists will be used to describe each of the cell types in the relative cell-type proportions (rCTPs) we are calculating for each cohort. 

There are marker lists calculated from all available brain regions in the Allen Brain Institute Cell Atlas, at the class and subclass level. This means at the the resolution of inhibitory/excitatory (or GABAergic/glutamatergic) neuronal cells and at the resolution of neuronal subclasses (SST cells, IT cells). 

```{r}
subclass_CgG_MTG = read_csv(url('https://raw.githubusercontent.com/sonnyc247/MarkerSelection/master/Data/Outputs/CSVs_and_Tables/Markers/MTG_and_CgG_lfct2/new_MTGnCgG_lfct2.5_results.csv'))
```

Now that we've loaded in all the marker lists we want to load in all our cohorts count matrices. The cohorts were defined in pre-processing and we will continue to use those definitions: the Religious Orders Study and Rush Memory and Aging Project cohort (delineated ROSMAP from here on) which has dorsolateral prefrontal cortex samples, the Mayo Clinic cohort, hereafter referred to as MAYO, which consists of temporal cortex samples and the Mount Sinai Brain Bank cohort (referred to as MSBB), which has expression data sampled from Brodmann area 10, 22, 36 and 44 (written BM10, BM22, BM36, BM44 hereafter).

```{r}
#load in new cohorts QC-ed data, final count matrices and keep the dataframes with ENSEMBL ids
cohorts <- c("ROSMAP", "MAYO", "MSBBM10", "MSBBM22", "MSBBM36", "MSBBM44")
for(cohort in cohorts){
  if(str_detect(cohort, "MSBB")){
    cohort <- gsub('MSB', '', cohort) 
  }
  print(cohort)
  matrix_name <- paste0(cohort, "_matrix")
  
  if(cohort=="ROSMAP"){
    filename <- "adjusted_counts.rds"
  }
  else{
    filename <- paste0(matrix_name, ".rds")
  }
  
  setwd('./finalCountMatrices')
  matrix <- readRDS(filename)
  setwd('../')
  assign(matrix_name, matrix)

  count_df <- as.data.frame(matrix)
  if(cohort!= "ROSMAP"){
     count_df <- tibble:: rownames_to_column(count_df, var="Gene")
  }
  else{
     count_df <- count_df %>% dplyr::rename(Gene = gene_ID)
     colnames(count_df) <- gsub(x = colnames(count_df), pattern = "X", replacement = "") 
  }
 
  df_name <- paste0(cohort, "_count_df")
  assign(df_name, count_df)
}
```


Now that we have our marker lists and our cohort data loaded, we want to get the list of top  marker genes that are found within each cohort. Each marker gene list defines a cell type through a list of marker genes that are highly expressed and uniquely expressed in that cell type (to an extent). However, not all of the marker genes that define cell type proportions are found in each of the cohorts, for example: *geneX* is a top marker and is used to define *celltypeX*. It is found in cohorts BM10, BM22, and ROSMAP. Since *geneX* is not found in MAYO, BM36 or BM44 we will not include it as a top marker in MAYO, BM36 and BM44, but we will in BM10, BM22 and ROSMAP. 

```{r}

cohorts <- c("ROSMAP", "MAYO", "MSBBBM10", "MSBBBM22", "MSBBBM36", "MSBBBM44")
marker_lists <- c("subclass_CgG_MTG")
for(markers in marker_lists){
  for(cohort in cohorts){
    if(str_detect(cohort, "MSBB")){
      cohort <- gsub('MSBB', '', cohort) 
    }
    df_name <- paste0(cohort, "_count_df")
    count_df <- get(df_name)
    bulk_gene_names <- count_df$Gene
    ## find overlapping markers between marker list and bulk expression data
    new_cell_types = get(markers) %>% pull(subclass) %>% unique
 
    new_marker_list_ensembl <- lapply(new_cell_types, function(cell_type){
      return(subclass_CgG_MTG) %>% 
               filter(subclass == cell_type,  ensembl_gene_id %in% 
                        count_df$Gene) %>% pull(ensembl_gene_id)
    })
    names(new_marker_list_ensembl) <- new_cell_types
    new_marker_list_ensembl <-
      new_marker_list_ensembl[lapply(new_marker_list_ensembl,length)>0]
    print(new_cell_types)
    marker_name <- gsub('_df', '', markers) 
    marker_list_name <- paste0(marker_name, "_", cohort)
    assign(marker_list_name, new_marker_list_ensembl)
    setwd('./commonMarkerLists')
    saveRDS(get(marker_list_name), paste0(marker_list_name, ".rds"))
    setwd('../')
  }
}

```



```{r}
#we define our QC algorithm such that it returns a dataframe with the cell type, markers_used (list of marker genes used ' per cell type), removed_marker_ratios (list of removed marker ratios per cell type) and percent_variance_PC1 (list of variance explained by the first PC per cell type)
mgpQCMetrics <-function(count_df, mgp_markers){
  
  if(!all(c("Gene") %in% colnames(count_df))){
    stop("The count_df argument must be a df with a column named Gene (HGNC gene symbols)")
  }
  mgp_est<- markerGeneProfile::mgpEstimate(exprData=count_df,
                                                  genes=mgp_markers,
                                                  geneColName="Gene",
                                                  outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
                                                  geneTransform =NULL, #function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
                                                  groups=NULL, #if there are experimental groups provide them here. if not desired set to NULL
                                                  seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
                                                  removeMinority = TRUE)
  i= 0
  for(cell in names(mgp_markers)){
    i = i + 1
    cells_df <- mgp_est$usedMarkerExpression[i] %>% as.data.frame()
    masterlist <- paste0(rownames(cells_df), collapse=', ')
    num_markers <- length(rownames(cells_df))
    rm_marker_ratios <- mgp_est$removedMarkerRatios[i]
    if(!is.null(mgp_est$trimmedPCAs[[i]])){
      percent_variance <- ((summary(mgp_est$trimmedPCAs[[i]]))[6]) %>% as.data.frame()
      percent_variance_PC1 <- percent_variance[2,1]
    }
    else{
      percent_variance_PC1 <- NA
    }
    if(i==1){
      master_df <- data.frame( "celltype" = cell,
                               "markers_used" = masterlist,
                               "removed_marker_ratios" = rm_marker_ratios[[1]],
                               "percent_variance_PC1" = percent_variance_PC1,
                               "num_markers" = num_markers)  
    }
    else{
      df <- data.frame( "celltype" = cell,
                               "markers_used" = masterlist,
                               "removed_marker_ratios" = rm_marker_ratios[[1]],
                               "percent_variance_PC1" = percent_variance_PC1,
                               "num_markers" = num_markers) 
      master_df <- rbind(master_df, df)
    }
  }
  rownames(master_df) <- NULL
  return(master_df)
}


#calculate QC metrics for each marker gene list for each cohort
cohorts <- c("ROSMAP", "MAYO", "MSBBBM10", "MSBBBM22", "MSBBBM36", "MSBBBM44")
marker_lists <- c("subclass_CgG_MTG")
for(markers in marker_lists){
  for(cohort in cohorts){
    if(str_detect(cohort, "MSBB")){
      cohort <- gsub('MSBB', '', cohort) 
    }
    final_markers <- paste0(markers, "_", cohort)
    print(final_markers)
    df_name <- paste0(cohort, "_count_df")
    mgpResult <- mgpQCMetrics(get(df_name), mgp_markers = get(final_markers))
    mgpResult$cohort <- cohort
    mgpName <- paste0("mgpQCResults",cohort, final_markers)
    assign(mgpName, mgpResult)
    setwd(paste0('./mgpQCperDataset_',markers))
    saveRDS(get(mgpName), paste0(mgpName, ".rds"))
    setwd('../')
    print(mgpName)
    assign(mgpName, mgpResult)
    if(cohort =="ROSMAP"){
      all_cohorts_QC <- mgpResult
    }
    else{
      all_cohorts_QC <-rbind(all_cohorts_QC, mgpResult)
    }
  }
  all_cohorts_QC$cohort <- factor(all_cohorts_QC$cohort, levels = c("ROSMAP", "BM10", "BM44", "BM22", "BM36", "MAYO"))
  all_cohorts_QC <- arrange(all_cohorts_QC, cohort)
  setwd(paste0('./mgpQCperDataset_',markers))
  saveRDS(all_cohorts_QC, "all_cohorts_QC.rds")
  setwd('../')
}


```
We now have a wide variety of QC results for each of the marker gene lists and each of the cohorts. Let's plot some of the data to see what our results look like.

```{r echo=FALSE, cache=FALSE, comment=FALSE, warning=FALSE}
marker_lists <- c("subclass_CgG_MTG")
for(markers in marker_lists){
    
    setwd(paste0('./mgpQCperDataset_',markers))
    all_cohorts_QC <- readRDS("all_cohorts_QC.rds")
    QC_name <- paste0('all_cohorts_QC_',markers)
    assign(QC_name, all_cohorts_QC)
    setwd('../')
    
    QCplot <- get(QC_name) %>%  ggplot(aes(x = celltype, y = num_markers)) +
      theme_minimal() +
      geom_bar(stat = "identity", fill = "#e0abf5") +
      geom_hline(yintercept = 4) + 
      facet_wrap(~cohort, scales = 'free_x',nrow=1) + 
      ggtitle(paste0("Plot of Number of Markers Per Celltype \n Available in Each Dataset for ", markers, " Marker List")) +
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
            axis.title.x = element_text('Cell Type'),
            axis.title.y = element_text('Markers Used')) +
      coord_flip()
    
     QCplot2 <- get(QC_name) %>%  ggplot(aes(x = celltype, y = removed_marker_ratios)) +
      theme_minimal() +
      geom_bar(stat = "identity", fill = "#bab86c") +
      geom_hline(yintercept = 0.4) + 
      facet_wrap(~cohort, scales = 'free_x',nrow=1) + 
      ggtitle(paste0("Plot Marker Ratio Removed Per Celltype \n in Each Dataset for ", markers, " Marker List")) +
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
            axis.title.x = element_text('Cell Type'),
            axis.title.y = element_text('Markers Used')) +
      coord_flip()
    
    
    QCplot_name <- paste0('QCplot_',markers)
    assign(QCplot_name, QCplot)
    print(QCplot)
    setwd(paste0('./mgpQCperDataset_',markers))
    ggsave(paste0("QC_", markers, "_plot", ".png"))
    setwd('../')
    
    assign(paste0(QCplot_name, "_rmr"), QCplot2)
    print(QCplot2)
    setwd(paste0('./mgpQCperDataset_',markers))
    ggsave(paste0("QC_", markers, "_plot_rmr", ".png"))
    setwd('../')
}
```

Let's run the MGP algorithm with our marker lists.
```{r}
marker_lists <- c("subclass_CgG_MTG")
mgpCalc<-function(count_df, meta_df, markers){
  # calculate MGPs per sample
  estimations_human_markers<- mgpEstimate(exprData=count_df,
                                          genes=markers,
                                          geneColName="Gene",
                                          outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
                                          geneTransform = NULL,
                                          #function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
                                          groups=NULL, #if there are experimental groups provide them here. if not desired set to NULL
                                          seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
                                          removeMinority = TRUE)
  mgp_info <- list("mgp_df"= count_df, "estimations_human_markers" = estimations_human_markers)
  estimations_human_markers <- mgp_info$estimations_human_markers
  
  # matrix of mgp estimates per cell type, column name SST stores SST MGPs
  mgp_est <- estimations_human_markers$estimates %>% as.data.frame() 
  colnames(mgp_est) <- names(markers)
  mgp_est <- mgp_est %>% tibble::rownames_to_column(var = 'projid')
  
  # merge mgp data frame with sample metadata data frame
  mgp_est = merge(meta_df, mgp_est, by = 'projid')
  return(mgp_est)
}


#cohort can be ROSMAP, MAYO, MSBBBM10, MSBBBM22, MSBBBM26, MSBBBM44
mgpsForAMPAD<-function(cohort, markers){
  if(cohort == "ROSMAP"){
    covars <- c("sex","age_death")
    #get metadata
    setwd('./cohortMetadata')
    pheno_df <- readRDS("ROSmasterCOMPLETE.rds")
    setwd('../')
    }
  else{
    covars <- c("msex","AgeAtDeath")
    setwd('./QCpipelineResults')
    v_name <-  paste0("v_", cohort)
    voom <- readRDS(paste0(v_name, ".rds"))
    assign(v_name, voom)
    setwd('../')
    pheno_df <- voom$targets
    if(cohort == "MAYO"){
      pheno_df<- pheno_df %>% 
        dplyr::rename(
          projid = SampleID,
        )
    }
    else{
      pheno_df<- pheno_df %>% 
       dplyr::rename(
          projid = sampleIdentifier,
        )
    }
  }
  count_df_name <- paste0(cohort, "_count_df")
  mgp_result <- mgpCalc(get(count_df_name), pheno_df, markers)
  model.data <- mgp_result
  
  return(list("model" = model.data, "covars" = covars))
}

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
marker_lists <- c("subclass_CgG_MTG")


for(markers in marker_lists){
  for(cohort in cohorts){
    final_markers <- paste0(markers, "_", cohort)
    print(final_markers)
    mgp_result <- mgpsForAMPAD(cohort, get(final_markers))
    mgp_name <- paste0("mgp_",cohort)
    print(mgp_name)
    assign(mgp_name, mgp_result)
    setwd(paste0('./mgpResultsPerDataset_',markers))
    saveRDS(get(mgp_name), paste0(mgp_name, ".rds"))
    setwd('../')
  }
}


```

Now we've calculated MGPs for each of the cohorts and each of the common marker gene lists. We want to Z-score these MGPs so we can compare them in a mega-analysis within each marker gene list.

```{r}
cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
marker_lists <- c("subclass_CgG_MTG")

#Z score mgps per subject
for(markers in marker_lists){
  for(cohort in cohorts){
    final_markers <- paste0(markers, "_", cohort)
    print(final_markers)
    mgp_name <- paste0("mgp_",cohort)
    setwd(paste0('./mgpResultsPerDataset_',markers))
    mgp <- readRDS(paste0(mgp_name, ".rds"))
    mgp_df <- mgp$model
    markers_mgp <- get(final_markers)
    for(cell in names(markers_mgp)){
      mgp_df[,cell] <- as.numeric(scale(mgp_df[,cell], center = TRUE, scale = TRUE))
    }
    mgp_ZScored_name <- paste0(mgp_name, "_ZScored")
    assign(mgp_ZScored_name, mgp_df)
    saveRDS(get(mgp_ZScored_name), paste0(mgp_ZScored_name, ".rds"))
    setwd('../')
  }
}

```

Now that we've Z scored all the MGPs, we can move onto our mega-analysis. Before we do though let's check something:

Let's use our common marker lists to determine how each of the genes that define the cell-types are behaving relative to AD diagnosis, i.e. is there a strong unidirectional association between all of the genes in a cell-type and AD diagnosis?

#need to change to HGNC symbols

```{r}
frenchFryPlot<-function(AD_coef, AD_pval, marker_list, cohort){
  
  
  pathology_df <- AD_pval
  
  marker_df <-  data.frame(unlist(marker_list, use.names=F),rep(names(marker_list),
                                                                lengths(marker_list)))
  colnames(marker_df) <- c("marker", "celltype")
  for (marker in names(marker_list)){
    df <- pathology_df %>% filter(marker == marker)
    if (marker == names(marker_list)[1]){
      result <- df
    }
    else{
      result <- rbind(result, df)
    }
  }
  
  final_result <- merge(result, AD_coef, by="gene")
  final_result$signedP <- -log10(final_result$pval) *(sign(final_result$coef))
  
  merge_markers<- marker_df %>%dplyr::rename( gene = marker)
  final_result <- merge(merge_markers, final_result)
  final_result$cohort <- c(cohort)
  return(unique(final_result))
  
}

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")
marker_lists <- c("subclass_CgG_MTG")
setwd('./geneAnno')
gene_anno <- readRDS("gene_anno.rds")
setwd('../')

for (markers in marker_lists){
  for(cohort in cohorts){
    markers_for_plot <- get(paste0(markers, "_", cohort))
    print(cohort)
    lmod_name <-  paste0("lmod_" ,cohort)
    
    filename <- paste0(lmod_name, ".rds")
    setwd('./cohortQCMods')
    lmod <- readRDS(filename)
    setwd('../')
    assign(lmod_name, lmod)
    eb <- eBayes(lmod,robust = T)
    
    if(cohort == "ROSMAP"){
      n=11
    }
    if(cohort == "MAYO"){
      n=24
    }
    if(cohort == "BM10"){
      n=18
    }
    if(cohort == "BM22"){
      n=21
    }
    else if(cohort == "BM36" || cohort == "BM44"){
      n=19
    }

    print(n)
    gene_v_AD  <-lmod$coefficients[,c(n)]
    gene_v_AD <- as.data.frame(gene_v_AD)
    gene_v_AD <- tibble:: rownames_to_column(gene_v_AD, var="Gene")
    
    final_df <- merge(gene_v_AD, gene_anno)
    final_df$new_gene<- final_df$Hgnc_Gene %>% make.names(unique = T) #this is ONE way of dealing with the duplicates, just making them into separate, unique names
    #add new_gene column w/ Hgnc_Gene values filtered to have no duplicates
    AD_coef <- final_df[,c(2,4)]
    colnames(AD_coef) <- c("coef", "gene")
    
    coef_df <- AD_coef #dataframe with coefficient column and HGNC gene name column 

    #get p values
    p_geneAD <- eb$p.value
    p_geneAD  <-p_geneAD[,c(n)]
    p_geneAD <- as.data.frame(p_geneAD)
    p_geneAD <- tibble:: rownames_to_column(p_geneAD, var="Gene")
    final_df <- merge(p_geneAD, gene_anno)
    new_gene <- final_df$Hgnc_Gene %>% make.names(unique = T) #this is ONE way of dealing with the duplicates, just making them into separate, unique names
    #add new_gene column w/ Hgnc_Gene values filtered to have no duplicates
    final_df$new_gene <- new_gene
    AD_pval <- final_df[,c(2,4)]
    colnames(AD_pval) <- c("pval", "gene")
    
    gene_pathology_association <- frenchFryPlot(AD_coef, AD_pval, cohort, marker_list =markers_for_plot)
    assign(paste0("fry_plot", cohort), gene_pathology_association)
  }
  fry_df <- rbind(fry_plotROSMAP, fry_plotMAYO, fry_plotBM10,
                  fry_plotBM22, fry_plotBM36, fry_plotBM44)
  fry_df$fdr <- p.adjust(fry_df$pval, method="fdr")
  fry_df$signedFDR <- -log10(fry_df$fdr) *(sign(fry_df$signedP))
  fry_df$sig <- ifelse(fry_df$pval < 0.05, "#94C973", "#808080")
  
  for (cell in names((markers_for_plot))){
    celltype_fry_df <- fry_df %>% filter(celltype == cell)
    celltype_heat_map <- ggplot(celltype_fry_df, aes(cohort, gene, fill= signedP))+
      theme_minimal() + geom_tile() + 
      scale_fill_gradient2(low="darkblue", high="darkgreen", guide="colorbar") +   
      ggtitle(paste0("Significance of \n" , cell, " ", markers , 
                     "\n marker genes per cohort")) 
    
    celltype_french <- ggplot(celltype_fry_df, aes(x= gene,y= signedP))+
      theme_minimal() + theme(axis.text.x = element_text(angle = 45)) +
      geom_bar(stat="identity", fill = celltype_fry_df$sig)+ 
      facet_wrap(~cohort, scales = 'free_x',nrow=6) + 
      ggtitle(paste0("Significance of \n" , cell, " ", markers , 
                     "\n marker genes per cohort"))
    
    celltype_heat_map
    celltype_french
    
    setwd('./frenchHeat')
    heat_name <- paste0("heat_map",make.names(cell))
    assign(heat_name, celltype_heat_map)
    print(get(heat_name))
    saveRDS(get(heat_name), paste0(heat_name, ".rds"))
    ggsave(paste0(heat_name, markers, "_plot", ".png"), width = 15, height = 12)
    
    french_name <- paste0("french_fry",make.names(cell))
    assign(french_name, celltype_french)
    print(get(french_name))
    saveRDS(get(french_name), paste0(french_name, ".rds"))
    ggsave(paste0(french_name, markers, "_plot", ".png"), width = 15, height = 20)
    setwd('../')
    
  }
}
```
There's one last thing we want to do before moving onto the mega-analysis. We want to use our marker list subclass_FINAL_common_final, to create a marker list that is identical save for one change: the SST cell type will be defined without the SST gene included, and we will rerun the MGP algorithm. This way we will proceed with the SST cell-type defined without SST to ensure that in the mega-analysis it is not just this one gene driving the association between the SST cell-type and AD. 

```{r}


setwd('./commonMarkerLists')
subclass_FINAL_common_final_noSST <- readRDS('subclass_FINAL_common_final_noSST.rds')
setwd('../')

mgpCalc<-function(count_df, meta_df, markers){
  # calculate MGPs per sample
  estimations_human_markers<- mgpEstimate(exprData=count_df,
                                          genes=markers,
                                          geneColName="Gene",
                                          outlierSampleRemove=F, # should outlier samples removed. This is done using boxplot stats.
                                          geneTransform = NULL,
                                          #function(x){homologene::mouse2human(x)$humanGene}, # this is the default option for geneTransform
                                          groups=NULL, #if there are experimental groups provide them here. if not desired set to NULL
                                          seekConsensus = FALSE, # ensures gene rotations are positive in both of the groups
                                          removeMinority = FALSE)
  mgp_info <- list("mgp_df"= count_df, "estimations_human_markers" = estimations_human_markers)
  estimations_human_markers <- mgp_info$estimations_human_markers
  
  # matrix of mgp estimates per cell type, column name SST stores SST MGPs
  mgp_est <- estimations_human_markers$estimates %>% as.data.frame() 
  colnames(mgp_est) <- names(markers)
  mgp_est <- mgp_est %>% tibble::rownames_to_column(var = 'projid')
  
  # merge mgp data frame with sample metadata data frame
  mgp_est = merge(meta_df, mgp_est, by = 'projid')
  return(mgp_est)
}


#cohort can be ROSMAP, MAYO, MSBBBM10, MSBBBM22, MSBBBM26, MSBBBM44
mgpsForAMPAD<-function(cohort, markers){
  if(cohort == "ROSMAP"){
    covars <- c("sex","age_death")
    }
  else{
    covars <- c("msex","AgeAtDeath")
  }
  setwd('./QCpipelineResults')
  v_name <-  paste0("v_", cohort)
  voom <- readRDS(paste0(v_name, ".rds"))
  assign(v_name, voom)
  setwd('../')
  pheno_df <- voom$targets
  if(cohort == "ROSMAP" || cohort == "MAYO"){
    pheno_df<- pheno_df %>% 
     dplyr::rename(
        projid = SampleID,
      )
  }
  else{
    pheno_df<- pheno_df %>% 
     dplyr::rename(
        projid = sampleIdentifier,
      )
  }
  count_df_name <- paste0(cohort, "_count_df")
  mgp_result <- mgpCalc(get(count_df_name), pheno_df, markers)
  model.data <- mgp_result
  
  return(list("model" = model.data, "covars" = covars))
}

cohorts <- c("ROSMAP", "MAYO", "BM10", "BM22", "BM36", "BM44")

for(cohort in cohorts){
  mgp_result <- mgpsForAMPAD(cohort, subclass_FINAL_common_final_noSST)
  mgp_name <- paste0("mgp_",cohort)
  print(mgp_name)
  assign(mgp_name, mgp_result)
  setwd(paste0('./mgpResults_subclass_FINAL_noSST'))
  saveRDS(get(mgp_name), paste0(mgp_name, ".rds"))
  mgp_df <- mgp_result$model
  
  for(cell in names(subclass_FINAL_common_final_noSST)){
    mgp_df[,cell] <- as.numeric(scale(mgp_df[,cell], center = TRUE, scale = TRUE))
  }
  mgp_ZScored_name <- paste0(mgp_name, "_ZScored")
  assign(mgp_ZScored_name, mgp_df)
  saveRDS(get(mgp_ZScored_name), paste0(mgp_ZScored_name, ".rds"))
  
  setwd('../')
}

```