#pkglibrary(tidyverse)library(readxl)library(WGCNA) library(ggcorrplot)library(moduleColor)library(genefilter) # for delet low variancelibrary(DESeq2)library(ggnewscale) # si tu veux plusieurs scales_fill plus tardlibrary(ggh4x)library(colorspace) library(pheatmap)library(WGCNA)library(doParallel)library(parallel)library(igraph)# srcsource(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot source(here::here("src/function/fig_export.R")) # This function saves a given plot (plot_x) as both a PDF and a high-resolution PNG file at specified dimensions.source(here::here("src/function/Evaluate_contrast.R"))# cosmeticssulfate_pallet=read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="sulfure_condition") %>% dplyr::select(color, treatment) %>%pull(color) %>%setNames(read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="sulfure_condition") %>%pull(treatment) )mutant_palette=read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="mutant") %>% dplyr::select(color, treatment) %>%pull(color) %>%setNames(read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="mutant") %>%pull(treatment) )my_color_palette <-read_excel(here::here("data/color_palette.xlsm"))# functionstatic_colors=function(hier1, h1=0.9,minsize1=50) {# here we define modules by using a height cut-off for the branches labelpred=cutree(hier1,h=h1) sort1=-sort(-table(labelpred)) modulename=as.numeric(names(sort1)) modulebranch= sort1>minsize1 no.modules=sum(modulebranch)# now we assume that there are fewer than a certain number of colors colorcode=GlobalStandardColors# "grey" means not in any module; colorhelp=rep("grey",length(labelpred))if ( no.modules==0| no.modules >length(colorcode)){ print(paste("The number of modules is problematic. Number of modules = ", as.character(no.modules)))} else { for (i inc(1:no.modules)) {colorhelp=ifelse(labelpred==modulename[i],colorcode[i],colorhelp)}; colorhelp=factor(colorhelp,levels=c(colorcode[1:no.modules],"grey")) }factor(colorhelp, levels=unique(colorhelp[hier1$order] ))}corFnc_i <- WGCNA::bicorclean_session <-function() {# Get all objects in the global environment objects <-ls(envir = .GlobalEnv)# Check if each object is a data frame or a matrix is_data <-vapply(objects, function(x) { obj <-get(x, envir = .GlobalEnv)is.data.frame(obj) ||is.matrix(obj) }, logical(1)) # Ensure logical vector output# Remove the identified objectsrm(list = objects[is_data], envir = .GlobalEnv)# Reset all plotswhile (dev.cur() >1) dev.off()cat_col("Session cleaned: data and plots removed.", "green")}add_prefix_except_plant_num <-function(type_of_date, df) { df_renamed <- df %>% dplyr::rename_with(~paste0(type_of_date, "_", .x),.cols =-plant_num )return(df_renamed)}correlation_summary <-function(df) { var_names <-colnames(df) results <-combn(var_names, 2, simplify =FALSE) |> purrr::map_df(function(vars) { x <- df[[vars[1]]] y <- df[[vars[2]]] ok <-complete.cases(x, y)if (sum(ok) >=3) { test <-cor.test(x[ok], y[ok], method ="pearson") tibble::tibble(V1 = vars[1],V2 = vars[2],corr =unname(test$estimate),pval = test$p.value,R2 = test$estimate^2 ) } else { tibble::tibble(V1 = vars[1],V2 = vars[2],corr =NA,pval =NA,R2 =NA ) } })# Add FDR correction results$fdr <-p.adjust(results$pval, method ="fdr")return(results)}
While WGCNA was initially developed for processing gene expression data, it can be adapted for other data types as well.
Weighted Gene Co-expression Network Analysis (WGCNA) is an approach dedicated to identifying correlation patterns between features of a dataset. In most cases, it is applied to gene expression data to find clusters -referred to as modules- of highly correlated genes. In transcriptomics, genes that exhibit strong correlations in their expression profiles are said co-expressed. One of the underlying principles of WGCNA is the concept of guilt-by-association. This principle suggests that genes sharing similar expression patterns are likely to be functionally related or involved in similar biological processes. By leveraging guilt-by-association through correlation, WGCNA provides a powerful tool for identifying potential functional relationships among genes and deciphering their roles in complex biological systems.
WGCNA builds a co-expression network where:
each feature (transcript/gene in transcriptomics data) is represented as a node
edges represent a correlation pattern between two features (nodes) and are weighted with the corresponding correlation coefficient, indicating the strength of their observed correlation.
16.2 Data importation
I keep only genes with variance > 0.1 across samples
Normalization of data according to the type of data and the type of analysis to perform. For RNA-seq data, WGCNA recommends to use the varianceStabilizingTransformation from the DESeq2 package, or to perform a log-transformation of normalized counts (RPFK/FPKM). I tested with and without normalization. May be it’s better if i rerun my script with VST normlisation. Some person say it’s RFPKM is a load methode (see reference below)
samples (e.g. patients, organisms …) are displayed on the rows
variables (e.g. genes, proteins …) are displayed on the columns
Code
load(file = here::here("data/multi_omic/WGCNA/expr_matrix.RData"))datExpr <-t(expr_matrix)#datExpr <- t(counts(dds, normalized = TRUE)) # perform the median of ratios method #datExpr_vst <- t(assay(vst)) # with varianceStabilizingTransformation# datExpr_fpkm <- t(fpkm) # with varianceStabilizingTransformation# Run WGCNA quality check## WGCNA provides a function goodSamplesGenes to check the quality of input data and remove features and samples with too many missing data as well as genes with zero variance.gsg =goodSamplesGenes(datExpr, verbose=3) # link zero variance samples, for example here#gsg_vst = goodSamplesGenes(datExpr_vst, verbose=3) #gsg_fpkm = goodSamplesGenes(datExpr_fpkm, verbose=3) # link zero variance samples, for example here# tchek if everithing is okif (!gsg$allOK){if (sum(!gsg$goodGenes)>0)printFlush(paste("Removing genes:", paste(names(datExpr)[!gsg$goodGenes], collapse =", ")));if (sum(!gsg$goodSamples)>0)printFlush(paste("Removing samples:", paste(rownames(datExpr)[!gsg$goodSamples], collapse =", "))); datExpr = datExpr[gsg$goodSamples, gsg$goodGenes]}# if (!gsg_vst$allOK){# if (sum(!gsg_vst$goodGenes)>0)# printFlush(paste("Removing genes:", paste(names(datExpr_vst)[!gsg_vst$goodGenes], collapse = ", ")));# if (sum(!gsg_vst$goodSamples)>0)# printFlush(paste("Removing samples:", paste(rownames(datExpr_vst)[!gsg_vst$goodSamples], collapse = ", ")));# datExpr = datExpr_vst[gsg_vst$goodSamples, gsg_vst$goodGenes]# }# # if (!gsg_fpkm$allOK){# if (sum(!gsg_fpkm$goodGenes)>0)# printFlush(paste("Removing genes:", paste(names(datExpr_fpkm)[!gsg_fpkm$goodGenes], collapse = ", ")));# if (sum(!gsg_vst$goodSamples)>0)# printFlush(paste("Removing samples:", paste(rownames(datExpr_fpkm)[!gsg_fpkm$goodSamples], collapse = ", ")));# datExpr = datExpr_fpkm[gsg_fpkm$goodSamples, gsg_fpkm$goodGenes]# }save(datExpr, file = here::here("data/multi_omic/WGCNA/datExpr.RData"))
Data verification according to WGCNA criteria
Code
load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))sampleTree =hclust(dist(datExpr), method ="average")#sampleTree_vst = hclust(dist(datExpr_vst), method = "average")#sampleTree_fpkm = hclust(dist(datExpr_fpkm), method = "average")# Plot the sample tree: Open a graphic output window of size 12 by 9 inches# The user should change the dimensions if the window is too large or too small.png(here::here(paste0("report/multi_omics/plot/WGCNA/sample_clustering_to_detect_outliers.png")),width =2500,height =2000,res =300)plot(sampleTree, main ="Sample clustering to detect outliers", sub="", xlab="", cex.lab =1.5,cex.axis =1.5, cex.main =2)dev.off()# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_vst.png")),width = 2500,height = 2000,res = 300)# plot(sampleTree_vst, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,# cex.axis = 1.5, cex.main = 2)# dev.off()# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_fpkm.png")),width = 2500,height = 2000,res = 300)# plot(sampleTree_fpkm, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,# cex.axis = 1.5, cex.main = 2)# dev.off()
Then you can use Pearson correlation or Spearman correlation.
The Pearson correlation measures linear associations between variables, whereas Spearman correlation measures monotonous associations between variables. Both can have interesting results, so it's completely up to you! By default, WGCNA uses the Pearson correlation. If you want to use the Spearman correlation, don't forget to specify it each time you notice a method parameter controlling the correlation method in one of the WGCNA functions!
16.4.1 Computing the correlation network (summary)
This first step of the analysis consists in creating the feature correlation network. This first step is composed of tree main substeps:
Power selection
Soft thresholding is a strategy introduced by WGCNA to transform correlation coefficients for building a weighted correlation network
The threshold should be set such that it makes the network scale-free
Similarity matrix construction
First, correlation are calculated between each pair of feature
Then, correlation or distance are transformed, using the selected soft-power
Topological overlap matrix construction (TOM)
The similarity matrix is transformed into a TOM to facilitate the module detection
16.4.2 Azrael code
Script use in server Azrael
You’ll probably be limited by the amount of RAM required to create the TOM. It’s best to use the following script on a server and then re-import the data.
Code
# Function of need in memoryBaseColors =c("turquoise","blue","brown","yellow","green","red","black","pink","magenta","purple","greenyellow","tan","salmon","cyan", "midnightblue", "lightcyan","grey60", "lightgreen", "lightyellow", "royalblue", "darkred", "darkgreen", "darkturquoise", "darkgrey", "orange", "darkorange","white", "skyblue", "saddlebrown", "steelblue", "paleturquoise", "violet", "darkolivegreen","darkmagenta", "white" );RColors =colors()[-grep("grey", colors())];RColors = RColors[-grep("gray", RColors)];InBase =match(BaseColors, RColors);ExtraColors = RColors[-c(InBase[!is.na(InBase)])];No.Extras =length(ExtraColors);# Here is the vector of colors that should be used by all functions:GlobalStandardColors =c(BaseColors, ExtraColors[rank(sin(13*c(1:No.Extras) +sin(13*c(1:No.Extras))) )] );rm(BaseColors, RColors, ExtraColors, No.Extras);############## Package and function (end) ################# Inputinput_name<-"datExpr"# based on previous script# Options WGCNAoptions(stringsAsFactors =FALSE);enableWGCNAThreads(8)# memory.limit(size=10000000)# Network parametersnetworkType_i <-"signed"#signed #unsigned# Module detection parameterminModuleSize <-30ch =0.99# for cutheight# Load WGCNA ready input dataload(file =paste0("data/multi_omic/WGCNA/input_",input_name,".RData")) # data previously created# Call the network topology analysis function.# You might want to change de minimal R^2 and the probed power value according# to your network topology.RsquaredCut =0.7powers =c(c(1:11), seq(from =12, to=40, by=2))cat("Begin of the script \n")# Correlation networks are by nature complete, which means every feature is associated to every other with a weight given by the correlation coefficient. To extract only meaningful associations, the correlation network must be thresholded.# # The most straight-forward way of thresholding a correlation dataset is to perform a hard-thresholding, i.e. two features are considered correlated if their absolute correlation coefficient exceeds a user-defined threshold.# # WGCNA propose another strategy, called soft-thresholding. It consists of raising the correlation coefficients to a certain power. This process transforms low correlation coefficients to values close to 0, while high correlation coefficients are transformed to values close to 1.# # Furthermore, WGCNA presents a function called pickSoftThreshold that aims to automatically determine the ideal power for your dataset. The optimal power is defined as the one that ensures your network exhibits a scale-free property. Let’s use it!# # First, we need to define a range of powers to try.cat("Power selection")# Sometimes, WGCNA can’t find the optimal soft-threshold. In such cases, the variable sft$powerEstimate contains NA. If the optimal soft-power was estimated but seems low, you might want to increase it. In both cases, you can plot the following figures to help you decide:# sft <- pickSoftThreshold(dataset, powerVector = powers, verbose = 5) # or allowWGCNAThreads(nThreads =12) # adapte à ton système # Sur mon macsft =pickSoftThreshold(get(input_name), #takes normalized data, samples data, creates an agacency matrix and compares with different agacency matricesRsquaredCut = RsquaredCut, powerVector = powers,corFnc = corFnc_i,networkType = networkType_i,verbose =5)# Which power should be used for your dataset? (don’t forget to look at the documentation of the pickSoftThreshold fucntion to known how).sprintf("Optimal soft-power = %d", sft$powerEstimate)# "Optimal soft-power = 11"# The optimal soft-power to use is defined by the topology of the soft-thresholded network, which has to be scale-free! Scale-free networks are networks in which the majority of nodes have a low degree, and a small number of nodes have a large degree. Multiple studies have shown that most biological networks are scale-free!# # Let’s visualize the difference between the raw correlations and the soft-thresholded correlations:corRaw =abs(cor(get(paste0(input_name))))corSoft =abs(corRaw**sft$powerEstimate)#export soft tresholdingpng(here::here(paste0("report/multi_omics/plot/WGCNA/","soft_tresholding_",input_name,"_",networkType_i,".png")),width =800, height =600)par(mfrow=c(1,2))hist(corRaw, main="Raw correlations")hist(corSoft, main="Soft-thresholded correlations")dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/", "soft_threshold_indices.png")),width =800, height =600)par(mfrow =c(1, 2)) plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab =paste("Scale Free Topology Model Fit,", networkType_i, "R^2"),type="n",main =paste("Scale independence", networkType_i));text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,col="red")top =0.9; mid =0.85; low =0.8; offset =0.02abline(h=top,col="green")text(y=top + offset, x=offset+0.5, col="green", "0.90")abline(h=mid, col="orange")text(y=mid + offset, x=offset+0.5, col="orange", "0.85")abline(h=low,col="red")text(y=low + offset, x=offset+0.5, col="red", "0.80")plot(sft$fitIndices[,1], sft$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main =paste("Mean connectivity (", networkType_i, " correlation)", sep =""))text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers,col="red")dev.off()# The Scale Independence plot (left) shows the fit indices for scale free topology for each power. The authors recommend to select a power with a R² > 0.8. The red line corresponds to using an R² cut-off of R²=0.80.# # The Mean Connectivity plot (right) shows the average of connectivity for each power. A balance needs to be identified between a fully connected network and too much disconnected edges.cat("Adjacency calculation (Similarity matrix construction)\n")# WGCNA provides the function adjacency to compute a features similarity matrix. This function operates in three steps:# # Compute features correlation coefficients# Transform the correlation coefficients depending on whether you are interested in absolute correlations (unsigned), signed correlations (signed) or positive correlations (signed hybrid). To have more details on this 3 types, check ?adjacency. By default, the constructed network is an unsigned network which means that sign are not reported into the network (absolute correlations).# Finally, the transformed correlations are raised using the previously estimated soft-threshold power.# Adjacency calculationsoftPower <- sft$powerEstimate #betaadjacency <-adjacency(datExpr =get(paste0(input_name)) , #normalise matrixtype = networkType_i,power = softPower, #betacorFnc = corFnc_i)head(adjacency[, c(1:5)])cat("Topological overlap matrix \n")# From the Similarity matrix, WGCNA proposes to extract correlated feature modules (co-expressed gene modules) based on Topological Overlap. The idea is that two features (genes) that share a high number of correlated (co-expressed) neighbors are likely to participate in the same correlation (co-expression) module.# # The Topological Overlap Matrix (TOM) is a concept used in network analysis to quantify the similarity or interconnectedness between nodes in a network.# # Transform the similarity matrix into a Topological Overlap Matrix that incorporates the level of overlapping between the neighborhoods of two nodes (features/genes), using the function TOMsimilarity. Notice that the output does not contains our row names and column names (i.e. node names). Make sure to rename it!# calculating topological superpositionTOM =TOMsimilarity(adjacency) # Because most clustering strategies are performed on distance (rather than similarities), we transform the TOM matrix into a dissimilarity matrix.dissTOM =1-TOMhead(dissTOM[, c(1:5)])cat("save disstom \n")save(dissTOM, file = here::here(paste0("data/multi_omic/",input_name,"_",networkType_i,"_dissTOM.RData)"))) #/!\ could be heavy################### Module detection #####################cat("Module detection \n")# Modules are groups of highly-connected nodes (features). This step is composed of two substeps:# Clustering: getting correlated feature modules# Use the hclust function to perform hierarchical clustering# Display the dendrogram# Use the Dynamic Tree Cut function to obtain the modules# Merging of modules: merging close modules# If they are too similar, modules should be merged# We compute module-module similarities by computing their eigenvectors (eigengene)# We merge two modules if their eigenvectors are correlated (the authors suggest a pearson corr r > 0.75)# export png genetree#png(file= paste0(input_name,"_",networkType_i,"_dissTOM_clustering-formation.png"),width = 3840 , height = 3840 ,res=300)# Call the hierarchical clustering functiongeneTree =hclust(as.dist(dissTOM), method ="average")# Plot the resulting clustering tree (dendrogram)plot(geneTree, xlab="", sub="", main ="Gene clustering on TOM-based dissimilarity",labels =FALSE, hang =0.04)dev.off()# The leaves represent the features (genes). The dendrogram branches group together densely interconnected and highly correlated features.# To identify modules from dendrograms, several methods exist. Here, we use the cutreeDynamic function.# First, we have to set the minimum module size minClusterSize (i.e. minimum number of genes inside a group/module). The default is 20, but you can try several values and choose according to you results. There are other parameters you can play with, such as deepSplit. Look at the documentation to select the parameters.########### from tutorial #############dynamicMods <-cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit =2, pamRespectsDendro =FALSE, minClusterSize = minModuleSize)table(dynamicMods)dynamicColors <-labels2colors(dynamicMods)table(dynamicColors)png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width =3840*2 , height =3840 ,res=600)plotDendroAndColors(dendro = geneTree, colors = dynamicColors, groupLabels ="Dynamic Tree Cut", dendroLabels =FALSE, hang =0.03, addGuide =TRUE, guideHang =0.05, main ="Gene dendrogram and module colors from tutorial")dev.off()# Clusturisation . there are lots of different ones. the best thing to do is to take the hybrid like everyone else.png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width =3840*2 , height =3840 ,res=600)colorh1=as.character(static_colors(geneTree,h1=ch,minsize1=minModuleSize))DynamicColor1 =labels2colors(cutreeDynamic(geneTree, cutHeight = ch,minClusterSize = minModuleSize, method ="tree", deepSplit =TRUE)); # deepSplit, regulates module coherence.# if we want to see what's inside table(DynamicColor1)DynamicColor2 =labels2colors(cutreeDynamic(geneTree, cutHeight = ch,minClusterSize = minModuleSize, method ="tree", deepSplit =FALSE));ClustColor1 =labels2colors(cutreeDynamic(dendro = geneTree, minClusterSize = minModuleSize,cutHeight = ch, method ="hybrid", maxCoreScatter =0.75, minGap =0.25, pamStage =TRUE, distM = dissTOM, useMedoids =FALSE, maxPamDist =0.9, respectSmallClusters =TRUE));Clusters1 =cutreeHybrid(dendro = geneTree, minClusterSize = minModuleSize,cutHeight = ch, maxCoreScatter =0.75, minGap =0.250, pamStage =TRUE, distM = dissTOM, useMedoids =FALSE, maxPamDist =0.9, respectSmallClusters =TRUE)CoreColor1 =labels2colors(Clusters1$cores)ClustColor2 =labels2colors(cutreeDynamic(dendro = geneTree, minClusterSize = minModuleSize,cutHeight = ch, method ="hybrid", maxCoreScatter =0.95, minGap =0.050, pamStage =TRUE, distM = dissTOM, useMedoids =FALSE, maxPamDist =0.90, respectSmallClusters =TRUE));Clusters2 =cutreeHybrid(dendro = geneTree, minClusterSize = minModuleSize,cutHeight = ch, maxCoreScatter =0.95, minGap =0.050, pamStage =TRUE, distM = dissTOM, useMedoids =FALSE, maxPamDist =0.90, respectSmallClusters =TRUE)CoreColor2 =labels2colors(Clusters2$cores)AutoColor =NULL;for (deepSplit in0:3){ AutoColor =cbind(AutoColor, labels2colors(cutreeDynamic(dendro = geneTree,minClusterSize =30-3*deepSplit,cutHeight = ch, method ="hybrid", deepSplit = deepSplit,distM = dissTOM)))AutoLabels =paste("Hybrid 'auto': dS =", c(0:3))par(mfrow=c(2,1))par(cex =1.4);par(mar =c(0,8.5,2,0));plot(geneTree,labels=F,main="Hierarchical dendrogram and module colors", sub="", xlab="")par(mar =c(1,8.5,0,0));plotHclustColors(geneTree, cbind(colorh1, DynamicColor1, DynamicColor2, ClustColor1, CoreColor1, ClustColor2, CoreColor2, AutoColor), c("Static", "Dynamic Tree (dS)", "Dynamic Tree (No dS)", "Dynamic Hybrid 1", "Cores in DHyb 1", "Dynamic Hybrid 2", "Cores in DHyb 2", AutoLabels), main ="")dev.off()}png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_hybride_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width =3840*2 , height =3840 ,res=600)# Same thing just with the hybrid function but with different settings# hybrid cutmColorh =NULLlabels =NULLpar(mfrow=c(2,1))par(cex =1.4);par(mar =c(0,5.5,2,0));plot(geneTree,labels=F,main="Hierarchical dendrogram and module colors with ", sub="", xlab="")par(mar =c(1,5.5,0,0));for (ds in0:3) for (mc inc(16,32,64)) { tree =cutreeHybrid(dendro = geneTree, pamStage=FALSE,minClusterSize = mc,#cutHeight = ch,deepSplit = ds,distM = dissTOM) mColorh=cbind(mColorh,labels2colors(tree$labels)) labels =cbind(labels,paste("dS =",ds,"; minSize =",mc))}#plotHclustColors(geneTree, staticColors, "tric", cex.rowLabels = 0.6)plotHclustColors(geneTree,mColorh,labels, cex.rowLabels =0.6)dev.off()deepSplit =2# The higher the parameters, the more coherent the modules created.mc=minModuleSizepng(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width =3840*2 , height =3840 ,res=600)tree =cutreeHybrid(dendro = geneTree,distM = dissTOM,#cutHeight = ch,deepSplit = deepSplit,minClusterSize = minModuleSize,pamStage=FALSE)dynamicMods <- tree$labels# Convert numeric lables into colorstable(dynamicMods)dynamicColors =labels2colors(dynamicMods)table(dynamicColors)write.csv(x=table(dynamicColors),file = here::here(paste0("data/multi_omic/WGCNA/module_detection_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))# Plot the dendrogram and colors underneathpar(mfrow =c(2, 1))plotDendroAndColors(geneTree, dynamicColors, "Tree Cut Hybrid",dendroLabels =FALSE, hang =0.03,addGuide =TRUE, guideHang =0.05,main ="Gene dendrogram and module colors selected")dev.off()cat("Merging modules \n")# Sometimes, especially when dealing with smaller clusters (for instance when the minClusterSize parameter is not set properly), merging modules with similar expression profiles might be usefull.# To do so, we first compute each module eigenvectors. A module eigenvector (or eigengene if you are dealing with gene data) represents the mean expression profile of the whole corresponding module. We can use the correlation of modules eigengenes to define a module-similarity metric that we can use to merge close modules.# Merging similar modules# Module-specific genesMEList = WGCNA::moduleEigengenes(get(paste0(input_name)), colors = dynamicColors)MEs = MEList$eigengenes # sample expression profile for each prob gene# Now let’s perform a hierarchical clustering using the correlation between module eigenvectors. Remember that for using hclust we need distances! Compute those distances and transform then with as.dist before running hclust.MEDiss <-1-cor(MEs)# merge certain modules# Classification of module genesMETree =hclust(as.dist(MEDiss), method ="average")# WGCNA’s authors recommend to merge cluster with a correlation greater than 0.75, hence a distance lower than 0.25. Let’s look at the tree with plotand abline, just as we did previously in this tutorial.MEDissThres =0.15# Disimilarite max thresholdpng(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width =3840*2 , height =3840 ,res=600)# This classification takes into account the type of networkplot(METree, main =paste("Clustering of module eigengenes pre-fusion (",networkType_i, " correlation)", sep =""),xlab ="", sub ="") # disimilarite tree# Similarity threshold for module mergingabline(h=MEDissThres, col ="red")dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select_verif_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width =3840*2 , height =3840 ,res=600)plotEigengeneNetworks(MEs, "",marDendro =c(0,4,1,2), marHeatmap =c(3,4,1,2),cex.lab =0.8, xLabelsAngle =90)dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select2_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width =3840*2 , height =3840 ,res=600)merge = WGCNA::mergeCloseModules(get(paste0(input_name)), dynamicColors, cutHeight = MEDissThres, verbose =3)mergedColors = merge$colors;mergedMEs = merge$newMEs;plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),c("Dynamic Tree Cut", "Merged dynamic"),dendroLabels =FALSE, hang =0.03,addGuide =TRUE, guideHang =0.05,main =paste("Merged clusters\nnetwork type = ", networkType_i,"\nmerging distance = ", MEDissThres, sep =""))dev.off()print("rename module")# Rename to moduleColorsmoduleColors = mergedColors# Construct numerical labels corresponding to the colorscolorOrder =c("grey", standardColors(50));moduleLabels =match(moduleColors, colorOrder)-1;MEs = mergedMEs;# Let’s save those clustering results! Choose your favorite export format (Rdata, csv, tsv, …) and save the results. Don’t forget to make sure you have all the information you need (in this case, we need the node names and their cluster labels).moduleColors = mergedColors # Or mergedColors, if some modules have been mergednames(moduleColors) =colnames(get(paste0(input_name)))head(moduleColors)unique(moduleColors)write.csv(x=moduleColors,file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))test =as.data.frame(moduleColors) %>%rownames_to_column("gene_id") %>%filter(gene_id =="PsCam042688")# Save module colors and labels for use in subsequent partssave(MEs, file = here::here(paste0("data/multi_omic/WGCNA/","modules_mrna_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".rda")))print("tomplot")# additional graphs# Transform dissTOM with a power to make moderately strong connections more visible in the heatmappng(file = here::here(paste0("report/multi_omics/plot/WGCNA/","plotTOM_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width =3840 , height =3840 ,res=600)plotTOM = TOM^0.6;# Set diagonal to NA for a nicer plotdiag(plotTOM) =NA#sizeGrWindow(9,9)TOMplot(plotTOM, geneTree, moduleColors, main ="Network heatmap plot, all genes")dev.off()# export importante valuesave(mergedMEs, file = here::here("data/multi_omic/WGCNA/mergedMEs.RData"))
16.4.3 Power selection
Correlation networks are by nature complete, which means every feature is associated to every other with a weight given by the correlation coefficient. To extract only meaningful associations, the correlation network must be thresholded.
The most straight-forward way of thresholding a correlation dataset is to perform a hard-thresholding, i.e. two features are considered correlated if their absolute correlation coefficient exceeds a user-defined threshold.
WGCNA propose another strategy, called soft-thresholding. It consists of raising the correlation coefficients to a certain power. This process transforms low correlation coefficients to values close to 0, while high correlation coefficients are transformed to values close to 1.
Furthermore, WGCNA presents a function called pickSoftThreshold that aims to automatically determine the ideal power for your dataset. The optimal power is defined as the one that ensures your network exhibits a scale-free property.
"Optimal soft-power = 10"
16.5 Show profile of each module
Code
load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))input_name ="datExpr"networkType_i ="signed"minModuleSize =30deepSplit =2moduleColors <-read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <-c("gene_id", "moduleColors")# specificity of the different modulesummarise_color <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(n_genes =n()) %>%arrange(desc(n_genes))# a. passer datExpr au format longexpr_long <-as.data.frame(datExpr) %>%rownames_to_column(var ="sample") %>%# garder l'ID échantillonpivot_longer(-sample,names_to ="gene_id",values_to ="expression")# b. ajouter la couleur de module pour chaque gèneexpr_long <- expr_long %>%left_join(., moduleColors, by ="gene_id") %>%left_join(., metadata %>%rownames_to_column("sample"), by ="sample") %>%mutate(condition =paste0(sulfur_condition, "_", genotype))## ─────────────────────────────## 2. Moyenne d’expression par échantillon & par module## ─────────────────────────────expr_module_sample <- expr_long %>% dplyr::group_by(moduleColors, sample, condition) %>% dplyr::summarise(mean_expr =mean(expression, na.rm =TRUE), # moyenne d’expressionn_genes = dplyr::n(), # nombre total de lignes (= gènes)# si tu préfères compter les gènes uniques :# n_genes = dplyr::n_distinct(gene_id),.groups ="drop" ) %>% dplyr::mutate(moduleColors_n =paste0(str_to_title(moduleColors), " (", n_genes, ")"), moduleColors_n = forcats::fct_reorder(moduleColors_n, n_genes, .fun = max, .desc =TRUE) ) %>%mutate(condition=fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")) strip_df <- expr_module_sample %>%distinct(moduleColors) %>%# une ligne par modulearrange(moduleColors) %>%mutate(strip_fill =as.character(moduleColors) )## On récupère les vecteurs (ordre = niveaux du facteur)strip_fill_vec <- strip_df$strip_fillnames(strip_fill_vec) <- strip_df$moduleColorsstrip_text_vec <- strip_df %>%mutate(couleur_texte_perso =sapply(strip_fill, evaluate_contrast)) %>%pull(couleur_texte_perso)names(strip_text_vec) <- strip_df$moduleColorsp1 <- expr_module_sample %>%mutate(genotype =sub("^[^_]+_", "", condition), sulfur_condition =sub("_.*$", "", condition),genotype =fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"),sulfur_condition =fct_relevel(sulfur_condition, "SS", "SD") ) %>%ggplot(.,aes(x = condition,y = mean_expr,fill = genotype, col = genotype)) +geom_boxplot(outlier.shape =NA, alpha = .4) +geom_jitter(width = .25, size = .6, alpha = .4) +# points individuellabs(x ="Module",y =expression("Mean expression level (log"["2"]*")"),colour ="Genotype", fill ="Genotype") +scale_color_manual(values = mutant_palette)+scale_fill_manual(values = mutant_palette)+# --> Ici tu peux ajouter new_scale_fill() si tu superposes d'autres couches ggh4x::facet_wrap2(~ moduleColors_n,scales ="free_y",nrow =7, # ajuste le layoutstrip = ggh4x::strip_themed(background_x = ggh4x::elem_list_rect(fill = strip_fill_vec),text_x = ggh4x::elem_list_text(colour = strip_text_vec) ) ) +theme_bw() +theme(panel.spacing =unit(0.5, "lines"),axis.ticks.y =element_blank(),axis.text.x =element_text(angle =90, vjust =0.5, hjust =1),legend.text =element_text(size =12),legend.title =element_text(size =12),legend.key.size =unit(10, "pt"),plot.caption =element_text(size =8),panel.grid.major =element_blank(),panel.grid.minor =element_blank() ) ;p1fig_export(path ="report/multi_omics/plot/WGCNA/average_module", plot_x = p1, height_i =9, width_i =12, res =600)expr_module_sample$condition# i think i need to use patchwork and change a litte bit my function stats. v_moduleColors =levels(as.factor(expr_module_sample$moduleColors_n))# add color to the table expr_module_sampleexpr_module_sample = expr_module_sample %>%mutate(couleur_texte_perso =sapply(moduleColors, evaluate_contrast)) # i = 4# moduleColors_i <- v_moduleColors[i]# # test <- stat_analyse(# data=expr_module_sample %>% # as.data.frame() %>% # filter(moduleColors_n == moduleColors_i) %>% # mutate(condition = as.factor(condition)),# # column_value = "mean_expr",# category_variables = "condition",# grp_var = "",# show_plot = T,# outlier_show = F, # label_outlier = "sample",# biologist_stats = T,# Ylab_i = "Chlorophyll content for SD and consistent pod \n number condition in N4_N5 at 40 DAP",# control_conditions = "",# strip_normale = T,# hex_pallet = as.character(c(mutant_palette, mutant_palette)), # strip = "moduleColors_n", #colone for texte# strip_fill_vec = "moduleColors", # color of the square # strip_text_vec = "couleur_texte_perso"# )# # test[["plot"]]######################### plots_list <-lapply(seq_along(v_moduleColors), function(i) { moduleColors_i <- v_moduleColors[i] p <-stat_analyse(data = expr_module_sample %>%as.data.frame() %>%filter(moduleColors_n == moduleColors_i) %>%mutate(condition =as.factor(condition)) %>%filter(condition %in%c("SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")), # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!column_value ="mean_expr",category_variables ="condition",grp_var ="",show_plot =TRUE,outlier_show =FALSE,label_outlier ="sample",biologist_stats =TRUE,Ylab_i =paste0("" ),control_conditions ="",strip_normale =TRUE,hex_pallet =as.character(rep(mutant_palette, 2)),strip ="moduleColors_n",strip_fill_vec ="moduleColors",strip_text_vec ="couleur_texte_perso" )[["plot"]]+theme_bw()+theme(axis.text.x =element_text(angle =90, hjust =1))## 1-a) On retire les titres d’axe des sous-graphiques p <- p +labs(x =NULL, y =NULL)## 1-b) Option : on masque aussi les graduations redondantes# (ici : on garde y seulement sur la 1re colonne# et x seulement sur la dernière ligne) ncol_layout <-5 nrow_layout <-ceiling(length(v_moduleColors) / ncol_layout) col_i <- ((i -1) %% ncol_layout) +1 row_i <-ceiling(i / ncol_layout)# if (col_i != 1) {# p# }if (row_i != nrow_layout) { p <- p +theme(axis.text.x =element_blank()) } p +labs(fill ="Genotype", colour ="Genotype")})# ------------------------------------------------------------------# 2) Mosaïque principale -------------------------------------------# ------------------------------------------------------------------panel <-wrap_plots(plots_list, ncol =5, guides ="collect") &theme(legend.position ="right", panel.grid.major =element_blank(),panel.grid.minor =element_blank() )## ------------## 2. Labels## ------------y_lab <-ggplot() +labs(y ="Mean expression") +theme_void() +theme(axis.title.y =element_text(angle =90, vjust =0.5,hjust =0.5))x_lab <-ggplot() +labs(x ="Condition") +theme_void() +theme(axis.title.x =element_text(vjust =-0.5, hjust =0.5))## ------------## 3. Assemblage 2 × 2## ------------final_plot_1 <- (y_lab | panel)+plot_layout(widths =c(0.02, 0.98), # 5 % pour la colonne Y, 95 % pour le panelguides ="collect" )final_plot_2 <- (plot_spacer() | x_lab)+plot_layout(widths =c(0.02, 0.98), # 5 % pour la colonne Y, 95 % pour le panelguides ="collect" )final_plot <- (final_plot_1/final_plot_2)+plot_layout(heights =c(0.99, 0.01), # 7 % pour le label Xguides ="collect" )final_plot#fig_export(path = "report/multi_omics/plot/WGCNA/average_module_stats", plot_x = final_plot, height_i = 9, width_i = 18, res = 600)fig_export(path ="report/multi_omics/plot/WGCNA/average_module_stats_SD", plot_x = final_plot, height_i =9, width_i =10, res =600)
16.6 Correlation with other data sets
Code
load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))load(file = here::here("data/multi_omic/WGCNA/mergedMEs.RData"))df_plant_info_2014_filter <-read_excel(here::here("data/plant_info_XP2014.xlsx")) %>%drop_na(code_run_transcripto) %>% dplyr::select(plant_num, code_run_transcripto, code_metabo_leaf_repro,code_metabo_leaf_vege, sulfate) %>% dplyr::rename(code_run_anion = sulfate)# load metabolomic resultsdf_metabo_global <-read_csv(here::here("data/metabolomic/output/df_metabolite_global_IBMP.csv"), show_col_types =FALSE) %>%mutate(sulfur_condition=fct_relevel(sulfur_condition, "SS", "SD"),genotype =fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"),type_genotype =fct_relevel(type_genotype, "WT", "Mut") ) %>%filter(plant_num %in% df_plant_info_2014_filter$plant_num, compartment =="vegetative leaf") %>%#() "reproductive leaf" "vegetative leaf" ) dplyr::select(variable, plant_num, value) %>%pivot_wider(names_from = variable, values_from = value)vect_metabo =colnames(df_metabo_global)[-1]################################## for more simplicity change here vect_metabo = vect_metabo[1:100]plant_measured <- df_metabo_global %>%drop_na(glucose) %>%pull(plant_num)metadata_transcripto_metabo_filter = metadata %>%rownames_to_column("sample_id") %>%mutate(sample_num =str_match(sample_id, "(Mut|WT)_[A-Z]+_(\\d+\\.\\d+)")[,3], Rep =paste0("Rep", str_match(sample_id, "_Rep(\\d+)")[,2])) %>%mutate(code_run_transcripto =gsub("\\.", "_", as.character(sample_num))) %>%left_join(., df_plant_info_2014_filter, by ="code_run_transcripto") %>%left_join(., df_metabo_global, by ="plant_num") %>%mutate(sample_name =paste(sep ="_", Mutant_type, sulfur_condition, sample_num, Rep, chip_color, chip_num)) %>%column_to_rownames("sample_name") %>%filter(plant_num %in% plant_measured) %>% dplyr::select(all_of(vect_metabo)) %>% dplyr::select(where(~any(. !=0)))# Add anionmetadata_transcripto_metabo_filterlist_sample <-rownames(metadata_transcripto_metabo_filter)mergedMEs_filter <- mergedMEs %>%rownames_to_column("sample_id") %>%filter(sample_id %in% list_sample) %>%column_to_rownames("sample_id")# filter one data i have for boath datasetstopifnot(identical(rownames(mergedMEs_filter), rownames(metadata_transcripto_metabo_filter)))moduleTraitCor <-cor(mergedMEs_filter, metadata_transcripto_metabo_filter, use ="pairwise.complete.obs", method ="pearson")moduleTraitPvalue <-corPvalueStudent(moduleTraitCor,nSamples =nrow(mergedMEs)) # WGCNA fonctiontextMatrix <-paste0(signif(moduleTraitCor, 2),"\n(",signif(moduleTraitPvalue, 1),")")dim(textMatrix) <-dim(moduleTraitCor)colnames(textMatrix) <-colnames(moduleTraitCor)rownames(textMatrix) <-rownames(moduleTraitCor)# 4. Masquer les r non significatifs (p > 0.05)for (r inrownames(moduleTraitCor)) {for (c incolnames(moduleTraitCor)) {if (moduleTraitPvalue[r, c] >0.05) { moduleTraitCor[r, c] <-0 textMatrix[r, c] <-"" } }}# 1. Récupérer noms + couleurs des modules -------------------------------module_names <-rownames(moduleTraitCor) # ex. "MEturquoise"module_colors <-sub("^ME", "", module_names) # ex. "turquoise"# 2. Data-frame d’annotation (une ligne = un module) ----------------------annotation_row <-data.frame(Module =factor(module_colors,levels =unique(module_colors)))rownames(annotation_row) <- module_names# 3. Palette : liste nommée (niveau du facteur -> couleur) ---------------annotation_colors <-list(Module =setNames(unique(module_colors), # valeurs uniques seulementunique(module_colors)))# Si tu veux des codes hex personnalisés, remplace la ligne ci-dessus par :# annotation_colors <- list(# Module = c(turquoise = "#00CED1", brown = "#8B4513", ...)# )# 4. Heatmap --------------------------------------------------------------my_col <-colorRampPalette(c("darkgreen", "white", "darkmagenta"))(50)png(here::here("report/multi_omics/plot/WGCNA/module_metabo_heatmap_clustered_test.png"),width =6000, height =4000, res =300)pheatmap(mat = moduleTraitCor,cluster_rows =TRUE,cluster_cols =TRUE,clustering_distance_rows ="euclidean",clustering_distance_cols ="euclidean",clustering_method ="complete",color = my_col,display_numbers = textMatrix,number_color ="black",fontsize_number =5,legend_breaks =seq(-1, 1, by =0.5),border_color =NA,main ="Module–metabolite relationships (clustered)",annotation_row = annotation_row, # ← la barre de couleur à gaucheannotation_colors = annotation_colors, # ← palette correspondantelegend =FALSE, # ← supprime la palette continueannotation_legend =FALSE# ← supprime la légende des modules)dev.off()
load(file = here::here("data/multi_omic/LMM/output/result_LMM_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/result_LMM_all_without_anion.RData"))df_result %>%filter(pval<0.05) %>%filter(V1 =="ANION_sulfate"|V2=="ANION_sulfate")standard_correlation %>%filter(fdr<0.05) %>%filter(V1 =="ANION_sulfate"|V2=="ANION_sulfate")ggplot(df_global_SD, aes(x = ANION_sulfate, y = MODULE_violet, col = genotype))+geom_point()ggplot(df_compile, aes(x = ANION_sulfate, y = METABO_glycolic_acid, col = genotype))+geom_point()df_result %>%filter(pval<0.05) %>%filter(V1 =="ANION_nitrate"|V2=="ANION_nitrate")standard_correlation %>%filter(fdr<0.05) %>%filter(V1 =="ANION_nitrate"|V2=="ANION_nitrate")
16.8 Selection based on pvalue and R²
R2m (Marginal R²): This represents the variance explained by the fixed effects in your model. It is similar in interpretation to the R² value in traditional linear regression models. R2c (Conditional R²): This represents the variance explained by both the fixed and random effects in your model. It considers both fixed and random effects, providing a more comprehensive measure of goodness-of-fit for mixed-effects models. Both Rm2 and R2c can be useful in understanding the proportion of variance explained by the predictors in your model.
Code
clean_session() # Remove data frames and matrices and Reset plots# data importationload(file = here::here("data/multi_omic/LMM/output/result_LMM_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/result_LMM_all_without_anion.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD_without_METABO.RData"))#load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_no_prefiltre.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_5.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_7.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_9.RData"))# parameterr2m_i =0.5r2c_i =0.85fdr_lim =0.05df_result = df_result %>%mutate(fdr =p.adjust(pval, method ="fdr"),log10fdr =-log10(fdr)) %>%mutate(color=ifelse(r2c*1.5- r2m*1.5>0.5, "red",ifelse(r2m<r2m_i, ifelse(r2c<r2c_i,"#DD6342","orange"),ifelse(fdr<fdr_lim, "green","forestgreen") ) ) )# for network related to mixed effect# test more restrictive for fixed effect# r2m_i = 1# r2c_i = 0.95# fdr_lim = 0.05# # df_result = df_result %>% # mutate(fdr = p.adjust(pval, method = "fdr"),# log10fdr = -log10(fdr)) %>% # # mutate(color=ifelse(r2c > r2c_i, "green","red"# )# )# if i whant to use normale correlationstandard_correlation %>%filter(V1 =="ANION_sulfate"|V2=="ANION_sulfate") %>%arrange (fdr)# Verifdf_result %>%filter(pval<0.05) %>%filter(V1 =="ANION_sulfate"|V2=="ANION_sulfate")# Verif# Visualization ## Use of randomized data to avoid displaying all pointsrandom_result <- df_result %>%sample_n(1000)random_result%>%ggplot(aes(x = r2m, y =-log10(pval), size = slope, col =-log10(pval))) +geom_point(alpha =0.7) +scale_color_viridis_c() +labs(title ="Dispersion diagramme", x ="R2m", y ="-log10(fdr)", color ="-log10(fdr)") +theme_minimal()random_result %>%ggplot(aes(x = r2c, y =-log10(pval), size = slope, col =-log10(pval))) +geom_point(alpha =0.7) +scale_color_viridis_c() +labs(title ="Dispersion diagramme", x ="R2c", y ="-log10(fdr)", color ="-log10(fdr)") +theme_minimal()random_result %>%ggplot(aes(x = r2m, y = r2c, col = color, size = slope)) +geom_point(aes(shape = slope <0), na.rm =TRUE) +scale_shape_manual(values =c(`TRUE`=17, `FALSE`=16)) +scale_color_manual(values =c("red"="red", "#DD6342"="#DD6342","orange"="orange", "green"="green", "forestgreen"="forestgreen")) +# Utilisation de couleurs fixesgeom_vline(xintercept = r2m_i, linetype ="dashed", color ="red", linewidth =1) +geom_hline(yintercept = r2c_i, linetype ="dashed", color ="red", linewidth =1)# filter the result to have variabledf_result_select <- df_result %>%filter(color=="green") %>%arrange(pval)cat_col(paste("Number of variable:",length(unique(c(df_result_select$V1,df_result_select$V2))), "\n"), "green")round(dim(df_result_select)[1]*100/dim(df_result)[1], 2)# for standard correlation ? #df_result_select_corr <- standard_correlation %>% df_result_select_corr <- standard_correlation %>%mutate(color =ifelse(fdr<0.05, "green", "red")) %>%filter(color=="green") %>%arrange(fdr)cat_col(paste("Number of variable:",length(unique(c(df_result_select_corr$V1,df_result_select_corr$V2))), "\n"), "green")# exportsave(df_result_select, file = here::here("data/multi_omic/LMM/output/result_LMM_filter.RData"))save(df_result_select_corr, file = here::here("data/multi_omic/LMM/output/result_corr_filter.RData"))
16.9 Network plot with Igraph
Code
clean_session() # Remove data frames and matrices and Reset plotsload(file = here::here("data/multi_omic/LMM/output/result_LMM_filter.RData"))load(file = here::here("data/multi_omic/LMM/output/result_corr_filter.RData"))############ warning i change that #################df_result_select = df_result %>%mutate(fdr =p.adjust(pval, method ="fdr"),log10fdr =-log10(fdr)) %>%mutate(color=ifelse(fdr<fdr_lim, "green", "red") )df_result_select = df_result_select_corr# pb <- progress::progress_bar$new(# format = " Calcul [:bar] :percent | Étape :current/:total | Temps écoulé :elapsed | Temps restant :eta",# total = nrow(df_result_select),# clear = FALSE,# width = 100# )# # # Fonction pour calculer la corrélation avec mise à jour de la barre de progression# calculate_corr <- function(var1, var2) {# pb$tick()# # # Vérifier si les variables existent dans df_compile# if (all(c(var1, var2) %in% colnames(df_compile))) {# # Calculer la corrélation# cor_value <- cor(df_compile[[var1]], df_compile[[var2]], use = "pairwise.complete.obs")# return(cor_value)# } else {# return(NA) # Retourne NA si une variable manque# }# }# Utiliser mapply pour vectoriser le calcul#df_result_select$corr <- mapply(calculate_corr, df_result_select$V1, df_result_select$V2)# filter the data previously generate# # words_to_remove <- c("sEUpE", "DW","Ba","As","Cd","Ti","Be","Fe","V","Cr","Tl")# words_to_remove <- c("sEUpE")# words_to_remove <- c("sEUpE_C")# words_to_remove <- c("sEUpE_U")# words_to_remove <- c("sEUpE_U", "IONO", "ECO", "RSA")#words_to_remove <- c("[", "]")df_select=df_result_select %>%filter(color=="green") %>%filter(between_IC05=="no") #%>% #filter(!str_detect(V1, paste(words_to_remove, collapse = "|"))) %>% #filter(!str_detect(V2, paste(words_to_remove, collapse = "|")))df_select <- df_result_select %>%filter(color=="green") %>%filter(!grepl("putative", V1, ignore.case =TRUE)) %>%filter(!grepl("putative", V2, ignore.case =TRUE))df_select <- df_select %>%filter(color=="green") %>%filter(!grepl("putative", V1, ignore.case =TRUE)) %>%filter(!grepl("putative", V2, ignore.case =TRUE))#or with correlation# df_select=df_result_select_corr %>% # #filter(select=="yes") %>% # filter(color=="green")#hist(df_select$slope_modif)# creation of the node nodes<-data.frame(variable =unique(c(df_select$V1,df_select$V2))) %>%mutate(id =paste0("s", 1:length(unique(c(df_select$V1,df_select$V2))))) %>%relocate(id, .before=variable) %>%#mutate(variable_type = sample(1:3, nrow(.), replace = TRUE)) %>% #generate for aleatoire color# mutate(organe= case_when(# grepl("leaf", variable) ~ "leaf",# grepl("stem", variable) ~ "stem",# grepl("root", variable) ~ "root",# grepl("surface", variable) ~ "root",# grepl("volume", variable) ~ "root",# grepl("density", variable) ~ "root",# grepl("sEUpE", variable) ~ "sEUpE",# grepl("TR", variable) ~ "leaf",# TRUE ~ "autre"# )) %>% mutate(organe_num= case_when(# grepl("leaf", variable) ~ 1,# grepl("stem", variable) ~ 2,# grepl("root", variable) ~ 3,# grepl("surface", variable) ~ 3,# grepl("volume", variable) ~ 3,# grepl("density", variable) ~ 3,# grepl("sEUpE", variable) ~ 5,# grepl("TR", variable) ~ 1,# TRUE ~ 4# )) %>% mutate(variable_type=case_when(grepl("METABO", variable) ~1,grepl("ANION", variable) ~4,grepl("MODULE", variable) ~2,#grepl("METABO.T", variable) ~ 5,TRUE~3 )) %>%mutate(# mutate(variable_cleaned = ifelse(variable_type==1,gsub("[0-9_]", "", variable),variable)) %>%# mutate(variable_cleaned = ifelse(variable_type==4,gsub("[0-9_]", "", variable_cleaned),variable_cleaned)) %>%# mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>%# mutate(variable_cleaned = gsub("(stem|leaf|root|concentration|sEUpE)", "", variable_cleaned)) %>% # dplyr::mutate(# variable_cleaned = gsub("weight", "DW", variable_cleaned), # variable_cleaned = gsub("sumtotalevapotranspiration", "ETtot", variable_cleaned),# variable_cleaned = gsub("Hydricpotential", " \u03C8 ", variable_cleaned),# variable_cleaned = gsub("Photo", "An", variable_cleaned),# variable_cleaned = gsub("volume", "Volume", variable_cleaned),# variable_cleaned = gsub("surface", "Surface", variable_cleaned),# variable_cleaned = gsub("shootratio", "S.R", variable_cleaned),# variable_cleaned = gsub("sumbiomass", "TotDW", variable_cleaned),# variable_cleaned = gsub("density", "Density", variable_cleaned),# variable_cleaned = gsub("area", "Area", variable_cleaned),# variable_cleaned = gsub("TRmmolms", "TR", variable_cleaned),# variable_cleaned = gsub("ConvexHull", "Hull", variable_cleaned),variable_cleaned =str_remove_all(variable, "METABO_|ANION_|MODULE_|\\[|\\]"),nb_char =nchar(variable_cleaned))# linkslinks<-df_select %>% dplyr::select(V1,V2, corr) %>% dplyr::mutate(weight =abs(corr)) %>% dplyr::mutate(weight_repulstion =1-corr) %>%inner_join(nodes %>% dplyr::rename(V1=variable) %>% dplyr::rename(X1=id),.,by="V1") %>%inner_join(nodes %>% dplyr::rename(V2=variable) %>% dplyr::rename(X2=id),.,by="V2") %>%relocate(X1, .before=X2) %>%relocate(V1, .before=V2) # Converting the data to an igraph object:net <-graph.data.frame(links, nodes, directed=T) save(nodes, links, net,file = here::here("data/multi_omic/LMM/output/result_LMM_graph.RData"))# Generate colors base on trend type:colrs2 <-c("#68a500", "#b8af83", "#ce7f50","black","#345995")V(net)$color <- colrs2[V(net)$variable_type]V(net)$color =ifelse(V(net)$variable_type==2, V(net)$variable_cleaned, V(net)$color )# Compute node degree (#links) and use it to set node size:deg <-degree(net, mode="all")# if(fc_or_pval=="pval"){# V(net)$size <- V(net)$log10_pval_test_t*6 #ixi modifie la taille des points# }else if(fc_or_pval=="fc"){# V(net)$size <- V(net)$abslog2fc*10 #ixi modifie la taille des points# }myrhombus <-function(coords, v =NULL, params) { vertex.color <-params("vertex", "color")if (length(vertex.color) !=1&&!is.null(v)) { vertex.color <- vertex.color[v] } vertex.size <-1/200*params("vertex", "size")if (length(vertex.size) !=1&&!is.null(v)) { vertex.size <- vertex.size[v] }symbols(x = coords[, 1], y = coords[, 2], bg = vertex.color,stars =cbind(1.2*vertex.size, vertex.size, 1.2*vertex.size, vertex.size),add =TRUE, inches =FALSE)}# clips as a circleadd_shape("rhombus", clip =shapes("circle")$clip,plot = myrhombus)## Function for plotting an elliptical nodemyellipse <-function(coords, v=NULL, params) { vertex.color <-params("vertex", "color")if (length(vertex.color) !=1&&!is.null(v)) { vertex.color <- vertex.color[v] } vertex.size <-1/70*params("vertex", "size") # largeur des elipseif (length(vertex.size) !=1&&!is.null(v)) { vertex.size <- vertex.size[v] } plotrix::draw.ellipse(x=coords[,1], y=coords[,2],a = vertex.size, b=0.024, col=vertex.color) #b = hauteur}## Register the shape with igraphadd_shape("ellipse", clip=shapes("circle")$clip,plot=myellipse)# Fonction pour dessiner un hexagonemyhexagon <-function(coords, v =NULL, params) { vertex.color <-params("vertex", "color")if(length(vertex.color) !=1&&!is.null(v)){ vertex.color <- vertex.color[v] } vertex.size <-1/200*params("vertex", "size")if(length(vertex.size) !=1&&!is.null(v)){ vertex.size <- vertex.size[v] }for(i inseq_len(nrow(coords))){ x <- coords[i, 1] y <- coords[i, 2] size <- vertex.size[i] theta <-seq(0, 2*pi, length.out =7)[-7] # Angles pour les sommets de l'hexagone xs <- x + size *cos(theta) ys <- y + size *sin(theta)polygon(x = xs, y = ys, col = vertex.color[i], border ="black") }}hexagon_clip <-function(coords, el, params) { vertex.size <-1/200*params("vertex", "size") angles <-seq(0, 2* pi, length.out =7)[-7] clip_coords <-sapply(1:nrow(coords), function(i) { size <- vertex.size[i] angle <-atan2(el[2, 2] - el[1, 2], el[2, 1] - el[1, 1]) # Edge direction idx <-which.min(abs(angles - angle)) # Find closest angle in hexagon x <- coords[i, 1] + size *cos(angles[idx]) y <- coords[i, 2] + size *sin(angles[idx])c(x, y) })return(t(clip_coords))}# Enregistrer la nouvelle forme avec igraphadd_shape("hexagon", clip =shapes("circle")$clip, plot = myhexagon)# V(net)$shape<-ifelse(V(net)$variable_type==2,"rectangle",ifelse(V(net)$variable_type==4,"ellipse",ifelse(V(net)$variable_type==3,"rhombus", ifelse(V(net)$variable_type==5, "hexagon", "hexagon"))))V(net)$size=1+deg*0.1#avant 10V(net)$size =ifelse(V(net)$shape =="rectangle", V(net)$nb_char *3.5+5,ifelse(V(net)$shape =="ellipse",ifelse(V(net)$nb_char >3, 10+V(net)$nb_char*0.16,1.7),ifelse(V(net)$shape =="hexagone", V(net)$nb_char *0.3+1, 7) # Taille ajustée pour les hexagones ))V(net)$size2=3.5# hauteurV(net)$label.cex <-0.8# Taille plus petite des étiquettes# The labels are currently node IDs.# Setting them to NA will render no labels:colrs1 <-c("black", "white", "black","white")V(net)$label.color <-colrs1[ifelse(V(net)$variable_type=="METABO",2,V(net)$variable_type)] #color of the labelV(net)$label.color <-sapply(V(net)$color, evaluate_contrast) #color of the labelV(net)$label.dist <-0#.7 #distance from the center of the vertexV(net)$label.family="sans"#type de police#V(net)$label.font=V(net)$variable_type # 2 for bolt 1 is plaint texte 3 for italicV(net)$label.font=ifelse(V(net)$variable_type==5,1,V(net)$variable_type) # 2 for bolt 1 is plaint texte 3 for italic#V(net)$label <- NA# Set edge width based on weight:#E(net)$width <- -log10(1-abs(E(net)$slope_modif))*1.5# avt /6E(net)$width <- (abs(E(net)$corr)-min(abs(E(net)$corr))+0.01)*5# avt /6#change arrow size and edge color:E(net)$arrow.size <- .2E(net)$edge.color <-"gray80"# Let's color the edges of the graph based on their source node color.# We'll get the starting node for each edge with "get.edges"edge.start <-get.edges(net, 1:ecount(net))[,1]edge.col <-V(net)$color[edge.start]vec_color=ifelse(E(net)$corr>0,"red","blue")edge.col <-vec_color #mettre red if is upper green if is lower 0l=layout_in_circlel=layout_with_kk # Kamada Kawai #Like Fruchterman Reingold, it attempts to minimize the energy in a spring system.l=layout_on_gridl=layout_with_fr #Fruchterman-Reingoldl <-layout_with_fr(net, niter =10000, grid ="nogrid")bound <-1# demi-largeur de la boîten <-vcount(net)for (sign inc("sign", "unsign")){if (sign =="unsign"){ scaling_factor_attraction_force <-3# A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly. l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$weight-min(igraph::E(net)$weight)+0.001), niter =500, grid ="nogrid") # without sign effect }elseif (sign =="sign"){ scaling_factor_attraction_force <-10# A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly. l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$corr+abs(min(igraph::E(net)$corr))+0.001), niter =500, grid ="nogrid")#, start.temp = bound, minx = rep(-bound, n), maxx = rep(bound, n),#miny = rep(-bound, n), maxy = rep(bound, n)) # withe sign effect }# exportsvg(width=8, height=8,filename = here::here(paste0("report/multi_omics/plot/LMM/network_SD_without_METABO_",sign,".svg")))set.seed(1)plot(net, layout = l, edge.lty =1,edge.arrow.size =0,vertex.label =V(net)$variable_cleaned,vertex.label.cex =0.5,vertex.size =V(net)$size*0.5, # Assurez-vous que la taille est passée iciedge.color = edge.col,edge.curved = .15)title(main ="Network plot based on LMM results show correlation between variables")mtext(paste0("Positive correlation: ",length(df_select %>%filter(corr>0) %>%pull(corr))),side =1, line =2, cex =0.8,col ="red")mtext(paste0("Negative correlation: ",length(df_select %>%filter(corr<0) %>%pull(corr))),side =1, line =1, cex =0.8,col ="blue")dev.off()}# save plot to analyse in cytoskapewrite_graph(net, file = here::here("data/multi_omics/output/B/network.graphml"), format ="graphml")## Degree Distribution Histogramdf_deg=tibble(nodes=names(deg),deg=as.vector(deg)) %>%full_join(nodes %>%as.data.frame() %>% dplyr::rename(nodes=id),.,by="nodes") %>%arrange(desc(deg)) %>%#mutate(variable_cleaned = ifelse(variable_type==1|variable_type==4,gsub("[0-9_]", "", variable),variable)) %>%# mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>%# mutate(variable_cleaned = gsub("(concentration)", "", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(stem)", "", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(root)", " ", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(leaf)", " ", variable_cleaned)) %>% mutate(organe=factor(organe, levels =c("leaf", "stem", "root","autre", "sEUpE"))) %>%#mutate(variable_cleaned = gsub("sEUpE", "sEUE ", variable_cleaned)) %>% mutate(variable_cleaned =paste0(variable_cleaned, " ", str_to_title(organe)))# mutate(variable_cleaned = gsub("(stem|leaf|root|concentration)", "", variable_cleaned))px=df_deg %>%filter(deg >4) %>%ggplot(aes(x = deg, y =reorder(variable_cleaned, -deg), fill = organe, label = deg)) +geom_col(position ='dodge') +geom_text(hjust =-0.4) +scale_fill_manual(values = colrs2, breaks =c("leaf","stem","root", "autre", "sEUpE")) +# scale_x_continuous(limits = c(0, 15),breaks = seq(0, 15, by = 1))+theme_minimal() +theme(axis.title.y =element_blank(),axis.text.x =element_blank(),axis.ticks.y =element_blank(),panel.grid =element_blank(),legend.position ="none")+labs(x="Degree", fill="Organe")ggsave(here::here("report/multi_omics/plot/LMM_metabo_iono/degree_all.svg"), px, height =15,width =10)############# module les plus correlele avec metabolitetop_corr_meta_module <- df_select %>%## 1. garder uniquement les lignes qui contiennent un METABO et un MODULEfilter( (str_starts(V1, "METABO_") &str_starts(V2, "MODULE_")) | (str_starts(V1, "MODULE_") &str_starts(V2, "METABO_")) ) %>%## 2. corrélation en valeur absoluemutate(abs_corr =abs(corr)) %>%## 3. trier du plus fort au plus faiblearrange(desc(abs_corr)) module_connection_counts <- df_select %>%## 1. ne garder que les couples METABO / MODULEfilter( (str_starts(V1, "METABO_") &str_starts(V2, "MODULE_")) | (str_starts(V1, "MODULE_") &str_starts(V2, "METABO_")) ) %>%## 2. isoler le nom du module et du métabolite dans deux nouvelles colonnesmutate(module =if_else(str_starts(V1, "MODULE_"), V1, V2),metabolite =if_else(str_starts(V1, "METABO_"), V1, V2) ) %>%## 3. compter le nombre de métabolites distincts par modulegroup_by(module) %>%summarise(n_metabolites =n_distinct(metabolite),.groups ="drop" ) %>%## 4. trier du plus connecté au moins connectéarrange(desc(n_metabolites))module_connection_counts
16.10 Analyse module
Code
input_name ="datExpr"networkType_i ="signed"minModuleSize =30deepSplit =2lfc_lim_i =0moduleColors <-read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <-c("PsCam", "moduleColors")# transformation into a listgene_list <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(genes =list(PsCam), .groups ="drop") %>%deframe() # transforme en liste nomméeload(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))# load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>%left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet ="Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by ="PsCam")# cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Down) ", comparison)), # purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Up) ", comparison))# ) %>% full_join(., df_info_gene, by = "PsCam") %>% # mutate(ID = Psat) %>% # drop_na(comparison, ID)col_order <-c("darkolivegreen", "orange", "yellowgreen", "purple","darkgreen", "midnightblue", "lightgreen", "blue","grey60", "steelblue")#for article, only for SD conditioncol_order <-c("saddlebrown", "black", "orange", "tan")col_order <-c("green", "steelblue", "orange")cluster_info <- purrr::imap_dfr(gene_list, ~tibble(PsCam = .x, cluster = .y)) %>%full_join(df_info_gene, by ="PsCam") %>%mutate(ID = Psat) %>%drop_na(ID) %>%filter(cluster %in% col_order) %>%# tu gardes les bonnes couleursmutate(# ① méthode dplyr de basecluster =factor(cluster, levels = col_order, ordered =TRUE)# ② méthode forcats, équivalente si tu préfères :# cluster = fct_relevel(cluster, !!!col_order) )#filter(cluster %in% c("tan", "saddlebrown"))source(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO termedf_GO =GO_on_different_group(functional_roles ="BP",group_info = cluster_info,group ="cluster",ID ="ID",top =10 )test = df_GO %>%filter(group =="tan")test = df_GO %>%filter(group =="saddlebrown")cyrcadian <-as.data.frame(gene_list[["orange"]]) ; colnames(cyrcadian) <-"PsCam"cyrcadian <- cyrcadian %>%left_join(., df_info_gene, by ="PsCam")# export for titouan#writexl::write_xlsx(cyrcadian, path = here::here("data/multi_omic/cyrcadian.xlsx"))# Merge this with your original data to fill missing combinationsCluster_selected_GO_filled <-df_GO %>%# on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>%mutate(cluster =factor(cluster,levels =unique(cluster))#, # comparison = str_trim(str_remove(cluster, "\\s*\\(.*?\\)")),# sign = str_extract(cluster, "(?<=\\().*?(?=\\))") ) #%>% #mutate(comparison = forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT"))#vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF")# Plot the enrichment by GOpx <-ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill =`|-log10(Pval)|`)) +geom_tile(color ="black") +scale_fill_gradient2(low ="white", high ="red", na.value ="white", limits =c(0, NA)) +theme_minimal() +theme(axis.text =element_text(size =8, colour ="black"),axis.text.x =element_text(angle =90, vjust =0.5, hjust =1),panel.grid.major =element_blank(),plot.background =element_rect(fill ="white", colour ="white")) +labs(fill ="-log10(FDR)",x="Biological process",y ="Cluster") +#scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+new_scale_fill() +scale_fill_manual(values = col_order,name ="Module" ) +geom_tile(aes(x = cluster,y =-1, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle thisfill = cluster,width =0.90,height =0.80 ),data = Cluster_selected_GO_filled,color ="black",alpha =1,inherit.aes =FALSE )# export #fig_export(path =paste0("report/multi_omics/plot/GO/darkolivegreen_orange_GO_SD_condition"), plot_x = px, height_i = 5.5, width_i = 4.5, res = 600)fig_export(path =paste0("report/multi_omics/plot/GO/some_module_GO_all_without_metabo"), plot_x = px, height_i =7, width_i =5.5, res =600)
## List of genes in each module
Code
input_name ="datExpr"networkType_i ="signed"minModuleSize =30deepSplit =2lfc_lim_i =0moduleColors <-read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <-c("PsCam", "moduleColors")# transformation into a listgene_list <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(genes =list(PsCam), .groups ="drop") %>%deframe() # transforme en liste nomméeload(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))# load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>%left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet ="Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by ="PsCam")# cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Down) ", comparison)), # purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Up) ", comparison))# ) %>% full_join(., df_info_gene, by = "PsCam") %>% # mutate(ID = Psat) %>% # drop_na(comparison, ID)col_order <-c("orange", "black", "tan", "saddlebrown")cluster_info <- purrr::imap_dfr(gene_list, ~tibble(PsCam = .x, cluster = .y)) %>%full_join(df_info_gene, by ="PsCam") %>%mutate(ID = Psat) %>%drop_na(ID)############# a revoirsource(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO termedf_GO =GO_on_different_group(functional_roles ="BP",group_info = cluster_info,group ="cluster",ID ="ID",top =10 )#info from excel info_from_excel =read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), skip =11) %>% dplyr::select("N°", "GENE_ID", "GENE_TYPE","DESCRIPTION", "...7","CAT") %>% dplyr::rename(PsCam ="N°",At ="...7")df_excel <- df_GO %>% dplyr::select(-c("GeneRatio", "BgRatio", "RichFactor", "FoldEnrichment", "zScore", "pvalue", "p.adjust", "qvalue", "Count","|-log10(Pval)|", "ID", "Description")) %>%drop_na(geneID) %>%separate_rows(geneID, sep ="/") %>% dplyr::group_by(geneID,group) %>%summarise(Description_GO =paste(Description_GO, collapse ="/"), .groups ="drop") %>% dplyr::rename(ModuleColor=group) %>% dplyr::rename(Psat = geneID) %>% dplyr::select(Psat, Description_GO) %>%full_join(.,cluster_info, by ="Psat") %>%left_join(., info_from_excel, by ="PsCam") %>%drop_na(cluster)colnames(df_excel)df_info <-as.data.frame(gene_list[["orange"]]) ; colnames(cyrcadian) <-"PsCam"cyrcadian <- cyrcadian %>%left_join(., df_info_gene, by ="PsCam")# export for titouanwritexl::write_xlsx(cyrcadian, path = here::here("data/multi_omic/cyrcadian.xlsx"))writexl::write_xlsx(df_excel, path = here::here("data/multi_omic/table_supp_all_genes.xlsx"))# Merge this with your original data to fill missing combinationsCluster_selected_GO_filled <-df_GO %>%# on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>%mutate(cluster =factor(cluster,levels =unique(cluster))#, # comparison = str_trim(str_remove(cluster, "\\s*\\(.*?\\)")),# sign = str_extract(cluster, "(?<=\\().*?(?=\\))") ) #%>% #mutate(comparison = forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT"))
16.11 Heatmap stats
The aim is to create heatmap portions with the stats inside, and then put them back on genes pathways.
Code
load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))df_info_gene_module_color <-read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) input_name ="datExpr"networkType_i ="signed"minModuleSize =30deepSplit =2moduleColors <-read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) colnames(moduleColors) <-c("gene_id", "moduleColors")moduleColors = moduleColors%>%left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet ="gene_symbol",col_types =c("text", "text", "text")) %>% dplyr::select(PsCam, gene_symbol) %>% dplyr::rename(gene_id = PsCam), by ="gene_id")moduleColors_select = moduleColors %>%drop_na("gene_symbol")#################### filter genes of interest# moduleColors_select = moduleColors %>% filter(moduleColors == "steelblue")#info for more simple name of genes# # input one metabolite# variable_i <- "Asp"# variable_i <- "OAS"# DAF_i <- 29v_variable <- moduleColors_select %>%pull(gene_id)expr_long <-as.data.frame(datExpr) %>%rownames_to_column(var ="sample") %>%# garder l'ID échantillonpivot_longer(-sample,names_to ="gene_id",values_to ="expression")# b. ajouter la couleur de module pour chaque gèneexpr_long_select <- expr_long %>%left_join(., moduleColors_select, by ="gene_id") %>%left_join(., metadata %>%rownames_to_column("sample"), by ="sample") %>%mutate(condition =paste0(sulfur_condition, "_", genotype)) %>%filter(gene_id %in% v_variable)#compartment_i <- "N6_N7"# measure min max of Zscore for color of the graduationdata_for_guides = expr_long_select %>% dplyr::group_by(gene_id) %>% dplyr::mutate(value_scaled =as.vector(scale(expression))) %>%ungroup() %>% dplyr::group_by(gene_id, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale =mean(value_scaled))# min(data_for_guides$mean_scale)# max(data_for_guides$mean_scale)for(variable_i in v_variable){# variable_i <- "OAS"# variable_i <- "Asp"# filter data df_select <- expr_long_select %>%filter( gene_id == variable_i#,#compartment == compartment_i ) %>%drop_na(expression) %>%as.data.frame()# make stats ylab_i <-paste0(variable_i, " in vegetative leaves") res <-stat_analyse(data = df_select,column_value ="expression",category_variables =c("condition"),grp_var ="",show_plot =TRUE,outlier_show =FALSE, label_outlier ="sample",biologist_stats =TRUE,Ylab_i = ylab_i,control_conditions ="",strip_normale =FALSE,hex_pallet =c("#003049", "#780000", "#7FACC7", "#EC323E", "#003049", "#780000", "#7FACC7", "#EC323E") ) stats_resum <- res$data_used %>%distinct(condition, gene_id, group) res_resum_raw <- res$summary_result %>%left_join(.,stats_resum, by ="condition") %>%mutate(column_value = variable_i) res_resum <- res$data_used %>%mutate(value_scaled =scale(expression)) %>%# scale based on z-score is beter than force in -1 ; 1 dplyr::group_by(genotype, sulfur_condition, condition) %>% dplyr::summarise(mean_scale =mean(value_scaled),.groups ='drop' ) %>%left_join(.,stats_resum, by ="condition")# combined_res[[variable_i]] <- res_resum final_res <- res_resum %>%mutate(color_text ="black") %>%mutate(condition =fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")) %>%#left_join(.,df_info_gene_module_color %>% dplyr::rename(gene_id = PsCam), by = "gene_id") %>% left_join(., moduleColors_select, by ="gene_id") min_max_zscore <-paste0("min:",round(min(final_res$mean_scale),2),";max:",round(max(final_res$mean_scale),2))# creation of the plot p_plot <-ggplot(final_res, aes(x = condition, y = gene_id, fill = mean_scale)) +geom_tile(color ="white") +# Pour dessiner les casesgeom_text(aes(label = group, color = color_text)) +# Ajouter les lettres statistiquesscale_fill_gradient2(low ="#1d4877", mid ="white", high ="#ee3e32", midpoint =0, limits =c(-2,2))+# Plage commune entre -2 et 2 pour tous les plots) + # Palette de couleurscale_color_manual(values =c("black"="black", "white"="white"))+theme_minimal() +theme(axis.title.y =element_blank(), # Supprime le titre de l'axe Yaxis.title.x =element_blank(), # Supprime le titre de l'axe Xaxis.text.y =element_text(size =10, face ="bold"), # Style du texte de la variableaxis.text.x =element_text(size =12, angle =90, vjust =+0.5, hjust =1.0),panel.grid.major =element_blank(), # Supprime les grillespanel.grid.minor =element_blank(),panel.border =element_rect(color = final_res %>%head(1) %>%pull(moduleColors),linewidth =2, fill =NA) ) +labs(title = final_res %>%head(1) %>%pull(gene_symbol))+#guides(fill = guide_colorbar(title = "Mean Value")) ; p_plotguides(fill ="none", color ="none") ; p_plot# export plot plot_name <-paste0(variable_i)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/", plot_name)), p_plot, height_i = .38*2+0.95, width_i =4, res_i =300,format ="svg")}#export legend# Création d'un ggplot minimal uniquement pour la légendep_legend <-ggplot(data_for_guides, aes(x = genotype, y = sulfur_condition, fill = mean_scale)) +geom_tile() +scale_fill_gradient2(low ="#1d4877", mid ="white", high ="#ee3e32", midpoint = (min(data_for_guides$mean_scale) +max(data_for_guides$mean_scale)) /2, limits =c(min(data_for_guides$mean_scale), max(data_for_guides$mean_scale)) ) +theme_void() +# Supprime les axes, titres, etc.guides(fill =guide_colorbar(title ="Mean Value") # Titre de la légende )# Extraire uniquement la légendelegend <- cowplot::get_legend(p_legend)# Affichage de la légende seulegrid::grid.newpage()grid::grid.draw(legend)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/legend")), legend, height_i =2, width_i =1, res_i =300,format ="svg")# test
Idem for metabolomic
Code
df_metabolite_global <-read_csv(here::here("data/metabolomic/output/df_metabolite_global.csv"), show_col_types =FALSE) %>%mutate(plant_num =paste(sep ="_", plant_num, condition)) %>%mutate(sulfur_condition=fct_relevel(sulfur_condition, "SS", "SD"),genotype =fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"),compartment =fct_relevel(compartment, "N4_N5", "N6_N7"),type_genotype =fct_relevel(type_genotype, "WT", "Mut"),condition =paste0(sulfur_condition, "_", genotype) ) %>%# filter(sulfur_condition == "SD") %>%filter(!(variable =="Met"& compartment =="N9_N10")) %>%# error of measure/ in WT2 = 0 in allfilter(compartment =="N6_N7") %>%filter(DAF ==29)df_metabolite_global_test <- df_metabolite_global %>%filter(variable =="Sulfate")# df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) # # input_name = "datExpr"# networkType_i = "signed"# minModuleSize = 30# deepSplit = 2# # moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) # colnames(moduleColors) <- c("gene_id", "moduleColors")# # moduleColors = moduleColors%>% # left_join(.,# read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol",# col_types = c("text", "text", "text")) %>% # dplyr::select(PsCam, gene_symbol) %>% # dplyr::rename(gene_id = PsCam), by = "gene_id")# # moduleColors_select = moduleColors %>% drop_na("gene_symbol")#################### filter genes of interest# moduleColors_select = moduleColors %>% filter(moduleColors == "steelblue")#info for more simple name of genes# # input one metabolite# variable_i <- "Asp"# variable_i <- "OAS"# DAF_i <- 29v_variable <- df_metabolite_global %>%distinct(variable) %>%pull(variable)# expr_long <- as.data.frame(datExpr) %>% # # rownames_to_column(var = "sample") %>% # garder l'ID échantillon# pivot_longer(-sample,# names_to = "gene_id",# values_to = "expression")# b. ajouter la couleur de module pour chaque gène# expr_long_select <- expr_long %>% # left_join(., moduleColors_select, by = "gene_id") %>% # left_join(., metadata %>% rownames_to_column("sample"), by = "sample") %>% # mutate(condition = paste0(sulfur_condition, "_", genotype)) %>% # filter(gene_id %in% v_variable)#compartment_i <- "N6_N7"# measure min max of Zscore for color of the graduationdata_for_guides = df_metabolite_global %>% dplyr::group_by(variable) %>% dplyr::mutate(value_scaled =as.vector(scale(value))) %>%ungroup() %>% dplyr::group_by(variable, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale =mean(value_scaled))# min(data_for_guides$mean_scale)# max(data_for_guides$mean_scale)for(variable_i in v_variable){# variable_i <- "OAS"# variable_i <- "Asp"# filter data df_select <- df_metabolite_global %>%filter( variable == variable_i#,#compartment == compartment_i ) %>%drop_na(value) %>%as.data.frame()# make stats ylab_i <-paste0(variable_i, " in vegetative leaves") res <-stat_analyse(data = df_select,column_value ="value",category_variables =c("condition"),grp_var ="",show_plot =TRUE,outlier_show =FALSE, label_outlier ="plant_num",biologist_stats =TRUE,Ylab_i = ylab_i,control_conditions ="",strip_normale =FALSE,hex_pallet =c("#003049", "#780000", "#7FACC7", "#EC323E", "#003049", "#780000", "#7FACC7", "#EC323E") ) stats_resum <- res$data_used %>%distinct(condition, variable, group) res_resum_raw <- res$summary_result %>%left_join(.,stats_resum, by ="condition") %>%mutate(column_value = variable_i) res_resum <- res$data_used %>%mutate(value_scaled =scale(value)) %>%# scale based on z-score is beter than force in -1 ; 1 dplyr::group_by(genotype, sulfur_condition, condition) %>% dplyr::summarise(mean_scale =mean(value_scaled),.groups ='drop' ) %>%left_join(.,stats_resum, by ="condition")# combined_res[[variable_i]] <- res_resum# # # final_res <- res_resum %>% # # mutate(color_text = ifelse(mean_scale<-2, "white", "black")) %>% # mutate(color_text = "black") %>% # mutate(condition = fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K"))# %>% # #left_join(.,df_info_gene_module_color %>% dplyr::rename(gene_id = PsCam), by = "gene_id") %>% # # left_join(., moduleColors_select, by = "gene_id") condition_levels <-c("SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K","SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")# ------------------------------------------------------------------# à INSÉRER juste avant la création du graphique, dans ta boucle# ------------------------------------------------------------------final_res <- res_resum %>% tidyr::complete( variable,condition =factor(condition_levels, levels = condition_levels),fill =list(mean_scale =NA_real_, group ="") ) %>%mutate(is_na =is.na(mean_scale), # <‑‑ nouvel indicateurcolor_text ="black",condition =fct_relevel(condition, !!!condition_levels) )# ------------------------------------------------------------------# 2e moitié inchangée, sauf na.value ajouté# ------------------------------------------------------------------p_plot <-ggplot(final_res, aes(x = condition, y = variable, fill = mean_scale)) +geom_tile(color ="white") +geom_text(aes(label = group, colour = color_text)) +geom_point(data =function(x) dplyr::filter(x, is_na),shape =4, # forme “×”stroke =1.2, # épaisseur des traits de la croixsize =5,colour ="white" ) +scale_fill_gradient2(low ="#55185D", mid ="white", high ="#ECB602",midpoint =0, limits =c(-2, 2), na.value ="gray60" ) +scale_colour_manual(values =c("black"="black", "white"="white")) +theme_minimal() +theme(axis.title =element_blank(),axis.text.y =element_text(size =10, face ="bold"),axis.text.x =element_text(size =12, angle =90, vjust =0.5, hjust =1),panel.grid =element_blank() ) +guides(fill ="none", colour ="none") min_max_zscore <-paste0("min:", round(min(final_res$mean_scale, na.rm =TRUE), 2),";max:", round(max(final_res$mean_scale, na.rm =TRUE), 2))print(min_max_zscore)# creation of the plot # p_plot <- ggplot(final_res, aes(x = condition, y = variable, fill = mean_scale)) +# geom_tile(color = "white") + # Pour dessiner les cases# geom_text(aes(label = group, color = color_text)) + # Ajouter les lettres statistiques# scale_fill_gradient2(low = "#55185D", mid = "white", high = "#ECB602", midpoint = 0, limits = c(-2,2))+# Plage commune entre -2 et 2 pour tous les plots) + # Palette de couleur# scale_color_manual(values = c("black" = "black", "white" = "white"))+# theme_minimal() +# theme(# axis.title.y = element_blank(), # Supprime le titre de l'axe Y# axis.title.x = element_blank(), # Supprime le titre de l'axe X# axis.text.y = element_text(size = 10, face = "bold"), # Style du texte de la variable# axis.text.x = element_text(size = 12, angle = 90, vjust = +0.5, hjust = 1.0),# panel.grid.major = element_blank(), # Supprime les grilles# panel.grid.minor = element_blank()# ) +# #labs(title = final_res %>% head(1) %>% pull(variable))+# #guides(fill = guide_colorbar(title = "Mean Value")) ; p_plot# guides(fill = "none", color = "none") ; p_plot# export plot plot_name <-paste0(variable_i)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/metabolomic/", plot_name)), p_plot, height_i = .38*1.4+0.95, width_i =4, res_i =300,format ="svg")}#export legend# Création d'un ggplot minimal uniquement pour la légendep_legend <-ggplot(data_for_guides, aes(x = genotype, y = sulfur_condition, fill = mean_scale)) +geom_tile() +scale_fill_gradient2(low ="#55185D", mid ="white", high ="#ECB602", midpoint = (min(data_for_guides$mean_scale) +max(data_for_guides$mean_scale)) /2, limits =c(min(data_for_guides$mean_scale), max(data_for_guides$mean_scale)) ) +theme_void() +# Supprime les axes, titres, etc.guides(fill =guide_colorbar(title ="Mean Value") # Titre de la légende )# Extraire uniquement la légendelegend <- cowplot::get_legend(p_legend)# Affichage de la légende seulegrid::grid.newpage()grid::grid.draw(legend)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/metabolomic/legend")), legend, height_i =2, width_i =1, res_i =300,format ="svg")
16.12 Heatmap cathegorie for microarray
16.12.1 If i use the normalized data (same as used for WGCNA)
Code
# data importation # Calculates the number of elements which satisfy the above condition (i.e. the number of elements which are neither NA nor 0) and checks whether this number is greater than 10% of the total length of the column (length(.)). If this is the case, it means that the column has more than 10% of values that are neither NA nor 0, and should therefore be kept.#group_pathway = compartment_iload(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))df_info_gene_module_color <-read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) input_name ="datExpr"networkType_i ="signed"minModuleSize =30deepSplit =2moduleColors <-read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) colnames(moduleColors) <-c("gene_id", "moduleColors")moduleColors = moduleColors%>%left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet ="gene_symbol",col_types =c("text", "text", "text", "text", "text", "text")) %>% dplyr::select(PsCam, gene_symbol, group_pathway, functional_class) %>% dplyr::rename(gene_id = PsCam), by ="gene_id")moduleColors_select = moduleColors %>%drop_na("gene_symbol")v_variable <- moduleColors_select %>%pull(gene_id)expr_long <-as.data.frame(datExpr) %>%rownames_to_column(var ="sample") %>%# garder l'ID échantillonpivot_longer(-sample,names_to ="gene_id",values_to ="expression")perform_pca <-function(data) {if (any(is.na(data))) { method <-"Method: missMDA (NA value in dataset)" nb <-estim_ncpPCA(data, quali.sup =1:4, method.cv ="Kfold", verbose =FALSE)plot(0:5, nb$criterion, xlab ="Number of dimensions", ylab ="MSEP") res.comp <-imputePCA(data, quali.sup =1:4, ncp = nb$ncp) PCA_result <-PCA(res.comp$completeObs, quali.sup =1:4, graph =TRUE) } else { method <-"Method: Normal PCA (no NA value in dataset)" PCA_result <-PCA(data, quali.sup =1:4, graph =TRUE) }return(list(PCA_result = PCA_result, method = method))}visualize_eigenvalues <-function(PCA_result) {fviz_eig(PCA_result, addlabels =TRUE, ylim =c(0, 50))}visualize_pca_variables <-function(PCA_result, group_pathway_i) { var <-get_pca_var(PCA_result)png(here::here(paste0("report/microarray/plot/PCA/PCA_seleced_genes_var_", group_pathway_i, ".png")), width =16, height =16, units ='cm', res =600)print(fviz_pca_var(PCA_result, col.var ="cos2",gradient.cols =c("#00AFBB", "#E7B800", "#FC4E07"),repel =TRUE))dev.off()}# b. ajouter la couleur de module pour chaque gèneexpr_long_select <- expr_long %>%left_join(., moduleColors_select, by ="gene_id") %>%left_join(., metadata %>%rownames_to_column("sample"), by ="sample") %>%mutate(condition =paste0(sulfur_condition, "_", genotype)) %>%filter(gene_id %in% v_variable)# measure min max of Zscore for color of the graduationdata_for_guides = expr_long_select %>% dplyr::group_by(gene_id) %>% dplyr::mutate(value_scaled =as.vector(scale(expression))) %>%ungroup() %>% dplyr::group_by(gene_id, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale =mean(value_scaled))prepare_data<-function(group_pathway_i){ expr_wide_select<- expr_long_select %>%filter(group_pathway == group_pathway_i) %>% dplyr::select(-c("gene_id", "moduleColors", "group_pathway", "Mutant_type","sulfur_condition","chip_color","chip_num","genotype", "functional_class")) %>%pivot_wider(names_from ="gene_symbol", values_from ="expression") %>%column_to_rownames("sample")}#compartment_i <- "N6_N7"# Function to calculate means by typecalculate_means_by_type <-function(data) { data %>% dplyr::group_by(condition) %>% dplyr::summarise(across(where(is.numeric), mean, na.rm =TRUE))}data = expr_wide_select # for teste to delcalculate_letter_signif <-function(data) {# data: a data frame where the first column is "climat_condition"# and the other columns represent metabolite values for each sample# Extract metabolite names (columns after "climat_condition") metabo_names <-colnames(data)[-1]# Replace hyphens and whitespace with underscores# Change colname of the dataframe# Initialize a list to store significance letters for each metabolite signif_letters <-list()# Loop through each metabolitefor (metabo in metabo_names) {#cat(metabo, "\n")# Build the ANOVA model# aov_model <- aov(as.formula(paste0("`", metabo, "` ~ condition")), data = data)# # # Apply the Tukey HSD test# tukey_result <- TukeyHSD(aov_model, "condition")# # # Extract the p-values for the condition comparisons# tukey_groups <- multcompLetters4(aov_model, tukey_result)# # # Retrieve significance letters for this metabolite# letters <- tukey_groups$`condition`$Letters####test new #### data2 <- data %>% dplyr::rename(response =all_of(metabo)) # on sécurise le nom aov_model <-aov(response ~ condition, data = data2) tukey <-TukeyHSD(aov_model, "condition") letters <-multcompLetters4(aov_model, tukey) letters <- letters$condition# Check if all letters are the sameif (length(unique(letters)) ==1) {# If all letters are the same, replace them with empty strings letters <-rep("", length(letters)) }# add the same name as before# Store the letters in the list signif_letters[[metabo]] <- letters }return(signif_letters)}# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional per‑cell numeric values.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; if TRUE (default) display numbers inside cells## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional numeric values and/or# significance letters inside each cell.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; show numeric means if TRUE (default FALSE)# show_letters : logical; show df_letter_signif letters if TRUE# (default TRUE)## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------create_heatmap <-function(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, type_colors, df_letter_signif, group_pathway_i, cutree_rows_i,show_values =FALSE,show_letters =TRUE) {# Replace non‑breaking spaces by standard spaces in 'functional_class' type_metabo_s$functional_class <- stringi::stri_replace_all_fixed( type_metabo_s$functional_class,"\u00A0", # non‑breaking space" ", # regular spacevectorize_all =FALSE ) type_metabo_s$functional_class <-factor(type_metabo_s$functional_class)# Order the mean table to match matrix columns mean_by_type$condition <-factor(mean_by_type$condition,levels =colnames(z_matrix_transposed_reordered)) mean_by_type_reorder <- mean_by_type[order(mean_by_type$condition), ]# -------------------------------------------------------------------------# Column annotations# ------------------------------------------------------------------------- col_labels <-colnames(z_matrix_transposed_reordered) col_annotation <-data.frame(Sulfur_Condition =ifelse(startsWith(col_labels, "SS_"), "SS", "SD"),Genotype =sub("^[A-Z]{2}_", "", col_labels),row.names = col_labels ) |> dplyr::mutate(Sulfur_Condition = forcats::fct_relevel(Sulfur_Condition, "SS", "SD"),Genotype = forcats::fct_relevel(Genotype, "WT1", "W78*", "WT2", "E568K") ) sulfur_pal <- sulfate_pallet # yellow / cyan genotype_pal <- mutant_palette # custom palette annotation_colors <-list(Sulfur_Condition = sulfur_pal,Genotype = genotype_pal,Condition = type_colors,moduleColors =setNames(moduleColors_select$moduleColors, moduleColors_select$moduleColors),functional_class =c(# Biosynthesis"Cysteine biosynthesis (O‑acetyl‑serine pathway)"="#7d092f","Glutathione biosynthesis & turnover"="#b32d45",# Methionine cycle / salvage"Methionine cycle & salvage"="#ef5659","Methionine cycle & salvage"="#D95F02",# Transcription factors"NAC transcription factors"="#ffb255","WRKY transcription factors"="#ffd673",# Regulation / signalling"Regulation / sulfur‑deficiency signalling"="#c5c35e",# Activation / reduction / transport / detox"Sulfate activation (entry step)"="#42a285","Sulfate reduction"="#577590","Sulfate uptake & distribution"="#6d597a","Sulfite detoxification & export"="#7D3C98" ) )# -------------------------------------------------------------------------# Heat‑map# ------------------------------------------------------------------------- ComplexHeatmap::pheatmap( z_matrix_transposed_reordered,cluster_rows =TRUE,cluster_cols =FALSE,show_rownames =TRUE,show_colnames =TRUE,cutree_rows = cutree_rows_i,labels_col = col_labels,annotation_col = col_annotation,annotation_colors = annotation_colors,annotation_row = type_metabo_s,main = stringr::str_to_title(group_pathway_i),color =colorRampPalette(c("#1d4877", "white", "#ee3e32"))(256),name ="Centering and Scaling",# Add text according to user’s choicescell_fun =function(j, i, x, y, width, height, fill) {if (isTRUE(show_values) ||isTRUE(show_letters)) {# Retrieve numeric value only if neededif (isTRUE(show_values)) { value_num <-round(apply(t(mean_by_type_reorder |> tibble::column_to_rownames("condition")), 2, as.numeric)[i, j],2 ) } else { value_num <-"" }# Retrieve significance letter letter <-if (isTRUE(show_letters)) df_letter_signif[i, j] else""# Decide what to print text_value <-""if (isTRUE(show_values) && value_num !="") { text_value <- value_num }if (isTRUE(show_letters) && letter !="") { text_value <-if (text_value !="") {paste0(text_value, " (", letter, ")") } else { letter } }# Draw only if something to drawif (text_value !="") { grid::grid.text(text_value, x, y, gp = grid::gpar(fontsize =8)) } } } )}# Main PCA functionmake_pca_condition <-function(group_pathway_i ="sulfate pathway", cutree_rows_i =7, show_values_i =TRUE) {# Step 1: Prepare the data all_data_select_h_x <-prepare_data(group_pathway_i)# Step 2: Perform PCA pca_results <-perform_pca(all_data_select_h_x) res_pca_targeted_metabolomic_r2_x <- pca_results$PCA_result method <- pca_results$method# Step 3: Visualize eigenvaluesvisualize_eigenvalues(res_pca_targeted_metabolomic_r2_x)# Step 4: Visualize PCA variablesvisualize_pca_variables(res_pca_targeted_metabolomic_r2_x, group_pathway_i)# Step 5: Calculate means by type and letter significant and prepare for heatmap type_metabo <- moduleColors_select %>%filter(group_pathway == group_pathway_i) %>% dplyr::select(moduleColors, gene_symbol, functional_class) %>%unique() %>% dplyr::rename(variable = gene_symbol) df_obs <- res_pca_targeted_metabolomic_r2_x$call$X #%>%#dplyr::select(-c(condition)) type_metabo_s <- type_metabo %>% dplyr::select(variable, functional_class, moduleColors) %>%# rename(Type = type) %>%filter(variable %in%colnames(df_obs)[-1]) %>%column_to_rownames("variable")# %>% # mutate(Type = recode(Type,# "AA" = "Amino acid",# "PHY" = "Phytohormone",# "AO" = "Organic acid",# "FLAVO" = "Flavonoid",# "CYT" = "Cytokinine",# "SUGAR" = "Sugar")) mean_by_type <-calculate_means_by_type(df_obs) letter_signif <-calculate_letter_signif(df_obs)# Step 6: Calculate Z-scores z_scores<-mean_by_type %>%column_to_rownames("condition") z_scores <-scale(z_scores) z_matrix <-as.matrix(z_scores) z_matrix_transposed <-t(z_matrix) z_matrix_transposed_reordered <- z_matrix_transposed[,c("SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")]# Step 7: Create table of letter significant (same as z_matrix_transposed) df_letter_signif <-as.data.frame(matrix(data ="", nrow =dim(z_matrix_transposed_reordered)[1], ncol =dim(z_matrix_transposed_reordered)[2]))colnames(df_letter_signif)<-colnames(z_matrix_transposed_reordered)rownames(df_letter_signif)<-rownames(z_matrix_transposed_reordered)for (metabolite_i inrownames(z_matrix_transposed_reordered)){for (condition_i incolnames(z_matrix_transposed_reordered)){ df_letter_signif[metabolite_i,condition_i] <- letter_signif[[metabolite_i]]$Letters[condition_i] } } df_letter_signif[is.na(df_letter_signif)] <-""# Step 8: Create heatmap type_colors <-c(setNames(as.character(mutant_palette), paste0("SS_", names(mutant_palette))),setNames(as.character(mutant_palette), paste0("SD_", names(mutant_palette))) ) # Define this variable as needed pheat<-create_heatmap(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, type_colors, df_letter_signif, group_pathway_i, cutree_rows_i, show_values = show_values_i )# Step 9: PCA biplot PCA_X <-fviz_pca_biplot(res_pca_targeted_metabolomic_r2_x,geom.ind ="point",pointshape =20,pointsize =3,col.ind = all_data_select_h_x$condition,col.quanti.sup ="red",palett = type_colors,addEllipses =TRUE,legend.title ="Treatment", ellipse.level =0.95,ellipse.alpha =0.3,ellipse.type ="convex",arrowsize =0.5,col.var ="black", axes =c(1, 2), repel = T ) +#ggtitle(label = paste0("Visualizing Individual PCA and Quantity Variable for ", group_pathway, " for Each Treatment"),# subtitle = method) +ggtitle(label =str_to_title(group_pathway_i))theme_classic()ggsave(filename = here::here(paste0("report/microarray/plot/PCA/PCA_seleced_genes_var_path_form", group_pathway_i, "_by_treatment.svg")), plot = PCA_X, width =18*1.5, height =15*1.5, units ="cm")dev.off()svg(filename = here::here(paste0("report/microarray/plot/pheatmap/pheatmap_", group_pathway_i, "_by_treatment.svg")), width =13 , height =10)print(pheat)dev.off()return(list(PCA_X = PCA_X, pheat = pheat)) # Adjust the return values as needed}# execution of the functionpca_heat_sulfate_pathway=make_pca_condition("sulfate pathway")pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i =2)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway")), p_sulfate_pathway, height_i =13, width_i =13, res_i =600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway")), p_scenescence_pathway, height_i =5, width_i =13, res_i =600)####### without value ###########pca_heat_sulfate_pathway=make_pca_condition("sulfate pathway", show_values_i =FALSE)pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i =2, show_values_i =FALSE)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway_without_value")), p_sulfate_pathway, height_i =10, width_i =9, res_i =600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway_without_value")), p_scenescence_pathway, height_i =4, width_i =8, res_i =600)
16.12.2 If i use the the data (and stats from the platforme)
Code
# data importation
# Calculates the number of elements which satisfy the above condition (i.e. the number of elements which are neither NA nor 0) and checks whether this number is greater than 10% of the total length of the column (length(.)). If this is the case, it means that the column has more than 10% of values that are neither NA nor 0, and should therefore be kept.
#group_pathway = compartment_i
# load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))
load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData"))
lfc_lim_i <- 0
list_ratio_stats
diff_exression <- list_ratio_stats %>%
map(~ select(.x, id_probe, logFC, BH)) %>% # garde seulement les 3 colonnes
bind_rows(.id = "source") %>%
mutate(
comparison = case_when(
source == "Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1" ~ "SD_WT2 vs SD_E568K",
source == "Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3" ~ "SS_WT2 vs SS_E568K",
source == "Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2" ~ "SD_WT1 vs SD_W78*",
source == "Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4" ~ "SS_WT1 vs SS_W78*",
source == "WT_SD_RepBio_5 vs WT_SS_RepBio_5" ~ "SS_WT vs SD_WT",
TRUE ~ source # si aucun des cas ne matche, on garde la valeur d'origine
)
) %>% dplyr::select(-c(source)) %>%
dplyr::rename(gene_id = id_probe)
# levels(as.factor(resultat$comparison))
load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))
df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx"))
input_name = "datExpr"
networkType_i = "signed"
minModuleSize = 30
deepSplit = 2
moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))
colnames(moduleColors) <- c("gene_id", "moduleColors")
moduleColors_select = moduleColors%>%
left_join(.,
read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol",
col_types = c("text", "text", "text", "text", "text", "text")) %>%
dplyr::select(PsCam, gene_symbol, group_pathway, functional_class) %>%
dplyr::rename(gene_id = PsCam), by = "gene_id") %>%
drop_na("gene_symbol")
v_variable <- moduleColors_select %>% pull(gene_id)
# b. ajouter la couleur de module pour chaque gène
expr_long_select <- diff_exression %>%
left_join(., moduleColors_select, by = "gene_id") %>%
drop_na(gene_symbol)
# measure min max of Zscore for color of the graduation
data_for_guides = expr_long_select %>%
dplyr::group_by(gene_id) %>%
dplyr::mutate(value_scaled = as.vector(scale(logFC))) %>%
ungroup() %>%
dplyr::group_by(gene_id, comparison) %>%
dplyr::summarise(mean_scale = mean(value_scaled))
prepare_data<-function(group_pathway_i){
expr_wide_select<- expr_long_select %>%
filter(group_pathway == group_pathway_i) %>%
dplyr::select(-c("gene_id", "moduleColors", "group_pathway","functional_class", "BH")) %>%
pivot_wider(names_from = "gene_symbol", values_from = "logFC")
return(expr_wide_select)
}
prepare_data_BH<-function(group_pathway_i){
expr_wide_select<- expr_long_select %>%
filter(group_pathway == group_pathway_i) %>%
dplyr::select(-c("gene_id", "moduleColors", "group_pathway","functional_class", "logFC")) %>%
pivot_wider(names_from = "gene_symbol", values_from = "BH")
return(expr_wide_select)
}
#compartment_i <- "N6_N7"
# Function to calculate means by type
calculate_means_by_type <- function(data) {
data %>%
dplyr::group_by(comparison) %>%
dplyr::summarise(across(where(is.numeric), mean, na.rm = TRUE))
}
# data = expr_wide_select # for teste to del
# calculate_letter_signif <- function(data) {
# # data: a data frame where the first column is "climat_condition"
# # and the other columns represent metabolite values for each sample
#
# # Extract metabolite names (columns after "climat_condition")
# metabo_names <- colnames(data)[-1]
#
# # Replace hyphens and whitespace with underscores
#
# # Change colname of the dataframe
#
# # Initialize a list to store significance letters for each metabolite
# signif_letters <- list()
#
# # Loop through each metabolite
# for (metabo in metabo_names) {
# #cat(metabo, "\n")
# # Build the ANOVA model
# # aov_model <- aov(as.formula(paste0("`", metabo, "` ~ condition")), data = data)
# #
# # # Apply the Tukey HSD test
# # tukey_result <- TukeyHSD(aov_model, "condition")
# #
# # # Extract the p-values for the condition comparisons
# # tukey_groups <- multcompLetters4(aov_model, tukey_result)
# #
# # # Retrieve significance letters for this metabolite
# # letters <- tukey_groups$`condition`$Letters
#
# ####test new ####
# data2 <- data %>%
# dplyr::rename(response = all_of(metabo)) # on sécurise le nom
#
# aov_model <- aov(response ~ condition, data = data2)
# tukey <- TukeyHSD(aov_model, "condition")
# letters <- multcompLetters4(aov_model, tukey)
#
# letters <- letters$condition
#
# # Check if all letters are the same
# if (length(unique(letters)) == 1) {
# # If all letters are the same, replace them with empty strings
# letters <- rep("", length(letters))
# }
#
# # add the same name as before
#
# # Store the letters in the list
# signif_letters[[metabo]] <- letters
# }
#
# return(signif_letters)
# }
calculate_star_signif <- function(data_BH) {
p_to_stars <- function(p) {
ifelse(p < 0.001, "***",
ifelse(p < 0.01, "**",
ifelse(p < 0.05, "*", "")))
}
# data: a data frame where the first column is "climat_condition"
# and the other columns represent metabolite values for each sample
# Extract metabolite names (columns after "climat_condition")
metabo_names <- colnames(data_BH[-1])
# Replace hyphens and whitespace with underscores
# Change colname of the dataframe
# Initialize a list to store significance letters for each metabolite
signif_stars <- vector("list", length(metabo_names))
names(signif_stars) <- metabo_names
# Loop through each metabolite
comp_labels <- as.character(data_BH[[1]]) # noms des comparaisons
for (metabo in metabo_names) {
stars_vec <- p_to_stars(data_BH[[metabo]])
names(stars_vec) <- comp_labels # on nomme le vecteur
signif_stars[[metabo]] <- stars_vec
}
return(signif_stars)
}
# ---------------------------------------------------------------------------
# create_heatmap
# ---------------------------------------------------------------------------
# Draws a ComplexHeatmap::pheatmap with optional per‑cell numeric values.
#
# Args:
# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)
# mean_by_type : data‑frame with mean values by condition
# type_metabo_s : data‑frame with row annotations
# type_colors : named vector of colours for 'Condition'
# df_letter_signif : matrix of significance letters, same dim as z_matrix_
# group_pathway_i : string used for the heat‑map title
# cutree_rows_i : integer; number of clusters for rows
# show_values : logical; if TRUE (default) display numbers inside cells
#
# Returns:
# A heat‑map plot (invisible; printed for its side‑effect).
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# create_heatmap
# ---------------------------------------------------------------------------
# Draws a ComplexHeatmap::pheatmap with optional numeric values and/or
# significance letters inside each cell.
#
# Args:
# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)
# mean_by_type : data‑frame with mean values by condition
# type_metabo_s : data‑frame with row annotations
# type_colors : named vector of colours for 'Condition'
# df_letter_signif : matrix of significance letters, same dim as z_matrix_
# group_pathway_i : string used for the heat‑map title
# cutree_rows_i : integer; number of clusters for rows
# show_values : logical; show numeric means if TRUE (default FALSE)
# show_letters : logical; show df_letter_signif letters if TRUE
# (default TRUE)
#
# Returns:
# A heat‑map plot (invisible; printed for its side‑effect).
# ---------------------------------------------------------------------------
create_heatmap_star <- function(z_matrix_transposed_reordered,
mean_by_type,
type_metabo_s,
type_colors,
df_star_signif,
group_pathway_i,
cutree_rows_i,
show_values = FALSE,
show_star = TRUE) {
# Replace non‑breaking spaces by standard spaces in 'functional_class'
type_metabo_s$functional_class <- stringi::stri_replace_all_fixed(
type_metabo_s$functional_class,
"\u00A0", # non‑breaking space
" ", # regular space
vectorize_all = FALSE
)
type_metabo_s <- type_metabo_s %>%
mutate(
functional_class = as.factor(functional_class) %>% # colonne d'origine
# stringi::stri_replace_all_fixed("\u00A0", " ", FALSE) %>% # remplace NBS
forcats::fct_relevel( # <-- ordre voulu ICI
"Sulfate uptake & distribution",
"Sulfate reduction",
"Sulfite detoxification & export",
"Cysteine biosynthesis",
"Glutathione biosynthesis & turnover",
"Methionine cycle & salvage",
"Methionine cycle & salvage",
"Regulation / sulfur deficiency signalling",
"NAC transcription factors",
"WRKY transcription factors",
"Sulfate activation (entry step)"
)
)
# Order the mean table to match matrix columns
mean_by_type$comparison <- factor(mean_by_type$comparison,
levels = colnames(z_matrix_transposed_reordered))
mean_by_type_reorder <- mean_by_type[order(mean_by_type$comparison), ]
# -------------------------------------------------------------------------
# Column annotations
# -------------------------------------------------------------------------
col_labels <- colnames(z_matrix_transposed_reordered)
col_annotation <- data.frame(
Sulfur_Condition = c("WT", "SS","SS","SD","SD"),
row.names = col_labels
) %>%
mutate(Sulfur_Condition = forcats::fct_relevel(Sulfur_Condition, "WT", "SS","SD"))
# sulfur_pal <- sulfate_pallet # yellow / cyan
sulfur_pal <- c(
"WT" = "#9381FF", # noir (SS_WT vs SD_WT)
"SS" = as.character(sulfate_pallet[1]), # vert (comparaisons SS only)
"SD" = as.character(sulfate_pallet[2]) # jaune (comparaisons SD only)
)
annotation_colors <- list(
Sulfur_Condition = sulfur_pal,
moduleColors = setNames(moduleColors_select$moduleColors,
moduleColors_select$moduleColors),
functional_class = c(
# Biosynthesis
"Sulfate uptake & distribution" = "#7d092f",
"Sulfate reduction" = "#ef5659",#"#b32d45",
"Sulfite detoxification & export" = "#D95F02",
"Cysteine biosynthesis" = "#ffd673",
"Glutathione biosynthesis & turnover" = "#42a285",#"#ffb255",
# Methionine cycle / salvage
"Methionine cycle & salvage" = "#6d597a",
"Methionine cycle & salvage" = "#6d597a",
# Transcription factors
# Regulation / signalling
"Regulation / sulfur deficiency signalling" = "#7D3C98",
"NAC transcription factors" = "#c5c35e",
"WRKY transcription factors" = "#ef5659"
# ,
# Activation / reduction / transport / detox
# "Sulfate activation (entry step)" = "#c5c35e",
# ,
# ,
# = "#7D3C98"
)#%>%
# mutate(functional_class = forcats::fct_relevel(functional_class,
# "Sulfate uptake & distribution" ,
# "Sulfate reduction" ,
# "Sulfite detoxification & export" ,
# "Cysteine biosynthesis" ,
# "Glutathione biosynthesis & turnover" ,
# # Methionine cycle / salvage
# "Methionine cycle & salvage" ,
# "Methionine cycle & salvage" ,
# "Regulation / sulfur deficiency signalling" ,
# # Transcription factors
# "NAC transcription factors" ,
# "WRKY transcription factors" ,
# # Regulation / signalling
# # Activation / reduction / transport / detox
# "Sulfate activation (entry step)"
#
#
#
# ))
)
# -------------------------------------------------------------------------
# Heat‑map
# -------------------------------------------------------------------------
ComplexHeatmap::pheatmap(
z_matrix_transposed_reordered,
cluster_rows = TRUE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
cutree_rows = cutree_rows_i,
labels_col = col_labels,
annotation_col = col_annotation,
annotation_colors = annotation_colors,
annotation_row = type_metabo_s,
main = stringr::str_to_title(group_pathway_i),
color = colorRampPalette(c("#1d4877", "white", "#ee3e32"))(256),
name = "Log2 ratio",
# Add text according to user’s choices
cell_fun = function(j, i, x, y, width, height, fill) {
if (isTRUE(show_values) || isTRUE(show_star)) {
# Retrieve numeric value only if needed
if (isTRUE(show_values)) {
value_num <- round(
apply(t(mean_by_type_reorder |>
tibble::column_to_rownames("comparison")), 2, as.numeric)[i, j],
2
)
} else {
value_num <- ""
}
# Retrieve significance letter
star <- if (isTRUE(show_star)) df_star_signif[i, j] else ""
# Decide what to print
text_value <- ""
if (isTRUE(show_values) && value_num != "") {
text_value <- value_num
}
if (isTRUE(show_star) && star != "") {
text_value <- if (text_value != "") {
paste0(text_value, " (", star, ")")
} else {
star
}
}
# Draw only if something to draw
if (text_value != "") {
grid::grid.text(text_value, x, y, gp = grid::gpar(fontsize = 8))
}
}
}
)
}
# Main PCA function
make_star_condition <- function(group_pathway_i = "sulfate pathway", cutree_rows_i = 7, show_values_i = TRUE) {
# Step 1: Prepare the data
all_data_select_h_x <- prepare_data(group_pathway_i)
all_data_select_h_x_BH <- prepare_data_BH(group_pathway_i)
# Step 5: Calculate means by type and letter significant and prepare for heatmap
type_metabo <- moduleColors_select %>%
filter(group_pathway == group_pathway_i) %>%
dplyr::select(moduleColors, gene_symbol, functional_class) %>%
unique() %>%
dplyr::rename(variable = gene_symbol)
type_metabo_s <- type_metabo %>%
dplyr::select(variable, functional_class, moduleColors) %>%
filter(variable %in% colnames(all_data_select_h_x)[-1]) %>%
column_to_rownames("variable")# %>%
# mutate(Type = recode(Type,
# "AA" = "Amino acid",
# "PHY" = "Phytohormone",
# "AO" = "Organic acid",
# "FLAVO" = "Flavonoid",
# "CYT" = "Cytokinine",
# "SUGAR" = "Sugar"))
mean_by_type <- calculate_means_by_type(all_data_select_h_x)
star_signif <- calculate_star_signif(all_data_select_h_x_BH)
# Step 6: Calculate Z-scores
z_scores<-mean_by_type %>% column_to_rownames("comparison")
# z_scores <- scale(z_scores)
z_matrix <- as.matrix(z_scores)
z_matrix_transposed <- t(z_matrix)
z_matrix_transposed_reordered <- z_matrix_transposed[,c("SS_WT vs SD_WT", "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K")]
# Step 7: Create table of letter significant (same as z_matrix_transposed)
df_star_signif <- as.data.frame(matrix(data = "", nrow = dim(z_matrix_transposed_reordered)[1], ncol = dim(z_matrix_transposed_reordered)[2]))
colnames(df_star_signif)<-colnames(z_matrix_transposed_reordered)
rownames(df_star_signif)<-rownames(z_matrix_transposed_reordered)
for (metabolite_i in rownames(z_matrix_transposed_reordered)){
for (condition_i in colnames(z_matrix_transposed_reordered)){
df_star_signif[metabolite_i,condition_i] <- star_signif[[metabolite_i]][condition_i]
}
}
df_star_signif[is.na(df_star_signif)] <- ""
# Step 8: Create heatmap
# type_colors <- c(
# setNames(as.character(mutant_palette), paste0("SS_", names(mutant_palette))),
# setNames(as.character(mutant_palette), paste0("SD_", names(mutant_palette)))
# ) # Define this variable as needed
pheat<-create_heatmap_star(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, df_star_signif = df_star_signif, group_pathway = group_pathway_i, cutree_rows = cutree_rows_i, show_values = show_values_i)
return(list(pheat = pheat)) # Adjust the return values as needed
}
# execution of the function
pca_heat_sulfate_pathway=make_star_condition("sulfate pathway")
pca_heat_scenescence_pathway=make_star_condition("scenescence pathway", cutree_rows_i = 2)
# Convertit le grob ComplexHeatmap en ggplot
p_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))
p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))
# export figure
fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway")), p_sulfate_pathway, height_i = 13, width_i = 13, res_i = 600)
fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway")), p_scenescence_pathway, height_i = 5, width_i = 13, res_i = 600)
####### without value ###########
pca_heat_sulfate_pathway=make_pca_condition("sulfate pathway", show_values_i = FALSE)
pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i = 2, show_values_i = FALSE)
# Convertit le grob ComplexHeatmap en ggplot
p_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))
p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))
# export figure
fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway_without_value_star")), p_sulfate_pathway, height_i = 10, width_i = 8, res_i = 600)
fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway_without_value_star")), p_scenescence_pathway, height_i = 4, width_i = 8, res_i = 600)
Source Code
---output: html_documenteditor_options: chunk_output_type: console---## WGCNA ```{r}#pkglibrary(tidyverse)library(readxl)library(WGCNA) library(ggcorrplot)library(moduleColor)library(genefilter) # for delet low variancelibrary(DESeq2)library(ggnewscale) # si tu veux plusieurs scales_fill plus tardlibrary(ggh4x)library(colorspace) library(pheatmap)library(WGCNA)library(doParallel)library(parallel)library(igraph)# srcsource(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot source(here::here("src/function/fig_export.R")) # This function saves a given plot (plot_x) as both a PDF and a high-resolution PNG file at specified dimensions.source(here::here("src/function/Evaluate_contrast.R"))# cosmeticssulfate_pallet=read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="sulfure_condition") %>% dplyr::select(color, treatment) %>%pull(color) %>%setNames(read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="sulfure_condition") %>%pull(treatment) )mutant_palette=read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="mutant") %>% dplyr::select(color, treatment) %>%pull(color) %>%setNames(read_excel(here::here("data/color_palette.xlsm")) %>%filter(set =="mutant") %>%pull(treatment) )my_color_palette <-read_excel(here::here("data/color_palette.xlsm"))# functionstatic_colors=function(hier1, h1=0.9,minsize1=50) {# here we define modules by using a height cut-off for the branches labelpred=cutree(hier1,h=h1) sort1=-sort(-table(labelpred)) modulename=as.numeric(names(sort1)) modulebranch= sort1>minsize1 no.modules=sum(modulebranch)# now we assume that there are fewer than a certain number of colors colorcode=GlobalStandardColors# "grey" means not in any module; colorhelp=rep("grey",length(labelpred))if ( no.modules==0| no.modules >length(colorcode)){ print(paste("The number of modules is problematic. Number of modules = ", as.character(no.modules)))} else { for (i inc(1:no.modules)) {colorhelp=ifelse(labelpred==modulename[i],colorcode[i],colorhelp)}; colorhelp=factor(colorhelp,levels=c(colorcode[1:no.modules],"grey")) }factor(colorhelp, levels=unique(colorhelp[hier1$order] ))}corFnc_i <- WGCNA::bicorclean_session <-function() {# Get all objects in the global environment objects <-ls(envir = .GlobalEnv)# Check if each object is a data frame or a matrix is_data <-vapply(objects, function(x) { obj <-get(x, envir = .GlobalEnv)is.data.frame(obj) ||is.matrix(obj) }, logical(1)) # Ensure logical vector output# Remove the identified objectsrm(list = objects[is_data], envir = .GlobalEnv)# Reset all plotswhile (dev.cur() >1) dev.off()cat_col("Session cleaned: data and plots removed.", "green")}add_prefix_except_plant_num <-function(type_of_date, df) { df_renamed <- df %>% dplyr::rename_with(~paste0(type_of_date, "_", .x),.cols =-plant_num )return(df_renamed)}correlation_summary <-function(df) { var_names <-colnames(df) results <-combn(var_names, 2, simplify =FALSE) |> purrr::map_df(function(vars) { x <- df[[vars[1]]] y <- df[[vars[2]]] ok <-complete.cases(x, y)if (sum(ok) >=3) { test <-cor.test(x[ok], y[ok], method ="pearson") tibble::tibble(V1 = vars[1],V2 = vars[2],corr =unname(test$estimate),pval = test$p.value,R2 = test$estimate^2 ) } else { tibble::tibble(V1 = vars[1],V2 = vars[2],corr =NA,pval =NA,R2 =NA ) } })# Add FDR correction results$fdr <-p.adjust(results$pval, method ="fdr")return(results)}```::: callout-tip## Helps for analyse- [Weighted Gene Co-expression Network Analysis (WGCNA)](https://moodle.france-bioinformatique.fr/pluginfile.php/703/mod_folder/content/0/WGCNA/WGCNA.html?forcedownload=1)- WGCNA: an R package for weighted correlation network analysis[@langfelder_wgcna_2008]:::## Weighted Gene Co-expression information- While WGCNA was initially developed for processing gene expression data, it can be adapted for other data types as well.- **Weighted Gene Co-expression Network Analysis (WGCNA)** is an approach dedicated to identifying **correlation** patterns between features of a dataset. In most cases, it is applied to gene expression data to find clusters -referred to as **modules**- of highly correlated genes. In transcriptomics, genes that exhibit strong correlations in their expression profiles are said **co-expressed**. One of the underlying principles of WGCNA is the concept of **guilt-by-association**. This principle suggests that genes sharing similar expression patterns are likely to be **functionally related** or involved in similar biological processes. By leveraging guilt-by-association through correlation, WGCNA provides a powerful tool for identifying potential functional relationships among genes and deciphering their roles in complex biological systems.- WGCNA builds a **co-expression** network where: - each feature (transcript/gene in transcriptomics data) is represented as a node - edges represent a correlation pattern between two features (nodes) and are weighted with the corresponding correlation coefficient, indicating the strength of their observed correlation.## Data importation ::: callout-warning## I keep only genes with variance > 0.1 across samples:::```{r, eval = F}load(file = here::here("data/microarray/output/df_RG_aquantil_log2.RData"))# Optional: filter out low-variability genesexpr_matrix <- df_RG_aquantil_log2[genefilter::rowVars(df_RG_aquantil_log2) > 0.1, ]# exportsave(expr_matrix, file = here::here("data/multi_omic/WGCNA/expr_matrix.RData"))```## Normalization ?**Normalization** of data according to the type of data and the type of analysis to perform. For RNA-seq data, WGCNA recommends to use the varianceStabilizingTransformation from the DESeq2 package, or to perform a log-transformation of normalized counts (RPFK/FPKM). I tested with and without normalization. May be it's better if i rerun my script with VST normlisation. Some person say it's RFPKM is a load methode (see reference below)```{r, eval=F}# dds <- DESeqDataSetFromMatrix(countData = expr_matrix, colData = col_data, design = ~ 1)# vst <- vst(dds)```Data must be **matrices** with specific shapes:- **samples** (e.g. patients, organisms ...) are displayed on the **rows**- **variables** (e.g. genes, proteins ...) are displayed on the **columns**```{r, eval=F}load(file = here::here("data/multi_omic/WGCNA/expr_matrix.RData"))datExpr <- t(expr_matrix)#datExpr <- t(counts(dds, normalized = TRUE)) # perform the median of ratios method #datExpr_vst <- t(assay(vst)) # with varianceStabilizingTransformation# datExpr_fpkm <- t(fpkm) # with varianceStabilizingTransformation# Run WGCNA quality check## WGCNA provides a function goodSamplesGenes to check the quality of input data and remove features and samples with too many missing data as well as genes with zero variance.gsg = goodSamplesGenes(datExpr, verbose=3) # link zero variance samples, for example here#gsg_vst = goodSamplesGenes(datExpr_vst, verbose=3) #gsg_fpkm = goodSamplesGenes(datExpr_fpkm, verbose=3) # link zero variance samples, for example here# tchek if everithing is okif (!gsg$allOK){ if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(datExpr)[!gsg$goodGenes], collapse = ", "))); if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(datExpr)[!gsg$goodSamples], collapse = ", "))); datExpr = datExpr[gsg$goodSamples, gsg$goodGenes]}# if (!gsg_vst$allOK){# if (sum(!gsg_vst$goodGenes)>0)# printFlush(paste("Removing genes:", paste(names(datExpr_vst)[!gsg_vst$goodGenes], collapse = ", ")));# if (sum(!gsg_vst$goodSamples)>0)# printFlush(paste("Removing samples:", paste(rownames(datExpr_vst)[!gsg_vst$goodSamples], collapse = ", ")));# datExpr = datExpr_vst[gsg_vst$goodSamples, gsg_vst$goodGenes]# }# # if (!gsg_fpkm$allOK){# if (sum(!gsg_fpkm$goodGenes)>0)# printFlush(paste("Removing genes:", paste(names(datExpr_fpkm)[!gsg_fpkm$goodGenes], collapse = ", ")));# if (sum(!gsg_vst$goodSamples)>0)# printFlush(paste("Removing samples:", paste(rownames(datExpr_fpkm)[!gsg_fpkm$goodSamples], collapse = ", ")));# datExpr = datExpr_fpkm[gsg_fpkm$goodSamples, gsg_fpkm$goodGenes]# }save(datExpr, file = here::here("data/multi_omic/WGCNA/datExpr.RData"))```**Data verification according to WGCNA criteria**```{r, eval=F}load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))sampleTree = hclust(dist(datExpr), method = "average")#sampleTree_vst = hclust(dist(datExpr_vst), method = "average")#sampleTree_fpkm = hclust(dist(datExpr_fpkm), method = "average")# Plot the sample tree: Open a graphic output window of size 12 by 9 inches# The user should change the dimensions if the window is too large or too small.png(here::here(paste0("report/multi_omics/plot/WGCNA/sample_clustering_to_detect_outliers.png")),width = 2500,height = 2000,res = 300)plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)dev.off()# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_vst.png")),width = 2500,height = 2000,res = 300)# plot(sampleTree_vst, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,# cex.axis = 1.5, cex.main = 2)# dev.off()# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_fpkm.png")),width = 2500,height = 2000,res = 300)# plot(sampleTree_fpkm, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,# cex.axis = 1.5, cex.main = 2)# dev.off()```**Metadata creation**```{r, eval=F}load(file = here::here("data/multi_omic/WGCNA/expr_matrix.RData"))load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))metadata <- colnames(expr_matrix) %>% as.data.frame() %>% dplyr::rename(sample_name = ".") %>% separate(sample_name, into = c("Mutant_type", "sulfur_condition", "sample_num", "Rep", "chip_color", "chip_num"), sep = "_", remove = FALSE) %>% mutate(genotype = ifelse(sample_num %in% c(1.1, 2.2, 3.3, 4.4, 17.1, 18.2, 19.3, 20.4), "WT2", sample_num),# Add genotype genotype = ifelse(sample_num %in% c(5.1, 6.2, 7.3, 8.4, 21.1, 22.2, 23.3, 24.4), "E568K", genotype),# Add genotype genotype = ifelse(sample_num %in% c(9.5, 10.6, 11.7, 12.8, 25.5, 26.6, 27.7, 28.8), "WT1", genotype),# Add genotype genotype = ifelse(sample_num %in% c(13.5, 14.6, 15.7, 16.8, 29.5, 30.6, 31.7, 32.8), "W78*", genotype), genotype = as.factor(genotype), genotype=fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")) %>% column_to_rownames("sample_name")metadata = metadata[rownames(datExpr),] %>% # We keep only samples that are present in our dataset dplyr::select(Mutant_type, sulfur_condition ,chip_color ,chip_num, genotype)save(metadata, file = here::here("data/multi_omic/WGCNA/metadata.RData"))```**Let's replot our sample dendrogram and visualize the metadata associated to each sample.**```{r, eval=F}load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))traitColors<-metadata %>% mutate(genotype = as.factor(genotype)) %>% mutate(Mutant_type=labels2colors(ifelse(Mutant_type=="Mut",0,1)), chip_num = labels2colors(as.numeric(chip_num)), genotype=mutant_palette[genotype], sulfur_condition=ifelse(sulfur_condition=="SS",sulfate_pallet[1],sulfate_pallet[2]), chip_color =ifelse(chip_color =="Green","green","red") )samplesTree <- hclust(d = dist(datExpr), method = "average")png(here::here(paste0("report/multi_omics/plot/WGCNA/sample_clustering_to_detect_outliers_color.png")),width = 2500,height = 2000,res = 300)plotDendroAndColors(dendro = samplesTree, colors = traitColors, groupLabels = names(metadata), main = "Samples clustering using hierarchical clustering method with \n metadata information and counts divided by sample-specific size factors \n determined by median ratio of gene counts relative to geometric mean per gene",cex.dendroLabels = 0.5, cex.colorLabels = 0.5) dev.off()# # WT1 W78* WT2 E568K # "#1B3764" "#931B1E" "#00A4CA" "#D91F26"# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_vst_color.png")),width = 2500,height = 2000,res = 300)# samplesTree <- hclust(d = dist(datExpr_vst), method = "average")# plotDendroAndColors(dendro = samplesTree, colors = traitColors, groupLabels = names(metadata), # main = "Samples clustering using hierarchical clustering \n method with metadata information and vst",cex.dendroLabels = 0.5, cex.colorLabels = 0.5) # dev.off()# # png(here::here(paste0("report/rnaseq/plot/WGCNA/sample_clustering_to_detect_outliers_fpkm_color.png")),width = 2500,height = 2000,res = 300)# samplesTree <- hclust(d = dist(datExpr_fpkm), method = "average")# plotDendroAndColors(dendro = samplesTree, colors = traitColors, groupLabels = names(metadata), # main = "Samples clustering using hierarchical clustering \n method with metadata information and fpkm",cex.dendroLabels = 0.5, cex.colorLabels = 0.5) # dev.off()```{fig-align="center"}**Saving genes after verification**```{r, eval=F}save(datExpr, file =here::here("data/multi_omic/WGCNA/input_datExpr.RData"))#save(datExpr_vst, file =here::here("data/rnaseq/output/WGCNA/input_vst.RData"))#save(datExpr_fpkm, file =here::here("data/rnaseq/output/WGCNA/input_fpkm.RData"))```## Construction of the correlation network**Then you can use Pearson correlation or Spearman correlation.**The Pearson correlation measures **linear** associations between variables, whereas Spearman correlation measures **monotonous** associations between variables. Both can have interesting results, so it\'s completely up to you! By default, WGCNA uses the Pearson correlation. If you want to use the Spearman correlation, don\'t forget to specify it each time you notice a **method** parameter controlling the correlation method in one of the WGCNA functions!### Computing the correlation network (summary)This first step of the analysis consists in creating the **feature correlation network**. This first step is composed of tree main substeps:1. **Power** selection - **Soft thresholding** is a strategy introduced by WGCNA to transform correlation coefficients for building a weighted correlation network - The threshold should be set such that it makes the network **scale-free**2. **Similarity matrix** construction - First, **correlation** are calculated between each **pair of feature** - Then, correlation or distance are **transformed**, using the selected soft-power3. **Topological overlap matrix** construction (TOM) - The similarity matrix is transformed into a TOM to facilitate the module detection### Azrael code ::: callout-caution## Script use in server AzraelYou'll probably be limited by the amount of RAM required to create the TOM. It's best to use the following script on a server and then re-import the data.:::```{r,eval=F}# Function of need in memoryBaseColors = c("turquoise","blue","brown","yellow","green","red","black","pink","magenta", "purple","greenyellow","tan","salmon","cyan", "midnightblue", "lightcyan","grey60", "lightgreen", "lightyellow", "royalblue", "darkred", "darkgreen", "darkturquoise", "darkgrey", "orange", "darkorange", "white", "skyblue", "saddlebrown", "steelblue", "paleturquoise", "violet", "darkolivegreen", "darkmagenta", "white" );RColors = colors()[-grep("grey", colors())];RColors = RColors[-grep("gray", RColors)];InBase = match(BaseColors, RColors);ExtraColors = RColors[-c(InBase[!is.na(InBase)])];No.Extras = length(ExtraColors);# Here is the vector of colors that should be used by all functions:GlobalStandardColors = c(BaseColors, ExtraColors[rank(sin(13*c(1:No.Extras) +sin(13*c(1:No.Extras))) )] );rm(BaseColors, RColors, ExtraColors, No.Extras);############## Package and function (end) ################# Inputinput_name<-"datExpr" # based on previous script# Options WGCNAoptions(stringsAsFactors = FALSE);enableWGCNAThreads(8)# memory.limit(size=10000000)# Network parametersnetworkType_i <- "signed" #signed #unsigned# Module detection parameterminModuleSize <- 30ch = 0.99 # for cutheight# Load WGCNA ready input dataload(file =paste0("data/multi_omic/WGCNA/input_",input_name,".RData")) # data previously created# Call the network topology analysis function.# You might want to change de minimal R^2 and the probed power value according# to your network topology.RsquaredCut = 0.7powers = c(c(1:11), seq(from = 12, to=40, by=2))cat("Begin of the script \n")# Correlation networks are by nature complete, which means every feature is associated to every other with a weight given by the correlation coefficient. To extract only meaningful associations, the correlation network must be thresholded.# # The most straight-forward way of thresholding a correlation dataset is to perform a hard-thresholding, i.e. two features are considered correlated if their absolute correlation coefficient exceeds a user-defined threshold.# # WGCNA propose another strategy, called soft-thresholding. It consists of raising the correlation coefficients to a certain power. This process transforms low correlation coefficients to values close to 0, while high correlation coefficients are transformed to values close to 1.# # Furthermore, WGCNA presents a function called pickSoftThreshold that aims to automatically determine the ideal power for your dataset. The optimal power is defined as the one that ensures your network exhibits a scale-free property. Let’s use it!# # First, we need to define a range of powers to try.cat("Power selection")# Sometimes, WGCNA can’t find the optimal soft-threshold. In such cases, the variable sft$powerEstimate contains NA. If the optimal soft-power was estimated but seems low, you might want to increase it. In both cases, you can plot the following figures to help you decide:# sft <- pickSoftThreshold(dataset, powerVector = powers, verbose = 5) # or allowWGCNAThreads(nThreads = 12) # adapte à ton système # Sur mon macsft = pickSoftThreshold(get(input_name), #takes normalized data, samples data, creates an agacency matrix and compares with different agacency matrices RsquaredCut = RsquaredCut, powerVector = powers, corFnc = corFnc_i, networkType = networkType_i, verbose = 5)# Which power should be used for your dataset? (don’t forget to look at the documentation of the pickSoftThreshold fucntion to known how).sprintf("Optimal soft-power = %d", sft$powerEstimate)# "Optimal soft-power = 11"# The optimal soft-power to use is defined by the topology of the soft-thresholded network, which has to be scale-free! Scale-free networks are networks in which the majority of nodes have a low degree, and a small number of nodes have a large degree. Multiple studies have shown that most biological networks are scale-free!# # Let’s visualize the difference between the raw correlations and the soft-thresholded correlations:corRaw = abs(cor(get(paste0(input_name))))corSoft = abs(corRaw**sft$powerEstimate)#export soft tresholdingpng(here::here(paste0("report/multi_omics/plot/WGCNA/","soft_tresholding_",input_name,"_",networkType_i,".png")),width = 800, height = 600)par(mfrow=c(1,2))hist(corRaw, main="Raw correlations")hist(corSoft, main="Soft-thresholded correlations")dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/", "soft_threshold_indices.png")),width = 800, height = 600)par(mfrow = c(1, 2)) plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)", ylab = paste("Scale Free Topology Model Fit,", networkType_i, "R^2"), type="n", main = paste("Scale independence", networkType_i));text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,col="red")top = 0.9; mid = 0.85; low = 0.8; offset = 0.02abline(h=top,col="green")text(y=top + offset, x=offset+0.5, col="green", "0.90")abline(h=mid, col="orange")text(y=mid + offset, x=offset+0.5, col="orange", "0.85")abline(h=low,col="red")text(y=low + offset, x=offset+0.5, col="red", "0.80")plot(sft$fitIndices[,1], sft$fitIndices[,5], xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n", main = paste("Mean connectivity (", networkType_i, " correlation)", sep = ""))text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers,col="red")dev.off()# The Scale Independence plot (left) shows the fit indices for scale free topology for each power. The authors recommend to select a power with a R² > 0.8. The red line corresponds to using an R² cut-off of R²=0.80.# # The Mean Connectivity plot (right) shows the average of connectivity for each power. A balance needs to be identified between a fully connected network and too much disconnected edges.cat("Adjacency calculation (Similarity matrix construction)\n")# WGCNA provides the function adjacency to compute a features similarity matrix. This function operates in three steps:# # Compute features correlation coefficients# Transform the correlation coefficients depending on whether you are interested in absolute correlations (unsigned), signed correlations (signed) or positive correlations (signed hybrid). To have more details on this 3 types, check ?adjacency. By default, the constructed network is an unsigned network which means that sign are not reported into the network (absolute correlations).# Finally, the transformed correlations are raised using the previously estimated soft-threshold power.# Adjacency calculationsoftPower <- sft$powerEstimate #betaadjacency <- adjacency(datExpr = get(paste0(input_name)) , #normalise matrix type = networkType_i, power = softPower, #beta corFnc = corFnc_i)head(adjacency[, c(1:5)])cat("Topological overlap matrix \n")# From the Similarity matrix, WGCNA proposes to extract correlated feature modules (co-expressed gene modules) based on Topological Overlap. The idea is that two features (genes) that share a high number of correlated (co-expressed) neighbors are likely to participate in the same correlation (co-expression) module.# # The Topological Overlap Matrix (TOM) is a concept used in network analysis to quantify the similarity or interconnectedness between nodes in a network.# # Transform the similarity matrix into a Topological Overlap Matrix that incorporates the level of overlapping between the neighborhoods of two nodes (features/genes), using the function TOMsimilarity. Notice that the output does not contains our row names and column names (i.e. node names). Make sure to rename it!# calculating topological superpositionTOM = TOMsimilarity(adjacency) # Because most clustering strategies are performed on distance (rather than similarities), we transform the TOM matrix into a dissimilarity matrix.dissTOM = 1-TOMhead(dissTOM[, c(1:5)])cat("save disstom \n")save(dissTOM, file = here::here(paste0("data/multi_omic/",input_name,"_",networkType_i,"_dissTOM.RData)"))) #/!\ could be heavy################### Module detection #####################cat("Module detection \n")# Modules are groups of highly-connected nodes (features). This step is composed of two substeps:# Clustering: getting correlated feature modules# Use the hclust function to perform hierarchical clustering# Display the dendrogram# Use the Dynamic Tree Cut function to obtain the modules# Merging of modules: merging close modules# If they are too similar, modules should be merged# We compute module-module similarities by computing their eigenvectors (eigengene)# We merge two modules if their eigenvectors are correlated (the authors suggest a pearson corr r > 0.75)# export png genetree#png(file= paste0(input_name,"_",networkType_i,"_dissTOM_clustering-formation.png"),width = 3840 , height = 3840 ,res=300)# Call the hierarchical clustering functiongeneTree = hclust(as.dist(dissTOM), method = "average")# Plot the resulting clustering tree (dendrogram)plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)dev.off()# The leaves represent the features (genes). The dendrogram branches group together densely interconnected and highly correlated features.# To identify modules from dendrograms, several methods exist. Here, we use the cutreeDynamic function.# First, we have to set the minimum module size minClusterSize (i.e. minimum number of genes inside a group/module). The default is 20, but you can try several values and choose according to you results. There are other parameters you can play with, such as deepSplit. Look at the documentation to select the parameters.########### from tutorial #############dynamicMods <- cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)table(dynamicMods)dynamicColors <- labels2colors(dynamicMods)table(dynamicColors)png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width = 3840*2 , height = 3840 ,res=600)plotDendroAndColors(dendro = geneTree, colors = dynamicColors, groupLabels = "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors from tutorial")dev.off()# Clusturisation . there are lots of different ones. the best thing to do is to take the hybrid like everyone else.png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width = 3840*2 , height = 3840 ,res=600)colorh1=as.character(static_colors(geneTree,h1=ch,minsize1=minModuleSize))DynamicColor1 = labels2colors(cutreeDynamic(geneTree, cutHeight = ch, minClusterSize = minModuleSize, method = "tree", deepSplit = TRUE)); # deepSplit, regulates module coherence.# if we want to see what's inside table(DynamicColor1)DynamicColor2 = labels2colors(cutreeDynamic(geneTree, cutHeight = ch, minClusterSize = minModuleSize, method = "tree", deepSplit = FALSE));ClustColor1 = labels2colors(cutreeDynamic(dendro = geneTree, minClusterSize = minModuleSize, cutHeight = ch, method = "hybrid", maxCoreScatter = 0.75, minGap = 0.25, pamStage = TRUE, distM = dissTOM, useMedoids = FALSE, maxPamDist = 0.9, respectSmallClusters = TRUE));Clusters1 = cutreeHybrid(dendro = geneTree, minClusterSize = minModuleSize, cutHeight = ch, maxCoreScatter = 0.75, minGap = 0.250, pamStage = TRUE, distM = dissTOM, useMedoids = FALSE, maxPamDist = 0.9, respectSmallClusters = TRUE)CoreColor1 = labels2colors(Clusters1$cores)ClustColor2 = labels2colors(cutreeDynamic(dendro = geneTree, minClusterSize = minModuleSize, cutHeight = ch, method = "hybrid", maxCoreScatter = 0.95, minGap = 0.050, pamStage = TRUE, distM = dissTOM, useMedoids = FALSE, maxPamDist = 0.90, respectSmallClusters = TRUE));Clusters2 = cutreeHybrid(dendro = geneTree, minClusterSize = minModuleSize, cutHeight = ch, maxCoreScatter = 0.95, minGap = 0.050, pamStage = TRUE, distM = dissTOM, useMedoids = FALSE, maxPamDist = 0.90, respectSmallClusters = TRUE)CoreColor2 = labels2colors(Clusters2$cores)AutoColor = NULL;for (deepSplit in 0:3){ AutoColor = cbind(AutoColor, labels2colors(cutreeDynamic(dendro = geneTree, minClusterSize = 30 - 3*deepSplit, cutHeight = ch, method = "hybrid", deepSplit = deepSplit, distM = dissTOM)))AutoLabels = paste("Hybrid 'auto': dS =", c(0:3))par(mfrow=c(2,1))par(cex = 1.4);par(mar = c(0,8.5,2,0));plot(geneTree,labels=F,main="Hierarchical dendrogram and module colors", sub="", xlab="")par(mar = c(1,8.5,0,0));plotHclustColors(geneTree, cbind(colorh1, DynamicColor1, DynamicColor2, ClustColor1, CoreColor1, ClustColor2, CoreColor2, AutoColor), c("Static", "Dynamic Tree (dS)", "Dynamic Tree (No dS)", "Dynamic Hybrid 1", "Cores in DHyb 1", "Dynamic Hybrid 2", "Cores in DHyb 2", AutoLabels), main = "")dev.off()}png(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_hybride_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,".png")),width = 3840*2 , height = 3840 ,res=600)# Same thing just with the hybrid function but with different settings# hybrid cutmColorh = NULLlabels = NULLpar(mfrow=c(2,1))par(cex = 1.4);par(mar = c(0,5.5,2,0));plot(geneTree,labels=F,main="Hierarchical dendrogram and module colors with ", sub="", xlab="")par(mar = c(1,5.5,0,0));for (ds in 0:3) for (mc in c(16,32,64)) { tree = cutreeHybrid(dendro = geneTree, pamStage=FALSE, minClusterSize = mc, #cutHeight = ch, deepSplit = ds, distM = dissTOM) mColorh=cbind(mColorh,labels2colors(tree$labels)) labels = cbind(labels,paste("dS =",ds,"; minSize =",mc))}#plotHclustColors(geneTree, staticColors, "tric", cex.rowLabels = 0.6)plotHclustColors(geneTree,mColorh,labels, cex.rowLabels = 0.6)dev.off()deepSplit = 2 # The higher the parameters, the more coherent the modules created.mc=minModuleSizepng(here::here(paste0("report/multi_omics/plot/WGCNA/module_detection_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width = 3840*2 , height = 3840 ,res=600)tree = cutreeHybrid(dendro = geneTree, distM = dissTOM, #cutHeight = ch, deepSplit = deepSplit, minClusterSize = minModuleSize, pamStage=FALSE)dynamicMods <- tree$labels# Convert numeric lables into colorstable(dynamicMods)dynamicColors = labels2colors(dynamicMods)table(dynamicColors)write.csv(x=table(dynamicColors),file = here::here(paste0("data/multi_omic/WGCNA/module_detection_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))# Plot the dendrogram and colors underneathpar(mfrow = c(2, 1))plotDendroAndColors(geneTree, dynamicColors, "Tree Cut Hybrid", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors selected")dev.off()cat("Merging modules \n")# Sometimes, especially when dealing with smaller clusters (for instance when the minClusterSize parameter is not set properly), merging modules with similar expression profiles might be usefull.# To do so, we first compute each module eigenvectors. A module eigenvector (or eigengene if you are dealing with gene data) represents the mean expression profile of the whole corresponding module. We can use the correlation of modules eigengenes to define a module-similarity metric that we can use to merge close modules.# Merging similar modules# Module-specific genesMEList = WGCNA::moduleEigengenes(get(paste0(input_name)), colors = dynamicColors)MEs = MEList$eigengenes # sample expression profile for each prob gene# Now let’s perform a hierarchical clustering using the correlation between module eigenvectors. Remember that for using hclust we need distances! Compute those distances and transform then with as.dist before running hclust.MEDiss <- 1-cor(MEs)# merge certain modules# Classification of module genesMETree = hclust(as.dist(MEDiss), method = "average")# WGCNA’s authors recommend to merge cluster with a correlation greater than 0.75, hence a distance lower than 0.25. Let’s look at the tree with plotand abline, just as we did previously in this tutorial.MEDissThres = 0.15 # Disimilarite max thresholdpng(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width = 3840*2 , height = 3840 ,res=600)# This classification takes into account the type of networkplot(METree, main = paste("Clustering of module eigengenes pre-fusion (",networkType_i, " correlation)", sep = ""), xlab = "", sub = "") # disimilarite tree# Similarity threshold for module mergingabline(h=MEDissThres, col = "red")dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select_verif_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width = 3840*2 , height = 3840 ,res=600)plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle = 90)dev.off()png(here::here(paste0("report/multi_omics/plot/WGCNA/module_merged_select2_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width = 3840*2 , height = 3840 ,res=600)merge = WGCNA::mergeCloseModules(get(paste0(input_name)), dynamicColors, cutHeight = MEDissThres, verbose = 3)mergedColors = merge$colors;mergedMEs = merge$newMEs;plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = paste("Merged clusters\nnetwork type = ", networkType_i,"\nmerging distance = ", MEDissThres, sep = ""))dev.off()print("rename module")# Rename to moduleColorsmoduleColors = mergedColors# Construct numerical labels corresponding to the colorscolorOrder = c("grey", standardColors(50));moduleLabels = match(moduleColors, colorOrder)-1;MEs = mergedMEs;# Let’s save those clustering results! Choose your favorite export format (Rdata, csv, tsv, …) and save the results. Don’t forget to make sure you have all the information you need (in this case, we need the node names and their cluster labels).moduleColors = mergedColors # Or mergedColors, if some modules have been mergednames(moduleColors) = colnames(get(paste0(input_name)))head(moduleColors)unique(moduleColors)write.csv(x=moduleColors,file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))test = as.data.frame(moduleColors) %>% rownames_to_column("gene_id") %>% filter(gene_id == "PsCam042688")# Save module colors and labels for use in subsequent partssave(MEs, file = here::here(paste0("data/multi_omic/WGCNA/","modules_mrna_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".rda")))print("tomplot")# additional graphs# Transform dissTOM with a power to make moderately strong connections more visible in the heatmappng(file = here::here(paste0("report/multi_omics/plot/WGCNA/","plotTOM_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".png")),width = 3840 , height = 3840 ,res=600)plotTOM = TOM^0.6;# Set diagonal to NA for a nicer plotdiag(plotTOM) = NA#sizeGrWindow(9,9)TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")dev.off()# export importante valuesave(mergedMEs, file = here::here("data/multi_omic/WGCNA/mergedMEs.RData"))```### Power selectionCorrelation networks are by nature complete, which means every feature is associated to every other with a weight given by the correlation coefficient. To extract only meaningful associations, the correlation network must be thresholded.The most straight-forward way of thresholding a correlation dataset is to perform a **hard-thresholding**, i.e. two features are considered correlated if their absolute correlation coefficient exceeds a user-defined threshold.WGCNA propose another strategy, called **soft-thresholding**. It consists of raising the correlation coefficients to a certain **power**. This process transforms low correlation coefficients to values close to 0, while high correlation coefficients are transformed to values close to 1.Furthermore, WGCNA presents a function called `pickSoftThreshold` that aims to automatically determine the ideal power for your dataset. The optimal power is defined as the one that ensures your network exhibits a **scale-free** property.`"Optimal soft-power = 10"`{fig-align="center"}{fig-align="center"}{fig-align="center"}{fig-align="center"}## Show profile of each module```{r, eval = F}load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <- c("gene_id", "moduleColors")# specificity of the different modulesummarise_color <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(n_genes = n()) %>% arrange(desc(n_genes))# a. passer datExpr au format longexpr_long <- as.data.frame(datExpr) %>% rownames_to_column(var = "sample") %>% # garder l'ID échantillon pivot_longer(-sample, names_to = "gene_id", values_to = "expression")# b. ajouter la couleur de module pour chaque gèneexpr_long <- expr_long %>% left_join(., moduleColors, by = "gene_id") %>% left_join(., metadata %>% rownames_to_column("sample"), by = "sample") %>% mutate(condition = paste0(sulfur_condition, "_", genotype))## ─────────────────────────────## 2. Moyenne d’expression par échantillon & par module## ─────────────────────────────expr_module_sample <- expr_long %>% dplyr::group_by(moduleColors, sample, condition) %>% dplyr::summarise( mean_expr = mean(expression, na.rm = TRUE), # moyenne d’expression n_genes = dplyr::n(), # nombre total de lignes (= gènes) # si tu préfères compter les gènes uniques : # n_genes = dplyr::n_distinct(gene_id), .groups = "drop" ) %>% dplyr::mutate( moduleColors_n = paste0(str_to_title(moduleColors), " (", n_genes, ")"), moduleColors_n = forcats::fct_reorder(moduleColors_n, n_genes, .fun = max, .desc = TRUE) ) %>% mutate(condition=fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")) strip_df <- expr_module_sample %>% distinct(moduleColors) %>% # une ligne par module arrange(moduleColors) %>% mutate( strip_fill = as.character(moduleColors) )## On récupère les vecteurs (ordre = niveaux du facteur)strip_fill_vec <- strip_df$strip_fillnames(strip_fill_vec) <- strip_df$moduleColorsstrip_text_vec <- strip_df %>% mutate(couleur_texte_perso = sapply(strip_fill, evaluate_contrast)) %>% pull(couleur_texte_perso)names(strip_text_vec) <- strip_df$moduleColorsp1 <- expr_module_sample %>% mutate(genotype = sub("^[^_]+_", "", condition), sulfur_condition = sub("_.*$", "", condition), genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), sulfur_condition = fct_relevel(sulfur_condition, "SS", "SD") ) %>% ggplot(., aes(x = condition, y = mean_expr, fill = genotype, col = genotype)) + geom_boxplot(outlier.shape = NA, alpha = .4) + geom_jitter(width = .25, size = .6, alpha = .4) + # points individuel labs(x = "Module", y = expression("Mean expression level (log"["2"]*")"), colour = "Genotype", fill = "Genotype") + scale_color_manual(values = mutant_palette)+ scale_fill_manual(values = mutant_palette)+ # --> Ici tu peux ajouter new_scale_fill() si tu superposes d'autres couches ggh4x::facet_wrap2( ~ moduleColors_n, scales = "free_y", nrow = 7, # ajuste le layout strip = ggh4x::strip_themed( background_x = ggh4x::elem_list_rect(fill = strip_fill_vec), text_x = ggh4x::elem_list_text(colour = strip_text_vec) ) ) + theme_bw() + theme( panel.spacing = unit(0.5, "lines"), axis.ticks.y = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), legend.text = element_text(size = 12), legend.title = element_text(size = 12), legend.key.size = unit(10, "pt"), plot.caption = element_text(size = 8), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) ;p1fig_export(path = "report/multi_omics/plot/WGCNA/average_module", plot_x = p1, height_i = 9, width_i = 12, res = 600)expr_module_sample$condition# i think i need to use patchwork and change a litte bit my function stats. v_moduleColors = levels(as.factor(expr_module_sample$moduleColors_n))# add color to the table expr_module_sampleexpr_module_sample = expr_module_sample %>% mutate(couleur_texte_perso = sapply(moduleColors, evaluate_contrast)) # i = 4# moduleColors_i <- v_moduleColors[i]# # test <- stat_analyse(# data=expr_module_sample %>% # as.data.frame() %>% # filter(moduleColors_n == moduleColors_i) %>% # mutate(condition = as.factor(condition)),# # column_value = "mean_expr",# category_variables = "condition",# grp_var = "",# show_plot = T,# outlier_show = F, # label_outlier = "sample",# biologist_stats = T,# Ylab_i = "Chlorophyll content for SD and consistent pod \n number condition in N4_N5 at 40 DAP",# control_conditions = "",# strip_normale = T,# hex_pallet = as.character(c(mutant_palette, mutant_palette)), # strip = "moduleColors_n", #colone for texte# strip_fill_vec = "moduleColors", # color of the square # strip_text_vec = "couleur_texte_perso"# )# # test[["plot"]]######################### plots_list <- lapply(seq_along(v_moduleColors), function(i) { moduleColors_i <- v_moduleColors[i] p <- stat_analyse( data = expr_module_sample %>% as.data.frame() %>% filter(moduleColors_n == moduleColors_i) %>% mutate(condition = as.factor(condition)) %>% filter(condition %in% c("SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")), # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! column_value = "mean_expr", category_variables = "condition", grp_var = "", show_plot = TRUE, outlier_show = FALSE, label_outlier = "sample", biologist_stats = TRUE, Ylab_i = paste0( "" ), control_conditions = "", strip_normale = TRUE, hex_pallet = as.character(rep(mutant_palette, 2)), strip = "moduleColors_n", strip_fill_vec = "moduleColors", strip_text_vec = "couleur_texte_perso" )[["plot"]]+ theme_bw()+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) ## 1-a) On retire les titres d’axe des sous-graphiques p <- p + labs(x = NULL, y = NULL) ## 1-b) Option : on masque aussi les graduations redondantes # (ici : on garde y seulement sur la 1re colonne # et x seulement sur la dernière ligne) ncol_layout <- 5 nrow_layout <- ceiling(length(v_moduleColors) / ncol_layout) col_i <- ((i - 1) %% ncol_layout) + 1 row_i <- ceiling(i / ncol_layout) # if (col_i != 1) { # p # } if (row_i != nrow_layout) { p <- p + theme(axis.text.x = element_blank()) } p +labs(fill = "Genotype", colour = "Genotype")})# ------------------------------------------------------------------# 2) Mosaïque principale -------------------------------------------# ------------------------------------------------------------------panel <- wrap_plots(plots_list, ncol = 5, guides = "collect") & theme(legend.position = "right", panel.grid.major = element_blank(), panel.grid.minor = element_blank() )## ------------## 2. Labels## ------------y_lab <- ggplot() + labs(y = "Mean expression") + theme_void() + theme(axis.title.y = element_text(angle = 90, vjust = 0.5, hjust = 0.5))x_lab <- ggplot() + labs(x = "Condition") + theme_void() + theme(axis.title.x = element_text(vjust = -0.5, hjust = 0.5))## ------------## 3. Assemblage 2 × 2## ------------final_plot_1 <- (y_lab | panel)+ plot_layout( widths = c(0.02, 0.98), # 5 % pour la colonne Y, 95 % pour le panel guides = "collect" )final_plot_2 <- (plot_spacer() | x_lab)+ plot_layout( widths = c(0.02, 0.98), # 5 % pour la colonne Y, 95 % pour le panel guides = "collect" )final_plot <- (final_plot_1/final_plot_2)+ plot_layout( heights = c(0.99, 0.01), # 7 % pour le label X guides = "collect" )final_plot#fig_export(path = "report/multi_omics/plot/WGCNA/average_module_stats", plot_x = final_plot, height_i = 9, width_i = 18, res = 600)fig_export(path = "report/multi_omics/plot/WGCNA/average_module_stats_SD", plot_x = final_plot, height_i = 9, width_i = 10, res = 600)```{fig-align="center"}{fig-align="center"}## Correlation with other data sets ```{r, eval = F}load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))load(file = here::here("data/multi_omic/WGCNA/mergedMEs.RData"))df_plant_info_2014_filter <- read_excel(here::here("data/plant_info_XP2014.xlsx")) %>% drop_na(code_run_transcripto) %>% dplyr::select(plant_num, code_run_transcripto, code_metabo_leaf_repro,code_metabo_leaf_vege, sulfate) %>% dplyr::rename(code_run_anion = sulfate)# load metabolomic resultsdf_metabo_global <- read_csv(here::here("data/metabolomic/output/df_metabolite_global_IBMP.csv"), show_col_types = FALSE) %>% mutate(sulfur_condition= fct_relevel(sulfur_condition, "SS", "SD"), genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), type_genotype = fct_relevel(type_genotype, "WT", "Mut") ) %>% filter(plant_num %in% df_plant_info_2014_filter$plant_num, compartment == "vegetative leaf") %>% #() "reproductive leaf" "vegetative leaf" ) dplyr::select(variable, plant_num, value) %>% pivot_wider(names_from = variable, values_from = value)vect_metabo = colnames(df_metabo_global)[-1]################################## for more simplicity change here vect_metabo = vect_metabo[1:100]plant_measured <- df_metabo_global %>% drop_na(glucose) %>% pull(plant_num)metadata_transcripto_metabo_filter = metadata %>% rownames_to_column("sample_id") %>% mutate(sample_num = str_match(sample_id, "(Mut|WT)_[A-Z]+_(\\d+\\.\\d+)")[,3], Rep = paste0("Rep", str_match(sample_id, "_Rep(\\d+)")[,2])) %>% mutate(code_run_transcripto = gsub("\\.", "_", as.character(sample_num))) %>% left_join(., df_plant_info_2014_filter, by = "code_run_transcripto") %>% left_join(., df_metabo_global, by = "plant_num") %>% mutate(sample_name = paste(sep = "_", Mutant_type, sulfur_condition, sample_num, Rep, chip_color, chip_num)) %>% column_to_rownames("sample_name") %>% filter(plant_num %in% plant_measured) %>% dplyr::select(all_of(vect_metabo)) %>% dplyr::select(where(~ any(. != 0)))# Add anionmetadata_transcripto_metabo_filterlist_sample <- rownames(metadata_transcripto_metabo_filter)mergedMEs_filter <- mergedMEs %>% rownames_to_column("sample_id") %>% filter(sample_id %in% list_sample) %>% column_to_rownames("sample_id")# filter one data i have for boath datasetstopifnot(identical(rownames(mergedMEs_filter), rownames(metadata_transcripto_metabo_filter)))moduleTraitCor <- cor(mergedMEs_filter, metadata_transcripto_metabo_filter, use = "pairwise.complete.obs", method = "pearson")moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples = nrow(mergedMEs)) # WGCNA fonctiontextMatrix <- paste0(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")")dim(textMatrix) <- dim(moduleTraitCor)colnames(textMatrix) <- colnames(moduleTraitCor)rownames(textMatrix) <- rownames(moduleTraitCor)# 4. Masquer les r non significatifs (p > 0.05)for (r in rownames(moduleTraitCor)) { for (c in colnames(moduleTraitCor)) { if (moduleTraitPvalue[r, c] > 0.05) { moduleTraitCor[r, c] <- 0 textMatrix[r, c] <- "" } }}# 1. Récupérer noms + couleurs des modules -------------------------------module_names <- rownames(moduleTraitCor) # ex. "MEturquoise"module_colors <- sub("^ME", "", module_names) # ex. "turquoise"# 2. Data-frame d’annotation (une ligne = un module) ----------------------annotation_row <- data.frame(Module = factor(module_colors, levels = unique(module_colors)))rownames(annotation_row) <- module_names# 3. Palette : liste nommée (niveau du facteur -> couleur) ---------------annotation_colors <- list( Module = setNames(unique(module_colors), # valeurs uniques seulement unique(module_colors)))# Si tu veux des codes hex personnalisés, remplace la ligne ci-dessus par :# annotation_colors <- list(# Module = c(turquoise = "#00CED1", brown = "#8B4513", ...)# )# 4. Heatmap --------------------------------------------------------------my_col <- colorRampPalette(c("darkgreen", "white", "darkmagenta"))(50)png(here::here("report/multi_omics/plot/WGCNA/module_metabo_heatmap_clustered_test.png"), width = 6000, height = 4000, res = 300)pheatmap( mat = moduleTraitCor, cluster_rows = TRUE, cluster_cols = TRUE, clustering_distance_rows = "euclidean", clustering_distance_cols = "euclidean", clustering_method = "complete", color = my_col, display_numbers = textMatrix, number_color = "black", fontsize_number = 5, legend_breaks = seq(-1, 1, by = 0.5), border_color = NA, main = "Module–metabolite relationships (clustered)", annotation_row = annotation_row, # ← la barre de couleur à gauche annotation_colors = annotation_colors, # ← palette correspondante legend = FALSE, # ← supprime la palette continue annotation_legend = FALSE # ← supprime la légende des modules)dev.off()```{fig-align="center"}## Find all correlation (with model mixte)```{r, eval = F}load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))load(file = here::here("data/multi_omic/WGCNA/mergedMEs.RData"))# data importation #####df_plant_info_2014 <- read_excel(here::here("data/plant_info_XP2014.xlsx")) %>% dplyr::select(plant_num, sulfur_condition, genotype, code_run_transcripto,code_metabo_leaf_vege, sulfate) %>% dplyr::rename(code_run_anion = sulfate)df_anion <- read_csv(file = here::here("data/metabolomic/output/anion_EVA.csv"), show_col_types = F) %>% dplyr::select(number_tube, anion,compartment, concentration) %>% filter(compartment == "vegetative_leaves") %>% pivot_wider(names_from = anion, values_from = concentration) %>% dplyr::rename(code_run_anion = number_tube) %>% left_join(., df_plant_info_2014 %>% dplyr::select(plant_num, code_run_anion), by = "code_run_anion") %>% dplyr::select(-c(code_run_anion, compartment))df_metabo = read_csv(here::here("data/metabolomic/output/df_metabolite_global_IBMP.csv"), show_col_types = FALSE) %>% filter(compartment == "vegetative leaf") %>% dplyr::select(variable, plant_num, value) %>% pivot_wider(names_from = variable, values_from = value) %>% dplyr::select(where(~ any(. != 0)))df_metadata <- metadata %>% rownames_to_column("sample_id") %>% mutate(sample_num = str_match(sample_id, "(Mut|WT)_[A-Z]+_(\\d+\\.\\d+)")[,3], Rep = paste0("Rep", str_match(sample_id, "_Rep(\\d+)")[,2])) %>% mutate(code_run_transcripto = gsub("\\.", "_", as.character(sample_num))) %>% left_join(., df_plant_info_2014, by = "code_run_transcripto") %>% dplyr::select(sample_id, plant_num)mergedMEs_input<-mergedMEs %>% rownames_to_column("sample_id") %>% left_join(., df_metadata, by = "sample_id") %>% dplyr::select(-sample_id) %>% dplyr::group_by(plant_num) %>% dplyr::summarise(across(everything(), mean))# add type of data before analyse and mergedf_global <- df_plant_info_2014 %>% dplyr::select(plant_num, sulfur_condition, genotype) %>% full_join(., add_prefix_except_plant_num("ANION", df_anion), by = "plant_num") %>% full_join(., add_prefix_except_plant_num("METABO", df_metabo), by = "plant_num") %>% full_join(., add_prefix_except_plant_num("MODULE", mergedMEs_input %>% dplyr::rename_with(~ gsub("ME", "", .x))), by = "plant_num")numeric_filtered <- df_global %>% dplyr::select(-c("plant_num", "sulfur_condition", "genotype")) %>% dplyr::select(where(~ sd(., na.rm = TRUE) > 0)) %>% select(where(~ any(!is.na(.) & . != 0))) # garde les colonnes ayant au moins une valeur ≠ 0 et ≠ NA# Combine avec les infos d'identitédf_global_clean <- df_global %>% dplyr::select(plant_num, sulfur_condition, genotype) %>% bind_cols(numeric_filtered)df_global = df_global_clean %>% dplyr::select(-c("METABO_Compound 22: glutamic acid","METABO_Compound 55: melibiose", "METABO_[threonic acid-1,4-lactone]", "METABO_maltose"))#### filter on SD beceau i havent sulfate for SSdf_global_SD = df_global_clean %>% filter(sulfur_condition == "SD") %>% dplyr::select(-c("METABO_Compound 22: glutamic acid","METABO_Compound 55: melibiose", "METABO_[threonic acid-1,4-lactone]", "METABO_maltose"))df_global_SD_without_metabo = df_global_clean %>% filter(sulfur_condition == "SD") %>% dplyr::select(-c("METABO_Compound 22: glutamic acid","METABO_Compound 55: melibiose", "METABO_[threonic acid-1,4-lactone]", "METABO_maltose")) %>% dplyr::select(-contains("METABO"))df_global_without_anion = df_global_clean %>% dplyr::select(-c("ANION_sulfate", "ANION_nitrate", "METABO_Compound 22: glutamic acid","METABO_Compound 55: melibiose", "METABO_[threonic acid-1,4-lactone]", "METABO_maltose"))### correlation standardcor_matrix <- cor(df_global[, -c(1:3)], use = "pairwise.complete.obs", method = "pearson")ggcorrplot(cor_matrix, type = "upper", lab = TRUE, colors = c("blue", "white", "red"), hc.order = TRUE) # active le tri par dendrogrammestandard_correlation = correlation_summary(df_global_SD[, -c(1:3)])standard_correlation = correlation_summary(df_global_without_anion[, -c(1:3)])standard_correlation = correlation_summary(df_global_SD_without_metabo[, -c(1:3)])####### creation of the LMM ####calculate_critical_value <- function(p_value) { # Calculate the corresponding cumulative probability for a two-way p-value prob_cumulative <- 1 - p_value / 2 # Find the corresponding z-value in the standard normal distribution critical_value <- qnorm(prob_cumulative) return(critical_value)}#here change datasetdf_compile = df_global_without_aniondf_compile = df_global_SD # change value of df_compile = df_compile %>% mutate(genotype = as.factor(genotype)) %>% dplyr:: rename_with(~ gsub("_+", "_", # remplace les __ par _ gsub("\\[|\\]", "", # supprime [ et ] gsub(":", "", gsub("-", "_", gsub(",", "_", gsub(" ", "_", .x)))))))#vec_variable=colnames(df_iono[8:length(df_iono)]) #data frame plant_num onlyvec_variable=colnames(df_compile[4:length(df_compile)])comb_vec <- combn(vec_variable, 2, simplify = F)comb_vec_df <- do.call(rbind, comb_vec) %>% as.data.frame()nb_cores <- detectCores()registerDoParallel(cores = nb_cores-2)tictoc::tic()df_result <- foreach(i = 1:length(comb_vec_df$V1), .combine = rbind) %dopar% { library(dplyr) library(tidyr) library(lme4) library(AICcmodavg) library(patchwork) library(performance) library(MuMIn) V1_x <- comb_vec_df[i, 1] V2_x <- comb_vec_df[i, 2] df_compile_select <- df_compile %>% drop_na(V1_x, V2_x) %>% as.data.frame() df_compile_select[,V1_x]<-(df_compile_select[,V1_x]-mean(df_compile_select[,V1_x]))/sd(df_compile_select[,V1_x]) df_compile_select[,V2_x]<-(df_compile_select[,V2_x]-mean(df_compile_select[,V2_x]))/sd(df_compile_select[,V2_x]) m1 <- lmer(formula(paste(V1_x, '~', V2_x ," + (1|genotype)+ (1|sulfur_condition) ")), data = df_compile_select, REML=T) #m1 <- lmer(formula(paste(V1_x, '~', V2_x ," + (1|genotype)")), data = df_compile_select, REML=T) AICc_val <- AICc(m1) df_coeff <- lmerTest:::get_coefmat(m1) |> as.data.frame() intercept <- df_coeff[1,1] slope <- df_coeff[2,1] pval <- df_coeff[2,5] result <- summary(m1) residuals <- result$coefficients[4] CI_sup05 <- slope + residuals * calculate_critical_value(0.05) CI_inf05 <- slope - residuals * calculate_critical_value(0.05) CI_sup01 <- slope + residuals * calculate_critical_value(0.01) CI_inf01 <- slope - residuals * calculate_critical_value(0.01) CI_sup001 <- slope + residuals * calculate_critical_value(0.001) CI_inf001 <- slope - residuals * calculate_critical_value(0.001) r2 <- r.squaredGLMM(m1) r2m <- r2[1] r2c <- r2[2] corr=cor(df_compile_select[,V1_x], df_compile_select[,V2_x], method = 'pearson') cat(i, "_", V1_x, "_", V2_x, "\n") return(data.frame(V1 = V1_x, V2 = V2_x, AICc = AICc_val, intercept = intercept, slope = slope, pval = pval, residuals = residuals, CI_sup05 = CI_sup05, CI_inf05 = CI_inf05, CI_sup01 = CI_sup01, CI_inf01 = CI_inf01, CI_sup001 = CI_sup001, CI_inf001 = CI_inf001, r2m = r2m, r2c = r2c, corr= corr) %>% mutate(between_IC05 = ifelse(CI_sup05 > CI_inf05, ifelse(0 >= CI_inf05 & 0 <= CI_sup05, "yes", "no"), "no")) %>% mutate(between_IC01 = ifelse(CI_sup01 > CI_inf01, ifelse(0 >= CI_inf01 & 0 <= CI_sup01, "yes", "no"), "no")) %>% mutate(between_IC001 = ifelse(CI_sup001 > CI_inf001, ifelse(0 >= CI_inf001 & 0 <= CI_sup001, "yes", "no"), "no")) )}stopImplicitCluster()tictoc::toc()save(df_result,file = here::here("data/multi_omic/LMM/output/result_LMM_all_without_anion.RData"))save(df_result,file = here::here("data/multi_omic/LMM/output/result_LMM_SD.RData"))save(standard_correlation,file = here::here("data/multi_omic/LMM/output/standard_correlation_SD.RData"))save(standard_correlation,file = here::here("data/multi_omic/LMM/output/standard_correlation_SD_without_METABO.RData"))```Verification ```{r, eval = F}load(file = here::here("data/multi_omic/LMM/output/result_LMM_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/result_LMM_all_without_anion.RData"))df_result %>% filter(pval<0.05) %>% filter(V1 == "ANION_sulfate"|V2== "ANION_sulfate")standard_correlation %>% filter(fdr<0.05) %>% filter(V1 == "ANION_sulfate"|V2== "ANION_sulfate")ggplot(df_global_SD, aes(x = ANION_sulfate, y = MODULE_violet, col = genotype))+ geom_point()ggplot(df_compile, aes(x = ANION_sulfate, y = METABO_glycolic_acid, col = genotype))+ geom_point()df_result %>% filter(pval<0.05) %>% filter(V1 == "ANION_nitrate"|V2== "ANION_nitrate")standard_correlation %>% filter(fdr<0.05) %>% filter(V1 == "ANION_nitrate"|V2== "ANION_nitrate")```## Selection based on pvalue and R²R2m (Marginal R²): This represents the variance explained by the fixed effects in your model. It is similar in interpretation to the R² value in traditional linear regression models. R2c (Conditional R²): This represents the variance explained by both the fixed and random effects in your model. It considers both fixed and random effects, providing a more comprehensive measure of goodness-of-fit for mixed-effects models. Both Rm2 and R2c can be useful in understanding the proportion of variance explained by the predictors in your model. ```{r, eval = F} clean_session() # Remove data frames and matrices and Reset plots# data importationload(file = here::here("data/multi_omic/LMM/output/result_LMM_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD.RData"))load(file = here::here("data/multi_omic/LMM/output/result_LMM_all_without_anion.RData"))load(file = here::here("data/multi_omic/LMM/output/standard_correlation_SD_without_METABO.RData"))#load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_no_prefiltre.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_5.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_7.RData"))# load(file = here::here("data/multi_omics/output/B/tmp1_LMM_test_prefiltre_correlation_9.RData"))# parameterr2m_i = 0.5r2c_i = 0.85fdr_lim = 0.05df_result = df_result %>% mutate(fdr = p.adjust(pval, method = "fdr"), log10fdr = -log10(fdr)) %>% mutate(color=ifelse(r2c*1.5 - r2m*1.5 > 0.5, "red", ifelse(r2m<r2m_i, ifelse(r2c<r2c_i,"#DD6342","orange"), ifelse(fdr<fdr_lim, "green","forestgreen") ) ) )# for network related to mixed effect# test more restrictive for fixed effect# r2m_i = 1# r2c_i = 0.95# fdr_lim = 0.05# # df_result = df_result %>% # mutate(fdr = p.adjust(pval, method = "fdr"),# log10fdr = -log10(fdr)) %>% # # mutate(color=ifelse(r2c > r2c_i, "green","red"# )# )# if i whant to use normale correlationstandard_correlation %>% filter(V1 == "ANION_sulfate"|V2== "ANION_sulfate") %>% arrange (fdr)# Verifdf_result %>% filter(pval<0.05) %>% filter(V1 == "ANION_sulfate"|V2== "ANION_sulfate")# Verif# Visualization ## Use of randomized data to avoid displaying all pointsrandom_result <- df_result %>% sample_n(1000)random_result%>% ggplot(aes(x = r2m, y = -log10(pval), size = slope, col = -log10(pval))) + geom_point(alpha = 0.7) + scale_color_viridis_c() + labs(title = "Dispersion diagramme", x = "R2m", y = "-log10(fdr)", color = "-log10(fdr)") + theme_minimal()random_result %>% ggplot(aes(x = r2c, y = -log10(pval), size = slope, col = -log10(pval))) + geom_point(alpha = 0.7) + scale_color_viridis_c() + labs(title = "Dispersion diagramme", x = "R2c", y = "-log10(fdr)", color = "-log10(fdr)") + theme_minimal()random_result %>% ggplot(aes(x = r2m, y = r2c, col = color, size = slope)) + geom_point(aes(shape = slope < 0), na.rm = TRUE) + scale_shape_manual(values = c(`TRUE` = 17, `FALSE` = 16)) + scale_color_manual(values = c("red" = "red", "#DD6342"="#DD6342","orange" = "orange", "green" = "green", "forestgreen" = "forestgreen")) + # Utilisation de couleurs fixes geom_vline(xintercept = r2m_i, linetype = "dashed", color = "red", linewidth = 1) + geom_hline(yintercept = r2c_i, linetype = "dashed", color = "red", linewidth = 1)# filter the result to have variabledf_result_select <- df_result %>% filter(color=="green") %>% arrange(pval)cat_col(paste("Number of variable:",length(unique(c(df_result_select$V1,df_result_select$V2))), "\n"), "green")round(dim(df_result_select)[1]*100/dim(df_result)[1], 2)# for standard correlation ? #df_result_select_corr <- standard_correlation %>% df_result_select_corr <- standard_correlation %>% mutate(color = ifelse(fdr<0.05, "green", "red")) %>% filter(color=="green") %>% arrange(fdr)cat_col(paste("Number of variable:",length(unique(c(df_result_select_corr$V1,df_result_select_corr$V2))), "\n"), "green")# exportsave(df_result_select, file = here::here("data/multi_omic/LMM/output/result_LMM_filter.RData"))save(df_result_select_corr, file = here::here("data/multi_omic/LMM/output/result_corr_filter.RData"))```## Network plot with Igraph```{r, eval = F}clean_session() # Remove data frames and matrices and Reset plotsload(file = here::here("data/multi_omic/LMM/output/result_LMM_filter.RData"))load(file = here::here("data/multi_omic/LMM/output/result_corr_filter.RData"))############ warning i change that #################df_result_select = df_result %>% mutate(fdr = p.adjust(pval, method = "fdr"), log10fdr = -log10(fdr)) %>% mutate(color=ifelse(fdr<fdr_lim, "green", "red") )df_result_select = df_result_select_corr# pb <- progress::progress_bar$new(# format = " Calcul [:bar] :percent | Étape :current/:total | Temps écoulé :elapsed | Temps restant :eta",# total = nrow(df_result_select),# clear = FALSE,# width = 100# )# # # Fonction pour calculer la corrélation avec mise à jour de la barre de progression# calculate_corr <- function(var1, var2) {# pb$tick()# # # Vérifier si les variables existent dans df_compile# if (all(c(var1, var2) %in% colnames(df_compile))) {# # Calculer la corrélation# cor_value <- cor(df_compile[[var1]], df_compile[[var2]], use = "pairwise.complete.obs")# return(cor_value)# } else {# return(NA) # Retourne NA si une variable manque# }# }# Utiliser mapply pour vectoriser le calcul#df_result_select$corr <- mapply(calculate_corr, df_result_select$V1, df_result_select$V2)# filter the data previously generate# # words_to_remove <- c("sEUpE", "DW","Ba","As","Cd","Ti","Be","Fe","V","Cr","Tl")# words_to_remove <- c("sEUpE")# words_to_remove <- c("sEUpE_C")# words_to_remove <- c("sEUpE_U")# words_to_remove <- c("sEUpE_U", "IONO", "ECO", "RSA")#words_to_remove <- c("[", "]")df_select=df_result_select %>% filter(color=="green") %>% filter(between_IC05=="no") #%>% #filter(!str_detect(V1, paste(words_to_remove, collapse = "|"))) %>% #filter(!str_detect(V2, paste(words_to_remove, collapse = "|")))df_select <- df_result_select %>% filter(color=="green") %>% filter(!grepl("putative", V1, ignore.case = TRUE)) %>% filter(!grepl("putative", V2, ignore.case = TRUE))df_select <- df_select %>% filter(color=="green") %>% filter(!grepl("putative", V1, ignore.case = TRUE)) %>% filter(!grepl("putative", V2, ignore.case = TRUE))#or with correlation# df_select=df_result_select_corr %>% # #filter(select=="yes") %>% # filter(color=="green")#hist(df_select$slope_modif)# creation of the node nodes<-data.frame(variable = unique(c(df_select$V1,df_select$V2))) %>% mutate(id = paste0("s", 1:length(unique(c(df_select$V1,df_select$V2))))) %>% relocate(id, .before=variable) %>% #mutate(variable_type = sample(1:3, nrow(.), replace = TRUE)) %>% #generate for aleatoire color # mutate(organe= case_when( # grepl("leaf", variable) ~ "leaf", # grepl("stem", variable) ~ "stem", # grepl("root", variable) ~ "root", # grepl("surface", variable) ~ "root", # grepl("volume", variable) ~ "root", # grepl("density", variable) ~ "root", # grepl("sEUpE", variable) ~ "sEUpE", # grepl("TR", variable) ~ "leaf", # TRUE ~ "autre" # )) %>% mutate(organe_num= case_when( # grepl("leaf", variable) ~ 1, # grepl("stem", variable) ~ 2, # grepl("root", variable) ~ 3, # grepl("surface", variable) ~ 3, # grepl("volume", variable) ~ 3, # grepl("density", variable) ~ 3, # grepl("sEUpE", variable) ~ 5, # grepl("TR", variable) ~ 1, # TRUE ~ 4 # )) %>% mutate(variable_type= case_when( grepl("METABO", variable) ~ 1, grepl("ANION", variable) ~ 4, grepl("MODULE", variable) ~ 2, #grepl("METABO.T", variable) ~ 5, TRUE ~ 3 )) %>% mutate( # mutate(variable_cleaned = ifelse(variable_type==1,gsub("[0-9_]", "", variable),variable)) %>% # mutate(variable_cleaned = ifelse(variable_type==4,gsub("[0-9_]", "", variable_cleaned),variable_cleaned)) %>% # mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(stem|leaf|root|concentration|sEUpE)", "", variable_cleaned)) %>% # dplyr::mutate( # variable_cleaned = gsub("weight", "DW", variable_cleaned), # variable_cleaned = gsub("sumtotalevapotranspiration", "ETtot", variable_cleaned), # variable_cleaned = gsub("Hydricpotential", " \u03C8 ", variable_cleaned), # variable_cleaned = gsub("Photo", "An", variable_cleaned), # variable_cleaned = gsub("volume", "Volume", variable_cleaned), # variable_cleaned = gsub("surface", "Surface", variable_cleaned), # variable_cleaned = gsub("shootratio", "S.R", variable_cleaned), # variable_cleaned = gsub("sumbiomass", "TotDW", variable_cleaned), # variable_cleaned = gsub("density", "Density", variable_cleaned), # variable_cleaned = gsub("area", "Area", variable_cleaned), # variable_cleaned = gsub("TRmmolms", "TR", variable_cleaned), # variable_cleaned = gsub("ConvexHull", "Hull", variable_cleaned), variable_cleaned = str_remove_all(variable, "METABO_|ANION_|MODULE_|\\[|\\]"), nb_char = nchar(variable_cleaned))# linkslinks<-df_select %>% dplyr::select(V1,V2, corr) %>% dplyr::mutate(weight = abs(corr)) %>% dplyr::mutate(weight_repulstion = 1-corr) %>% inner_join(nodes %>% dplyr::rename(V1=variable) %>% dplyr::rename(X1=id),.,by="V1") %>% inner_join(nodes %>% dplyr::rename(V2=variable) %>% dplyr::rename(X2=id),.,by="V2") %>% relocate(X1, .before=X2) %>% relocate(V1, .before=V2) # Converting the data to an igraph object:net <- graph.data.frame(links, nodes, directed=T) save(nodes, links, net,file = here::here("data/multi_omic/LMM/output/result_LMM_graph.RData"))# Generate colors base on trend type:colrs2 <- c("#68a500", "#b8af83", "#ce7f50","black","#345995")V(net)$color <- colrs2[V(net)$variable_type]V(net)$color = ifelse( V(net)$variable_type==2, V(net)$variable_cleaned, V(net)$color )# Compute node degree (#links) and use it to set node size:deg <- degree(net, mode="all")# if(fc_or_pval=="pval"){# V(net)$size <- V(net)$log10_pval_test_t*6 #ixi modifie la taille des points# }else if(fc_or_pval=="fc"){# V(net)$size <- V(net)$abslog2fc*10 #ixi modifie la taille des points# }myrhombus <- function(coords, v = NULL, params) { vertex.color <- params("vertex", "color") if (length(vertex.color) != 1 && !is.null(v)) { vertex.color <- vertex.color[v] } vertex.size <- 1/200 * params("vertex", "size") if (length(vertex.size) != 1 && !is.null(v)) { vertex.size <- vertex.size[v] } symbols(x = coords[, 1], y = coords[, 2], bg = vertex.color, stars = cbind(1.2*vertex.size, vertex.size, 1.2*vertex.size, vertex.size), add = TRUE, inches = FALSE)}# clips as a circleadd_shape("rhombus", clip = shapes("circle")$clip, plot = myrhombus)## Function for plotting an elliptical nodemyellipse <- function(coords, v=NULL, params) { vertex.color <- params("vertex", "color") if (length(vertex.color) != 1 && !is.null(v)) { vertex.color <- vertex.color[v] } vertex.size <- 1/70 * params("vertex", "size") # largeur des elipse if (length(vertex.size) != 1 && !is.null(v)) { vertex.size <- vertex.size[v] } plotrix::draw.ellipse(x=coords[,1], y=coords[,2], a = vertex.size, b=0.024, col=vertex.color) #b = hauteur}## Register the shape with igraphadd_shape("ellipse", clip=shapes("circle")$clip, plot=myellipse)# Fonction pour dessiner un hexagonemyhexagon <- function(coords, v = NULL, params) { vertex.color <- params("vertex", "color") if(length(vertex.color) != 1 && !is.null(v)){ vertex.color <- vertex.color[v] } vertex.size <- 1/200 * params("vertex", "size") if(length(vertex.size) !=1 && !is.null(v)){ vertex.size <- vertex.size[v] } for(i in seq_len(nrow(coords))){ x <- coords[i, 1] y <- coords[i, 2] size <- vertex.size[i] theta <- seq(0, 2*pi, length.out = 7)[-7] # Angles pour les sommets de l'hexagone xs <- x + size * cos(theta) ys <- y + size * sin(theta) polygon(x = xs, y = ys, col = vertex.color[i], border = "black") }}hexagon_clip <- function(coords, el, params) { vertex.size <- 1/200 * params("vertex", "size") angles <- seq(0, 2 * pi, length.out = 7)[-7] clip_coords <- sapply(1:nrow(coords), function(i) { size <- vertex.size[i] angle <- atan2(el[2, 2] - el[1, 2], el[2, 1] - el[1, 1]) # Edge direction idx <- which.min(abs(angles - angle)) # Find closest angle in hexagon x <- coords[i, 1] + size * cos(angles[idx]) y <- coords[i, 2] + size * sin(angles[idx]) c(x, y) }) return(t(clip_coords))}# Enregistrer la nouvelle forme avec igraphadd_shape("hexagon", clip = shapes("circle")$clip, plot = myhexagon)# V(net)$shape<-ifelse(V(net)$variable_type==2,"rectangle",ifelse(V(net)$variable_type==4,"ellipse",ifelse(V(net)$variable_type==3,"rhombus", ifelse(V(net)$variable_type==5, "hexagon", "hexagon"))))V(net)$size=1+deg*0.1 #avant 10V(net)$size = ifelse( V(net)$shape == "rectangle", V(net)$nb_char * 3.5 + 5, ifelse( V(net)$shape == "ellipse", ifelse(V(net)$nb_char > 3, 10+V(net)$nb_char*0.16,1.7), ifelse(V(net)$shape == "hexagone", V(net)$nb_char * 0.3 + 1, 7) # Taille ajustée pour les hexagones ))V(net)$size2=3.5 # hauteurV(net)$label.cex <- 0.8 # Taille plus petite des étiquettes# The labels are currently node IDs.# Setting them to NA will render no labels:colrs1 <- c("black", "white", "black","white")V(net)$label.color <-colrs1[ifelse(V(net)$variable_type=="METABO",2,V(net)$variable_type)] #color of the labelV(net)$label.color <- sapply(V(net)$color, evaluate_contrast) #color of the labelV(net)$label.dist <-0#.7 #distance from the center of the vertexV(net)$label.family="sans" #type de police#V(net)$label.font=V(net)$variable_type # 2 for bolt 1 is plaint texte 3 for italicV(net)$label.font=ifelse(V(net)$variable_type==5,1,V(net)$variable_type) # 2 for bolt 1 is plaint texte 3 for italic#V(net)$label <- NA# Set edge width based on weight:#E(net)$width <- -log10(1-abs(E(net)$slope_modif))*1.5# avt /6E(net)$width <- (abs(E(net)$corr)-min(abs(E(net)$corr))+0.01)*5# avt /6#change arrow size and edge color:E(net)$arrow.size <- .2E(net)$edge.color <- "gray80"# Let's color the edges of the graph based on their source node color.# We'll get the starting node for each edge with "get.edges"edge.start <- get.edges(net, 1:ecount(net))[,1]edge.col <- V(net)$color[edge.start]vec_color=ifelse(E(net)$corr>0,"red","blue")edge.col <-vec_color #mettre red if is upper green if is lower 0l=layout_in_circlel=layout_with_kk # Kamada Kawai #Like Fruchterman Reingold, it attempts to minimize the energy in a spring system.l=layout_on_gridl=layout_with_fr #Fruchterman-Reingoldl <-layout_with_fr(net, niter = 10000, grid = "nogrid")bound <- 1 # demi-largeur de la boîten <- vcount(net)for (sign in c("sign", "unsign")){ if (sign =="unsign"){ scaling_factor_attraction_force <- 3 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly. l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$weight-min(igraph::E(net)$weight)+0.001), niter = 500, grid = "nogrid") # without sign effect }else if (sign == "sign"){ scaling_factor_attraction_force <- 10 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly. l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$corr+abs(min(igraph::E(net)$corr))+0.001), niter = 500, grid = "nogrid")#, start.temp = bound, minx = rep(-bound, n), maxx = rep(bound, n), #miny = rep(-bound, n), maxy = rep(bound, n)) # withe sign effect } # export svg(width=8, height=8,filename = here::here(paste0("report/multi_omics/plot/LMM/network_SD_without_METABO_",sign,".svg"))) set.seed(1) plot(net, layout = l, edge.lty = 1, edge.arrow.size = 0, vertex.label = V(net)$variable_cleaned, vertex.label.cex = 0.5, vertex.size = V(net)$size*0.5, # Assurez-vous que la taille est passée ici edge.color = edge.col, edge.curved = .15) title(main = "Network plot based on LMM results show correlation between variables") mtext(paste0("Positive correlation: ",length(df_select %>% filter(corr>0) %>% pull(corr))), side = 1, line = 2, cex = 0.8,col = "red") mtext(paste0("Negative correlation: ",length(df_select %>% filter(corr<0) %>% pull(corr))), side = 1, line = 1, cex = 0.8,col = "blue") dev.off()}# save plot to analyse in cytoskapewrite_graph(net, file = here::here("data/multi_omics/output/B/network.graphml"), format = "graphml")## Degree Distribution Histogramdf_deg=tibble(nodes=names(deg),deg=as.vector(deg)) %>% full_join(nodes %>% as.data.frame() %>% dplyr::rename(nodes=id),.,by="nodes") %>% arrange(desc(deg)) %>% #mutate(variable_cleaned = ifelse(variable_type==1|variable_type==4,gsub("[0-9_]", "", variable),variable)) %>% # mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(concentration)", "", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(stem)", "", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(root)", " ", variable_cleaned)) %>% # mutate(variable_cleaned = gsub("(leaf)", " ", variable_cleaned)) %>% mutate(organe=factor(organe, levels = c("leaf", "stem", "root","autre", "sEUpE"))) %>% #mutate(variable_cleaned = gsub("sEUpE", "sEUE ", variable_cleaned)) %>% mutate(variable_cleaned = paste0(variable_cleaned, " ", str_to_title(organe))) # mutate(variable_cleaned = gsub("(stem|leaf|root|concentration)", "", variable_cleaned))px=df_deg %>% filter(deg > 4) %>% ggplot(aes(x = deg, y = reorder(variable_cleaned, -deg), fill = organe, label = deg)) + geom_col(position = 'dodge') + geom_text(hjust = -0.4) + scale_fill_manual(values = colrs2, breaks = c("leaf","stem","root", "autre", "sEUpE")) + # scale_x_continuous(limits = c(0, 15),breaks = seq(0, 15, by = 1))+ theme_minimal() + theme(axis.title.y = element_blank(), axis.text.x = element_blank(), axis.ticks.y = element_blank(), panel.grid = element_blank(), legend.position = "none")+ labs(x="Degree", fill="Organe")ggsave(here::here("report/multi_omics/plot/LMM_metabo_iono/degree_all.svg"), px, height = 15,width = 10)############# module les plus correlele avec metabolitetop_corr_meta_module <- df_select %>% ## 1. garder uniquement les lignes qui contiennent un METABO et un MODULE filter( (str_starts(V1, "METABO_") & str_starts(V2, "MODULE_")) | (str_starts(V1, "MODULE_") & str_starts(V2, "METABO_")) ) %>% ## 2. corrélation en valeur absolue mutate(abs_corr = abs(corr)) %>% ## 3. trier du plus fort au plus faible arrange(desc(abs_corr)) module_connection_counts <- df_select %>% ## 1. ne garder que les couples METABO / MODULE filter( (str_starts(V1, "METABO_") & str_starts(V2, "MODULE_")) | (str_starts(V1, "MODULE_") & str_starts(V2, "METABO_")) ) %>% ## 2. isoler le nom du module et du métabolite dans deux nouvelles colonnes mutate( module = if_else(str_starts(V1, "MODULE_"), V1, V2), metabolite = if_else(str_starts(V1, "METABO_"), V1, V2) ) %>% ## 3. compter le nombre de métabolites distincts par module group_by(module) %>% summarise( n_metabolites = n_distinct(metabolite), .groups = "drop" ) %>% ## 4. trier du plus connecté au moins connecté arrange(desc(n_metabolites))module_connection_counts```{fig-align="center"}## Analyse module ```{r, eval = F}input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2lfc_lim_i = 0moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <- c("PsCam", "moduleColors")# transformation into a listgene_list <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(genes = list(PsCam), .groups = "drop") %>% deframe() # transforme en liste nomméeload(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))# load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>% left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam")# cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Down) ", comparison)), # purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Up) ", comparison))# ) %>% full_join(., df_info_gene, by = "PsCam") %>% # mutate(ID = Psat) %>% # drop_na(comparison, ID)col_order <- c( "darkolivegreen", "orange", "yellowgreen", "purple", "darkgreen", "midnightblue", "lightgreen", "blue", "grey60", "steelblue")#for article, only for SD conditioncol_order <- c( "saddlebrown", "black", "orange", "tan")col_order <- c( "green", "steelblue", "orange")cluster_info <- purrr::imap_dfr(gene_list, ~ tibble(PsCam = .x, cluster = .y)) %>% full_join(df_info_gene, by = "PsCam") %>% mutate(ID = Psat) %>% drop_na(ID) %>% filter(cluster %in% col_order) %>% # tu gardes les bonnes couleurs mutate( # ① méthode dplyr de base cluster = factor(cluster, levels = col_order, ordered = TRUE) # ② méthode forcats, équivalente si tu préfères : # cluster = fct_relevel(cluster, !!!col_order) ) #filter(cluster %in% c("tan", "saddlebrown"))source(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO termedf_GO = GO_on_different_group(functional_roles = "BP", group_info = cluster_info, group = "cluster", ID = "ID", top = 10 )test = df_GO %>% filter(group == "tan")test = df_GO %>% filter(group == "saddlebrown")cyrcadian <- as.data.frame(gene_list[["orange"]]) ; colnames(cyrcadian) <- "PsCam"cyrcadian <- cyrcadian %>% left_join(., df_info_gene, by = "PsCam")# export for titouan#writexl::write_xlsx(cyrcadian, path = here::here("data/multi_omic/cyrcadian.xlsx"))# Merge this with your original data to fill missing combinationsCluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>% mutate(cluster = factor(cluster, levels = unique(cluster))#, # comparison = str_trim(str_remove(cluster, "\\s*\\(.*?\\)")), # sign = str_extract(cluster, "(?<=\\().*?(?=\\))") ) #%>% #mutate(comparison = forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT"))#vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF")# Plot the enrichment by GOpx <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + geom_tile(color = "black") + scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) + theme_minimal() + theme(axis.text = element_text(size = 8, colour = "black"), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid.major = element_blank(), plot.background = element_rect(fill = "white", colour = "white")) + labs(fill = "-log10(FDR)", x= "Biological process", y = "Cluster") + #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+ new_scale_fill() + scale_fill_manual(values = col_order, name = "Module" ) +geom_tile( aes( x = cluster, y = -1, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = cluster, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE )# export #fig_export(path =paste0("report/multi_omics/plot/GO/darkolivegreen_orange_GO_SD_condition"), plot_x = px, height_i = 5.5, width_i = 4.5, res = 600)fig_export(path =paste0("report/multi_omics/plot/GO/some_module_GO_all_without_metabo"), plot_x = px, height_i = 7, width_i = 5.5, res = 600)```{fig-align="center"}{fig-align="center"}## List of genes in each module ```{r, eval = F}input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2lfc_lim_i = 0moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv")))colnames(moduleColors) <- c("PsCam", "moduleColors")# transformation into a listgene_list <- moduleColors %>% dplyr::group_by(moduleColors) %>% dplyr::summarise(genes = list(PsCam), .groups = "drop") %>% deframe() # transforme en liste nomméeload(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))# load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>% left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam")# cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Down) ", comparison)), # purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% # mutate(cluster = paste0("(Up) ", comparison))# ) %>% full_join(., df_info_gene, by = "PsCam") %>% # mutate(ID = Psat) %>% # drop_na(comparison, ID)col_order <- c( "orange", "black", "tan", "saddlebrown")cluster_info <- purrr::imap_dfr(gene_list, ~ tibble(PsCam = .x, cluster = .y)) %>% full_join(df_info_gene, by = "PsCam") %>% mutate(ID = Psat) %>% drop_na(ID)############# a revoirsource(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO termedf_GO = GO_on_different_group(functional_roles = "BP", group_info = cluster_info, group = "cluster", ID = "ID", top = 10 )#info from excel info_from_excel = read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), skip = 11) %>% dplyr::select("N°", "GENE_ID", "GENE_TYPE","DESCRIPTION", "...7","CAT") %>% dplyr::rename(PsCam = "N°", At = "...7")df_excel <- df_GO %>% dplyr::select(-c("GeneRatio", "BgRatio", "RichFactor", "FoldEnrichment", "zScore", "pvalue", "p.adjust", "qvalue", "Count","|-log10(Pval)|", "ID", "Description")) %>% drop_na(geneID) %>% separate_rows(geneID, sep = "/") %>% dplyr::group_by(geneID,group) %>% summarise(Description_GO = paste(Description_GO, collapse = "/"), .groups = "drop") %>% dplyr::rename(ModuleColor=group) %>% dplyr::rename(Psat = geneID) %>% dplyr::select(Psat, Description_GO) %>% full_join(.,cluster_info, by ="Psat") %>% left_join(., info_from_excel, by ="PsCam") %>% drop_na(cluster)colnames(df_excel)df_info <- as.data.frame(gene_list[["orange"]]) ; colnames(cyrcadian) <- "PsCam"cyrcadian <- cyrcadian %>% left_join(., df_info_gene, by = "PsCam")# export for titouanwritexl::write_xlsx(cyrcadian, path = here::here("data/multi_omic/cyrcadian.xlsx"))writexl::write_xlsx(df_excel, path = here::here("data/multi_omic/table_supp_all_genes.xlsx"))# Merge this with your original data to fill missing combinationsCluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>% mutate(cluster = factor(cluster, levels = unique(cluster))#, # comparison = str_trim(str_remove(cluster, "\\s*\\(.*?\\)")), # sign = str_extract(cluster, "(?<=\\().*?(?=\\))") ) #%>% #mutate(comparison = forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT"))```## Heatmap statsThe aim is to create heatmap portions with the stats inside, and then put them back on genes pathways.```{r, eval = F}load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) colnames(moduleColors) <- c("gene_id", "moduleColors")moduleColors = moduleColors%>% left_join(., read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol", col_types = c("text", "text", "text")) %>% dplyr::select(PsCam, gene_symbol) %>% dplyr::rename(gene_id = PsCam), by = "gene_id")moduleColors_select = moduleColors %>% drop_na("gene_symbol")#################### filter genes of interest# moduleColors_select = moduleColors %>% filter(moduleColors == "steelblue")#info for more simple name of genes# # input one metabolite# variable_i <- "Asp"# variable_i <- "OAS"# DAF_i <- 29v_variable <- moduleColors_select %>% pull(gene_id)expr_long <- as.data.frame(datExpr) %>% rownames_to_column(var = "sample") %>% # garder l'ID échantillon pivot_longer(-sample, names_to = "gene_id", values_to = "expression")# b. ajouter la couleur de module pour chaque gèneexpr_long_select <- expr_long %>% left_join(., moduleColors_select, by = "gene_id") %>% left_join(., metadata %>% rownames_to_column("sample"), by = "sample") %>% mutate(condition = paste0(sulfur_condition, "_", genotype)) %>% filter(gene_id %in% v_variable)#compartment_i <- "N6_N7"# measure min max of Zscore for color of the graduationdata_for_guides = expr_long_select %>% dplyr::group_by(gene_id) %>% dplyr::mutate(value_scaled = as.vector(scale(expression))) %>% ungroup() %>% dplyr::group_by(gene_id, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale = mean(value_scaled))# min(data_for_guides$mean_scale)# max(data_for_guides$mean_scale)for(variable_i in v_variable){ # variable_i <- "OAS" # variable_i <- "Asp" # filter data df_select <- expr_long_select %>% filter( gene_id == variable_i#, #compartment == compartment_i ) %>% drop_na(expression) %>% as.data.frame() # make stats ylab_i <- paste0(variable_i, " in vegetative leaves") res <- stat_analyse( data = df_select, column_value = "expression", category_variables = c("condition"), grp_var = "", show_plot = TRUE, outlier_show = FALSE, label_outlier = "sample", biologist_stats = TRUE, Ylab_i = ylab_i, control_conditions = "", strip_normale = FALSE, hex_pallet = c("#003049", "#780000", "#7FACC7", "#EC323E", "#003049", "#780000", "#7FACC7", "#EC323E") ) stats_resum <- res$data_used %>% distinct(condition, gene_id, group) res_resum_raw <- res$summary_result %>% left_join(.,stats_resum, by = "condition") %>% mutate(column_value = variable_i) res_resum <- res$data_used %>% mutate(value_scaled = scale(expression)) %>% # scale based on z-score is beter than force in -1 ; 1 dplyr::group_by(genotype, sulfur_condition, condition) %>% dplyr::summarise( mean_scale = mean(value_scaled), .groups = 'drop' ) %>% left_join(.,stats_resum, by = "condition") # combined_res[[variable_i]] <- res_resum final_res <- res_resum %>% mutate(color_text = "black") %>% mutate(condition = fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")) %>% #left_join(.,df_info_gene_module_color %>% dplyr::rename(gene_id = PsCam), by = "gene_id") %>% left_join(., moduleColors_select, by = "gene_id") min_max_zscore <- paste0("min:",round(min(final_res$mean_scale),2),";max:",round(max(final_res$mean_scale),2)) # creation of the plot p_plot <- ggplot(final_res, aes(x = condition, y = gene_id, fill = mean_scale)) + geom_tile(color = "white") + # Pour dessiner les cases geom_text(aes(label = group, color = color_text)) + # Ajouter les lettres statistiques scale_fill_gradient2(low = "#1d4877", mid = "white", high = "#ee3e32", midpoint = 0, limits = c(-2,2))+# Plage commune entre -2 et 2 pour tous les plots) + # Palette de couleur scale_color_manual(values = c("black" = "black", "white" = "white"))+ theme_minimal() + theme( axis.title.y = element_blank(), # Supprime le titre de l'axe Y axis.title.x = element_blank(), # Supprime le titre de l'axe X axis.text.y = element_text(size = 10, face = "bold"), # Style du texte de la variable axis.text.x = element_text(size = 12, angle = 90, vjust = +0.5, hjust = 1.0), panel.grid.major = element_blank(), # Supprime les grilles panel.grid.minor = element_blank(), panel.border = element_rect(color = final_res %>% head(1) %>% pull(moduleColors), linewidth = 2, fill = NA) ) + labs(title = final_res %>% head(1) %>% pull(gene_symbol))+ #guides(fill = guide_colorbar(title = "Mean Value")) ; p_plot guides(fill = "none", color = "none") ; p_plot # export plot plot_name <- paste0(variable_i) fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/", plot_name)), p_plot, height_i = .38*2+0.95, width_i = 4, res_i = 300,format = "svg")}#export legend# Création d'un ggplot minimal uniquement pour la légendep_legend <- ggplot(data_for_guides, aes(x = genotype, y = sulfur_condition, fill = mean_scale)) + geom_tile() + scale_fill_gradient2( low = "#1d4877", mid = "white", high = "#ee3e32", midpoint = (min(data_for_guides$mean_scale) + max(data_for_guides$mean_scale)) / 2, limits = c(min(data_for_guides$mean_scale), max(data_for_guides$mean_scale)) ) + theme_void() + # Supprime les axes, titres, etc. guides( fill = guide_colorbar(title = "Mean Value") # Titre de la légende )# Extraire uniquement la légendelegend <- cowplot::get_legend(p_legend)# Affichage de la légende seulegrid::grid.newpage()grid::grid.draw(legend)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/legend")), legend, height_i = 2, width_i = 1, res_i = 300,format = "svg")# test```Idem for metabolomic```{r, eval = F}df_metabolite_global <- read_csv(here::here("data/metabolomic/output/df_metabolite_global.csv"), show_col_types = FALSE) %>% mutate(plant_num = paste(sep = "_", plant_num, condition)) %>% mutate(sulfur_condition= fct_relevel(sulfur_condition, "SS", "SD"), genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), compartment = fct_relevel(compartment, "N4_N5", "N6_N7"), type_genotype = fct_relevel(type_genotype, "WT", "Mut"), condition = paste0(sulfur_condition, "_", genotype) ) %>% # filter(sulfur_condition == "SD") %>% filter(!(variable == "Met" & compartment == "N9_N10")) %>% # error of measure/ in WT2 = 0 in all filter(compartment == "N6_N7") %>% filter(DAF == 29)df_metabolite_global_test <- df_metabolite_global %>% filter(variable == "Sulfate")# df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) # # input_name = "datExpr"# networkType_i = "signed"# minModuleSize = 30# deepSplit = 2# # moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) # colnames(moduleColors) <- c("gene_id", "moduleColors")# # moduleColors = moduleColors%>% # left_join(.,# read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol",# col_types = c("text", "text", "text")) %>% # dplyr::select(PsCam, gene_symbol) %>% # dplyr::rename(gene_id = PsCam), by = "gene_id")# # moduleColors_select = moduleColors %>% drop_na("gene_symbol")#################### filter genes of interest# moduleColors_select = moduleColors %>% filter(moduleColors == "steelblue")#info for more simple name of genes# # input one metabolite# variable_i <- "Asp"# variable_i <- "OAS"# DAF_i <- 29v_variable <- df_metabolite_global %>% distinct(variable) %>% pull(variable)# expr_long <- as.data.frame(datExpr) %>% # # rownames_to_column(var = "sample") %>% # garder l'ID échantillon# pivot_longer(-sample,# names_to = "gene_id",# values_to = "expression")# b. ajouter la couleur de module pour chaque gène# expr_long_select <- expr_long %>% # left_join(., moduleColors_select, by = "gene_id") %>% # left_join(., metadata %>% rownames_to_column("sample"), by = "sample") %>% # mutate(condition = paste0(sulfur_condition, "_", genotype)) %>% # filter(gene_id %in% v_variable)#compartment_i <- "N6_N7"# measure min max of Zscore for color of the graduationdata_for_guides = df_metabolite_global %>% dplyr::group_by(variable) %>% dplyr::mutate(value_scaled = as.vector(scale(value))) %>% ungroup() %>% dplyr::group_by(variable, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale = mean(value_scaled))# min(data_for_guides$mean_scale)# max(data_for_guides$mean_scale)for(variable_i in v_variable){ # variable_i <- "OAS" # variable_i <- "Asp" # filter data df_select <- df_metabolite_global %>% filter( variable == variable_i#, #compartment == compartment_i ) %>% drop_na(value) %>% as.data.frame() # make stats ylab_i <- paste0(variable_i, " in vegetative leaves") res <- stat_analyse( data = df_select, column_value = "value", category_variables = c("condition"), grp_var = "", show_plot = TRUE, outlier_show = FALSE, label_outlier = "plant_num", biologist_stats = TRUE, Ylab_i = ylab_i, control_conditions = "", strip_normale = FALSE, hex_pallet = c("#003049", "#780000", "#7FACC7", "#EC323E", "#003049", "#780000", "#7FACC7", "#EC323E") ) stats_resum <- res$data_used %>% distinct(condition, variable, group) res_resum_raw <- res$summary_result %>% left_join(.,stats_resum, by = "condition") %>% mutate(column_value = variable_i) res_resum <- res$data_used %>% mutate(value_scaled = scale(value)) %>% # scale based on z-score is beter than force in -1 ; 1 dplyr::group_by(genotype, sulfur_condition, condition) %>% dplyr::summarise( mean_scale = mean(value_scaled), .groups = 'drop' ) %>% left_join(.,stats_resum, by = "condition") # combined_res[[variable_i]] <- res_resum # # # final_res <- res_resum %>% # # mutate(color_text = ifelse(mean_scale<-2, "white", "black")) %>% # mutate(color_text = "black") %>% # mutate(condition = fct_relevel(condition, "SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K"))# %>% # #left_join(.,df_info_gene_module_color %>% dplyr::rename(gene_id = PsCam), by = "gene_id") %>% # # left_join(., moduleColors_select, by = "gene_id") condition_levels <- c("SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")# ------------------------------------------------------------------# à INSÉRER juste avant la création du graphique, dans ta boucle# ------------------------------------------------------------------final_res <- res_resum %>% tidyr::complete( variable, condition = factor(condition_levels, levels = condition_levels), fill = list(mean_scale = NA_real_, group = "") ) %>% mutate( is_na = is.na(mean_scale), # <‑‑ nouvel indicateur color_text = "black", condition = fct_relevel(condition, !!!condition_levels) )# ------------------------------------------------------------------# 2e moitié inchangée, sauf na.value ajouté# ------------------------------------------------------------------p_plot <- ggplot(final_res, aes(x = condition, y = variable, fill = mean_scale)) + geom_tile(color = "white") + geom_text(aes(label = group, colour = color_text)) + geom_point( data = function(x) dplyr::filter(x, is_na), shape = 4, # forme “×” stroke = 1.2, # épaisseur des traits de la croix size = 5, colour = "white" ) + scale_fill_gradient2( low = "#55185D", mid = "white", high = "#ECB602", midpoint = 0, limits = c(-2, 2), na.value = "gray60" ) + scale_colour_manual(values = c("black" = "black", "white" = "white")) + theme_minimal() + theme( axis.title = element_blank(), axis.text.y = element_text(size = 10, face = "bold"), axis.text.x = element_text(size = 12, angle = 90, vjust = 0.5, hjust = 1), panel.grid = element_blank() ) + guides(fill = "none", colour = "none") min_max_zscore <- paste0( "min:", round(min(final_res$mean_scale, na.rm = TRUE), 2), ";max:", round(max(final_res$mean_scale, na.rm = TRUE), 2)) print(min_max_zscore) # creation of the plot # p_plot <- ggplot(final_res, aes(x = condition, y = variable, fill = mean_scale)) + # geom_tile(color = "white") + # Pour dessiner les cases # geom_text(aes(label = group, color = color_text)) + # Ajouter les lettres statistiques # scale_fill_gradient2(low = "#55185D", mid = "white", high = "#ECB602", midpoint = 0, limits = c(-2,2))+# Plage commune entre -2 et 2 pour tous les plots) + # Palette de couleur # scale_color_manual(values = c("black" = "black", "white" = "white"))+ # theme_minimal() + # theme( # axis.title.y = element_blank(), # Supprime le titre de l'axe Y # axis.title.x = element_blank(), # Supprime le titre de l'axe X # axis.text.y = element_text(size = 10, face = "bold"), # Style du texte de la variable # axis.text.x = element_text(size = 12, angle = 90, vjust = +0.5, hjust = 1.0), # panel.grid.major = element_blank(), # Supprime les grilles # panel.grid.minor = element_blank() # ) + # #labs(title = final_res %>% head(1) %>% pull(variable))+ # #guides(fill = guide_colorbar(title = "Mean Value")) ; p_plot # guides(fill = "none", color = "none") ; p_plot # export plot plot_name <- paste0(variable_i) fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/metabolomic/", plot_name)), p_plot, height_i = .38*1.4+0.95, width_i = 4, res_i = 300,format = "svg")}#export legend# Création d'un ggplot minimal uniquement pour la légendep_legend <- ggplot(data_for_guides, aes(x = genotype, y = sulfur_condition, fill = mean_scale)) + geom_tile() + scale_fill_gradient2( low = "#55185D", mid = "white", high = "#ECB602", midpoint = (min(data_for_guides$mean_scale) + max(data_for_guides$mean_scale)) / 2, limits = c(min(data_for_guides$mean_scale), max(data_for_guides$mean_scale)) ) + theme_void() + # Supprime les axes, titres, etc. guides( fill = guide_colorbar(title = "Mean Value") # Titre de la légende )# Extraire uniquement la légendelegend <- cowplot::get_legend(p_legend)# Affichage de la légende seulegrid::grid.newpage()grid::grid.draw(legend)fig_export(here::here(paste0("report/multi_omics/plot/for_pathway/metabolomic/legend")), legend, height_i = 2, width_i = 1, res_i = 300,format = "svg")```## Heatmap cathegorie for microarray ### If i use the normalized data (same as used for WGCNA) ```{r, eval = F}# data importation # Calculates the number of elements which satisfy the above condition (i.e. the number of elements which are neither NA nor 0) and checks whether this number is greater than 10% of the total length of the column (length(.)). If this is the case, it means that the column has more than 10% of values that are neither NA nor 0, and should therefore be kept.#group_pathway = compartment_iload(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx")) input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) colnames(moduleColors) <- c("gene_id", "moduleColors")moduleColors = moduleColors%>% left_join(., read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol", col_types = c("text", "text", "text", "text", "text", "text")) %>% dplyr::select(PsCam, gene_symbol, group_pathway, functional_class) %>% dplyr::rename(gene_id = PsCam), by = "gene_id")moduleColors_select = moduleColors %>% drop_na("gene_symbol")v_variable <- moduleColors_select %>% pull(gene_id)expr_long <- as.data.frame(datExpr) %>% rownames_to_column(var = "sample") %>% # garder l'ID échantillon pivot_longer(-sample, names_to = "gene_id", values_to = "expression")perform_pca <- function(data) { if (any(is.na(data))) { method <- "Method: missMDA (NA value in dataset)" nb <- estim_ncpPCA(data, quali.sup = 1:4, method.cv = "Kfold", verbose = FALSE) plot(0:5, nb$criterion, xlab = "Number of dimensions", ylab = "MSEP") res.comp <- imputePCA(data, quali.sup = 1:4, ncp = nb$ncp) PCA_result <- PCA(res.comp$completeObs, quali.sup = 1:4, graph = TRUE) } else { method <- "Method: Normal PCA (no NA value in dataset)" PCA_result <- PCA(data, quali.sup = 1:4, graph = TRUE) } return(list(PCA_result = PCA_result, method = method))}visualize_eigenvalues <- function(PCA_result) { fviz_eig(PCA_result, addlabels = TRUE, ylim = c(0, 50))}visualize_pca_variables <- function(PCA_result, group_pathway_i) { var <- get_pca_var(PCA_result) png(here::here(paste0("report/microarray/plot/PCA/PCA_seleced_genes_var_", group_pathway_i, ".png")), width = 16, height = 16, units = 'cm', res = 600) print(fviz_pca_var(PCA_result, col.var = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE)) dev.off()}# b. ajouter la couleur de module pour chaque gèneexpr_long_select <- expr_long %>% left_join(., moduleColors_select, by = "gene_id") %>% left_join(., metadata %>% rownames_to_column("sample"), by = "sample") %>% mutate(condition = paste0(sulfur_condition, "_", genotype)) %>% filter(gene_id %in% v_variable)# measure min max of Zscore for color of the graduationdata_for_guides = expr_long_select %>% dplyr::group_by(gene_id) %>% dplyr::mutate(value_scaled = as.vector(scale(expression))) %>% ungroup() %>% dplyr::group_by(gene_id, genotype, sulfur_condition) %>% dplyr::summarise(mean_scale = mean(value_scaled))prepare_data<-function(group_pathway_i){ expr_wide_select<- expr_long_select %>% filter(group_pathway == group_pathway_i) %>% dplyr::select(-c("gene_id", "moduleColors", "group_pathway", "Mutant_type","sulfur_condition","chip_color","chip_num","genotype", "functional_class")) %>% pivot_wider(names_from = "gene_symbol", values_from = "expression") %>% column_to_rownames("sample")}#compartment_i <- "N6_N7"# Function to calculate means by typecalculate_means_by_type <- function(data) { data %>% dplyr::group_by(condition) %>% dplyr::summarise(across(where(is.numeric), mean, na.rm = TRUE))}data = expr_wide_select # for teste to delcalculate_letter_signif <- function(data) { # data: a data frame where the first column is "climat_condition" # and the other columns represent metabolite values for each sample # Extract metabolite names (columns after "climat_condition") metabo_names <- colnames(data)[-1] # Replace hyphens and whitespace with underscores # Change colname of the dataframe # Initialize a list to store significance letters for each metabolite signif_letters <- list() # Loop through each metabolite for (metabo in metabo_names) { #cat(metabo, "\n") # Build the ANOVA model # aov_model <- aov(as.formula(paste0("`", metabo, "` ~ condition")), data = data) # # # Apply the Tukey HSD test # tukey_result <- TukeyHSD(aov_model, "condition") # # # Extract the p-values for the condition comparisons # tukey_groups <- multcompLetters4(aov_model, tukey_result) # # # Retrieve significance letters for this metabolite # letters <- tukey_groups$`condition`$Letters ####test new #### data2 <- data %>% dplyr::rename(response = all_of(metabo)) # on sécurise le nom aov_model <- aov(response ~ condition, data = data2) tukey <- TukeyHSD(aov_model, "condition") letters <- multcompLetters4(aov_model, tukey) letters <- letters$condition # Check if all letters are the same if (length(unique(letters)) == 1) { # If all letters are the same, replace them with empty strings letters <- rep("", length(letters)) } # add the same name as before # Store the letters in the list signif_letters[[metabo]] <- letters } return(signif_letters)}# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional per‑cell numeric values.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; if TRUE (default) display numbers inside cells## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional numeric values and/or# significance letters inside each cell.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; show numeric means if TRUE (default FALSE)# show_letters : logical; show df_letter_signif letters if TRUE# (default TRUE)## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------create_heatmap <- function(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, type_colors, df_letter_signif, group_pathway_i, cutree_rows_i, show_values = FALSE, show_letters = TRUE) { # Replace non‑breaking spaces by standard spaces in 'functional_class' type_metabo_s$functional_class <- stringi::stri_replace_all_fixed( type_metabo_s$functional_class, "\u00A0", # non‑breaking space " ", # regular space vectorize_all = FALSE ) type_metabo_s$functional_class <- factor(type_metabo_s$functional_class) # Order the mean table to match matrix columns mean_by_type$condition <- factor(mean_by_type$condition, levels = colnames(z_matrix_transposed_reordered)) mean_by_type_reorder <- mean_by_type[order(mean_by_type$condition), ] # ------------------------------------------------------------------------- # Column annotations # ------------------------------------------------------------------------- col_labels <- colnames(z_matrix_transposed_reordered) col_annotation <- data.frame( Sulfur_Condition = ifelse(startsWith(col_labels, "SS_"), "SS", "SD"), Genotype = sub("^[A-Z]{2}_", "", col_labels), row.names = col_labels ) |> dplyr::mutate( Sulfur_Condition = forcats::fct_relevel(Sulfur_Condition, "SS", "SD"), Genotype = forcats::fct_relevel(Genotype, "WT1", "W78*", "WT2", "E568K") ) sulfur_pal <- sulfate_pallet # yellow / cyan genotype_pal <- mutant_palette # custom palette annotation_colors <- list( Sulfur_Condition = sulfur_pal, Genotype = genotype_pal, Condition = type_colors, moduleColors = setNames(moduleColors_select$moduleColors, moduleColors_select$moduleColors), functional_class = c( # Biosynthesis "Cysteine biosynthesis (O‑acetyl‑serine pathway)" = "#7d092f", "Glutathione biosynthesis & turnover" = "#b32d45", # Methionine cycle / salvage "Methionine cycle & salvage" = "#ef5659", "Methionine cycle & salvage" = "#D95F02", # Transcription factors "NAC transcription factors" = "#ffb255", "WRKY transcription factors" = "#ffd673", # Regulation / signalling "Regulation / sulfur‑deficiency signalling" = "#c5c35e", # Activation / reduction / transport / detox "Sulfate activation (entry step)" = "#42a285", "Sulfate reduction" = "#577590", "Sulfate uptake & distribution" = "#6d597a", "Sulfite detoxification & export" = "#7D3C98" ) ) # ------------------------------------------------------------------------- # Heat‑map # ------------------------------------------------------------------------- ComplexHeatmap::pheatmap( z_matrix_transposed_reordered, cluster_rows = TRUE, cluster_cols = FALSE, show_rownames = TRUE, show_colnames = TRUE, cutree_rows = cutree_rows_i, labels_col = col_labels, annotation_col = col_annotation, annotation_colors = annotation_colors, annotation_row = type_metabo_s, main = stringr::str_to_title(group_pathway_i), color = colorRampPalette(c("#1d4877", "white", "#ee3e32"))(256), name = "Centering and Scaling", # Add text according to user’s choices cell_fun = function(j, i, x, y, width, height, fill) { if (isTRUE(show_values) || isTRUE(show_letters)) { # Retrieve numeric value only if needed if (isTRUE(show_values)) { value_num <- round( apply(t(mean_by_type_reorder |> tibble::column_to_rownames("condition")), 2, as.numeric)[i, j], 2 ) } else { value_num <- "" } # Retrieve significance letter letter <- if (isTRUE(show_letters)) df_letter_signif[i, j] else "" # Decide what to print text_value <- "" if (isTRUE(show_values) && value_num != "") { text_value <- value_num } if (isTRUE(show_letters) && letter != "") { text_value <- if (text_value != "") { paste0(text_value, " (", letter, ")") } else { letter } } # Draw only if something to draw if (text_value != "") { grid::grid.text(text_value, x, y, gp = grid::gpar(fontsize = 8)) } } } )}# Main PCA functionmake_pca_condition <- function(group_pathway_i = "sulfate pathway", cutree_rows_i = 7, show_values_i = TRUE) { # Step 1: Prepare the data all_data_select_h_x <- prepare_data(group_pathway_i) # Step 2: Perform PCA pca_results <- perform_pca(all_data_select_h_x) res_pca_targeted_metabolomic_r2_x <- pca_results$PCA_result method <- pca_results$method # Step 3: Visualize eigenvalues visualize_eigenvalues(res_pca_targeted_metabolomic_r2_x) # Step 4: Visualize PCA variables visualize_pca_variables(res_pca_targeted_metabolomic_r2_x, group_pathway_i) # Step 5: Calculate means by type and letter significant and prepare for heatmap type_metabo <- moduleColors_select %>% filter(group_pathway == group_pathway_i) %>% dplyr::select(moduleColors, gene_symbol, functional_class) %>% unique() %>% dplyr::rename(variable = gene_symbol) df_obs <- res_pca_targeted_metabolomic_r2_x$call$X #%>% #dplyr::select(-c(condition)) type_metabo_s <- type_metabo %>% dplyr::select(variable, functional_class, moduleColors) %>% # rename(Type = type) %>% filter(variable %in% colnames(df_obs)[-1]) %>% column_to_rownames("variable")# %>% # mutate(Type = recode(Type, # "AA" = "Amino acid", # "PHY" = "Phytohormone", # "AO" = "Organic acid", # "FLAVO" = "Flavonoid", # "CYT" = "Cytokinine", # "SUGAR" = "Sugar")) mean_by_type <- calculate_means_by_type(df_obs) letter_signif <- calculate_letter_signif(df_obs) # Step 6: Calculate Z-scores z_scores<-mean_by_type %>% column_to_rownames("condition") z_scores <- scale(z_scores) z_matrix <- as.matrix(z_scores) z_matrix_transposed <- t(z_matrix) z_matrix_transposed_reordered <- z_matrix_transposed[,c("SS_WT1", "SS_W78*", "SS_WT2", "SS_E568K", "SD_WT1", "SD_W78*", "SD_WT2", "SD_E568K")] # Step 7: Create table of letter significant (same as z_matrix_transposed) df_letter_signif <- as.data.frame(matrix(data = "", nrow = dim(z_matrix_transposed_reordered)[1], ncol = dim(z_matrix_transposed_reordered)[2])) colnames(df_letter_signif)<-colnames(z_matrix_transposed_reordered) rownames(df_letter_signif)<-rownames(z_matrix_transposed_reordered) for (metabolite_i in rownames(z_matrix_transposed_reordered)){ for (condition_i in colnames(z_matrix_transposed_reordered)){ df_letter_signif[metabolite_i,condition_i] <- letter_signif[[metabolite_i]]$Letters[condition_i] } } df_letter_signif[is.na(df_letter_signif)] <- "" # Step 8: Create heatmap type_colors <- c( setNames(as.character(mutant_palette), paste0("SS_", names(mutant_palette))), setNames(as.character(mutant_palette), paste0("SD_", names(mutant_palette))) ) # Define this variable as needed pheat<-create_heatmap(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, type_colors, df_letter_signif, group_pathway_i, cutree_rows_i, show_values = show_values_i ) # Step 9: PCA biplot PCA_X <- fviz_pca_biplot(res_pca_targeted_metabolomic_r2_x, geom.ind = "point", pointshape = 20, pointsize = 3, col.ind = all_data_select_h_x$condition, col.quanti.sup = "red", palett = type_colors, addEllipses = TRUE, legend.title = "Treatment", ellipse.level = 0.95, ellipse.alpha = 0.3, ellipse.type = "convex", arrowsize = 0.5, col.var = "black", axes = c(1, 2), repel = T ) + #ggtitle(label = paste0("Visualizing Individual PCA and Quantity Variable for ", group_pathway, " for Each Treatment"), # subtitle = method) + ggtitle(label = str_to_title(group_pathway_i)) theme_classic() ggsave(filename = here::here(paste0("report/microarray/plot/PCA/PCA_seleced_genes_var_path_form", group_pathway_i, "_by_treatment.svg")), plot = PCA_X, width = 18 * 1.5, height = 15 * 1.5, units = "cm") dev.off() svg(filename = here::here(paste0("report/microarray/plot/pheatmap/pheatmap_", group_pathway_i, "_by_treatment.svg")), width = 13 , height = 10) print(pheat) dev.off() return(list(PCA_X = PCA_X, pheat = pheat)) # Adjust the return values as needed}# execution of the functionpca_heat_sulfate_pathway=make_pca_condition("sulfate pathway")pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i = 2)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway")), p_sulfate_pathway, height_i = 13, width_i = 13, res_i = 600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway")), p_scenescence_pathway, height_i = 5, width_i = 13, res_i = 600)####### without value ###########pca_heat_sulfate_pathway=make_pca_condition("sulfate pathway", show_values_i = FALSE)pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i = 2, show_values_i = FALSE)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway_without_value")), p_sulfate_pathway, height_i = 10, width_i = 9, res_i = 600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway_without_value")), p_scenescence_pathway, height_i = 4, width_i = 8, res_i = 600)```{fig-align="center"}{fig-align="center"}### If i use the the data (and stats from the platforme) ```{eval = F}# data importation # Calculates the number of elements which satisfy the above condition (i.e. the number of elements which are neither NA nor 0) and checks whether this number is greater than 10% of the total length of the column (length(.)). If this is the case, it means that the column has more than 10% of values that are neither NA nor 0, and should therefore be kept.#group_pathway = compartment_i# load(file = here::here("data/multi_omic/WGCNA/datExpr.RData"))load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData"))lfc_lim_i <- 0list_ratio_statsdiff_exression <- list_ratio_stats %>% map(~ select(.x, id_probe, logFC, BH)) %>% # garde seulement les 3 colonnes bind_rows(.id = "source") %>% mutate( comparison = case_when( source == "Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1" ~ "SD_WT2 vs SD_E568K", source == "Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3" ~ "SS_WT2 vs SS_E568K", source == "Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2" ~ "SD_WT1 vs SD_W78*", source == "Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4" ~ "SS_WT1 vs SS_W78*", source == "WT_SD_RepBio_5 vs WT_SS_RepBio_5" ~ "SS_WT vs SD_WT", TRUE ~ source # si aucun des cas ne matche, on garde la valeur d'origine ) ) %>% dplyr::select(-c(source)) %>% dplyr::rename(gene_id = id_probe)# levels(as.factor(resultat$comparison))load(file = here::here("data/multi_omic/WGCNA/metadata.RData"))df_info_gene_module_color <- read_excel(here::here("data/multi_omic/table_supp_all_genes.xlsx"))input_name = "datExpr"networkType_i = "signed"minModuleSize = 30deepSplit = 2moduleColors <- read_csv(file = here::here(paste0("data/multi_omic/WGCNA/","gene_color_end_",input_name,"_sign",networkType_i,"_minsize",minModuleSize,"_dp",deepSplit,".csv"))) colnames(moduleColors) <- c("gene_id", "moduleColors")moduleColors_select = moduleColors%>% left_join(., read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "gene_symbol", col_types = c("text", "text", "text", "text", "text", "text")) %>% dplyr::select(PsCam, gene_symbol, group_pathway, functional_class) %>% dplyr::rename(gene_id = PsCam), by = "gene_id") %>% drop_na("gene_symbol")v_variable <- moduleColors_select %>% pull(gene_id)# b. ajouter la couleur de module pour chaque gèneexpr_long_select <- diff_exression %>% left_join(., moduleColors_select, by = "gene_id") %>% drop_na(gene_symbol)# measure min max of Zscore for color of the graduationdata_for_guides = expr_long_select %>% dplyr::group_by(gene_id) %>% dplyr::mutate(value_scaled = as.vector(scale(logFC))) %>% ungroup() %>% dplyr::group_by(gene_id, comparison) %>% dplyr::summarise(mean_scale = mean(value_scaled))prepare_data<-function(group_pathway_i){ expr_wide_select<- expr_long_select %>% filter(group_pathway == group_pathway_i) %>% dplyr::select(-c("gene_id", "moduleColors", "group_pathway","functional_class", "BH")) %>% pivot_wider(names_from = "gene_symbol", values_from = "logFC") return(expr_wide_select)}prepare_data_BH<-function(group_pathway_i){ expr_wide_select<- expr_long_select %>% filter(group_pathway == group_pathway_i) %>% dplyr::select(-c("gene_id", "moduleColors", "group_pathway","functional_class", "logFC")) %>% pivot_wider(names_from = "gene_symbol", values_from = "BH") return(expr_wide_select)}#compartment_i <- "N6_N7"# Function to calculate means by typecalculate_means_by_type <- function(data) { data %>% dplyr::group_by(comparison) %>% dplyr::summarise(across(where(is.numeric), mean, na.rm = TRUE))}# data = expr_wide_select # for teste to del# calculate_letter_signif <- function(data) {# # data: a data frame where the first column is "climat_condition"# # and the other columns represent metabolite values for each sample# # # Extract metabolite names (columns after "climat_condition")# metabo_names <- colnames(data)[-1]# # # Replace hyphens and whitespace with underscores# # # Change colname of the dataframe# # # Initialize a list to store significance letters for each metabolite# signif_letters <- list()# # # Loop through each metabolite# for (metabo in metabo_names) {# #cat(metabo, "\n")# # Build the ANOVA model# # aov_model <- aov(as.formula(paste0("`", metabo, "` ~ condition")), data = data)# # # # # Apply the Tukey HSD test# # tukey_result <- TukeyHSD(aov_model, "condition")# # # # # Extract the p-values for the condition comparisons# # tukey_groups <- multcompLetters4(aov_model, tukey_result)# # # # # Retrieve significance letters for this metabolite# # letters <- tukey_groups$`condition`$Letters# # ####test new ##### data2 <- data %>% # dplyr::rename(response = all_of(metabo)) # on sécurise le nom# # aov_model <- aov(response ~ condition, data = data2)# tukey <- TukeyHSD(aov_model, "condition")# letters <- multcompLetters4(aov_model, tukey)# # letters <- letters$condition# # # Check if all letters are the same# if (length(unique(letters)) == 1) {# # If all letters are the same, replace them with empty strings# letters <- rep("", length(letters))# }# # # add the same name as before# # # Store the letters in the list# signif_letters[[metabo]] <- letters# }# # return(signif_letters)# }calculate_star_signif <- function(data_BH) { p_to_stars <- function(p) { ifelse(p < 0.001, "***", ifelse(p < 0.01, "**", ifelse(p < 0.05, "*", ""))) } # data: a data frame where the first column is "climat_condition" # and the other columns represent metabolite values for each sample # Extract metabolite names (columns after "climat_condition") metabo_names <- colnames(data_BH[-1]) # Replace hyphens and whitespace with underscores # Change colname of the dataframe # Initialize a list to store significance letters for each metabolite signif_stars <- vector("list", length(metabo_names)) names(signif_stars) <- metabo_names # Loop through each metabolite comp_labels <- as.character(data_BH[[1]]) # noms des comparaisons for (metabo in metabo_names) { stars_vec <- p_to_stars(data_BH[[metabo]]) names(stars_vec) <- comp_labels # on nomme le vecteur signif_stars[[metabo]] <- stars_vec } return(signif_stars)}# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional per‑cell numeric values.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; if TRUE (default) display numbers inside cells## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------# ---------------------------------------------------------------------------# create_heatmap# ---------------------------------------------------------------------------# Draws a ComplexHeatmap::pheatmap with optional numeric values and/or# significance letters inside each cell.## Args:# z_matrix_transposed_reordered : numeric matrix (rows = metabolites, cols = samples)# mean_by_type : data‑frame with mean values by condition# type_metabo_s : data‑frame with row annotations# type_colors : named vector of colours for 'Condition'# df_letter_signif : matrix of significance letters, same dim as z_matrix_# group_pathway_i : string used for the heat‑map title# cutree_rows_i : integer; number of clusters for rows# show_values : logical; show numeric means if TRUE (default FALSE)# show_letters : logical; show df_letter_signif letters if TRUE# (default TRUE)## Returns:# A heat‑map plot (invisible; printed for its side‑effect).# ---------------------------------------------------------------------------create_heatmap_star <- function(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, type_colors, df_star_signif, group_pathway_i, cutree_rows_i, show_values = FALSE, show_star = TRUE) { # Replace non‑breaking spaces by standard spaces in 'functional_class' type_metabo_s$functional_class <- stringi::stri_replace_all_fixed( type_metabo_s$functional_class, "\u00A0", # non‑breaking space " ", # regular space vectorize_all = FALSE ) type_metabo_s <- type_metabo_s %>% mutate( functional_class = as.factor(functional_class) %>% # colonne d'origine # stringi::stri_replace_all_fixed("\u00A0", " ", FALSE) %>% # remplace NBS forcats::fct_relevel( # <-- ordre voulu ICI "Sulfate uptake & distribution", "Sulfate reduction", "Sulfite detoxification & export", "Cysteine biosynthesis", "Glutathione biosynthesis & turnover", "Methionine cycle & salvage", "Methionine cycle & salvage", "Regulation / sulfur deficiency signalling", "NAC transcription factors", "WRKY transcription factors", "Sulfate activation (entry step)" ) ) # Order the mean table to match matrix columns mean_by_type$comparison <- factor(mean_by_type$comparison, levels = colnames(z_matrix_transposed_reordered)) mean_by_type_reorder <- mean_by_type[order(mean_by_type$comparison), ] # ------------------------------------------------------------------------- # Column annotations # ------------------------------------------------------------------------- col_labels <- colnames(z_matrix_transposed_reordered) col_annotation <- data.frame( Sulfur_Condition = c("WT", "SS","SS","SD","SD"), row.names = col_labels ) %>% mutate(Sulfur_Condition = forcats::fct_relevel(Sulfur_Condition, "WT", "SS","SD")) # sulfur_pal <- sulfate_pallet # yellow / cyan sulfur_pal <- c( "WT" = "#9381FF", # noir (SS_WT vs SD_WT) "SS" = as.character(sulfate_pallet[1]), # vert (comparaisons SS only) "SD" = as.character(sulfate_pallet[2]) # jaune (comparaisons SD only))annotation_colors <- list( Sulfur_Condition = sulfur_pal, moduleColors = setNames(moduleColors_select$moduleColors, moduleColors_select$moduleColors), functional_class = c( # Biosynthesis "Sulfate uptake & distribution" = "#7d092f", "Sulfate reduction" = "#ef5659",#"#b32d45", "Sulfite detoxification & export" = "#D95F02", "Cysteine biosynthesis" = "#ffd673", "Glutathione biosynthesis & turnover" = "#42a285",#"#ffb255", # Methionine cycle / salvage "Methionine cycle & salvage" = "#6d597a", "Methionine cycle & salvage" = "#6d597a", # Transcription factors # Regulation / signalling "Regulation / sulfur deficiency signalling" = "#7D3C98", "NAC transcription factors" = "#c5c35e", "WRKY transcription factors" = "#ef5659" # , # Activation / reduction / transport / detox # "Sulfate activation (entry step)" = "#c5c35e", # , # , # = "#7D3C98" )#%>% # mutate(functional_class = forcats::fct_relevel(functional_class, # "Sulfate uptake & distribution" , # "Sulfate reduction" , # "Sulfite detoxification & export" , # "Cysteine biosynthesis" , # "Glutathione biosynthesis & turnover" , # # Methionine cycle / salvage # "Methionine cycle & salvage" , # "Methionine cycle & salvage" , # "Regulation / sulfur deficiency signalling" , # # Transcription factors # "NAC transcription factors" , # "WRKY transcription factors" , # # Regulation / signalling # # Activation / reduction / transport / detox # "Sulfate activation (entry step)" # # # # )) ) # ------------------------------------------------------------------------- # Heat‑map # ------------------------------------------------------------------------- ComplexHeatmap::pheatmap( z_matrix_transposed_reordered, cluster_rows = TRUE, cluster_cols = FALSE, show_rownames = TRUE, show_colnames = TRUE, cutree_rows = cutree_rows_i, labels_col = col_labels, annotation_col = col_annotation, annotation_colors = annotation_colors, annotation_row = type_metabo_s, main = stringr::str_to_title(group_pathway_i), color = colorRampPalette(c("#1d4877", "white", "#ee3e32"))(256), name = "Log2 ratio", # Add text according to user’s choices cell_fun = function(j, i, x, y, width, height, fill) { if (isTRUE(show_values) || isTRUE(show_star)) { # Retrieve numeric value only if needed if (isTRUE(show_values)) { value_num <- round( apply(t(mean_by_type_reorder |> tibble::column_to_rownames("comparison")), 2, as.numeric)[i, j], 2 ) } else { value_num <- "" } # Retrieve significance letter star <- if (isTRUE(show_star)) df_star_signif[i, j] else "" # Decide what to print text_value <- "" if (isTRUE(show_values) && value_num != "") { text_value <- value_num } if (isTRUE(show_star) && star != "") { text_value <- if (text_value != "") { paste0(text_value, " (", star, ")") } else { star } } # Draw only if something to draw if (text_value != "") { grid::grid.text(text_value, x, y, gp = grid::gpar(fontsize = 8)) } } } )}# Main PCA functionmake_star_condition <- function(group_pathway_i = "sulfate pathway", cutree_rows_i = 7, show_values_i = TRUE) { # Step 1: Prepare the data all_data_select_h_x <- prepare_data(group_pathway_i) all_data_select_h_x_BH <- prepare_data_BH(group_pathway_i) # Step 5: Calculate means by type and letter significant and prepare for heatmap type_metabo <- moduleColors_select %>% filter(group_pathway == group_pathway_i) %>% dplyr::select(moduleColors, gene_symbol, functional_class) %>% unique() %>% dplyr::rename(variable = gene_symbol) type_metabo_s <- type_metabo %>% dplyr::select(variable, functional_class, moduleColors) %>% filter(variable %in% colnames(all_data_select_h_x)[-1]) %>% column_to_rownames("variable")# %>% # mutate(Type = recode(Type, # "AA" = "Amino acid", # "PHY" = "Phytohormone", # "AO" = "Organic acid", # "FLAVO" = "Flavonoid", # "CYT" = "Cytokinine", # "SUGAR" = "Sugar")) mean_by_type <- calculate_means_by_type(all_data_select_h_x) star_signif <- calculate_star_signif(all_data_select_h_x_BH) # Step 6: Calculate Z-scores z_scores<-mean_by_type %>% column_to_rownames("comparison") # z_scores <- scale(z_scores) z_matrix <- as.matrix(z_scores) z_matrix_transposed <- t(z_matrix) z_matrix_transposed_reordered <- z_matrix_transposed[,c("SS_WT vs SD_WT", "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K")] # Step 7: Create table of letter significant (same as z_matrix_transposed) df_star_signif <- as.data.frame(matrix(data = "", nrow = dim(z_matrix_transposed_reordered)[1], ncol = dim(z_matrix_transposed_reordered)[2])) colnames(df_star_signif)<-colnames(z_matrix_transposed_reordered) rownames(df_star_signif)<-rownames(z_matrix_transposed_reordered) for (metabolite_i in rownames(z_matrix_transposed_reordered)){ for (condition_i in colnames(z_matrix_transposed_reordered)){ df_star_signif[metabolite_i,condition_i] <- star_signif[[metabolite_i]][condition_i] } } df_star_signif[is.na(df_star_signif)] <- "" # Step 8: Create heatmap # type_colors <- c( # setNames(as.character(mutant_palette), paste0("SS_", names(mutant_palette))), # setNames(as.character(mutant_palette), paste0("SD_", names(mutant_palette))) # ) # Define this variable as needed pheat<-create_heatmap_star(z_matrix_transposed_reordered, mean_by_type, type_metabo_s, df_star_signif = df_star_signif, group_pathway = group_pathway_i, cutree_rows = cutree_rows_i, show_values = show_values_i) return(list(pheat = pheat)) # Adjust the return values as needed}# execution of the functionpca_heat_sulfate_pathway=make_star_condition("sulfate pathway")pca_heat_scenescence_pathway=make_star_condition("scenescence pathway", cutree_rows_i = 2)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway")), p_sulfate_pathway, height_i = 13, width_i = 13, res_i = 600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway")), p_scenescence_pathway, height_i = 5, width_i = 13, res_i = 600)####### without value ###########pca_heat_sulfate_pathway=make_pca_condition("sulfate pathway", show_values_i = FALSE)pca_heat_scenescence_pathway=make_pca_condition("scenescence pathway", cutree_rows_i = 2, show_values_i = FALSE)# Convertit le grob ComplexHeatmap en ggplotp_sulfate_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_sulfate_pathway[["pheat"]])))p_scenescence_pathway <- ggplotify::as.ggplot(grid::grid.grabExpr(ComplexHeatmap::draw(pca_heat_scenescence_pathway[["pheat"]])))# export figure fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_sulfat_pathway_without_value_star")), p_sulfate_pathway, height_i = 10, width_i = 8, res_i = 600)fig_export(here::here(paste0("report/microarray/plot/pheatmap/fig_x_pheatmap_scenescence_pathway_without_value_star")), p_scenescence_pathway, height_i = 4, width_i = 8, res_i = 600)```
## List of genes in each module
