主题
08 · 单细胞定位与伪 bulk
⚠ 代码目录与 01–07 不同。01–07 讲
/home/research/multiomics-mr/(多组学 MR); 08–10 讲/home/research/dr-amd-targets/analysis/(功能验证),是独立项目目录。本页回答组会最常被问的一句:"你单细胞到底怎么做的、凭什么说 MERTK 定位在小胶质?" 所有参数以实际脚本为准,不是通用教程值。
对应脚本:11_dr_scrna.R(DR)、12_amd_scrna_rpe.R(AMD)、13_amd_pseudobulk.R(严格确认)
数据与样本
| DR | AMD | |
|---|---|---|
| GEO | GSE165784 纤维血管膜(PDR/PVR) | GSE230348 RPE/脉络膜 |
| 样本 | 6 例(4 例 dense tsv + 2 例 10X tar) | mRPE/tRPE 多例,读 amd_manifest.csv 定分组 |
| 分组 | PDR(含 F47) vs PVR_RRD | AMD vs Normal(由 disease 字段 grepl "AMD") |
| 额外 | — | region = Macula / Peripheral(由样本名首字母 m/t 判) |
QC 门槛(两套数据不同,别记混)
r
# DR (11_dr_scrna.R)
nFeature_RNA >= 200 & nFeature_RNA <= 6000 & percent.mt < 20
# AMD (12_amd_scrna_rpe.R)
nFeature_RNA >= 200 & nFeature_RNA <= 7000 & percent.mt < 25建对象时两边都是 min.cells = 3, min.features = 200。
AMD 侧有个必须主动交代的处理:CAP <- 2500,每样本随机下采样到 2500 个细胞(内存限制)。 另外 QC 后 样本细胞数 < 50 的样本整个丢弃。下采样对"定位在哪种细胞"稳健,但任何涉及绝对细胞数/比例的结论都要标注。
整合与聚类(两边一致)
r
NormalizeData → FindVariableFeatures(nfeatures = 2000) → ScaleData
→ RunPCA(npcs = 30) → RunHarmony(group.by.vars = "sample")
→ RunUMAP(reduction = "harmony", dims = 1:30)
→ FindNeighbors(dims = 1:30) → FindClusters(resolution = 0.5)set.seed(1)。批次变量是 sample(不是 disease,避免把疾病效应整合掉)。
★ 细胞注释:不是手工点的,是模块评分 argmax
这是最容易被追问的一点。没有人工指认细胞类型,用的是标准 marker panel 打模块分、再按 cluster 均分取最大值自动指派:
r
markers <- list(
Microglia_Myeloid = c("PTPRC","AIF1","CX3CR1","P2RY12","CSF1R","C1QA","C1QB","ITGAM","CD68","LYZ"),
Endothelial = c("PECAM1","VWF","CLDN5","CDH5","FLT1"),
Mural_Fibroblast = c("PDGFRB","ACTA2","RGS5","COL1A1","DCN","LUM"),
T_NK = c("CD3D","CD3E","NKG7","GNLY"),
RPE = c("RPE65","BEST1","RLBP1","TYR"),
Photoreceptor = c("RCVRN","RHO","PDE6A"),
Muller_Glia = c("GLUL","VIM","SLC1A3","SOX9"))
for (lin in names(markers)) merged <- AddModuleScore(merged, features = list(g), name = paste0("sc_", lin), seed = 1)
cl_mean <- aggregate(...[, sc_cols], by = list(cl = merged$seurat_clusters), FUN = mean)
merged$celltype <- cl2lin[...] # 每个 cluster 取模块分均值最大的那一系优点:可复现、无主观。 局限(要如实讲):argmax 是强制指派——每个 cluster 必被分到某一系,没有 "Unknown" 出口。若某 cluster 是双胞体或未覆盖的细胞类型,会被硬塞进最近的一系。正式发表建议补一张 cluster × 模块分热图作为注释可信度证据。
靶点定位的两个读出
r
avg <- AverageExpression(merged, features = targets, group.by = "celltype", assays = "RNA")$RNA
pct <- tapply(FetchData(merged, vars = g)[,1] > 0, merged$celltype, function(x) 100*mean(x))即 平均表达 + 阳性率,两者一起看。PPT 里"MERTK 平均表达 0.71、阳性率 25% vs 其他细胞类型 7%"就出自这里 (target_avgexpr_by_celltype.csv / target_pctexpr_by_celltype.csv)。
★ 2026-07-26 独立复算查出两处表述错误(数值本身无误)
tools/verify_singlecell.R 脱离流水线重算后发现:
- 0.71 不是"模块评分",是平均表达。
AddModuleScore在本脚本里只用于细胞类型 marker panel 的自动注释,从未用于 MERTK。原表述会让人以为 MERTK 有一个模块分。 - 对照组标签写错:原写"阳性率 25% vs 整体 8%"。实测全细胞池阳性率是 22.4%,而 7.0% 才是「非小胶质/髓系细胞」的阳性率。所以正确表述是 25% vs 其他细胞类型 7%。
⚠️ 口径提醒:AverageExpression 对 log 归一化数据是先 expm1 再求均值(0.7126);若直接对 log 值求均值会得到 0.3063,两者不是同一个量,复核时容易误判为"算错了"。
查的靶点:MERTK, SPRY2, AGER, TNFRSF10A, CFH, TGFB1(取与矩阵行名的交集,缺的自动跳过)。
★ 伪 bulk:为什么不用细胞级 Wilcoxon
13_amd_pseudobulk.R 开头的注释就写明了理由:
avoiding the anticonservative cell-level Wilcoxon (Squair et al. 2021 Nat Commun)
细胞级检验把同一个体的上万个细胞当独立样本,p 值被严重高估。所以 AMD vs Normal 的差异走每样本伪 bulk:
r
cnt <- GetAssayData(rpe, assay = "RNA", layer = "counts") # 原始 counts,不是 normalized
pb <- sapply(levels(factor(samp)), function(s) Matrix::rowSums(cnt[, samp == s]))
mertk_cpm <- 1e6 * pb["MERTK", ] / colSums(pb) # CPM
df <- df[df$n_cells >= 20, ] # 每样本至少 20 个 RPE 细胞
w <- wilcox.test(a, n); t <- t.test(log1p(a), log1p(n)) # 两种检验都报判定:样本级 Wilcoxon + Welch t(log1p)两个都报,不挑好看的。 PPT 结果图 7 的"AMD vs 正常无差异"就是这里出来的——这是阴性结果,如实保留。 ⚠ 2026-07-27 更新:样本分组补全后(22→34 份),Wilcoxon p=0.076、log2FC +0.71, 仍未达显著,但方向为 AMD 侧偏高(与功能丧失假说相反)。旧的 p=0.82 是在 22 份、 且 40.4% 细胞被误判为「分组未知」的情况下算出的,已作废。
这一层能证明什么、不能证明什么
| 能 | 不能 |
|---|---|
| 靶点表达在哪种细胞 | ❌ 不能证明因果 |
| 疾病 vs 正常在该细胞里有无表达差异(伪 bulk) | ❌ 不能替代共定位——MERTK 在这里定位清晰,但 coloc H4=0.01 |
| 为选细胞系提供依据(定位到 RPE → 选 ARPE-19) | ❌ 样本量小(DR 6 例、AMD 下采样),效应量不可外推 |
★ 写作红线:单细胞定位是功能层证据,与因果遗传层(MR+共定位)分开陈述,不可混为"多层收敛因此因果"。
输出物
results/scrna_dr/ dr_fvm_seurat.rds, celltype_by_sample.csv,
target_avgexpr_by_celltype.csv, target_pctexpr_by_celltype.csv
figures/ umap_celltype.png, umap_sample.png, dotplot_targets.png, feature_*.png
results/scrna_amd/ amd_rpe_seurat.rds, MERTK_RPE_pseudobulk_by_sample.csv
figures/ MERTK_RPE_pseudobulk_box.png全部图 dpi = 600、英文标签,符合结果图统一作图规范。
下一页 → 09 · 细胞通讯与活化轨迹