Python 3中如何用for循环对两组基因表达数据批量执行t检验？

Got it, let's walk through how to automate t-tests across all your gene columns without writing repetitive code. I'll cover both Python and R since they're the most common tools for this kind of bioinformatics task.

Python Implementation (Using Pandas & SciPy) #

First, assuming you've already split your 50-row dataset into two 25-row groups (long_term and short_term) as pandas DataFrames:

import pandas as pd
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

# Example: Split your original DataFrame (replace with your actual data)
# original_df = pd.read_csv("your_data.csv")
# long_term = original_df.iloc[:25, :]  # First 25 rows = long-term survivors
# short_term = original_df.iloc[25:, :] # Last 25 rows = short-term survivors

# Initialize a list to store results
t_test_results = []

# Loop through each gene column
for gene in long_term.columns:
    # Extract expression values for the current gene from both groups
    long_expr = long_term[gene]
    short_expr = short_term[gene]
    
    # Run independent samples t-test
    # Use equal_var=False for Welch's t-test if variances are unequal
    t_stat, p_val = ttest_ind(long_expr, short_expr, equal_var=True)
    
    # Append results to the list
    t_test_results.append({
        "Gene": gene,
        "t-statistic": round(t_stat, 4),
        "p-value": round(p_val, 4)
    })

# Convert results to a DataFrame for readability
results_df = pd.DataFrame(t_test_results)

# Add multiple testing correction (FDR-BH) to avoid false positives
results_df["p-value_corrected"] = round(
    multipletests(results_df["p-value"], alpha=0.05, method="fdr_bh")[1],
    4
)
results_df["Significant (α=0.05)"] = results_df["p-value_corrected"] < 0.05

print(results_df)

Notes for Python: #

• If your data is a numpy array instead of a DataFrame, loop over column indices instead of column names:

  for i in range(long_term_arr.shape[1]):
      long_expr = long_term_arr[:, i]
      short_expr = short_term_arr[:, i]
      t_stat, p_val = ttest_ind(long_expr, short_expr)
      # Store results as before
  

• Adjust equal_var based on whether your groups have equal variance (test with Levene's test if unsure).

R Implementation #

If you're working in R, here's an equivalent workflow:

# Example: Split your original data frame
# original_df <- read.csv("your_data.csv")
# long_term <- original_df[1:25, ]
# short_term <- original_df[26:50, ]

# Initialize a list to hold results
results_list <- list()

# Loop through each gene column
for (gene in colnames(long_term)) {
    # Extract expression values
    long_expr <- long_term[[gene]]
    short_expr <- short_term[[gene]]
    
    # Run t-test (use var.equal=FALSE for Welch's test)
    t_test_output <- t.test(long_expr, short_expr, var.equal = TRUE)
    
    # Store results in the list
    results_list[[gene]] <- data.frame(
        Gene = gene,
        t_statistic = round(t_test_output$statistic, 4),
        p_value = round(t_test_output$p.value, 4),
        degrees_of_freedom = t_test_output$parameter,
        stringsAsFactors = FALSE
    )
}

# Combine list into a single data frame
results_df <- do.call(rbind, results_list)

# Add FDR correction
results_df$p_value_corrected <- round(p.adjust(results_df$p_value, method = "fdr"), 4)
results_df$Significant <- results_df$p_value_corrected < 0.05

print(results_df)

Key Considerations: #

• Assumptions: Independent t-tests require normality and (for the standard test) equal variance. If your data violates these, consider Welch's t-test (as noted above) or non-parametric alternatives like the Mann-Whitney U test.
• Multiple Testing: Since you're testing 30 genes, correcting for multiple comparisons is critical to avoid overinterpreting false positives. FDR-BH is a common choice for this kind of genomic data.

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