WGCNA

Last update: 24-10-16

Last markdown compiled: 24-10-16

knitr::opts_chunk$set(warning = FALSE, message = FALSE, results = FALSE) 

# dirs
## Sebastian @Linux
dir_moa <- "/home/isar/Dropbox/scn_moa/"
# dir_moa <- "/Users/home_sh/Library/CloudStorage/Dropbox/scn_moa"  # Sebastian @Mac

# data, revised 2023.10.18
# exp
load(file = paste0(dir_moa, "/data/final/", "gnp.pro.RData"))
load(file = paste0(dir_moa, "/data/final/", "gnp.mrn.RData"))
load(file = paste0(dir_moa, "/data/final/", "gnp.mic.RData"))

load(file = paste0(dir_moa, "/data/final/", "gnp.Isets.RData"))

# results
load(file = paste0(dir_moa, "/data/final/", "dex.all.RData"))

## packages
source(paste0(dir_moa, "/scripts/", "scn_moa_fct.R"))


# load packages
sapply(packages_general, require, character.only = TRUE)
sapply(packages_WGCNA, require, character.only = TRUE)
#sapply(packages_cluster, require, character.only = TRUE)
#sapply(packages_heatmap, require, character.only = TRUE)


library(gt) # for tables

WGCNA

fct

wgcna

# calc 0.2 : 0.5
get_blockwiseModules <- function(exp, soft_power, mergeCutHeight_val){
  blockwiseModules(exp,
      maxBlockSize = ncol(exp),
      networkType = "signed", # network type
      TOMType = 'signed', # redundant
      power = soft_power,
      mergeCutHeight = mergeCutHeight_val,
      numericLabels = T,
      randomSeed = 1234,
      verbose = 3)
}

# plotter Dendrogram
plot_WGCNA_dendro <- function(bwnet){
  plotDendroAndColors(bwnet$dendrograms[[1]], cbind(bwnet$unmergedColors, bwnet$colors),
                    c("unmerged", "merged"),
                    dendroLabels = F,
                    addGuide = T,
                    hang = 0.03,
                    guideHang = 0.05)
}


# extract eigengene to module correlation
extract_Eigengene_cor2Module <- function(module_eigengenes, traits, nSamples){
  module.popGen.corr <- cor(module_eigengenes, traits, use = 'p')
  module.popGen.corr.pval <- corPvalueStudent(module.popGen.corr, nSamples)
  list <- list(
    module.popGen.corr = module.popGen.corr,
    module.popGen.corr.pval = module.popGen.corr.pval)
  return(list)
}


# nr of genes/module as bar graph
plotter_genesPERmodule <- function(bwnet, title){
  # get number
  df_genesPerModule <- as.data.frame(table(bwnet$colors)) 
  # plot
  df_genesPerModule %>%
    ggplot(aes(Var1, Freq)) +
    geom_col()+
    geom_text(aes(label=Freq, vjust = -sign(Freq))) +
    ggtitle(title)
}


## prep data for heatmap of eigengene to module correlation
prep_HM_EigG2Mod <- function(module_eigengenes, traits){
  # prep
  heatmap.data <- merge(module_eigengenes, traits, by = 'row.names')
  heatmap.data %<>% column_to_rownames('Row.names')
  return(heatmap.data)
}


# module_eigengenes

# extract eigengenes
extract_EigenG_per_Module <- function(module_eigengenes, san, levels){
  # eigengenes / module
  moduleDF <- module_eigengenes %>%
    tibble::rownames_to_column("id") %>%
    # Here we are performing an inner join with a subset of metadata
    dplyr::inner_join(san %>%
      dplyr::select(id, stage),
    by = c("id")
    )
  # add levels
  moduleDF$stage <- factor(moduleDF$stage, levels = levels)
  # return
  return(moduleDF)
}


# Eigengenes in modules
plotter_eigengenesINmodules <- function(moduleDF, module, title){
  moduleDF %>%
    select(id, val = module, stage) %>%
    ggplot(aes( x = stage, y = val, color = stage)) +
    geom_boxplot(width = 0.2) +
    ggforce::geom_sina(maxwidth = 0.3) +
    #theme_classic() +
    
    ggtitle(title) +
    scale_color_manual(values = colors_stage)+
    theme_gray() +
    theme(legend.position="none") 
}


# old name: module_gene_mapping_mrn
extract_featuresOFmodules <- function(bwnet){
  fOm <-  as.data.frame(bwnet$colors)
  names(fOm)[1] <- "module"
  fOm %<>% rownames_to_column('ENSG')
  fOm$module <- paste0("ME", fOm$module)
  return(fOm)
}


# get most significant genes correlating with stages
extract_features_of_traits <- function(exp_4wgcna, traits, stage, nSamp){
  trait <- traits %>% select(stage)
  gene.signf.corr <- cor(exp_4wgcna, trait, use = "p")
  gene.signf.corr.pval <- corPvalueStudent(gene.signf.corr, nSamp)
  gene.signf.corr.pval <- as.data.frame(gene.signf.corr.pval) 
  names(gene.signf.corr.pval)[1] <- "pval"
  gene.signf.corr.pval %<>% rownames_to_column("SYMBOL") 
  gene.signf.corr.pval %<>% arrange(pval)
  return(gene.signf.corr.pval)
}


plotter_eigengenesINmodules.white <- function(moduleDF, module){
  moduleDF %>%
    select(id, val = module, stage) %>%
    ggplot(aes( x = stage, y = val, color = stage)) +
    geom_boxplot(width = 1) +
    ggforce::geom_sina(maxwidth = 0.3, size = 4) +
    scale_color_manual(values = colors_stage)+
    theme_minimal(base_size=22) +
    theme(legend.position="none") 
}


plotter_eigengenesINmodules <- function(moduleDF, module, title){
  moduleDF %>%
    select(id, val = module, stage) %>%
    ggplot(aes( x = stage, y = val, color = stage)) +
    geom_boxplot(width = 0.2) +
    ggforce::geom_sina(maxwidth = 0.3) +
    #theme_classic() +
    
    ggtitle(title) +
    scale_color_manual(values = colors_stage)+
    theme_gray() +
    theme(legend.position="none") 
}

heatmap

Creates heatmap directly from features in modules. Also use own heatmap for single features (later)

make_module_heatmap <- function(module_name,
                                expression_mat,
                                metadata_df,
                                gene_module_key_df,
                                module_eigengenes_df) {
  # Create a summary heatmap of a given module.
  #
  # Args:
  # module_name: a character indicating what module should be plotted, e.g. "ME19"
  # expression_mat: The full gene expression matrix. Default is `normalized_counts`.
  # metadata_df: a data frame with refinebio_accession_code and time_point
  #              as columns. Default is `metadata`.
  # gene_module_key: a data.frame indicating what genes are a part of what modules. Default is `gene_module_key`.
  # module_eigengenes: a sample x eigengene data.frame with samples as row names. Default is `module_eigengenes`.
  #
  # Returns:
  # A heatmap of expression matrix for a module's genes, with a barplot of the
  # eigengene expression for that module.

  # Set up the module eigengene with its refinebio_accession_code
  module_eigengene <- module_eigengenes_df %>%
    dplyr::select(all_of(module_name), id)
    #tibble::rownames_to_column("id")

  # Set up column annotation from metadata
  col_annot_df <- metadata_df %>%
    # Only select the treatment and sample ID columns
    dplyr::select(id, stage) %>%
    # Add on the eigengene expression by joining with sample IDs
    dplyr::inner_join(module_eigengene, by = "id") %>%
    # Arrange by patient and time point
    dplyr::arrange(stage) %>%
    # Store sample
    tibble::column_to_rownames("id")
  col_annot_df$stage <- factor(col_annot_df$stage, levels = c('MB', 'PM', 'MC', 'MM', 'B', 'S', 'PMN'))

  # Create the ComplexHeatmap column annotation object
  col_annot <- ComplexHeatmap::HeatmapAnnotation(
    # Supply treatment labels
    stage = col_annot_df$stage,
    # Add annotation barplot
    module_eigengene = ComplexHeatmap::anno_barplot(dplyr::select(col_annot_df, module_name)),
    # Pick colors for each experimental group in time_point
    col = list(stage = colors_stage)
  )

  # Get a vector of the Ensembl gene IDs that correspond to this module
  module_genes <- gene_module_key_df %>%
    #rename(module = bwnet_mrn$colors) %>%
    dplyr::filter(module == module_name) %>%
    dplyr::pull(ENSG)

  # Set up the gene expression data frame
  mod_mat <- expression_mat %>%
    t() %>%
    as.data.frame() %>%
    # Only keep genes from this module
    dplyr::filter(rownames(.) %in% module_genes) %>%
    # Order the samples to match col_annot_df
    dplyr::select(rownames(col_annot_df)) %>%
    # Data needs to be a matrix
    as.matrix()

  # Normalize the gene expression values
  mod_mat <- mod_mat %>%
    # Scale can work on matrices, but it does it by column so we will need to
    # transpose first
    t() %>%
    scale() %>%
    # And now we need to transpose back
    t()

  # Create a color function based on standardized scale
  color_func <- circlize::colorRamp2(
    c(-2, 0, 2),
    c("#67a9cf", "#f7f7f7", "#ef8a62")
  )

  # Plot on a heatmap
  heatmap <- ComplexHeatmap::Heatmap(mod_mat,
    name = module_name,
    # Supply color function
    col = color_func,
    # Supply column annotation
    bottom_annotation = col_annot,
    # We don't want to cluster samples
    cluster_columns = F,
    # We don't need to show sample or gene labels
    show_row_names = FALSE,
    show_column_names = FALSE,
    show_row_dend = FALSE,
  )

  # Return heatmap
  return(heatmap)
}

options

# allow parallel processing (Sebastian at LINUX)
enableWGCNAThreads(nThreads = 10)

# vector to test soft threshold powers
power_single <- c(c(1:10), seq(from = 12, to = 50, by =2)) #


set.seed(1409)

graphics

border_up <- border(color = color_up, linetype = "solid", size = 1)
border_down <- border(color = color_down, linetype = "solid", size = 1)

# , out.height = "200%"

MRN

prep

# prep data
exp_mrn_4wgcna <- gnp.mrn$exp_mrn_tpm_vsn_NIPALS
dim(exp_mrn_4wgcna)

# QC
qc_mrn <- goodSamplesGenes(exp_mrn_4wgcna)
qc_mrn$allOK # all oK!

analyse soft power

• choice of the soft thresholding power β to which co-expression similarity is raised to calculate adjacency
• raising the correlation to a power will reduce the noise of the correlations in the adjacency matrix
• You want to a power that gives you a big enough R2 but is not excessively large.

WGCNA’s authors recommend using a power that has an signed R2 above 0.80, otherwise they warn your results may be too noisy to be meaningful. Online, I also read the R2 should be 0.9. No idea whats correct.

If you find you have all very low R2 values this may be because there are too many genes with low expression values that are cluttering up the calculations. (not my problem)

What is R2? The mean connectivity???

sft_mrn <- pickSoftThreshold(exp_mrn_4wgcna, powerVector = power_single, networkType = "signed", verbose = 5)

# visualize
sft.data_mrn <- sft_mrn$fitIndices

a1_mrn <- sft.data_mrn %>% ggplot(aes(Power, SFT.R.sq, label = Power)) +
  geom_point() +
  geom_text(nudge_y = 0.3, size = 2) +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Power", y = "Scale free topology model fit, signed R^2") +
  theme_classic() +
  theme(text=element_text(size=6))

a2_mrn <- sft.data_mrn %>% ggplot((aes(Power, mean.k., label = Power))) +
  geom_point() +
  geom_text(nudge_y = 500, size = 2) +
  labs(x = "Power", y = "Mean connectivity") +
  theme_classic() +
  theme(text=element_text(size=6))

fig_fst_mrn_pre <- ggarrange(
    a1_mrn, 
    a2_mrn, 
    nrow = 2)

fig_fst_mrn <- annotate_figure(fig_fst_mrn_pre,
               top = text_grob("Softpower mRNA", color = "black", face = "bold", size = 12),)

fig_fst_mrn # determine: sft_mrn = 9

constants

soft_power_mrn <- 9

## numbers for pval calc
nSamples_mrn <- nrow(exp_mrn_4wgcna)
nGenes_mrn <- ncol(exp_mrn_4wgcna)

levels_mrn = c('MB', 'PM', 'MC', 'MM', 'B', 'S', 'PMN')

network

What is the difference between networkType and TOMtype?

Explanations:

networkType controls the adjacency matrix with: signed=[(1 + cor)/2]^pow unsigned=|cor|^pow signed hybrid = ifelse(cor < 0, 0, (1 + cor)/2)^pow

The TOMType determines whether the signs of the correlation are used to ‘reinforce’ or ‘cancel out’ edge projects of the topological overlap (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/TechnicalReports/signedTOM.pdf). Note that for signed and signed-hybrid networks, the adjacency for negative correlations are 0 or close to 0, so the impact of using a signed topological overlap is negligible.

The standard, or “unsigned” TOM assumes that neighbor-mediated connections can always be considered as “reinforcing” the direct connection. This may not always be the case, and the signed TOM is an attempt to take this into account.

The take-home message from these notes is this: signed TOM takes into account possible anti-reinforcing connection strengths that may occur in unsigned networks. Since the anti-reinforcing connection strengths (practically) cannot occur in signed networks, in signed networks the signed and unsigned TOM are (practically) identical.

# make matrix numeric
exp_mrn_4wgcna[] <- sapply(exp_mrn_4wgcna, as.numeric)

# get blockModules from 0.2 - 0.5
# res.bwnet_mrn <- lapply(seq(0.2, 0.5, by=0.1), get_blockwiseModules, exp = exp_mrn_4wgcna, soft_power = soft_power_mrn)
# names(res.bwnet_mrn) <- c('point2', 'point3', 'point4', 'point5')

# save(sft.data_mrn,
#     res.bwnet_mrn,
#     file = paste0(dir_moa, "/data/final/preps/wgcna/", "mrn_bwnet.RData"))

load(file = paste0(dir_moa, "/data/final/preps/wgcna/", "mrn_bwnet.RData"))

dendro

plot_WGCNA_dendro(res.bwnet_mrn$point2)

plot_WGCNA_dendro(res.bwnet_mrn$point3)

plot_WGCNA_dendro(res.bwnet_mrn$point4)

plot_WGCNA_dendro(res.bwnet_mrn$point5)

traits to associate

Apparently, there are more trait options: genes of defined pathways, groups, also coded as binary (via tool) ?

But how do I do this? Bc these are feature annos, not sample annos????

## define phenotypes to associate with modules: stages
stage_mrn <- gnp.mrn$san_mrn %>% select(id, stage)
stage_mrn$stage <- factor(stage_mrn$stage, levels = c('MB', 'PM', 'MC', 'MM', 'B', 'S', 'PMN'))
stage_mrn %<>% column_to_rownames('id')

trait_MB_mrn <- stage_mrn %>% mutate(MB = ifelse(stage == "MB", 1, 0)) %>% select(MB)
trait_rest_mrn <- binarizeCategoricalColumns(stage_mrn, 
                                         includePairwise = F,
                                         includeLevelVsAll = T,
                                         minCount = 3)

rownames(trait_rest_mrn) <- rownames(stage_mrn)
colnames(trait_rest_mrn) <- gsub("stage.", "", colnames(trait_rest_mrn))
colnames(trait_rest_mrn) <- gsub(".vs.all", "", colnames(trait_rest_mrn))

traits_mrn <- cbind(trait_MB_mrn, trait_rest_mrn)

Cor: modules & traits

## module eigengenes
MEs_mrn_0.2 <- res.bwnet_mrn$point2$MEs
MEs_mrn_0.3 <- res.bwnet_mrn$point3$MEs
MEs_mrn_0.4 <- res.bwnet_mrn$point4$MEs
MEs_mrn_0.5 <- res.bwnet_mrn$point5$MEs

## correlation of eigengenes to modules
EG2M_mrn_0.2 <- extract_Eigengene_cor2Module(MEs_mrn_0.2, traits_mrn, nSamples_mrn)
EG2M_mrn_0.3 <- extract_Eigengene_cor2Module(MEs_mrn_0.3, traits_mrn, nSamples_mrn)
EG2M_mrn_0.4 <- extract_Eigengene_cor2Module(MEs_mrn_0.4, traits_mrn, nSamples_mrn)
EG2M_mrn_0.5 <- extract_Eigengene_cor2Module(MEs_mrn_0.5, traits_mrn, nSamples_mrn)


## prep for heatmap
datHM_mrn_0.2 <- prep_HM_EigG2Mod(MEs_mrn_0.2, traits_mrn)
datHM_mrn_0.3 <- prep_HM_EigG2Mod(MEs_mrn_0.3, traits_mrn)
datHM_mrn_0.4 <- prep_HM_EigG2Mod(MEs_mrn_0.4, traits_mrn)
datHM_mrn_0.5 <- prep_HM_EigG2Mod(MEs_mrn_0.5, traits_mrn)

# plot heatmap (need to pick cols!)
colnames(datHM_mrn_0.2)
CorLevelPlot(datHM_mrn_0.2,
           x = names(datHM_mrn_0.2)[14:20], # picks traits
           y = names(datHM_mrn_0.2)[1:13], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mRNA (cut 0.2)")

colnames(datHM_mrn_0.3)
CorLevelPlot(datHM_mrn_0.3,
           x = names(datHM_mrn_0.3)[11:16], # picks traits
           y = names(datHM_mrn_0.3)[1:10], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mRNA (cut 0.3)")

CorLevelPlot(datHM_mrn_0.4,
           x = names(datHM_mrn_0.4)[8:14], # picks traits
           y = names(datHM_mrn_0.4)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mRNA (cut 0.4)")

CorLevelPlot(datHM_mrn_0.5,
           x = names(datHM_mrn_0.5)[7:13], # picks traits
           y = names(datHM_mrn_0.5)[1:6], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mRNA (cut 0.5)")

Genes per Module

# requires san
san_mrn = gnp.mrn$san_mrn

EGpM_mrn_0.2 <- extract_EigenG_per_Module(MEs_mrn_0.2, san_mrn, levels_mrn)
EGpM_mrn_0.3 <- extract_EigenG_per_Module(MEs_mrn_0.3, san_mrn, levels_mrn)
EGpM_mrn_0.4 <- extract_EigenG_per_Module(MEs_mrn_0.4, san_mrn, levels_mrn)
EGpM_mrn_0.5 <- extract_EigenG_per_Module(MEs_mrn_0.5, san_mrn, levels_mrn)

# total number of genes per module
plotter_genesPERmodule(res.bwnet_mrn$point2, "mRNA 0.2") + ylim(0, 3500)

plotter_genesPERmodule(res.bwnet_mrn$point3, "mRNA 0.3") + ylim(0, 3500)

plotter_genesPERmodule(res.bwnet_mrn$point4, "mRNA 0.4") + ylim(0, 3600)

plotter_genesPERmodule(res.bwnet_mrn$point5, "mRNA 0.5") + ylim(0, 3800)

Modules

0.2

plot_ME1_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME1', "ME1: 3161 features")
plot_ME2_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME2', "ME2: 1714 features")
plot_ME3_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME3', "ME3: 1145 features")
plot_ME4_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME4', "ME4: 407 features")
plot_ME5_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME5', "ME5: 306 features")
plot_ME6_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME6', "ME6: 283 features")
plot_ME7_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME7', "ME7: 278 features")
plot_ME8_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME8', "ME8: 277 features")
plot_ME9_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME9', "ME9: 164 features")
plot_ME10_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME10', "ME10: 145 features")
plot_ME11_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME11', "ME11: 114 features")
plot_ME12_mrn_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_0.2, 'ME12', "ME12: 43 features")


fig_EigGinModules_mrn_0.2.pt1<- ggarrange(
  plot_ME1_mrn_0.2 + border_down,
  plot_ME2_mrn_0.2 + border_up,
  plot_ME3_mrn_0.2 + border_up,
  plot_ME4_mrn_0.2,
  plot_ME5_mrn_0.2,
  plot_ME6_mrn_0.2 + border_up,
  ncol = 3,
  nrow = 2)
fig_EigGinModules_mrn_0.2.pt1 <- annotate_figure(fig_EigGinModules_mrn_0.2.pt1,
            top = text_grob("Modules mRNA (cut 0.2), pt1", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.2.pt1

fig_EigGinModules_mrn_0.2.pt2 <- ggarrange(
  plot_ME7_mrn_0.2,
  plot_ME8_mrn_0.2,
  plot_ME9_mrn_0.2,
  plot_ME10_mrn_0.2,
  plot_ME11_mrn_0.2,
  plot_ME12_mrn_0.2,
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.2.pt2 <- annotate_figure(fig_EigGinModules_mrn_0.2.pt2,
            top = text_grob("Modules mRNA (cut 0.2), pt2", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.2.pt2

0.3

plot_ME1_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME1', "ME1: 3161 features")
plot_ME2_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME2', "ME2: 2121 features")
plot_ME3_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME3', "ME3: 1145 features")
plot_ME4_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME4', "ME4: 706 features")
plot_ME5_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME5', "ME5: 306 features")
plot_ME6_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME6', "ME6: 277 features")
plot_ME7_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME7', "ME7: 164 features")
plot_ME8_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME8', "ME8: 114 features")
plot_ME9_mrn_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_0.3, 'ME9', "ME9: 43 features")

# first test: to look at them
fig_EigGinModules_mrn_0.3.pt1 <- ggarrange(
  plot_ME1_mrn_0.3,
  plot_ME2_mrn_0.3,
  plot_ME3_mrn_0.3,
  plot_ME4_mrn_0.3,
  plot_ME5_mrn_0.3,
  plot_ME6_mrn_0.3,
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.3.pt1 <- annotate_figure(fig_EigGinModules_mrn_0.3.pt1,
            top = text_grob("Modules mRNA (cut 0.3), pt1", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.3.pt1

fig_EigGinModules_mrn_0.3.pt2 <- ggarrange(
  plot_ME7_mrn_0.3,
  plot_ME8_mrn_0.3,
  plot_ME9_mrn_0.3,
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.3.pt2 <- annotate_figure(fig_EigGinModules_mrn_0.3.pt2,
            top = text_grob("Modules mRNA (cut 0.3), pt2", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.3.pt2

0.4 = final

plot_ME1_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME1', "ME1: 3467 features")
plot_ME2_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME2', "ME2: 3266 features")
plot_ME3_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME3', "ME3: 706 features")
plot_ME4_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME4', "ME4: 277 features")
plot_ME5_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME5', "ME5: 164 features")
plot_ME6_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME6', "ME6: 114 features")
plot_ME7_mrn_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_0.4, 'ME7', "ME7: 43 features")

fig_EigGinModules_mrn_0.4.pt1 <- ggarrange(
  plot_ME1_mrn_0.4,
  plot_ME2_mrn_0.4,
  plot_ME3_mrn_0.4,
  plot_ME4_mrn_0.4,
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.4.pt1 <- annotate_figure(fig_EigGinModules_mrn_0.4.pt1,
            top = text_grob("Modules mRNA (cut 0.4), pt1", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.4.pt1

fig_EigGinModules_mrn_0.4.pt2 <- ggarrange(
  plot_ME5_mrn_0.4,
  plot_ME6_mrn_0.4, 
  plot_ME7_mrn_0.4,
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.4.pt2 <- annotate_figure(fig_EigGinModules_mrn_0.4.pt2,
            top = text_grob("Modules mRNA (cut 0.4), pt2", 
                color = "black", face = "bold", size = 18))
fig_EigGinModules_mrn_0.4.pt2

### final
# singles

plot_ME2_mrn_0.4 + ggtitle("mRNA ME2: 3266 ft. cont. up")

plot_ME1_mrn_0.4 + ggtitle("mRNA ME1: 3467 ft. cont. down")

plot_ME7_mrn_0.4 + ggtitle("mRNA ME7: 43 ft. down in MB")

plot_ME3_mrn_0.4 + ggtitle("mRNA ME3: 706 ft. up in MC,MM,B")

plot_ME4_mrn_0.4 + ggtitle("mRNA ME4: 277 ft. down in MC,MM,B")

plot_ME6_mrn_0.4 + ggtitle("mRNA ME6: 144 ft. up in PMN")

fig_Modules_mrn_0.4.fin <- ggarrange(
  plot_ME2_mrn_0.4,
  plot_ME1_mrn_0.4,
  plot_ME7_mrn_0.4,
  plot_ME3_mrn_0.4,
  plot_ME4_mrn_0.4,
  plot_ME6_mrn_0.4,
  ncol = 3,
  nrow = 2)

fig_Modules_mrn_0.4.fin <- annotate_figure(fig_Modules_mrn_0.4.fin,
            top = text_grob("WGCNA modules mRNA", 
                color = "black", face = "bold", size = 18))
fig_Modules_mrn_0.4.fin

0.5

plot_ME1_mrn_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_0.5, 'ME1', "ME1: 3467 features")
plot_ME2_mrn_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_0.5, 'ME2', "ME2: 3266 features")
plot_ME3_mrn_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_0.5, 'ME3', "ME3: 749 features")
plot_ME4_mrn_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_0.5, 'ME4', "ME4: 391 features")
plot_ME5_mrn_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_0.5, 'ME5', "ME5: 164 features")


  
fig_EigGinModules_mrn_0.5.pre <- ggarrange(
  plot_ME1_mrn_0.5,
  plot_ME2_mrn_0.5,
  plot_ME3_mrn_0.5,
  plot_ME4_mrn_0.5,
  plot_ME5_mrn_0.5,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_0.5 <- annotate_figure(
  fig_EigGinModules_mrn_0.5.pre,
  top = text_grob("Modules mRNA (cut 0.5)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_0.5

Features of modules

fom_mrn_0.2 <- extract_featuresOFmodules(res.bwnet_mrn$point2)
fom_mrn_0.3 <- extract_featuresOFmodules(res.bwnet_mrn$point3)
fom_mrn_0.4 <- extract_featuresOFmodules(res.bwnet_mrn$point4)
fom_mrn_0.5 <- extract_featuresOFmodules(res.bwnet_mrn$point5)

fom_mrn <- list(
  fom_mrn_0.2 = fom_mrn_0.2,
  fom_mrn_0.3 = fom_mrn_0.3,
  fom_mrn_0.4 = fom_mrn_0.4,
  fom_mrn_0.5 = fom_mrn_0.5
)

## now you can export features of modules to eg clipboard
# featuresOFmodules_mrn_0.2 %>% filter(module == "ME2") %>% pull(ENSG) %>% write_clip()

Features of traits

Can give the most important feature of each trait.

Cave: this is NOT module specific but only trait specific. Tables have to be filtered by pValue.

Not 100% sure how to use this info.

fot_MB_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'MB', nSamples_mrn)
fot_PM_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'PM', nSamples_mrn)
fot_MC_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'MC', nSamples_mrn)
fot_MM_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'MM', nSamples_mrn)
fot_B_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'B', nSamples_mrn)
fot_S_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'S', nSamples_mrn)
fot_PMN_mrn <- extract_features_of_traits(exp_mrn_4wgcna, traits_mrn, 'PMN', nSamples_mrn)

fot_mrn <- list(
  fot_MB_mrn = fot_MB_mrn,
  fot_PM_mrn = fot_PM_mrn,
  fot_MC_mrn = fot_MC_mrn,
  fot_MM_mrn = fot_MM_mrn,
  fot_B_mrn = fot_B_mrn,
  fot_S_mrn = fot_S_mrn,
  fot_PMN_mrn = fot_PMN_mrn
)

heatmaps

0.2

# plot_hm_0.2_mrn_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_mrn_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_mrn_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_mrn_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_mrn_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_mrn_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)

0.4

# plot_hm_0.4_mrn_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_mrn_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_mrn_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_mrn_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_mrn_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_mrn_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)

0.5

# plot_hm_0.5_mrn_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.5, EGpM_mrn_0.5)
# plot_hm_0.5_mrn_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.5, EGpM_mrn_0.5)
# plot_hm_0.5_mrn_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.5, EGpM_mrn_0.5)
# plot_hm_0.5_mrn_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.5, EGpM_mrn_0.5)
# plot_hm_0.5_mrn_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.5, EGpM_mrn_0.5)

export

res.mrn <- list(
  fom_mrn = fom_mrn,
  fot_mrn = fot_mrn
)

PRO

prep

# prep data
exp_pro_4wgcna <- gnp.pro$exp_pro_vsn # using the TNP data
dim(exp_pro_4wgcna) # 27, 3156

# QC
qc_pro <- goodSamplesGenes(exp_pro_4wgcna)
qc_pro$allOK # all oK!

analyse soft power

sft_pro <- pickSoftThreshold(exp_pro_4wgcna, powerVector = power_single, networkType = "signed", verbose = 5)

# visualize
sft.data_pro <- sft_pro$fitIndices

a1_pro <- sft.data_pro %>% ggplot(aes(Power, SFT.R.sq, label = Power)) +
  geom_point() +
  geom_text(nudge_y = 0.3, size = 2) +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Power", y = "Scale free topology model fit, signed R^2") +
  theme_classic() +
  theme(text=element_text(size=6))

a2_pro <- sft.data_pro %>% ggplot((aes(Power, mean.k., label = Power))) +
  geom_point() +
  geom_text(nudge_y = 500, size = 2) +
  labs(x = "Power", y = "Mean connectivity") +
  theme_classic() +
  theme(text=element_text(size=6))

fig_fst_pro_pre <- ggarrange(
    a1_pro, 
    a2_pro, 
    nrow = 2)

fig_fst_pro <- annotate_figure(fig_fst_pro_pre,
               top = text_grob("Softpower Proteome", color = "black", face = "bold", size = 12),)

fig_fst_pro # determine: sft_pro = 7

constants (data specific)

soft_power_pro <- 7

## numbers for pval calc
nSamples_pro <- nrow(exp_pro_4wgcna)
nGenes_pro <- ncol(exp_pro_4wgcna)

levels_pro = c('MB', 'PM', 'MC', 'MM', 'B', 'PMN')

network

# make matrix numeric
exp_pro_4wgcna[] <- sapply(exp_pro_4wgcna, as.numeric)

# get blockModules from 0.2 - 0.5
# res.bwnet_pro <- lapply(seq(0.2, 0.5, by=0.1), get_blockwiseModules, exp = exp_pro_4wgcna, soft_power = soft_power_pro)
# names(res.bwnet_pro) <- c('point2', 'point3', 'point4', 'point5')
# 
# save(sft.data_pro,
#      res.bwnet_pro,
#      file = paste0(dir_moa, "/data/final/preps/wgcna/", "pro_bwnet.RData"))

load(file = paste0(dir_moa, "/data/final/preps/wgcna/", "pro_bwnet.RData"))

dendro

plot_WGCNA_dendro(res.bwnet_mrn$point2)

plot_WGCNA_dendro(res.bwnet_mrn$point3)

plot_WGCNA_dendro(res.bwnet_mrn$point4)

plot_WGCNA_dendro(res.bwnet_mrn$point5)

traits to associate

Apparently, there are more trait options: genes of defined pathways, groups, also coded as binary (via tool) ?

But how do I do this? Bc these are feature annos, not sample annos????

## define phenotypes to associate with modules: stages
stage_pro <- gnp.pro$san_pro %>% select(id, stage)
stage_pro$stage <- factor(stage_pro$stage, levels = c('MB', 'PM', 'MC', 'MM', 'B', 'PMN'))
stage_pro %<>% column_to_rownames('id')

trait_MB_pro <- stage_pro %>% mutate(MB = ifelse(stage == "MB", 1, 0)) %>% select(MB)
trait_rest_pro <- binarizeCategoricalColumns(stage_pro, 
                                         includePairwise = F,
                                         includeLevelVsAll = T,
                                         minCount = 2)

rownames(trait_rest_pro) <- rownames(stage_pro)
colnames(trait_rest_pro) <- gsub("stage.", "", colnames(trait_rest_pro))
colnames(trait_rest_pro) <- gsub(".vs.all", "", colnames(trait_rest_pro))

traits_pro <- cbind(trait_MB_pro, trait_rest_pro)

Cor: modules & traits

## module eigengenes
MEs_pro_0.2 <- res.bwnet_pro$point2$MEs
MEs_pro_0.3 <- res.bwnet_pro$point3$MEs
MEs_pro_0.4 <- res.bwnet_pro$point4$MEs
MEs_pro_0.5 <- res.bwnet_pro$point5$MEs

## correlation of eigengenes to modules
EG2M_pro_0.2 <- extract_Eigengene_cor2Module(MEs_pro_0.2, traits_pro, nSamples_pro)
EG2M_pro_0.3 <- extract_Eigengene_cor2Module(MEs_pro_0.3, traits_pro, nSamples_pro)
EG2M_pro_0.4 <- extract_Eigengene_cor2Module(MEs_pro_0.4, traits_pro, nSamples_pro)
EG2M_pro_0.5 <- extract_Eigengene_cor2Module(MEs_pro_0.5, traits_pro, nSamples_pro)


## prep for heatmap
datHM_pro_0.2 <- prep_HM_EigG2Mod(MEs_pro_0.2, traits_pro)
datHM_pro_0.3 <- prep_HM_EigG2Mod(MEs_pro_0.3, traits_pro)
datHM_pro_0.4 <- prep_HM_EigG2Mod(MEs_pro_0.4, traits_pro)
datHM_pro_0.5 <- prep_HM_EigG2Mod(MEs_pro_0.5, traits_pro)

# plot heatmap (need to pick cols!)
CorLevelPlot(datHM_pro_0.2,
           x = names(datHM_pro_0.2)[15:20], # picks traits
           y = names(datHM_pro_0.2)[1:14], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits proA (cut 0.2)")

names(datHM_pro_0.3)
CorLevelPlot(datHM_pro_0.3,
           x = names(datHM_pro_0.3)[12:17], # picks traits
           y = names(datHM_pro_0.3)[1:11], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits proA (cut 0.3)")

names(datHM_pro_0.4)
CorLevelPlot(datHM_pro_0.4,
           x = names(datHM_pro_0.4)[8:13], # picks traits
           y = names(datHM_pro_0.4)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits proA (cut 0.4)")

CorLevelPlot(datHM_pro_0.5,
           x = names(datHM_pro_0.5)[7:12], # picks traits
           y = names(datHM_pro_0.5)[1:6], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits proA (cut 0.5)")

Genes per Module

# requires san
san_pro = gnp.pro$san_pro

EGpM_pro_0.2 <- extract_EigenG_per_Module(MEs_pro_0.2, san_pro, levels_pro)
EGpM_pro_0.3 <- extract_EigenG_per_Module(MEs_pro_0.3, san_pro, levels_pro)
EGpM_pro_0.4 <- extract_EigenG_per_Module(MEs_pro_0.4, san_pro, levels_pro)
EGpM_pro_0.5 <- extract_EigenG_per_Module(MEs_pro_0.5, san_pro, levels_pro)

# total number of genes per module
plotter_genesPERmodule(res.bwnet_pro$point2, "proA 0.2")+ ylim(0, 750)

plotter_genesPERmodule(res.bwnet_pro$point3, "proA 0.3")+ ylim(0, 750)

plotter_genesPERmodule(res.bwnet_pro$point4, "proA 0.4") + ylim(0, 1250)

plotter_genesPERmodule(res.bwnet_pro$point5, "proA 0.5") + ylim(0, 1250)

Modules

0.2

plot_ME1_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME1', "ME1: 558 features")
plot_ME2_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME2', "ME2: 461 features")
plot_ME3_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME3', "ME3: 351 features")
plot_ME4_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME4', "ME4: 316 features")
plot_ME5_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME5', "ME5: 245 features")
plot_ME6_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME6', "ME6: 184 features")
plot_ME7_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME7', "ME7: 184 features")
plot_ME8_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME8', "ME8: 132 features")
plot_ME9_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME9', "ME9: 116 features")
plot_ME10_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME10', "ME10: 90 features")
plot_ME11_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME11', "ME11: 70 features")
plot_ME12_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME12', "ME12: 59 features")
plot_ME13_pro_0.2 <- plotter_eigengenesINmodules(EGpM_pro_0.2, 'ME13', "ME13: 49 features")


fig_EigGinModules_pro_0.2.pt1 <- ggarrange(
  plot_ME1_pro_0.2,
  plot_ME2_pro_0.2,
  plot_ME3_pro_0.2,
  plot_ME4_pro_0.2,
  plot_ME5_pro_0.2,
  plot_ME6_pro_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_pro_0.2.pt1 <- annotate_figure(
  fig_EigGinModules_pro_0.2.pt1,
  top = text_grob("Modules pro (cut 0.2), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.2.pt1

fig_EigGinModules_pro_0.2.pt2 <- ggarrange(
 plot_ME7_pro_0.2,
  plot_ME8_pro_0.2,
  plot_ME9_pro_0.2,
  plot_ME10_pro_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_0.2.pt2 <- annotate_figure(
  fig_EigGinModules_pro_0.2.pt2,
  top = text_grob("Modules pro (cut 0.2), pt2", color = "black", face = "bold", size = 18))

fig_EigGinModules_pro_0.2.pt2

fig_EigGinModules_pro_0.2.pt3 <- ggarrange(
  plot_ME11_pro_0.2,
  plot_ME12_pro_0.2,
  plot_ME13_pro_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_0.2.pt3 <- annotate_figure(
  fig_EigGinModules_pro_0.2.pt3,
  top = text_grob("Modules pro (cut 0.2), pt3", color = "black", face = "bold", size = 18))

fig_EigGinModules_pro_0.2.pt3

0.3

plot_ME1_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME1', "ME1: 558 features")
plot_ME2_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME2', "ME2: 467 features")
plot_ME3_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME3', "ME3: 461 features")
plot_ME4_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME4', "ME4: 429 features")
plot_ME5_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME5', "ME5: 316 features")
plot_ME6_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME6', "ME6: 233 features")
plot_ME7_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME7', "ME7: 132 features")
plot_ME8_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME8', "ME8: 90 features")
plot_ME9_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME9', "ME9: 70 features")
plot_ME10_pro_0.3 <- plotter_eigengenesINmodules(EGpM_pro_0.3, 'ME10', "ME10: 59 features")

# pt1
fig_EigGinModules_pro_0.3.pt1 <- ggarrange(
  plot_ME1_pro_0.3,
  plot_ME2_pro_0.3,
  plot_ME3_pro_0.3,
  plot_ME4_pro_0.3,
  plot_ME5_pro_0.3,
  plot_ME6_pro_0.3,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_pro_0.3.pt1 <- annotate_figure(
  fig_EigGinModules_pro_0.3.pt1,
  top = text_grob("Modules pro (cut 0.3), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.3.pt1

# first test: to look at them
fig_EigGinModules_pro_0.3.pt2 <- ggarrange(
  plot_ME7_pro_0.3,
  plot_ME8_pro_0.3,
  plot_ME9_pro_0.3,
  plot_ME10_pro_0.3,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_pro_0.3.pt2 <- annotate_figure(
  fig_EigGinModules_pro_0.3.pt2,
  top = text_grob("Modules pro (cut 0.3), pt2", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.3.pt2

0.4

plot_ME1_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME1', "ME1: 1025 features")
plot_ME2_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME2', "ME2: 461 features")
plot_ME3_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME3', "ME3: 448 features")
plot_ME4_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME4', "ME4: 429 features")
plot_ME5_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME5', "ME5: 303 features")
plot_ME6_pro_0.4 <- plotter_eigengenesINmodules(EGpM_pro_0.4, 'ME6', "ME6: 149 features")


fig_EigGinModules_pro_0.4 <- ggarrange(
  plot_ME1_pro_0.4,
  plot_ME2_pro_0.4,
  plot_ME3_pro_0.4,
  plot_ME4_pro_0.4,
  plot_ME5_pro_0.4,
  plot_ME6_pro_0.4,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_0.4 <- annotate_figure(
  fig_EigGinModules_pro_0.4,
  top = text_grob("Modules pro (cut 0.4)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.4

# ### fin
# fig_EigGinModules_pro_0.4_fin <- ggarrange(
#   plot_ME1_pro_0.4,
#   plot_ME5_pro_0.4,
#   plot_ME2_pro_0.4,
#   plot_ME3_pro_0.4,
#   plot_ME4_pro_0.4,
# 
#   NULL,
#   common.legend = T,
#   legend = "none",
#   ncol = 3,
#   nrow = 2)
# 
# fig_EigGinModules_pro_0.4_fin <- annotate_figure(
#   fig_EigGinModules_pro_0.4_fin,
#   top = text_grob("WGCNA modules proteome", color = "black", face = "bold", size = 18))
#                 
# fig_EigGinModules_pro_0.4_fin

0.5 = final

plot_ME1_pro_0.5 <- plotter_eigengenesINmodules(EGpM_pro_0.5, 'ME1', "ME1: 1025 features")
plot_ME2_pro_0.5 <- plotter_eigengenesINmodules(EGpM_pro_0.5, 'ME2', "ME2: 732 features")
plot_ME3_pro_0.5 <- plotter_eigengenesINmodules(EGpM_pro_0.5, 'ME3', "ME3: 461 features")
plot_ME4_pro_0.5 <- plotter_eigengenesINmodules(EGpM_pro_0.5, 'ME4', "ME4: 448 features")
plot_ME5_pro_0.5 <- plotter_eigengenesINmodules(EGpM_pro_0.5, 'ME5', "ME5: 149 features")

fig_EigGinModules_pro_0.5 <- ggarrange(
  plot_ME1_pro_0.5,
  plot_ME2_pro_0.5,
  plot_ME3_pro_0.5,
  plot_ME4_pro_0.5,
  plot_ME5_pro_0.5,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_0.5 <- annotate_figure(
  fig_EigGinModules_pro_0.5,
  top = text_grob("Modules pro (cut 0.5)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.5

#### fin
## singles
plot_ME1_pro_0.5 + ggtitle("Prot ME1: 1025 ft. cont. up")

plot_ME2_pro_0.5  + ggtitle("Prot ME2: 732 ft. cont. down")

plot_ME4_pro_0.5 + ggtitle("Prot ME4: 448 up in MB")

plot_ME3_pro_0.5 + ggtitle("Prot ME3: 461 up in PM & MC")

fig_EigGinModules_pro_0.5_fin <- ggarrange(
  plot_ME1_pro_0.5,
  plot_ME2_pro_0.5,
  plot_ME4_pro_0.5,
  plot_ME3_pro_0.5,
  common.legend = T,
  legend = "none",
  ncol = 2,
  nrow = 2)

fig_EigGinModules_pro_0.5_fin <- annotate_figure(
  fig_EigGinModules_pro_0.5_fin,
  top = text_grob("Modules proteome", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_0.5_fin

Features of modules

fom_pro_0.2 <- extract_featuresOFmodules(res.bwnet_pro$point2)
fom_pro_0.3 <- extract_featuresOFmodules(res.bwnet_pro$point3)
fom_pro_0.4 <- extract_featuresOFmodules(res.bwnet_pro$point4)
fom_pro_0.5 <- extract_featuresOFmodules(res.bwnet_pro$point5)

fom_pro <- list(
  fom_pro_0.2 = fom_pro_0.2,
  fom_pro_0.3 = fom_pro_0.3,
  fom_pro_0.4 = fom_pro_0.4,
  fom_pro_0.5 = fom_pro_0.5
)

## now you can export features of modules to eg clipboard
# featuresOFmodules_mrn_0.2 %>% filter(module == "ME2") %>% pull(ENSG) %>% write_clip()

Features of traits

fot_MB_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'MB', nSamples_pro)
fot_PM_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'PM', nSamples_pro)
fot_MC_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'MC', nSamples_pro)
fot_MM_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'MM', nSamples_pro)
fot_B_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'B', nSamples_pro)
fot_PMN_pro <- extract_features_of_traits(exp_pro_4wgcna, traits_pro, 'PMN', nSamples_pro)

fot_pro <- list(
  fot_MB_pro = fot_MB_pro,
  fot_PM_pro = fot_PM_pro,
  fot_MC_pro = fot_MC_pro,
  fot_MM_pro = fot_MM_pro,
  fot_B_pro = fot_B_pro,
  fot_PMN_pro = fot_PMN_pro
)

heatmaps

0.2

# plot_hm_0.2_pro_ME1 <- make_module_heatmap("ME1", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)
# plot_hm_0.2_pro_ME2 <- make_module_heatmap("ME2", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)
# plot_hm_0.2_pro_ME3 <- make_module_heatmap("ME3", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)
# plot_hm_0.2_pro_ME4 <- make_module_heatmap("ME4", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)
# plot_hm_0.2_pro_ME5 <- make_module_heatmap("ME5", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)
# plot_hm_0.2_pro_ME6 <- make_module_heatmap("ME6", exp_pro_4wgcna, san_pro, fom_pro_0.2, EGpM_pro_0.2)

0.4

# plot_hm_0.4_pro_ME1 <- make_module_heatmap("ME1", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)
# plot_hm_0.4_pro_ME2 <- make_module_heatmap("ME2", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)
# plot_hm_0.4_pro_ME3 <- make_module_heatmap("ME3", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)
# plot_hm_0.4_pro_ME4 <- make_module_heatmap("ME4", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)
# plot_hm_0.4_pro_ME5 <- make_module_heatmap("ME5", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)
# plot_hm_0.4_pro_ME6 <- make_module_heatmap("ME6", exp_pro_4wgcna, san_pro, fom_pro_0.4, EGpM_pro_0.4)

export

res.pro <- list(
  fom_pro = fom_pro,
  fot_pro = fot_pro
)

MIC

prep

# prep data
exp_mic_4wgcna <- gnp.mic$exp_mic_vsn_NIPALS # nipals imputed
dim(exp_mic_4wgcna) # 24, 283

# QC
qc_mic <- goodSamplesGenes(exp_mic_4wgcna)
qc_mic$allOK # all oK!

analyse soft power

sft_mic <- pickSoftThreshold(exp_mic_4wgcna, powerVector = power_single, networkType = "signed", verbose = 5)

# visualize
sft.data_mic <- sft_mic$fitIndices

a1_mic <- sft.data_mic %>% ggplot(aes(Power, SFT.R.sq, label = Power)) +
  geom_point() +
  geom_text(nudge_y = 0.3, size = 2) +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Power", y = "Scale free topology model fit, signed R^2") +
  theme_classic() +
  theme(text=element_text(size=6))

a2_mic <- sft.data_mic %>% ggplot((aes(Power, mean.k., label = Power))) +
  geom_point() +
  geom_text(nudge_y = 500, size = 2) +
  labs(x = "Power", y = "Mean connectivity") +
  theme_classic() +
  theme(text=element_text(size=6))

fig_fst_mic_pre <- ggarrange(
    a1_mic, 
    a2_mic, 
    nrow = 2)

fig_fst_mic <- annotate_figure(fig_fst_mic_pre,
               top = text_grob("Softpower mic", color = "black", face = "bold", size = 12),)

fig_fst_mic # determine: sft_mic = 12

constants

soft_power_mic <- 11

## numbers for pval calc
nSamples_mic <- nrow(exp_mic_4wgcna)
nSamples_mic
nGenes_mic <- ncol(exp_mic_4wgcna)
nGenes_mic

levels_mic = c('MB', 'MC', 'MM', 'B', 'S', 'PMN')

network

# make matrix numeric
exp_mic_4wgcna[] <- sapply(exp_mic_4wgcna, as.numeric)

# # get blockModules from 0.2 - 0.5
# res.bwnet_mic <- lapply(seq(0.2, 0.5, by=0.1), get_blockwiseModules, exp = exp_mic_4wgcna, soft_power = soft_power_mic)
# names(res.bwnet_mic) <- c('point2', 'point3', 'point4', 'point5')
# 
# save(sft.data_mic,
#      res.bwnet_mic,
#      file = paste0(dir_moa, "/data/final/preps/wgcna/", "mic_bwnet.RData"))

load(file = paste0(dir_moa, "/data/final/preps/wgcna/", "mic_bwnet.RData"))

dendro

plot_WGCNA_dendro(res.bwnet_mic$point2) # same! 
plot_WGCNA_dendro(res.bwnet_mic$point3) # same! 

plot_WGCNA_dendro(res.bwnet_mic$point4) # same! 
plot_WGCNA_dendro(res.bwnet_mic$point5) # different

Seems that 0.2, 0.3 and 0.4 all deliver the same results!

traits to associate

Apparently, there are more trait options: genes of defined pathways, groups, also coded as binary (via tool) ?

But how do I do this? Bc these are feature annos, not sample annos????

## define phenotypes to associate with modules: stages
stage_mic <- gnp.mic$san_mic %>% select(id, stage)
stage_mic$stage <- factor(stage_mic$stage, levels = c('MB', 'MC', 'MM', 'B', 'S', 'PMN'))
stage_mic %<>% column_to_rownames('id')

trait_MB_mic <- stage_mic %>% mutate(MB = ifelse(stage == "MB", 1, 0)) %>% select(MB)
trait_rest_mic <- binarizeCategoricalColumns(stage_mic, 
                                         includePairwise = F,
                                         includeLevelVsAll = T,
                                         minCount = 2)

rownames(trait_rest_mic) <- rownames(stage_mic)
colnames(trait_rest_mic) <- gsub("stage.", "", colnames(trait_rest_mic))
colnames(trait_rest_mic) <- gsub(".vs.all", "", colnames(trait_rest_mic))

traits_mic <- cbind(trait_MB_mic, trait_rest_mic)

Cor: modules & traits

## module eigengenes
MEs_mic_0.2 <- res.bwnet_mic$point2$MEs
MEs_mic_0.3 <- res.bwnet_mic$point3$MEs
MEs_mic_0.4 <- res.bwnet_mic$point4$MEs
MEs_mic_0.5 <- res.bwnet_mic$point5$MEs

## correlation of eigengenes to modules
EG2M_mic_0.2 <- extract_Eigengene_cor2Module(MEs_mic_0.2, traits_mic, nSamples_mic)
EG2M_mic_0.3 <- extract_Eigengene_cor2Module(MEs_mic_0.3, traits_mic, nSamples_mic)
EG2M_mic_0.4 <- extract_Eigengene_cor2Module(MEs_mic_0.4, traits_mic, nSamples_mic)
EG2M_mic_0.5 <- extract_Eigengene_cor2Module(MEs_mic_0.5, traits_mic, nSamples_mic)


## prep for heatmap
datHM_mic_0.2 <- prep_HM_EigG2Mod(MEs_mic_0.2, traits_mic)
datHM_mic_0.3 <- prep_HM_EigG2Mod(MEs_mic_0.3, traits_mic)
datHM_mic_0.4 <- prep_HM_EigG2Mod(MEs_mic_0.4, traits_mic)
datHM_mic_0.5 <- prep_HM_EigG2Mod(MEs_mic_0.5, traits_mic)

# plot heatmap (need to pick cols!)
names(datHM_mic_0.2) 
CorLevelPlot(datHM_mic_0.2,
           x = names(datHM_mic_0.2)[8:13], # picks traits
           y = names(datHM_mic_0.2)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mic (cut 0.2)")

names(datHM_mic_0.3)
CorLevelPlot(datHM_mic_0.3,
           x = names(datHM_mic_0.3)[8:13], # picks traits
           y = names(datHM_mic_0.3)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mic (cut 0.3)")

names(datHM_mic_0.4)
CorLevelPlot(datHM_mic_0.4,
           x = names(datHM_mic_0.4)[8:13], # picks traits
           y = names(datHM_mic_0.4)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mic (cut 0.4)")

names(datHM_mic_0.5)
CorLevelPlot(datHM_mic_0.5,
           x = names(datHM_mic_0.5)[7:12], # picks traits
           y = names(datHM_mic_0.5)[1:6], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mic (cut 0.5)")

0.2-0.4: all the same!

Genes per Module

# requires san
san_mic = gnp.mic$san_mic

EGpM_mic_0.2 <- extract_EigenG_per_Module(MEs_mic_0.2, san_mic, levels_mic)
EGpM_mic_0.3 <- extract_EigenG_per_Module(MEs_mic_0.3, san_mic, levels_mic)
EGpM_mic_0.4 <- extract_EigenG_per_Module(MEs_mic_0.4, san_mic, levels_mic)
EGpM_mic_0.5 <- extract_EigenG_per_Module(MEs_mic_0.5, san_mic, levels_mic)

# total number of genes per module
plotter_genesPERmodule(res.bwnet_mic$point2, "micA 0.2")+ ylim(0, 130)

plotter_genesPERmodule(res.bwnet_mic$point3, "micA 0.3")+ ylim(0, 130)

plotter_genesPERmodule(res.bwnet_mic$point4, "micA 0.4")+ ylim(0, 130)

plotter_genesPERmodule(res.bwnet_mic$point5, "micA 0.5")+ ylim(0, 130)

Modules

0.2 = final

plot_ME1_mic_0.2 <- plotter_eigengenesINmodules(EGpM_mic_0.2, 'ME1', "ME1: 101 features")
plot_ME2_mic_0.2 <- plotter_eigengenesINmodules(EGpM_mic_0.2, 'ME2', "ME2: 84 features")
plot_ME3_mic_0.2 <- plotter_eigengenesINmodules(EGpM_mic_0.2, 'ME3', "ME3: 28 features")
plot_ME4_mic_0.2 <- plotter_eigengenesINmodules(EGpM_mic_0.2, 'ME4', "ME4: 25 features")
plot_ME5_mic_0.2 <- plotter_eigengenesINmodules(EGpM_mic_0.2, 'ME5', "ME5: 20 features")

fig_EigGinModules_mic_0.2 <- ggarrange(
  plot_ME1_mic_0.2,
  plot_ME2_mic_0.2,
  plot_ME3_mic_0.2,
  plot_ME4_mic_0.2,
  plot_ME5_mic_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mic_0.2 <- annotate_figure(
  fig_EigGinModules_mic_0.2,
  top = text_grob("Modules mic (cut 0.2)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mic_0.2

### fin
fig_EigGinModules_mic_0.2_fin <- ggarrange(
    plot_ME1_mic_0.2,
  plot_ME3_mic_0.2,
  plot_ME2_mic_0.2,
  plot_ME5_mic_0.2,

  common.legend = T,
  legend = "none",
  ncol = 2,
  nrow = 2)

fig_EigGinModules_mic_0.2_fin <- annotate_figure(
  fig_EigGinModules_mic_0.2_fin,
  top = text_grob("WGCNA modules microRNA", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mic_0.2_fin

figure

# ME1
p.ME1.mic_0.2 <- plotter_eigengenesINmodules.white(EGpM_mic_0.2, 'ME1')
p.ME1.mic_0.2_fin <- p.ME1.mic_0.2 + theme_noX + ylab("Intensity")
p.ME1.mic_0.2_fin

ggsave(
  filename = paste0(dir_moa, "/figures/fig2/", "p.ME1.mic_0.2_fin", ".png"),
  plot = p.ME1.mic_0.2_fin,
  device = png,
  width = 5.5,
  height = 5)

# ME2
p.ME2.mic_0.2 <- plotter_eigengenesINmodules.white(EGpM_mic_0.2, 'ME2')
p.ME2.mic_0.2_fin <- p.ME2.mic_0.2 + theme_nada
p.ME2.mic_0.2_fin

ggsave(
  filename = paste0(dir_moa, "/figures/fig2/", "p.ME2.mic_0.2_fin", ".png"),
  plot = p.ME2.mic_0.2_fin,
  device = png,
  width = 5.5,
  height = 5)


# ME3
p.ME3.mic_0.2 <- plotter_eigengenesINmodules.white(EGpM_mic_0.2, 'ME3')
p.ME3.mic_0.2_fin <- p.ME3.mic_0.2 + theme_nada
p.ME3.mic_0.2_fin

ggsave(
  filename = paste0(dir_moa, "/figures/fig2/", "p.ME3.mic_0.2_fin", ".png"),
  plot = p.ME3.mic_0.2_fin,
  device = png,
  width = 5.5,
  height = 5)


# ME4
p.ME4.mic_0.2 <- plotter_eigengenesINmodules.white(EGpM_mic_0.2, 'ME4')
p.ME4.mic_0.2_fin <- p.ME4.mic_0.2 + ylab("Intensity") + xlab("Stage")
p.ME4.mic_0.2_fin

ggsave(
  filename = paste0(dir_moa, "/figures/fig2/", "p.ME4.mic_0.2_fin", ".png"),
  plot = p.ME4.mic_0.2_fin,
  device = png,
  width = 5.5,
  height = 5)

# ME5
p.ME5.mic_0.2 <- plotter_eigengenesINmodules.white(EGpM_mic_0.2, 'ME5')
p.ME5.mic_0.2_fin <- p.ME5.mic_0.2 + theme_noY + xlab("Stage")
#p.ME5.mic_0.2_fin

ggsave(
  filename = paste0(dir_moa, "/figures/fig2/", "p.ME5.mic_0.2_fin", ".png"),
  plot = p.ME5.mic_0.2_fin,
  device = png,
  width = 5.5,
  height = 5)

MIC numbers

# total
ncol(gnp.mic$exp_mic_lg2) # total = 283

# not assigned
# 283 - (101 + 84 + 28 + 25 + 20)

# create
mic.numbers <- data.frame(
  Module = c('ME1: up', 'ME2: dn', 'ME3: MC up','ME4: B/S up', 'ME5: up->dn', 'ME0: NA'),
  number = c(101, 84, 28, 25, 20, 25)
)

# ratio
mic.numbers$"of total" <- round(mic.numbers$number / 283 * 100, 0)
mic.numbers$"of total" <- paste0("(", mic.numbers$"of total", "%)" )

# t.mic.wgcna <- mic.numbers %>% gt()

# gtsave(t.mic.wgcna, 
#        filename = paste0(dir_moa, "/figures/fig2/", "t.mic.wgcna", ".png")) ,
#        path = paste0(dir_moa, "/figures/fig2/"

0.5

plot_ME1_mic_0.5 <- plotter_eigengenesINmodules(EGpM_mic_0.5, 'ME1', "ME1: 101 features")
plot_ME2_mic_0.5 <- plotter_eigengenesINmodules(EGpM_mic_0.5, 'ME2', "ME2: 84 features")
plot_ME3_mic_0.5 <- plotter_eigengenesINmodules(EGpM_mic_0.5, 'ME3', "ME3: 28 features")
plot_ME4_mic_0.5 <- plotter_eigengenesINmodules(EGpM_mic_0.5, 'ME4', "ME4: 25 features")
plot_ME5_mic_0.5 <- plotter_eigengenesINmodules(EGpM_mic_0.5, 'ME5', "ME5: 20 features")

fig_EigGinModules_mic_0.5 <- ggarrange(
  plot_ME1_mic_0.5,
  plot_ME2_mic_0.5,
  plot_ME3_mic_0.5,
  plot_ME4_mic_0.5,
  plot_ME5_mic_0.5,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_mic_0.5 <- annotate_figure(
  fig_EigGinModules_mic_0.5,
  top = text_grob("Modules mic (cut 0.5)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mic_0.5

# fig_EigGinModules_mic_0.5 <- annotate_figure(fig_EigGinModules_mic_0.5.pre,
#             top = text_grob("Modules micA (cut 0.5)", 
#                 color = "black", face = "bold", size = 18))
# fig_EigGinModules_mic_0.5

Features of modules

fom_mic_0.2 <- extract_featuresOFmodules(res.bwnet_mic$point2)
fom_mic_0.3 <- extract_featuresOFmodules(res.bwnet_mic$point3)
fom_mic_0.4 <- extract_featuresOFmodules(res.bwnet_mic$point4)
fom_mic_0.5 <- extract_featuresOFmodules(res.bwnet_mic$point5)

identical(fom_mic_0.2$module, fom_mic_0.4$module) # same! 
identical(fom_mic_0.3$module, fom_mic_0.4$module) # same! 

fom_mic <- list(
  fom_mic_0.2 = fom_mic_0.2,
  fom_mic_0.5 = fom_mic_0.5
)


## now you can export features of modules to eg clipboard
# featuresOFmodules_mrn_0.2 %>% filter(module == "ME2") %>% pull(ENSG) %>% write_clip()

Features of traits

Can give the most important feature of each trait.

Cave: this is NOT module specific but only trait specific. Tables have to be filtered by pValue.

Not 100% sure how to use this info.

fot_MB_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'MB', nSamples_mic)
fot_MC_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'MC', nSamples_mic)
fot_MM_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'MM', nSamples_mic)
fot_B_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'B', nSamples_mic)
fot_S_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'S', nSamples_mic)
fot_PMN_mic <- extract_features_of_traits(exp_mic_4wgcna, traits_mic, 'PMN', nSamples_mic)

fot_mic <- list(
  fot_MB_mic = fot_MB_mic,
  fot_MC_mic = fot_MC_mic,
  fot_MM_mic = fot_MM_mic,
  fot_B_mic = fot_B_mic,
  fot_S_mic = fot_S_mic,
  fot_PMN_mic = fot_PMN_mic
)

heatmaps

0.2

plot_hm_0.2_mic_ME1 <- make_module_heatmap("ME1", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)
plot_hm_0.2_mic_ME2 <- make_module_heatmap("ME2", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)
plot_hm_0.2_mic_ME3 <- make_module_heatmap("ME3", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)
plot_hm_0.2_mic_ME4 <- make_module_heatmap("ME4", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)
plot_hm_0.2_mic_ME5 <- make_module_heatmap("ME5", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)
plot_hm_0.2_mic_ME6 <- make_module_heatmap("ME6", exp_mic_4wgcna, san_mic, fom_mic_0.2, EGpM_mic_0.2)

0.5

# plot_hm_0.5_mic_ME1 <- make_module_heatmap("ME1", exp_mic_4wgcna, san_mic, fom_mic_0.5, EGpM_mic_0.5)
# plot_hm_0.5_mic_ME2 <- make_module_heatmap("ME2", exp_mic_4wgcna, san_mic, fom_mic_0.5, EGpM_mic_0.5)
# plot_hm_0.5_mic_ME3 <- make_module_heatmap("ME3", exp_mic_4wgcna, san_mic, fom_mic_0.5, EGpM_mic_0.5)
# plot_hm_0.5_mic_ME4 <- make_module_heatmap("ME4", exp_mic_4wgcna, san_mic, fom_mic_0.5, EGpM_mic_0.5)
# plot_hm_0.5_mic_ME5 <- make_module_heatmap("ME5", exp_mic_4wgcna, san_mic, fom_mic_0.5, EGpM_mic_0.5)

export

res.mic <- list(
  fom_mic = fom_mic,
  fot_mic = fot_mic
)


save(res.mic,
     file = paste0(dir_moa, "/data/final/", "res_wgcna_mic.RData"))

MRN_if

prep

# prep data
exp_mrn_4wgcna_if <- gnp.Isets$gnp_if1$exp_If1_mrn # 

# QC
qc_mrn_if <- goodSamplesGenes(exp_mrn_4wgcna_if)
qc_mrn_if$allOK # all oK!

analyse soft power

sft_mrn_if <- pickSoftThreshold(exp_mrn_4wgcna_if, powerVector = power_single, networkType = "signed", verbose = 5)

# visualize
sft.data_mrn_if <- sft_mrn_if$fitIndices

a1_mrn_if <- sft.data_mrn_if %>% ggplot(aes(Power, SFT.R.sq, label = Power)) +
  geom_point() +
  geom_text(nudge_y = 0.3, size = 2) +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Power", y = "Scale free topology model fit, signed R^2") +
  theme_classic() +
  theme(text=element_text(size=6))

a2_mrn_if <- sft.data_mrn_if %>% ggplot((aes(Power, mean.k., label = Power))) +
  geom_point() +
  geom_text(nudge_y = 500, size = 2) +
  labs(x = "Power", y = "Mean connectivity") +
  theme_classic() +
  theme(text=element_text(size=6))

fig_fst_mrn_if_pre <- ggarrange(
    a1_mrn_if, 
    a2_mrn_if, 
    nrow = 2)

fig_fst_mrn_if <- annotate_figure(fig_fst_mrn_if_pre,
               top = text_grob("Softpower mrn_if", color = "black", face = "bold", size = 12),)

fig_fst_mrn_if # determine: sft_mrn_if = 10

constants

soft_power_mrn_if <- 10

## numbers for pval calc
nSamples_mrn_if <- nrow(exp_mrn_4wgcna_if)
nSamples_mrn_if
nGenes_mrn_if <- ncol(exp_mrn_4wgcna_if)
nGenes_mrn_if

network

# make matrix numeric
exp_mrn_4wgcna_if[] <- sapply(exp_mrn_4wgcna_if, as.numeric)

# # get blockModules from 0.2 - 0.5
# res.bwnet_mrn_if <- lapply(seq(0.2, 0.5, by=0.1), get_blockwiseModules, exp = exp_mrn_4wgcna_if, soft_power = soft_power_mrn_if)
# names(res.bwnet_mrn_if) <- c('point2', 'point3', 'point4', 'point5')
# 
# save(sft.data_mrn_if,
#      res.bwnet_mrn_if,
#      file = paste0(dir_moa, "/data/final/preps/wgcna/", "mrn_if_bwnet.RData"))

load(file = paste0(dir_moa, "/data/final/preps/wgcna/", "mrn_if_bwnet.RData"))

dendro

plot_WGCNA_dendro(res.bwnet_mrn_if$point2)

plot_WGCNA_dendro(res.bwnet_mrn_if$point3)

plot_WGCNA_dendro(res.bwnet_mrn_if$point4)

plot_WGCNA_dendro(res.bwnet_mrn_if$point5)

traits to associate:

Same as mrn!

Cor: modules & traits

## module eigengenes
MEs_mrn_if_0.2 <- res.bwnet_mrn_if$point2$MEs
MEs_mrn_if_0.3 <- res.bwnet_mrn_if$point3$MEs
MEs_mrn_if_0.4 <- res.bwnet_mrn_if$point4$MEs
MEs_mrn_if_0.5 <- res.bwnet_mrn_if$point5$MEs

## correlation of eigengenes to modules
EG2M_mrn_if_0.2 <- extract_Eigengene_cor2Module(MEs_mrn_if_0.2, traits_mrn, nSamples_mrn_if)
EG2M_mrn_if_0.3 <- extract_Eigengene_cor2Module(MEs_mrn_if_0.3, traits_mrn, nSamples_mrn_if)
EG2M_mrn_if_0.4 <- extract_Eigengene_cor2Module(MEs_mrn_if_0.4, traits_mrn, nSamples_mrn_if)
EG2M_mrn_if_0.5 <- extract_Eigengene_cor2Module(MEs_mrn_if_0.5, traits_mrn, nSamples_mrn_if)


## prep for heatmap
datHM_mrn_if_0.2 <- prep_HM_EigG2Mod(MEs_mrn_if_0.2, traits_mrn)
datHM_mrn_if_0.3 <- prep_HM_EigG2Mod(MEs_mrn_if_0.3, traits_mrn)
datHM_mrn_if_0.4 <- prep_HM_EigG2Mod(MEs_mrn_if_0.4, traits_mrn)
datHM_mrn_if_0.5 <- prep_HM_EigG2Mod(MEs_mrn_if_0.5, traits_mrn)

# plot heatmap (need to pick cols!)
names(datHM_mrn_if_0.2)
CorLevelPlot(datHM_mrn_if_0.2,
           x = names(datHM_mrn_if_0.2)[12:18], # picks traits
           y = names(datHM_mrn_if_0.2)[1:11], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mrn_if (cut 0.2)")

names(datHM_mrn_if_0.3)
CorLevelPlot(datHM_mrn_if_0.3,
           x = names(datHM_mrn_if_0.3)[8:14], # picks traits
           y = names(datHM_mrn_if_0.3)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mrn_ifA (cut 0.3)")

names(datHM_mrn_if_0.4)
CorLevelPlot(datHM_mrn_if_0.4,
           x = names(datHM_mrn_if_0.4)[7:13], # picks traits
           y = names(datHM_mrn_if_0.4)[1:6], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mrn_ifA (cut 0.4)")

names(datHM_mrn_if_0.5)
CorLevelPlot(datHM_mrn_if_0.5,
           x = names(datHM_mrn_if_0.5)[7:13], # picks traits
           y = names(datHM_mrn_if_0.5)[1:6], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits mrn_ifA (cut 0.5)")

Genes per Module

# requires san: both defined above
# san_mrn
# levels_mrn

EGpM_mrn_if_0.2 <- extract_EigenG_per_Module(MEs_mrn_if_0.2, san_mrn, levels_mrn)
EGpM_mrn_if_0.3 <- extract_EigenG_per_Module(MEs_mrn_if_0.3, san_mrn, levels_mrn)
EGpM_mrn_if_0.4 <- extract_EigenG_per_Module(MEs_mrn_if_0.4, san_mrn, levels_mrn)
EGpM_mrn_if_0.5 <- extract_EigenG_per_Module(MEs_mrn_if_0.5, san_mrn, levels_mrn)

# total number of genes per module
plotter_genesPERmodule(res.bwnet_mrn_if$point2, "mrn_if 0.2")+ ylim(0, 1500)

plotter_genesPERmodule(res.bwnet_mrn_if$point3, "mrn_if 0.3")+ ylim(0, 1500)

plotter_genesPERmodule(res.bwnet_mrn_if$point4, "mrn_if 0.4")+ ylim(0, 1500)

plotter_genesPERmodule(res.bwnet_mrn_if$point5, "mrn_if 0.5")+ ylim(0, 1500)

Modules

0.2

plot_ME1_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME1', "ME1: 1242 features")
plot_ME2_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME2', "ME2: 308 features")
plot_ME3_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME3', "ME3: 259 features")
plot_ME4_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME4', "ME4: 170 features")
plot_ME5_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME5', "ME5: 119 features")
plot_ME6_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME6', "ME6: 94 features")
plot_ME7_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME7', "ME7: 78 features")
plot_ME8_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME8', "ME8: 72 features")
plot_ME9_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME9', "ME9: 48 features")
plot_ME10_mrn_if_0.2 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.2, 'ME10', "ME10: 36 features")

fig_EigGinModules_mrn_if_0.2.pt1 <- ggarrange(
  plot_ME1_mrn_if_0.2,
  plot_ME2_mrn_if_0.2,
  plot_ME3_mrn_if_0.2,
  plot_ME4_mrn_if_0.2,
  plot_ME5_mrn_if_0.2,
  plot_ME6_mrn_if_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_if_0.2.pt1 <- annotate_figure(
  fig_EigGinModules_mrn_if_0.2.pt1,
  top = text_grob("Modules mrn IF (cut 0.2), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.2.pt1

fig_EigGinModules_mrn_if_0.2.pt2 <- ggarrange(
  plot_ME6_mrn_if_0.2,
  plot_ME7_mrn_if_0.2,
  plot_ME8_mrn_if_0.2,
  plot_ME9_mrn_if_0.2,
  plot_ME10_mrn_if_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_mrn_if_0.2.pt2 <- annotate_figure(
  fig_EigGinModules_mrn_if_0.2.pt2,
  top = text_grob("Modules mrn IF (cut 0.2), pt2", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.2.pt2

0.3

plot_ME1_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME1', "ME1: 1242 features")
plot_ME2_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME2', "ME2: 734 features")
plot_ME3_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME3', "ME3: 170 features")
plot_ME4_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME4', "ME4: 166 features")
plot_ME5_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME5', "ME5: 78 features")
plot_ME6_mrn_if_0.3 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.3, 'ME6', "ME6: 36 features")

# first test: to look at them
fig_EigGinModules_mrn_if_0.3 <- ggarrange(
  plot_ME1_mrn_if_0.3,
  plot_ME2_mrn_if_0.3,
  plot_ME3_mrn_if_0.3,
  plot_ME4_mrn_if_0.3,
  plot_ME5_mrn_if_0.3,
  plot_ME6_mrn_if_0.3,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_if_0.3 <- annotate_figure(
  fig_EigGinModules_mrn_if_0.3,
  top = text_grob("Modules mrn IF (cut 0.3)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.3

0.4

plot_ME1_mrn_if_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.4, 'ME1', "ME1: 1242 features")
plot_ME2_mrn_if_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.4, 'ME2', "ME2: 734 features")
plot_ME3_mrn_if_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.4, 'ME3', "ME3: 206 features")
plot_ME4_mrn_if_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.4, 'ME4', "ME4: 166 features")
plot_ME5_mrn_if_0.4 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.4, 'ME5', "ME5: 78 features")

fig_EigGinModules_mrn_if_0.4 <- ggarrange(
  plot_ME1_mrn_if_0.4 + border_up,
  plot_ME2_mrn_if_0.4 + border_down,
  plot_ME3_mrn_if_0.4 + border(color = 'black'),
  plot_ME4_mrn_if_0.4 + border_down,
  plot_ME5_mrn_if_0.4 + border_down,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_if_0.4 <- annotate_figure(
  fig_EigGinModules_mrn_if_0.4,
  top = text_grob("Modules mrn IF (cut 0.4)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.4

## final
fig_EigGinModules_mrn_if_0.4_fin <- ggarrange(
  plot_ME1_mrn_if_0.4 + border_up,
  NULL,
  NULL,
  plot_ME2_mrn_if_0.4 + border_down,
  plot_ME4_mrn_if_0.4 + border_down,
  plot_ME5_mrn_if_0.4 + border_down,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_if_0.4_fin <- annotate_figure(
  fig_EigGinModules_mrn_if_0.4_fin,
  top = text_grob("Modules mrn IF (cut 0.4)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.4_fin

0.5

plot_ME1_mrn_if_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.5, 'ME1', "ME1: 1242 features")
plot_ME2_mrn_if_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.5, 'ME2', "ME2: 734 features")
plot_ME3_mrn_if_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.5, 'ME3', "ME3: 206 features")
plot_ME4_mrn_if_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.5, 'ME4', "ME4: 166 features")
plot_ME5_mrn_if_0.5 <- plotter_eigengenesINmodules(EGpM_mrn_if_0.5, 'ME5', "ME5: 78 features")


fig_EigGinModules_mrn_if_0.5 <- ggarrange(
  plot_ME1_mrn_if_0.5,
  plot_ME2_mrn_if_0.5,
  plot_ME3_mrn_if_0.5,
  plot_ME4_mrn_if_0.5,
  plot_ME5_mrn_if_0.5,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_mrn_if_0.5 <- annotate_figure(
  fig_EigGinModules_mrn_if_0.5,
  top = text_grob("Modules mrn IF (cut 0.5)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_mrn_if_0.5

Features of modules

fom_mrn_if_0.2 <- extract_featuresOFmodules(res.bwnet_mrn_if$point2)
fom_mrn_if_0.3 <- extract_featuresOFmodules(res.bwnet_mrn_if$point3)
fom_mrn_if_0.4 <- extract_featuresOFmodules(res.bwnet_mrn_if$point4)
fom_mrn_if_0.5 <- extract_featuresOFmodules(res.bwnet_mrn_if$point5)

fom_mrn_if <- list(
  fom_mrn_if_0.2 = fom_mrn_if_0.2,
  fom_mrn_if_0.3 = fom_mrn_if_0.3,
  fom_mrn_if_0.4 = fom_mrn_if_0.4,
  fom_mrn_if_0.5 = fom_mrn_if_0.5
)

Features of traits

Can give the most important feature of each trait.

Cave: this is NOT module specific but only trait specific. Tables have to be filtered by pValue.

Not 100% sure how to use this info.

fot_MB_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'MB', nSamples_mrn_if)
fot_PM_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'PM', nSamples_mrn_if)
fot_MC_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'MC', nSamples_mrn_if)
fot_MM_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'MM', nSamples_mrn_if)
fot_B_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'B', nSamples_mrn_if)
fot_S_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'S', nSamples_mrn_if)
fot_PMN_mrn_if <- extract_features_of_traits(exp_mrn_4wgcna_if, traits_mrn, 'PMN', nSamples_mrn_if)

fot_mrn_if <- list(
  fot_MB_mrn_if = fot_MB_mrn_if,
  fot_PM_mrn_if = fot_PM_mrn_if,
  fot_MC_mrn_if = fot_MC_mrn_if,
  fot_MM_mrn_if = fot_MM_mrn_if,
  fot_B_mrn_if = fot_B_mrn_if,
  fot_S_mrn_if = fot_S_mrn_if,
  fot_PMN_mrn_if = fot_PMN_mrn_if
)

heatmaps

0.2

# 
# plot_hm_0.2_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)

0.4

# plot_hm_0.4_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)

export

res.mrn_if <- list(
  fom_mrn_if = fom_mrn_if,
  fot_mrn_if = fot_mrn_if
)

pro_if

prep

# prep data

exp_pro_4wgcna_if<- gnp.Isets$gnp_if1$exp_If1_pro # using the TNP data

colnames(exp_pro_4wgcna_if) # SYMBOLS = ok
rownames(exp_pro_4wgcna_if) # samples = ok

# QC
qc_pro_if <- goodSamplesGenes(exp_pro_4wgcna_if)
qc_pro_if$allOK # all oK!

analyse soft power

sft_pro_if <- pickSoftThreshold(exp_pro_4wgcna_if, powerVector = power_single, networkType = "signed", verbose = 5)

# visualize
sft.data_pro_if <- sft_pro_if$fitIndices

a1_pro_if <- sft.data_pro_if %>% ggplot(aes(Power, SFT.R.sq, label = Power)) +
  geom_point() +
  geom_text(nudge_y = 0.3, size = 2) +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(x = "Power", y = "Scale free topology model fit, signed R^2") +
  theme_classic() +
  theme(text=element_text(size=6))

a2_pro_if <- sft.data_pro_if %>% ggplot((aes(Power, mean.k., label = Power))) +
  geom_point() +
  geom_text(nudge_y = 500, size = 2) +
  labs(x = "Power", y = "Mean connectivity") +
  theme_classic() +
  theme(text=element_text(size=6))

fig_fst_pro_if_pre <- ggarrange(
    a1_pro_if, 
    a2_pro_if, 
    nrow = 2)

fig_fst_pro_if <- annotate_figure(fig_fst_pro_if_pre,
               top = text_grob("Softpower pro_if", color = "black", face = "bold", size = 12),)

fig_fst_pro_if # determine: sft_pro_if = 5

constants

soft_power_pro_if <- 5

## numbers for pval calc
nSamples_pro_if <- nrow(exp_pro_4wgcna_if)
nSamples_pro_if
nGenes_pro_if <- ncol(exp_pro_4wgcna_if)
nGenes_pro_if

network

# make matrix numeric
exp_pro_4wgcna_if[] <- sapply(exp_pro_4wgcna_if, as.numeric)

# # get blockModules from 0.2 - 0.5
# res.bwnet_pro_if <- lapply(seq(0.2, 0.5, by=0.1), get_blockwiseModules, exp = exp_pro_4wgcna_if, soft_power = soft_power_pro_if)
# names(res.bwnet_pro_if) <- c('point2', 'point3', 'point4', 'point5')
# 
# save(sft.data_pro_if,
#      res.bwnet_pro_if,
#      file = paste0(dir_moa, "/data/final/preps/wgcna/", "pro_if_bwnet.RData"))

load(file = paste0(dir_moa, "/data/final/preps/wgcna/", "pro_if_bwnet.RData"))

dendro

plot_WGCNA_dendro(res.bwnet_pro_if$point2)

plot_WGCNA_dendro(res.bwnet_pro_if$point3)

plot_WGCNA_dendro(res.bwnet_pro_if$point4)

plot_WGCNA_dendro(res.bwnet_pro_if$point5)

traits to associate:

Same as mrn!

Cor: modules & traits

## module eigengenes
MEs_pro_if_0.2 <- res.bwnet_pro_if$point2$MEs
MEs_pro_if_0.3 <- res.bwnet_pro_if$point3$MEs
MEs_pro_if_0.4 <- res.bwnet_pro_if$point4$MEs
MEs_pro_if_0.5 <- res.bwnet_pro_if$point5$MEs

## correlation of eigengenes to modules
EG2M_pro_if_0.2 <- extract_Eigengene_cor2Module(MEs_pro_if_0.2, traits_pro, nSamples_pro_if)
EG2M_pro_if_0.3 <- extract_Eigengene_cor2Module(MEs_pro_if_0.3, traits_pro, nSamples_pro_if)
EG2M_pro_if_0.4 <- extract_Eigengene_cor2Module(MEs_pro_if_0.4, traits_pro, nSamples_pro_if)
EG2M_pro_if_0.5 <- extract_Eigengene_cor2Module(MEs_pro_if_0.5, traits_pro, nSamples_pro_if)


## prep for heatmap
datHM_pro_if_0.2 <- prep_HM_EigG2Mod(MEs_pro_if_0.2, traits_pro)
datHM_pro_if_0.3 <- prep_HM_EigG2Mod(MEs_pro_if_0.3, traits_pro)
datHM_pro_if_0.4 <- prep_HM_EigG2Mod(MEs_pro_if_0.4, traits_pro)
datHM_pro_if_0.5 <- prep_HM_EigG2Mod(MEs_pro_if_0.5, traits_pro)

# plot heatmap (need to pick cols!)
names(datHM_pro_if_0.2)
CorLevelPlot(datHM_pro_if_0.2,
           x = names(datHM_pro_if_0.2)[14:19], # picks traits
           y = names(datHM_pro_if_0.2)[1:13], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits pro_if (cut 0.2)")

names(datHM_pro_if_0.3)
CorLevelPlot(datHM_pro_if_0.3,
           x = names(datHM_pro_if_0.3)[12:17], # picks traits
           y = names(datHM_pro_if_0.3)[1:11], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits pro_ifA (cut 0.3)")

names(datHM_pro_if_0.4)
CorLevelPlot(datHM_pro_if_0.4,
           x = names(datHM_pro_if_0.4)[9:14], # picks traits
           y = names(datHM_pro_if_0.4)[1:8], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits pro_ifA (cut 0.4)")

names(datHM_pro_if_0.5)
CorLevelPlot(datHM_pro_if_0.5,
           x = names(datHM_pro_if_0.5)[8:13], # picks traits
           y = names(datHM_pro_if_0.5)[1:7], # picks eigengene values
           col = c('blue3', 'skyblue', 'white', 'pink', 'red'),
           main = "Cor. modules/traits pro_ifA (cut 0.5)")

Genes per Module

# requires san & levels from above

EGpM_pro_if_0.2 <- extract_EigenG_per_Module(MEs_pro_if_0.2, san_pro, levels_pro)
EGpM_pro_if_0.3 <- extract_EigenG_per_Module(MEs_pro_if_0.3, san_pro, levels_pro)
EGpM_pro_if_0.4 <- extract_EigenG_per_Module(MEs_pro_if_0.4, san_pro, levels_pro)
EGpM_pro_if_0.5 <- extract_EigenG_per_Module(MEs_pro_if_0.5, san_pro, levels_pro)

# total number of genes per module
plotter_genesPERmodule(res.bwnet_pro_if$point2, "pro_if 0.2")+ ylim(0, 600)

plotter_genesPERmodule(res.bwnet_pro_if$point3, "pro_if 0.3")+ ylim(0, 600)

plotter_genesPERmodule(res.bwnet_pro_if$point4, "pro_if 0.4")+ ylim(0, 1000)

plotter_genesPERmodule(res.bwnet_pro_if$point5, "pro_if 0.5")+ ylim(0, 1000)

Modules

0.2

plot_ME1_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME1', "ME1: 471 features")
plot_ME2_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME2', "ME2: 414 features")
plot_ME3_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME3', "ME3: 374 features")
plot_ME4_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME4', "ME4: 246 features")
plot_ME5_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME5', "ME5: 164 features")
plot_ME6_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME6', "ME6: 162 features")
plot_ME7_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME7', "ME7: 116 features")
plot_ME8_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME8', "ME8: 93 features")
plot_ME9_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME9', "ME9: 93 features")
plot_ME10_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME10', "ME10: 63 features")
plot_ME11_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME11', "ME11: 35 features")
plot_ME12_pro_if_0.2 <- plotter_eigengenesINmodules(EGpM_pro_if_0.2, 'ME12', "ME12: 30 features")

fig_EigGinModules_pro_if_0.2.pt1 <- ggarrange(
  plot_ME1_pro_if_0.2,
  plot_ME2_pro_if_0.2,
  plot_ME3_pro_if_0.2,
  plot_ME4_pro_if_0.2,
  plot_ME5_pro_if_0.2,
  plot_ME6_pro_if_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)
  
fig_EigGinModules_pro_if_0.2.pt1 <- annotate_figure(
  fig_EigGinModules_pro_if_0.2.pt1,
  top = text_grob("Modules pro IF (cut 0.2), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.2.pt1

fig_EigGinModules_pro_if_0.2.pt2 <- ggarrange(
  plot_ME7_pro_if_0.2,
  plot_ME8_pro_if_0.2,
  plot_ME9_pro_if_0.2,
  plot_ME10_pro_if_0.2,
  plot_ME11_pro_if_0.2,
  plot_ME12_pro_if_0.2,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_if_0.2.pt2 <- annotate_figure(
  fig_EigGinModules_pro_if_0.2.pt2,
  top = text_grob("Modules pro IF (cut 0.2), pt2", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.2.pt2

# # fig FC level comparison
# fig_EigGinModules_pro_0.2_pre <- ggarrange(
#   plot_ME2_pro_0.2,
#   plot_ME3_pro_0.2,
#   plot_ME7_pro_0.2,
#   plot_ME4_pro_0.2,
#   plot_ME5_pro_0.2,
#   plot_ME9_pro_0.2,
#   plot_ME10_pro_0.2,
#   plot_ME11_pro_0.2,
#   plot_ME13_pro_0.2,
#   plot_ME1_pro_0.2,
#   NULL,
#   NULL,
#   plot_ME6_pro_0.2,
#   plot_ME8_pro_0.2,
#   plot_ME12_pro_0.2,
#   ncol = 3,
#   nrow = 5
# )
# fig_EigGinModules_pro_0.2_pre

# fig_EigGinModules_pro_if_0.2.test <- annotate_figure(fig_EigGinModules_pro_if_0.2.test,
#             top = text_grob("Eigengene expression in modules: pro if (cut 0.2)", 
#                 color = "black", face = "bold", size = 18))
# fig_EigGinModules_pro_if_0.2.test

0.3

plot_ME1_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME1', "ME1: 471 features")
plot_ME2_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME2', "ME2: 414 features")
plot_ME3_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME3', "ME3: 408 features")
plot_ME4_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME4', "ME4: 374 features")
plot_ME5_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME5', "ME5: 164 features")
plot_ME6_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME6', "ME6: 156 features")
plot_ME7_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME7', "ME7: 116 features")
plot_ME8_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME8', "ME8: 93 features")
plot_ME9_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME9', "ME9: 35 features")
plot_ME10_pro_if_0.3 <- plotter_eigengenesINmodules(EGpM_pro_if_0.3, 'ME10', "ME10: 30 features")

# first test: to look at them
fig_EigGinModules_pro_if_0.3.pt1 <- ggarrange(
  plot_ME1_pro_if_0.3,
  plot_ME2_pro_if_0.3,
  plot_ME3_pro_if_0.3,
  plot_ME4_pro_if_0.3,
  plot_ME5_pro_if_0.3,
  plot_ME6_pro_if_0.3,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_if_0.3.pt1 <- annotate_figure(
  fig_EigGinModules_pro_if_0.3.pt1,
  top = text_grob("Modules pro IF (cut 0.3), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.3.pt1

fig_EigGinModules_pro_if_0.3.pt2 <- ggarrange(
  plot_ME7_pro_if_0.3,
  plot_ME8_pro_if_0.3,
  plot_ME9_pro_if_0.3,
  plot_ME10_pro_if_0.3,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_if_0.3.pt2 <- annotate_figure(
  fig_EigGinModules_pro_if_0.3.pt2,
  top = text_grob("Modules pro IF (cut 0.3), pt2", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.3.pt2

0.4

plot_ME1_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME1', "ME1: 885 features")
plot_ME2_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME2', "ME2: 408 features")
plot_ME3_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME3', "ME3: 374 features")
plot_ME4_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME4', "ME4: 280 features")
plot_ME5_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME5', "ME5: 156 features")
plot_ME6_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME6', "ME6: 128 features")
plot_ME7_pro_if_0.4 <- plotter_eigengenesINmodules(EGpM_pro_if_0.4, 'ME7', "ME7: 30 features")

fig_EigGinModules_pro_if_0.4.pt1 <- ggarrange(
  plot_ME1_pro_if_0.4,
  plot_ME2_pro_if_0.4,
  plot_ME3_pro_if_0.4,
  plot_ME4_pro_if_0.4,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_if_0.4.pt1 <- annotate_figure(
  fig_EigGinModules_pro_if_0.4.pt1,
  top = text_grob("Modules pro IF (cut 0.4), pt1", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.4.pt1

fig_EigGinModules_pro_if_0.4.pt2 <- ggarrange(
  plot_ME5_pro_if_0.4,
  plot_ME6_pro_if_0.4,
  plot_ME7_pro_if_0.4,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)

fig_EigGinModules_pro_if_0.4.pt2 <- annotate_figure(
  fig_EigGinModules_pro_if_0.4.pt2,
  top = text_grob("Modules pro IF (cut 0.4), pt2", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.4.pt2

0.5

plot_ME1_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME1', "ME1: 885 features")
plot_ME2_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME2', "ME2: 689 features")
plot_ME3_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME3', "ME3: 408 features")
plot_ME4_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME4', "ME4: 156 features")
plot_ME5_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME5', "ME5: 93 features")
plot_ME6_pro_if_0.5 <- plotter_eigengenesINmodules(EGpM_pro_if_0.5, 'ME6', "ME6: 30 features")


fig_EigGinModules_pro_if_0.5 <- ggarrange(
  plot_ME1_pro_if_0.5,
  plot_ME2_pro_if_0.5,
  plot_ME3_pro_if_0.5,
  plot_ME4_pro_if_0.5,
  plot_ME5_pro_if_0.5,
  plot_ME6_pro_if_0.5,
  common.legend = T,
  legend = "none",
  ncol = 3,
  nrow = 2)


fig_EigGinModules_pro_if_0.5 <- annotate_figure(
  fig_EigGinModules_pro_if_0.5,
  top = text_grob("Modules pro IF (cut 0.5)", color = "black", face = "bold", size = 18))
                
fig_EigGinModules_pro_if_0.5

Features of modules

fom_pro_if_0.2 <- extract_featuresOFmodules(res.bwnet_pro_if$point2)
fom_pro_if_0.3 <- extract_featuresOFmodules(res.bwnet_pro_if$point3)
fom_pro_if_0.4 <- extract_featuresOFmodules(res.bwnet_pro_if$point4)
fom_pro_if_0.5 <- extract_featuresOFmodules(res.bwnet_pro_if$point5)

fom_pro_if <- list(
  fom_pro_if_0.2 = fom_pro_if_0.2,
  fom_pro_if_0.3 = fom_pro_if_0.3,
  fom_pro_if_0.4 = fom_pro_if_0.4,
  fom_pro_if_0.5 = fom_pro_if_0.5
)

Features of traits

Can give the most important feature of each trait.

Cave: this is NOT module specific but only trait specific. Tables have to be filtered by pValue.

Not 100% sure how to use this info.

fot_MB_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'MB', nSamples_pro_if)
fot_PM_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'PM', nSamples_pro_if)
fot_MC_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'MC', nSamples_pro_if)
fot_MM_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'MM', nSamples_pro_if)
fot_B_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'B', nSamples_pro_if)
fot_PMN_pro_if <- extract_features_of_traits(exp_pro_4wgcna_if, traits_pro, 'PMN', nSamples_pro_if)

fot_pro_if <- list(
  fot_MB_pro_if = fot_MB_pro_if,
  fot_PM_pro_if = fot_PM_pro_if,
  fot_MC_pro_if = fot_MC_pro_if,
  fot_MM_pro_if = fot_MM_pro_if,
  fot_B_pro_if = fot_B_pro_if,
  fot_PMN_pro_if = fot_PMN_pro_if
)

heatmaps

0.2

# plot_hm_0.2_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)
# plot_hm_0.2_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.2, EGpM_mrn_0.2)

0.4

# plot_hm_0.4_ME1 <- make_module_heatmap("ME1", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME2 <- make_module_heatmap("ME2", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME3 <- make_module_heatmap("ME3", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME4 <- make_module_heatmap("ME4", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME5 <- make_module_heatmap("ME5", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)
# plot_hm_0.4_ME6 <- make_module_heatmap("ME6", exp_mrn_4wgcna, san_mrn, fom_mrn_0.4, EGpM_mrn_0.4)

export

res.pro_if <- list(
  fom_pro_if = fom_pro_if,
  fot_pro_if = fot_pro_if
)

export

# res.mrn
# res.mic
# res.pro
# res.mrn_if
# res.pro_if

res.wgcna <- list(
  res.mrn = res.mrn,
  res.pro = res.pro,
  res.mic = res.mic,
  res.mrn_if = res.mrn_if,
  res.pro_if = res.pro_if
)

save(res.wgcna,
     file = paste0(dir_moa, "/data/final/", "res.wgcna.RData"))