蛋白ID为行、观测为列的数据框组间统计分析R代码咨询

Hey Eva, let's break down the steps to analyze your protein data (3 groups, 20 samples each) in R. I'll walk you through data reshaping, statistical tests, batch processing, and visualization—all with reproducible code that you can adapt to your actual dataset.

1. First, Simulate Example Data (Match Your Structure) #

Let's create a dummy dataset that mirrors your setup (protein IDs as rows, samples as columns) so you can test the code right away:

# Load core tidyverse packages (for data wrangling + visualization)
library(tidyverse)

# Set seed for reproducible random data
set.seed(123)

# Create wide-format data: 10 proteins, 3 groups × 20 samples each
protein_data <- data.frame(
  Protein_ID = paste0("P", 1:10),
  # Group 1 samples (baseline expression)
  paste0("G1_", 1:20) %>% set_names() %>% map_dbl(~rnorm(1, mean = 20, sd = 3)),
  # Group 2 samples (higher average expression)
  paste0("G2_", 1:20) %>% set_names() %>% map_dbl(~rnorm(1, mean = 25, sd = 3)),
  # Group 3 samples (intermediate expression)
  paste0("G3_", 1:20) %>% set_names() %>% map_dbl(~rnorm(1, mean = 22, sd = 3))
)

# Check the first few rows to confirm structure
head(protein_data)

2. Reshape Data from Wide to Long Format #

Most statistical tools in R work better with long-format data (one row per protein-sample pair). We'll also extract group labels from your sample names:

long_data <- protein_data %>%
  pivot_longer(
    cols = -Protein_ID,  # Keep Protein_ID as our identifier column
    names_to = "Sample",
    values_to = "Expression"
  ) %>%
  # Pull group name from sample IDs (e.g., "G1_1" → "G1")
  mutate(Group = str_extract(Sample, "^G\\d"))

# Verify the reshaped data
head(long_data)

3. Check Statistical Assumptions #

Before running tests like ANOVA, we need to confirm two key assumptions:

• Normality: Expression values within each group follow a normal distribution
• Homogeneity of variances: Variances across groups are roughly equal

Here's how to check these for a single protein (e.g., P1):

# Filter data for protein P1
p1_data <- long_data %>% filter(Protein_ID == "P1")

# Test normality per group (Shapiro-Wilk test)
p1_data %>%
  group_by(Group) %>%
  summarize(shapiro_p_value = shapiro.test(Expression)$p.value)

# Test variance homogeneity (Bartlett test)
bartlett.test(Expression ~ Group, data = p1_data)

• If p-values > 0.05, assumptions are met → use one-way ANOVA
• If assumptions are violated → use Kruskal-Wallis test (non-parametric alternative)

4. Run Batch Tests for All Proteins #

We'll use dplyr to automate tests for every protein in your dataset:

Option 1: One-Way ANOVA + Tukey Post-Hoc Test #

Use this if your data meets ANOVA assumptions:

# Run ANOVA for each protein
anova_results <- long_data %>%
  group_by(Protein_ID) %>%
  summarize(
    anova_model = list(aov(Expression ~ Group, data = cur_data())),
    anova_p_value = summary(anova_model[[1]])[[1]][["Pr(>F)"]][1],
    .groups = "drop"
  )

# Add Tukey post-hoc results for proteins with significant ANOVA findings
tukey_posthoc <- anova_results %>%
  filter(anova_p_value < 0.05) %>%
  mutate(
    tukey_test = map(anova_model, ~TukeyHSD(.x)),
    tukey_p_values = map(tukey_test, ~.x$Group[, "p adj"])
  )

# View the first few results
head(anova_results)

Option 2: Kruskal-Wallis Test (Non-Parametric) #

Use this if your data violates ANOVA assumptions:

kruskal_results <- long_data %>%
  group_by(Protein_ID) %>%
  summarize(
    kruskal_p_value = kruskal.test(Expression ~ Group, data = cur_data())$p.value,
    .groups = "drop"
  )

# View results
head(kruskal_results)

5. Visualize Group Differences #

Let's make a boxplot to highlight expression differences for a significant protein:

# Plot protein P2 (which should have a low ANOVA p-value in our dummy data)
long_data %>%
  filter(Protein_ID == "P2") %>%
  ggplot(aes(x = Group, y = Expression, fill = Group)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5, size = 1.5) +
  labs(title = "Protein P2 Expression Across Groups",
       x = "Experimental Group", y = "Expression Level") +
  theme_minimal()

Bonus: Adjust for Multiple Comparisons #

Since we're testing dozens/hundreds of proteins, always adjust p-values to avoid false positives:

# Add FDR-adjusted p-values to ANOVA results
anova_results <- anova_results %>%
  mutate(adj_p_value = p.adjust(anova_p_value, method = "fdr"))

head(anova_results)

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