如何使用R分析三组（处理前、后1/2个月）的200个基因表达数据？

Hey there! Let's walk through how to analyze your 3-group gene expression dataset step by step. You've got pre-treatment (BEF, 16 samples), 1 month post-treatment (AFT-1, 16 samples), and 2 months post-treatment (AFT-2, 9 samples) with 200 genes total. Here's a structured approach tailored to your data:

Step 1: Data Import & Preprocessing #

First, let's get your data into a usable format and clean up missing values (I noticed NaN in your example—R uses NA for missing values, so we'll adjust that):

# Create a data frame from your sample data (I filled in incomplete values for demo)
data <- data.frame(
  Group = c("AFT-1", "AFT-2", "BEF", "AFT-1", "AFT-2", "BEF"),
  Sample_code = c(10, 10, 10, 11, 11, 11),
  Sex = c("F", "F", "M", "F", "M", "F"),
  GENE_a = c(22.88006, 23.46812, 24.23213, 24.10917, 24.55109, 23.87654),
  GENE_b = c(23.24577, 23.08145, 23.19317, NA, NA, 22.98765),
  GENE_c = c(20.00031, 19.57906, 20.27364, 19.87654, 20.12345, 19.98765)
)

# Check for missing values across genes/samples
colSums(is.na(data))

# Handle missing values (choose one approach):
# Option 1: Remove samples with any missing values (good if missingness is low)
data_clean <- na.omit(data)

# Option 2: Impute missing gene values (better for preserving sample size)
library(mice)
gene_cols <- data[,4:ncol(data)]
imputed_genes <- complete(mice(gene_cols, m=1, method="mean")) # Mean imputation
data_imputed <- cbind(data[,1:3], imputed_genes)

Step 2: Exploratory Data Analysis (EDA) #

Before diving into stats, let's visualize data distribution and sample clustering:

library(ggplot2)
library(tidyr)
library(patchwork)

# Reshape data to long format for plotting
data_long <- pivot_longer(data_imputed, cols = starts_with("GENE"), 
                          names_to = "Gene", values_to = "Expression")

# Boxplot for a single gene (e.g., GENE_a) across groups
boxplot_plot <- ggplot(data_long[data_long$Gene == "GENE_a",], 
                       aes(x=Group, y=Expression, fill=Group)) +
  geom_boxplot(alpha=0.7) +
  labs(title="GENE_a Expression by Treatment Group", 
       x="Group", y="Expression Level") +
  theme_minimal()

# PCA to check sample clustering
expr_matrix <- as.matrix(data_imputed[,4:ncol(data_imputed)])
pca_result <- prcomp(t(expr_matrix), scale. = TRUE) # Transpose for sample-wise PCA

# Prepare PCA data for plotting
pca_df <- as.data.frame(pca_result$x[,1:2])
pca_df$Group <- data_imputed$Group
pca_df$Sex <- data_imputed$Sex

pca_plot <- ggplot(pca_df, aes(x=PC1, y=PC2, color=Group, shape=Sex)) +
  geom_point(size=3) +
  labs(title="PCA of Gene Expression", 
       x=paste0("PC1 (", round(pca_result$sdev[1]^2/sum(pca_result$sdev^2)*100,1), "%)"),
       y=paste0("PC2 (", round(pca_result$sdev[2]^2/sum(pca_result$sdev^2)*100,1), "%)")) +
  theme_minimal()

# Combine plots
boxplot_plot + pca_plot

Step 3: Differential Expression Analysis #

For multi-group comparisons, we have two reliable methods:

Method 1: ANOVA + Tukey's Post-Hoc Test (Single Gene) #

Great for testing individual genes:

# Test GENE_a across groups
anova_result <- aov(Expression ~ Group, data = data_long[data_long$Gene == "GENE_a",])
summary(anova_result) # Check overall group difference

# Post-hoc test to find which groups differ
tukey_result <- TukeyHSD(anova_result)
print(tukey_result)

Method 2: Limma (Batch Analysis for All Genes) #

The gold standard for transcriptome data—handles unbalanced sample sizes (like your AFT-2 group) and can include covariates (e.g., Sex):

library(limma)

# Build design matrix (include Sex as a covariate if it affects expression)
design <- model.matrix(~ 0 + Group + Sex, data = data_imputed)
colnames(design) <- gsub("Group", "", colnames(design)) # Clean column names

# Fit linear model
fit <- lmFit(expr_matrix, design)

# Define contrasts (all pairwise group comparisons)
contrast_matrix <- makeContrasts(
  BEF_vs_AFT1 = BEF - AFT-1,
  BEF_vs_AFT2 = BEF - AFT-2,
  AFT1_vs_AFT2 = AFT-1 - AFT-2,
  levels = design
)

# Calculate differential expression
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

# Get results for BEF vs AFT-1 (adjust number=200 to include all genes)
de_results <- topTable(fit2, coef="BEF_vs_AFT1", number=200, adjust="fdr")

# Filter significant genes (adjust thresholds as needed)
sig_de_genes <- de_results[de_results$adj.P.Val < 0.05 & abs(de_results$logFC) > 1,]

Step 4: Visualize Differential Expression #

A volcano plot helps highlight significant genes:

# Add significance label
de_results$Significance <- ifelse(de_results$adj.P.Val < 0.05 & abs(de_results$logFC) > 1,
                                  "Significant", "Not Significant")

# Plot volcano
ggplot(de_results, aes(x=logFC, y=-log10(adj.P.Val), color=Significance)) +
  geom_point(alpha=0.7) +
  scale_color_manual(values=c("gray50", "firebrick")) +
  geom_vline(xintercept=c(-1, 1), linetype="dashed", color="black") +
  geom_hline(yintercept=-log10(0.05), linetype="dashed", color="black") +
  labs(title="Volcano Plot: Pre-Treatment vs 1 Month Post-Treatment",
       x="Log2 Fold Change", y="-Log10 Adjusted P-Value") +
  theme_minimal()

Key Notes: #

• If your data is from RNA-seq, make sure it's normalized (e.g., TMM) before analysis.
• Adjust significance thresholds (p-value, logFC) based on your study's needs.
• For small sample sizes (like AFT-2), consider using empirical Bayes methods (which limma does) to improve power.

内容的提问来源于stack exchange，提问作者Reza