fix order of DGE table columns

Author: torres-alexis

RNAseq/Pipeline_GL-DPPD-7101_Versions/GL-DPPD-7101-G.md

@@ -4,7 +4,7 @@
 
 ---
 
-**Date:** January 21, 2025  
+**Date:** January 24, 2025  
 **Revision:** G  
 **Document Number:** GL-DPPD-7101-G  
 
@@ -23,7 +23,7 @@ Lauren Sanders (GeneLab Project Scientist)
 
 ## Updates from previous version  
 
-Added separate pipeline document: [Microbes.md](Microbes.md) to document the pipeline steps for Bowtie2 alignment, used when the `--microbes` parameter is specified. In short, reads are aligned to a reference genome using Bowtie2 rather than STAR, gene counts are quantified using FeatureCounts rather than RSEM. Other steps remain unchanged.
+Added separate pipeline document: [GL-DPPD-7XXX.md](../Pipeline_GL-DPPD-7XXX_Versions/GL-DPPD-7XXX.md) to document the pipeline steps for Bowtie2 alignment, used when the `--microbes` parameter is specified. In short, reads are aligned to a reference genome using Bowtie2 rather than STAR, gene counts are quantified using FeatureCounts rather than RSEM. Other steps remain unchanged.
 
 Added "_GLbulkRNAseq" suffix to output files to prevent naming conflicts with files relevant to other assays.
 
@@ -1104,7 +1104,7 @@ echo "*: ${rRNA_count} rRNA entries removed." > *_rRNA_counts.txt
 
 ## 9. Normalize Read Counts and Perform Differential Gene Expression Analysis
 
-> **Note:** Steps 9a-9c are performed once as they set up the analysis framework. Steps 9d-9f are performed twice with different input files:
+> **Note:** DGE Analysis is performed twice with different sets of input files:
 > 1. Using raw RSEM genes.results files
 > 2. Using rRNA-removed raw RSEM genes.results files
 
@@ -1444,23 +1444,27 @@ output_table <- cbind(output_table, res_df)
 ### Add summary statistics ###
 output_table$All.mean <- rowMeans(normCounts, na.rm = TRUE)
 output_table$All.stdev <- rowSds(as.matrix(normCounts), na.rm = TRUE)
+output_table$LRT.p.value <- res_1_lrt@listData$padj
 
 ### Add group-wise statistics ###
 tcounts <- as.data.frame(t(normCounts))
 tcounts$group <- names(group)
 
-# Calculate group means
+# Calculate group means and standard deviations
 group_means <- as.data.frame(t(aggregate(. ~ group, data = tcounts, mean)))
-colnames(group_means) <- paste0("Group.Mean_", group_means['group',])
-
-# Calculate group standard deviations
 group_stdev <- as.data.frame(t(aggregate(. ~ group, data = tcounts, sd)))
-colnames(group_stdev) <- paste0("Group.Stdev_", group_stdev['group',])
 
-# Remove group name rows and combine statistics
+# Remove group name rows
 group_means <- group_means[-1,]
 group_stdev <- group_stdev[-1,]
-output_table <- cbind(output_table, group_means, group_stdev)
+
+# For each group, add mean and stdev columns
+for (group_name in names(group)) {
+    mean_col <- paste0("Group.Mean_(", group_name, ")")
+    stdev_col <- paste0("Group.Stdev_(", group_name, ")")
+    output_table[[mean_col]] <- group_means[, paste0("Group.Mean_", group_means['group',])]
+    output_table[[stdev_col]] <- group_stdev[, paste0("Group.Stdev_", group_stdev['group',])]
+}
 
 ### Read in GeneLab annotation table for the organism of interest ###
 annot <- read.table(annotations_link, 
@@ -1480,7 +1484,7 @@ output_table <- output_table %>%
 
 ```
 
-### 9h. Export DGE Tables
+### 9g. Export DGE Tables
 
 ```R
 ### Export unnormalized and normalized counts tables ###