RNA-seq 复现（一）

母文献：Labonté B, et al. Sex-specific transcriptional signatures in human depression. Nat Med. 2017
公开数据库：GSE102556

一、复现背景和目的（AI撰写）

抑郁症（Major Depressive Disorder, MDD）在流行病学、临床表现及治疗反应上均表现出显著的性别差异。既有转录组学研究提示，大脑皮层中与神经活动和突触可塑性相关的基因在抑郁症中可能发生系统性改变，但这些改变是否在不同性别中具有不同的分子特征仍不清楚。

本复现研究基于公开 RNA-seq 数据集 GSE102556，系统复现并验证原研究中关于性别分层差异表达（sex-stratified differential expression） 的核心结论，重点关注活动依赖型即刻早期基因（Immediate Early Genes, IEGs） ，具体包括NPAS4、CCN1、EGR4、FOSB、ARC、EGR1、EGR2、FOS、JUN、JUNB九个基因。

二、复现所用样本及材料

数据集：GEO数据库的 GSE102556 数据集，具体包括 GSE102556_raw_counts_GRCh38.p13 和 GSE102556-GPL11154_series_matrix 两套文件

GSE102556_raw_counts_GRCh38.p13：基因表达矩阵，行名为基因名称（使用Entrez ID进行编号），列名为GSM样本编号，数值为每一个基因在对应样本体内的表达量。
- 具体格式：为表格形式。
- 基因表达矩阵不提供样本相关的统计学数据（如性别、诊断、脑区等），如需要进一步分组，应该结合series数据来进行分析。（这也就是为什么之后需要读取两个文件并做数据对齐的原因）。

GSE102556-GPL11154_series_matrix：GSM样本的元数据，为每一个样本对应的相关信息的标注，存储了每一个样本的属性（性别、诊断、脑区等相关数据）、实验信息、平台信息、处理流程等信息

具体格式：为 series matrix 格式，其形式为 “! Sampel_characteristics_ch1 = xxx : xxx…” 字样的重复，如：

!Sample_characteristics_ch1 = tissue: Orbitofrontal (OFC; BA11)
!Sample_characteristics_ch1 = gender: Male
!Sample_characteristics_ch1 = phenotype: CTRL
!Sample_characteristics_ch1 = age: 47
!Sample_characteristics_ch1 = rin: 7.9
...

其特点为，同一个字段重复很多行；并且样本特征是以 "xxx : xxx"的形式而非表格形式展现，这意味着我们在进行分析之前，需要将series matrix格式整理成可分析的形式，并将series matrix和基因表达矩阵的数据对齐。

研究样本类型：人类死后脑组织 RNA-seq，样本中所含包括六个脑区：OFC（眶额皮层，Orbitofrontal Cortex）、aINS （前岛叶 / 前脑岛，Anterior Insula）、dlPFC（背外侧前额叶皮层Dorsolateral Prefrontal Cortex）、vmPFC（腹内侧前额叶皮层，Ventromedial Prefrontal Cortex）、NAc（伏隔核，Nucleus Accumbens）、vSUB（腹侧下托，Ventral Subiculum）
数据基本概括
- Count Martix : 39157 genes x 261 samples
- 样本情况：男 139 ：女 122；对照 122 ：确诊（MDD）139

三、具体数据分析过程

文件读取与数据预处理：
a. 读取 raw_count 文件，进行初步清洗并大致预览样本概况
b. 读取 series 文件（这里最好确保找对文件）
c. 对 series 文件的处理：提取样本ID（GSM编号）、提取数据并将其转化为矩阵 （难点一） ，转化为数据框，并将样本ID与GSM编号对应
d. 进一步处理 series 数据框：提取本次分析所需要的所有变量并将其转化为可用于分析的形式：性别（Male/Female）、诊断（CTRL/MDD）、年龄、pmi、rin（以上三个转化为数值作为分析协变量）、脑区来源（6个）、饮酒状况和用药状况（yes/no/NA）
e. 因子转换：将饮酒和用药状况数值化（0/1/NA）并转化为因子形式（性别及诊断状况同理）：此为还原论文变量操作，真实研究可能不同。
数据对齐：将处理后的series数据框和表达数据raw count文件对齐，检查有无多余GSM编号，去除多于编号的数据（可列出被去除的样本特征方便复盘）。
构造模型并进行统计学分析：
a. 复现论文的数据筛选：筛选出在至少80%的样本里表达量高于5的基因，以避免低表达基因造成噪声（问题：可能会造成某些性别特异表达基因的误筛）
b. 因子设置（实际上前面已有因子化设置）：定义参考基线：性别变量以 Male 作为reference line，诊断变量以 CTRL 作为reference line，后面的模型构建的系数全部 Male + CTRL 作为“零点”。
c. 矩阵设计：每个基因的线性模型为：

$Y_i = \beta_0 + \beta_{\text{sex}}\cdot \mathbb{I}(\text{sex}=\text{Female}) + \beta_{\text{dx}}\cdot \mathbb{I}(\text{dx}=\text{MDD}) + \beta_{\text{int}}\cdot \mathbb{I}(\text{sex}=\text{Female})\,\mathbb{I}(\text{dx}=\text{MDD}) + \beta_{\text{age}}\cdot \text{age} + \beta_{\text{rin}}\cdot \text{rin} + \beta_{\text{alcool}}\cdot \text{alcool} + \beta_{\text{med}}\cdot \text{med} + \varepsilon$

其中，$Y_i$ 是是 voom 处理后的带权重的基因表达量。
在这一模型下，各个系数的直觉解释为：在以 **Male + CTRL** 组别为参照的情况下，$β_0$ 表示 MDD+CTRL 组的平均表达（即在协变量为0或者“协变量中心未做”时的基线）；$β_{sex}$ 表示在 CTRL 条件下，Female vs Male 的差异；$β_{sex}$ 表示在 Male 条件下，MDD vs CTRL 的差异（即男性内部差异）；$β_{int}$ 表示 MDD 效益在女性条件下比在男性条件下额外多增加了多少（也就是本次复现主要的目标效应值），即：$(Female MDD - CTRL) - (Male MDD - CTRL)$ 
d. 使用 voom + limma 建模：通过 voom ，将 raw count 数据转换为近似正态的表达尺度即 logCPM 化，并通过**估计均值-方差关系**，给每个观测一个加权：即**低表达、方差大**的基因的权重更低，因此，后续的 liFit 实际上是**加权线性回归**。此外，在 liFit 分析结束后，对每个基因的残差方差做经验贝叶斯收缩，提升结果小样本/多检验下的稳定性，这也是 limma 的经典步骤 。
e. 进行三组跨组检验：$β_{dx}$ 检验男性内部 MDD vs CTRL 的差异；$β_{dx} + β_{int}$ 检验女性内部 MDD vs CTRL 的差异；$β_{int}$ 则检验**性别是否调剂了MDD的效应**。该检验显著意味着对应的基因在 MDD 的差异存在性别差异性（方向/大小不同）。
f. 对来自六个脑区的样本分别展开检验并输出初步结果：各脑区男性内部、女性内部以及跨性别表达差异显著基因的数量

对特定基因组结果（IEGs）的筛选
a. 指定 IEGs 基因组的基因（此为复现需要，在实际研究中会根据结果做相应调整）
b. 在六组结果中筛选出相应的基因及对应的三组检验的结果（男性、女性、跨性别），先选择前两项性别内检验结果，将总脑区结果制成一张长表（男女分列）+一张总表（男女并列）
c. 单独筛选出第三组的跨性别检验的 int 结果，每个脑区各输出一张表。
d. 各脑区跨性别检验组中 p value < 0.05 的结果即为我们的目标基因（注：更严格的条件应选择 padj value 作为检验结果指标）
数据可视化（火山图/热图等）

四、代码文件与结果解读

代码文件分为两部分，一部分为数据处理代码，另一部分为数据可视化代码。

重点展示数据处理代码：

# ============================================================
# GSE102556 (GPL11154) RNA-seq raw counts: voom+limma differential expression
# - Overall: MDD vs CTRL (adjusting covariates)
# - Sex-stratified: Male-only and Female-only analyses
# Output: CSV table for selected immediate early genes (IEGs)
# ============================================================

# ---------------------------
# 0) Packages and setup
# ---------------------------

suppressPackageStartupMessages({
  library(data.table)
  library(stringr)
  library(DESeq2)
  library(edgeR)
  library(limma)
  library(AnnotationDbi)
  library(org.Hs.eg.db)
})
  
setwd("D:/GEO")
  
# ---------------------------
# 1) File paths (EDIT THESE)
# ---------------------------

counts_file <- "D:/GEO/data/GSE102556_raw_counts_GRCh38.p13_NCBI.tsv.gz"
series_file <- "D:/GEO/data/GSE102556-GPL11154_series_matrix.txt.gz"
  
# ---------------------------
# 2) Helpers
# ---------------------------
  
# Convert strings like "age: 42" -> "42"
strip_prefix <- function(x) {
  x <- as.character(x)
  x <- sub("^[^:]+:\\s*", "", x)
  trimws(x)
}
  
# Safe SYMBOL <-> ENTREZ mapping wrappers
symbol_to_entrez <- function(symbols) {
  out <- mapIds(
    org.Hs.eg.db,
    keys = symbols,
    keytype = "SYMBOL",
    column = "ENTREZID",
    multiVals = "first"
  )
  unname(out)
}
  
entrez_to_symbol <- function(entrez_ids) {
  out <- mapIds(
    org.Hs.eg.db,
    keys = entrez_ids,
    keytype = "ENTREZID",
    column = "SYMBOL",
    multiVals = "first"
  )
  unname(out)
}
  
# Build limma result with basic filtering + ordering
run_voom_limma <- function(count_mat, coldata) {
  stopifnot(identical(colnames(count_mat), rownames(coldata)))
  
  # Filter lowly expressed genes
  keep <- rowSums(count_mat > 5) >= ceiling(0.8 * ncol(count_mat))
  y <- DGEList(counts = count_mat[keep, , drop=FALSE])
  y <- calcNormFactors(y, method = "TMM")
  
  # Set factor reference levels (EDIT if your data uses different labels)
  coldata$dx  <- factor(coldata$dx,  levels = c("CTRL", "MDD"))
  coldata$sex <- factor(coldata$sex, levels = c("Male", "Female"))
  
  # model matrix with interaction
  design <- model.matrix(
    ~ sex + dx + sex:dx + age_num + rin_num + alcool_f + med_f,
    data = coldata
  )
  colnames(design) <- make.names(colnames(design))
  # print(colnames(design))
  v <- voom(y, design, plot = FALSE)
  fit <- lmFit(v, design)
  fit <- eBayes(fit)
  ct <- makeContrasts(
  male        = dxMDD,
  female      = dxMDD + sexFemale.dxMDD,
  interaction = sexFemale.dxMDD,
  levels = design
  )
  
fit2 <- contrasts.fit(fit, ct)
fit2 <- eBayes(fit2)
  
return(list(
  res_male = topTable(fit2, coef="male", number=Inf, sort.by="P"),
  res_female = topTable(fit2, coef="female", number=Inf, sort.by="P"),
  res_interaction = topTable(fit2, coef="interaction", number=Inf, sort.by="P"),
  keep_genes = rownames(y),
  design = design
  ))
}
  
# ---------------------------
# 3) Load counts
# ---------------------------

expr <- read.delim(
  counts_file,
  header = TRUE,
  row.names = 1,
  stringsAsFactors = FALSE,
  check.names = FALSE
)
  
# Basic sanity + remove all-zero genes
expr <- expr[rowSums(expr) > 0, ]
stopifnot(sum(is.na(expr)) == 0)
  
message("Counts matrix: ", nrow(expr), " genes x ", ncol(expr), " samples")
  
# ---------------------------
# 4) Parse sample metadata from series matrix
# ---------------------------

x <- readLines(gzfile(series_file))
  
# sample IDs (GSM)
gsm <- strsplit(x[grepl("^!Sample_geo_accession", x)], "\t", fixed = TRUE)[[1]][-1]
gsm <- gsub('^"|"$', "", gsm)
  
# all characteristics lines
char_lines <- x[grepl("^!Sample_characteristics_ch1", x)]
  
# convert to matrix: rows = feature lines, cols = samples
char_mat <- do.call(rbind, lapply(char_lines, function(line) {
  v <- strsplit(line, "\t", fixed = TRUE)[[1]][-1]
  gsub('^"|"$', "", v)
}))
  
# derive rownames from the label segment (before ':')
rownames(char_mat) <- sapply(char_lines, function(line) {
  # second field looks like: !Sample_characteristics_ch1  "gender: Male"
  lab <- strsplit(line, "\t", fixed = TRUE)[[1]][2]
  lab <- gsub('^"|"$', "", lab)
  sub(":.*$", "", lab)
})
  
meta <- as.data.frame(t(char_mat), stringsAsFactors = FALSE)
meta$GSM <- gsm
  
# clean column names
colnames(meta) <- gsub('"', "", colnames(meta))
colnames(meta) <- trimws(colnames(meta))
colnames(meta) <- gsub("\\s+", "_", colnames(meta))
  
# pull key covariates
meta$sex <- strip_prefix(meta$gender)
meta$dx  <- strip_prefix(meta$phenotype)
meta$age_num <- as.numeric(strip_prefix(meta$age))
meta$pmi_num <- as.numeric(strip_prefix(meta$pmi))
meta$rin_num <- as.numeric(strip_prefix(meta$rin))
meta$tissue <- strip_prefix(meta$tissue)
meta$alcool_raw <- strip_prefix(meta$alcool)
meta$med_raw <- strip_prefix(meta$medication)
  
# standardize alcohol and medication to 0/1/NA
to01 <- function(x) {
  x <- tolower(trimws(as.character(x)))
  x[x %in% c("", "na", "n/a", "unknown", "unk", "not available")] <- NA
  x[x %in% c("no", "non", "none", "0", "false")] <- "0"
  x[x %in% c("yes", "oui", "1", "true")] <- "1"
  suppressWarnings(as.integer(x))
}
meta$alcool <- to01(meta$alcool_raw)
meta$medication_status <- to01(meta$med_raw)
to_yn_unk <- function(x){
  x <- tolower(trimws(as.character(x)))
  x[x %in% c("", "na", "n/a", "unknown", "unk", "not available")] <- NA
  out <- ifelse(is.na(x), "Unknown",
                ifelse(x %in% c("yes","oui","1","true"), "Yes",
                       ifelse(x %in% c("no","non","0","false"), "No", "Unknown")))
  factor(out, levels=c("No","Yes","Unknown"))
}
meta$alcool_f <- to_yn_unk(meta$alcool_raw)
meta$med_f    <- to_yn_unk(meta$med_raw)
  
# handle the tissue naming inconsistency
map_region <- function(tissue) {
  tissue <- trimws(tissue)
  dplyr::case_when(
    tissue == "Orbitofrontal (OFC; BA11)" ~ "OFC",
    tissue == "Dorsolateral prefrontal cortex (dlPFC; BA8/9)" ~ "dlPFC",
    tissue == "Cingulate gyrus 25 (Cg25)" ~ "vmPFC",   # BA25
    tissue == "Anterior Insula (aINS)" ~ "aINS",
    tissue == "Nucleus Accumbens (Nac)" ~ "NAc",
    tissue == "Subiculum (Sub)" ~ "vSUB",
    TRUE ~ NA_character_
  )
}
meta$region <- map_region(meta$tissue)
stopifnot(!any(is.na(meta$region)))
  
# quick checks
message("Sex counts:"); print(table(meta$sex, useNA = "ifany"))
message("Dx counts:");  print(table(meta$dx,  useNA = "ifany"))
message("Covariate summaries:"); print(summary(meta[, c("age_num","pmi_num","rin_num")]))
message("Alcohol (0/1/NA):"); print(table(meta$alcool, useNA="ifany"))
message("Medication (0/1/NA):"); print(table(meta$medication_status, useNA="ifany"))
message("Region counts:"); print(table(meta$region))
  
# ---------------------------
# 5) Match samples between expr and meta
# ---------------------------

keep <- meta$GSM %in% colnames(expr)
if (!all(keep)) {
  warning("Some meta samples are missing from counts: ", paste(meta$GSM[!keep], collapse = ", "))
}
meta2 <- meta[keep, ]
expr2 <- expr[, meta2$GSM]

# ---- DIAG: which GSM are missing from counts, by region ----
meta$in_counts <- meta$GSM %in% colnames(expr)
message("Total meta samples: ", nrow(meta))
message("Present in counts:  ", sum(meta$in_counts))
message("Missing in counts:  ", sum(!meta$in_counts))
message("Missing-by-region:")
print(table(meta$region[!meta$in_counts]))
if (sum(!meta$in_counts) > 0) {
  missing_df <- meta[!meta$in_counts, c("GSM","region","tissue","sex","dx")]
  fwrite(missing_df, file="results/missing_samples_not_in_counts.tsv", sep="\t")
  message("Wrote: results/missing_samples_not_in_counts.tsv")
}
  
stopifnot(identical(colnames(expr2), meta2$GSM))
stopifnot(!any(duplicated(meta2$GSM)))
  
# Minimal colData for modeling
colData <- meta2[, c("GSM","region","dx","sex","age_num","rin_num","pmi_num",
                     "alcool_f","med_f")]
rownames(colData) <- colData$GSM
  
# Set factor reference levels (EDIT if your data uses different labels)
colData$dx  <- factor(colData$dx,  levels = c("CTRL", "MDD"))
colData$sex <- factor(colData$sex, levels = c("Male", "Female"))
  
# ---------------------------
# 6) Overall limma (all samples)
# ---------------------------
  
dir.create("results", showWarnings = FALSE)
regions <- c("vmPFC","OFC","dlPFC","aINS","NAc","vSUB")
for (rg in regions) {
  message("\n=== Running region: ", rg, " ===")
  cd <- colData[colData$region==rg, , drop=FALSE]
  # 关键：先不要盲目 complete.cases，把缺失情况打印出来
  message("Missingness:");  print(sapply(cd[,c("age_num","rin_num","alcool_f","med_f")], \(v) sum(is.na(v))))
  # 暂时：只剔除 age/rin 缺失
  cd <- cd[!is.na(cd$age_num) & !is.na(cd$rin_num), , drop = FALSE]
  print(table(cd$sex, cd$dx))
  cm <- expr2[, rownames(cd), drop = FALSE]
  fit <- run_voom_limma(cm, cd)
  out_m <- data.frame(entrez_id = rownames(fit$res_male), fit$res_male, row.names = NULL)
  out_f <- data.frame(entrez_id = rownames(fit$res_female), fit$res_female, row.names = NULL)
  out_m$symbol <- entrez_to_symbol(out_m$entrez_id)
  out_f$symbol <- entrez_to_symbol(out_f$entrez_id)
  fwrite(out_m, file = file.path("results", paste0("DE_", rg, "_male.tsv")), sep = "\t")
  fwrite(out_f, file = file.path("results", paste0("DE_", rg, "_female.tsv")), sep = "\t")
  message("Male DEGs p<0.05: ", sum(out_m$P.Value < 0.05, na.rm = TRUE))
  message("Male DEGs FDR<0.05: ", sum(out_m$adj.P.Val < 0.05, na.rm = TRUE))
  message("Female DEGs p<0.05: ", sum(out_f$P.Value < 0.05, na.rm = TRUE))
  message("Female DEGs FDR<0.05: ", sum(out_f$adj.P.Val < 0.05, na.rm = TRUE))
  out_i <- data.frame(entrez_id = rownames(fit$res_interaction), fit$res_interaction, row.names = NULL)
  
  # 加 symbol 注释，保持三个文件一致
  out_m$symbol <- entrez_to_symbol(out_m$entrez_id)
  out_f$symbol <- entrez_to_symbol(out_f$entrez_id)
  out_i$symbol <- entrez_to_symbol(out_i$entrez_id)
  
  fwrite(out_m, file = file.path("results", paste0("DE_", rg, "_male.tsv")), sep = "\t")
  fwrite(out_f, file = file.path("results", paste0("DE_", rg, "_female.tsv")), sep = "\t")
  fwrite(out_i, file = file.path("results", paste0("DE_", rg, "_interaction.tsv")), sep = "\t")
  
  message("Interaction DEGs p<0.05: ", sum(fit$res_interaction$P.Value < 0.05, na.rm = TRUE))
}
  
# ---------------------------
# 7) Summarize IEGs across regions (paper-like table)
# ---------------------------
  
ieg_symbols <- c("NPAS4","CCN1","EGR4","FOSB","ARC","EGR1","EGR2","FOS","JUN","JUNB")  # 可补充或替换
ieg_entrez  <- symbol_to_entrez(ieg_symbols)
ieg_map <- data.frame(symbol = ieg_symbols, entrez_id = ieg_entrez, stringsAsFactors = FALSE)
  
read_res <- function(rg, sex){
  fn <- file.path("results", paste0("DE_", rg, "_", sex, ".tsv"))
  dt <- fread(fn)
  # dt 里你已经有 entrez_id / logFC / P.Value / adj.P.Val / symbol（如果你按我上面加了）
  dt[, region := rg]
  dt[, sex := sex]
  dt
}
  
all_dt <- rbindlist(lapply(regions, function(rg){
  rbind(read_res(rg,"male"), read_res(rg,"female"))
}), fill=TRUE)
  
all_dt2 <- copy(all_dt)
all_dt2[, symbol := NULL]   # 删掉 DE 注释列，避免 merge 产生 symbol.x/y
  
ieg_dt <- merge(ieg_map, all_dt2, by="entrez_id", all.x=TRUE)
  
ieg_dt2 <- ieg_dt[, c("symbol","entrez_id","region","sex","logFC","P.Value","adj.P.Val")]
  
library(data.table)
ieg_dt2 <- as.data.table(ieg_dt2)
  
fwrite(ieg_dt2, file = "results/IEG_summary_long.tsv", sep = "\t")
  
# 宽表：每个 region 一张（male/female并列）
ieg_wide <- dcast(
  ieg_dt2,
  symbol + entrez_id + region ~ sex,
  value.var = c("logFC","P.Value","adj.P.Val")
)
fwrite(ieg_wide, file="results/IEG_summary_wide.tsv", sep="\t")
  
message("Wrote: results/IEG_summary_long.tsv and results/IEG_summary_wide.tsv")
setDT(ieg_wide)
  
#---------------------------
# 8) IEG counts ———— 总览各脑区 IEG 差异表达情况
# ---------------------------
  
# 论文口径：nominal p<0.05
ieg_counts_p <- ieg_wide[, .(
  n_female_p = sum(P.Value_female < 0.05, na.rm=TRUE),
  n_male_p   = sum(P.Value_male   < 0.05, na.rm=TRUE)
), by=region][order(region)]
  
ieg_counts_p
  
# 严格口径：FDR<0.05
ieg_counts_fdr <- ieg_wide[, .(
  n_female_fdr = sum(adj.P.Val_female < 0.05, na.rm=TRUE),
  n_male_fdr   = sum(adj.P.Val_male   < 0.05, na.rm=TRUE)
), by=region][order(region)]

ieg_counts_fdr