Hi all, I am working on WGCNA of Arabidopsis data. I noticed an odd observation about inconsistent results from my analysis when using the same expression dataset but with differently ordered genes. I have my regular expression dataset (ordered by chromosome and gene ID) and the same expression data with the genes in a randomized order. When I run the analyses I get the exact same trait-module relationships (correlation and p-value) and the same numbers of genes per module, but when I look at the genes that make up the module the genes themselves are different. When using the expression data ordered by chromosome/gene ID my module is made up of all genes on Chromosome 1 (which seems biologically highly unlikely); when I run the random ordered expression data I get the same module construction (all values are identical; mod-trait relationship, kWithin values, etc) but different genes (in a seemingly random order). To clarify, the expression data is identical in both instances, but genes (rows) are ordered differently.

Here is the top of my script up to the network construction. When looking at the same module (ie. "green") all values are the same regardless of which version of the expression data I use, but the genes that make up the module are different. Any help or suggestions would be so very much appreciated.

```# Load expression data
    exp_csv_file <- "/Path/to/expression/data.csv"
    expression.data <- read.csv(exp_csv_file, header = TRUE, row.names = 1, check.names = FALSE) %>% as.data.frame()

    #transforming the data.frame so columns now represent genes and rows represent samples
    expression.data <- as.data.frame(t(expression.data))

    gsg <-goodSamplesGenes(expression.data)
    summary(gsg)
    gsg$allOK        #should return TRUE, if not, use filtering step below 

    #gsg filtering
    if (!gsg$allOK)
    {
      if (sum(!gsg$goodGenes)>0) 
        printFlush(paste("Removing genes:", paste(names(expression.data)[!gsg$goodGenes], collapse = ", "))); #Identifies and prints outlier genes
      if (sum(!gsg$goodSamples)>0)
        printFlush(paste("Removing samples:", paste(rownames(expression.data)[!gsg$goodSamples], collapse = ", "))); #Identifies and prints oulier samples
      expression.data <- expression.data[gsg$goodSamples == TRUE, gsg$goodGenes == TRUE] # Removes the offending genes and samples from the data
    }

    #after filtering gsg, rerun to ensure gsg$allOK returns TRUE
    gsg <-goodSamplesGenes(expression.data)
    summary(gsg)
    gsg$allOK

    sampleTree <- hclust(dist(expression.data), method = "average") #Clustering samples based on distance 

    #Setting the graphical parameters
    par(cex = 0.6);
    par(mar = c(0,4,2,0))

    #Plotting the cluster dendrogram
    plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
         cex.axis = 1.5, cex.main = 2)

    #if WGCNAThreads was enable, disable before running pickSoftThreshold
    disableWGCNAThreads()

    #Determining the Soft Power Threshold ...be patient
    spt <- pickSoftThreshold(expression.data) 
    spt

    par(mar=c(4,4,4,4))
    plot(spt$fitIndices[,1],spt$fitIndices[,2],
         xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
         main = paste("Scale independence"))
    text(spt$fitIndices[,1],spt$fitIndices[,2],col="firebrick3")
    abline(h=0.80,col="firebrick3")

    par(mar=c(5,5,5,5))
    plot(spt$fitIndices[,1], spt$fitIndices[,5],
         xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
         main = paste("Mean connectivity"))
    text(spt$fitIndices[,1], spt$fitIndices[,5], labels= spt$fitIndices[,1],col="firebrick3")

    # Load WGCNA
    library(WGCNA)
    options(stringsAsFactors = FALSE)

    # Set a fixed seed for reproducibility
    set.seed(1)

    # Prepare expression data (samples as rows, genes as columns)
    datExpr <- expression.data  # Replace with your data frame

    # Choose soft-threshold power (from pickSoftThreshold)
    softPower <- 8

    net <- blockwiseModules(
          datExpr,
          power = softPower,
          networkType = "signed",
          corType = "bicor",
          maxPOutliers = 0.1,
          quickCor = 0,
          TOMType = "signed",
          minModuleSize = 20,
          minKMEtoStay = 0.9,
          maxBlockSize = 35000,
          reassignThreshold = 0.05,
          mergeCutHeight = 0.05,
          numericLabels = FALSE,
          pamRespectsDendro = FALSE,
          saveTOMs = TRUE,                         # Save TOM for reproducibility
          saveTOMFileBase = "WGCNA_TOM_clubroot",  # Base name for TOM files
          randomSeed = 1,                          # Locks in module assignment
          verbose = 3
        )

# Extract module colors and eigengenes
moduleColors <- net$colors
MEs <- net$MEs

#calculate connectivity and retain gene IDs as rows
connectivity_all <- intramodularConnectivity.fromExpr(datExpr, moduleColors, power = softPower)
rownames(connectivity_all) <- colnames(datExpr)

# Plot dendrogram with module colors
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module Colors", dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

# Save module assignments for later use
write.csv(data.frame(Gene = colnames(datExpr), Module = moduleColors))
# Load phenotype data
pheno_csv_file <- "/Path/to/phenotype/data.csv"

# Load phenotype data
traitData <- read.csv(pheno_csv_file, header = TRUE, check.names = FALSE, stringsAsFactors = FALSE) %>% as.data.frame()

# Match samples using the 'Sample' column
Samples <- rownames(expression.data)
traitRows <- match(Samples, traitData$Sample)
datTraits <- traitData[traitRows, "Status", drop = FALSE]  # Keep as data frame

# Set rownames to match expression data
rownames(datTraits) <- Samples

# Convert Status to numeric (e.g., Control = 1, Infected = 2)
datTraits$Status <- as.numeric(as.factor(datTraits$Status))

# Check for missing values
anyNA(datTraits)  # Should be FALSE
str(datTraits)    # Should show Status as numeric

# Run correlation
module.trait.correlation <- cor(MEs, datTraits, use = "p")
nSamples <- nrow(expression.data)  # number of samples
module.trait.Pvalue <- corPvalueStudent(module.trait.correlation, nSamples)
```