I am conducting WGCNA analysis on my data and would like to know the best way to adjust for covariates. I have first normalized the data, and then adjusting for covariates (Age and Gender), followed by WGCNA analysis, as seen below:

covariates <- phenotype[, c("Age_yr", "Gender_M0F1")]
covariates$Gender_M0F1 <- as.factor(covariates$Gender_M0F1)

# Apply empirical Bayes linear model to adjust for Age and Gender
adjustedData <- empiricalBayesLM(
  data = exposures_Normalized,
  removedCovariates = covariates,  # Specify Age and Gender as covariates
  automaticWeights = "none",       # No automatic weights; adjust if needed
  getEBadjustedData = TRUE,        # Get the covariate-adjusted data
  verbose = 1                      # Show output for clarity
)
exposures_Adjusted <- adjustedData$adjustedData

# Perform WGCNA on the adjusted data
modules.omics <- blockwiseModules(exposures_Adjusted,
                                  power = 12, 
                                  networkType = "signed", 
                                  TOMType = "signed",
                                  corType = "bicor",
                                  deepSplit = 4, 
                                  minModuleSize = 15,  
                                  numericLabels = TRUE,
                                  verbose = 1,
                                  maxBlockSize = 8000, 
                                  nThreads = 6)

Does this approach seem reasonable? Any opinions would be appreciated.