实操：GSE163973单细胞数据处理流程

一、写在前面

在处理单细胞转录组数据时，如何精准地剔除双细胞（Doublets）并整合多样本数据，是后续所有高级分析的基础。本篇教程基于GSE163973数据集，手把手带大家走一遍单细胞分析的标准流程：从基础的质控、利用 DoubletFinder 识别并去除伪细胞，到使用Harmony算法进行多样本批次校正。

更进一步，我们不仅会演示如何对全细胞类型进行注释，还会聚焦于成纤维细胞（FIB）子集，展示如何通过亚聚类分析（Sub-clustering）挖掘组织内部的细微异质性。最后，教程还整合了 CellPhoneDB 的跨语言调用（Python 运算 + R 可视化），并绘制了细胞通讯热图。

如果需要单细胞数据分析指导、生信热点全文复现、自测数据个性化分析辅导、常态化实验学习，欢迎联系[Biomamba_zhushou]。

二、实操流程

说明：

本Rmarkdown文件包含两部分代码：

1 运行环境准备（R）

1.1导入必要的 R 包

suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(remotes)
library(Seurat)
library(data.table)
library(SeuratDisk)
library(patchwork)
library(ggtext)
library(gghalves)
library(cowplot)
library(reshape2)
library(stringr)
library(reticulate)
library(DoubletFinder)
library(harmony)
library(clustree)
library(tibble)
library(RColorBrewer)
library(parallel)
library(pbapply)
library(future)
library(future.apply)
})

3 读取原始数据与合并（R）

## 读取数据
data_path <– "/gpfs/hpc/home/GSE163973_RAW"
# 获取所有子文件夹路径
folders <– list.dirs(data_path, recursive = FALSE, full.names = TRUE)
# 读取并合并数据
seurat_list <– lapply(folders, function(folder) {
orig_id <– strsplit(basename(folder), "_")[[1]][1]
data <– Read10X(data.dir = folder)
CreateSeuratObject(counts = data, project = orig_id, min.cells = 3, min.features = 200)
})
merged_seurat <– merge(seurat_list[[1]], seurat_list[-1])
rm(seurat_list)
# 计算线粒体比例
merged_seurat[["percent.mt"]] <– PercentageFeatureSet(merged_seurat, pattern = "^MT-")
VlnPlot(merged_seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

4 以样本为单位进行低质量细胞过滤（R）

# 过滤数据阈值（按样本）
thresholds <– list(
'KL1' = list(nFeature_RNA = c(200, 5000), percent.mt = 15),
'KL2' = list(nFeature_RNA = c(200, 5000), percent.mt = 15),
'KL3' = list(nFeature_RNA = c(200, 5000), percent.mt = 15),
'NS1' = list(nFeature_RNA = c(200, 4000), percent.mt = 15),
'NS2' = list(nFeature_RNA = c(200, 4000), percent.mt = 15),
'NS3' = list(nFeature_RNA = c(200, 4000), percent.mt = 15)
)
thresholds <– bind_rows(thresholds, .id = "orig.ident") %>%
group_by(orig.ident) %>%
summarise(
nFeature_RNA_min = min(nFeature_RNA),
nFeature_RNA_max = max(nFeature_RNA),
percent_mt_max = max(percent.mt),
.groups = "drop"
)
meta_data <– merged_seurat@meta.data %>% left_join(thresholds, by = "orig.ident")
keep_cells <– with(meta_data,
nFeature_RNA >= nFeature_RNA_min &
nFeature_RNA <= nFeature_RNA_max &
percent.mt <= percent_mt_max)
merged_seurat <– merged_seurat[, keep_cells]

5 使用DoubletFinder分样本处理双联体（R）

# 常规预处理
merged_seurat <– NormalizeData(merged_seurat)
merged_seurat <– FindVariableFeatures(merged_seurat, selection.method = "vst", nfeatures = 2000)
all.genes <– rownames(merged_seurat)
merged_seurat <– ScaleData(merged_seurat, features = all.genes)
merged_seurat <– RunPCA(merged_seurat, features = VariableFeatures(object = merged_seurat))
# 拆分为每个 orig.ident
seu_list <– SplitObject(merged_seurat, split.by = "orig.ident")
names(seu_list) <– paste0("seu_", names(seu_list))
# 定义每个样本的处理任务：UMAP/Neighbors/Clusters + DoubletFinder
task <– function(seu) {
library(Seurat)
dim.usage <– 15
seu <– RunUMAP(seu, dims = 1:dim.usage)
seu <– FindNeighbors(seu, dims = 1:dim.usage)
seu <– FindClusters(seu, resolution = 0.8)
sweep.res.list <– paramSweep(seu, PCs = 1:dim.usage, sct = FALSE)
sweep.stats <– summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <– find.pK(sweep.stats)
# 选择 MeanBC 最大的 pK
p <– as.numeric(as.vector(bcmvn[bcmvn$MeanBC==max(bcmvn$MeanBC),]$pK))
homotypic.prop <– modelHomotypic(seu$seurat_clusters)
nExp_poi <– round(0.08*ncol(seu))
nExp_poi.adj <– round(nExp_poi*(1–homotypic.prop)) # 异源双细胞数估计，供参考
seu <– doubletFinder(seu, PCs = 1:dim.usage, pN = 0.25, pK = p, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
colnames(seu@meta.data)[ncol(seu@meta.data)] <– "doublet_info"
seu
}
# 并行处理
results <– future_lapply(seu_list, task, future.seed = TRUE)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 8251

## Number of edges: 273151

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8812

## Number of communities: 20

## Elapsed time: 0 seconds

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 7632

## Number of edges: 251945

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8829

## Number of communities: 19

## Elapsed time: 0 seconds

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 7332

## Number of edges: 244452

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8771

## Number of communities: 19

## Elapsed time: 0 seconds

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 7122

## Number of edges: 236543

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8734

## Number of communities: 16

## Elapsed time: 0 seconds

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 7646

## Number of edges: 245902

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8886

## Number of communities: 16

## Elapsed time: 0 seconds

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 5157

## Number of edges: 164980

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8852

## Number of communities: 17

## Elapsed time: 0 seconds

merged_seurat <– Reduce(function(x, y) merge(x, y), results)
rm(results)
# 清理 pANN 列，仅保留 doublet_info
del_columns <– grep("pANN_", colnames(merged_seurat@meta.data))
if (length(del_columns)) merged_seurat@meta.data <– merged_seurat@meta.data[, –del_columns, drop = FALSE]
# 仅保留 Singlet
merged_seurat <– subset(merged_seurat, subset = doublet_info=="Singlet")

6 Harmony 去批次 + 降维聚类（R）

merged_seurat <- NormalizeData(merged_seurat)
merged_seurat <- FindVariableFeatures(merged_seurat, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(merged_seurat)
merged_seurat <- ScaleData(merged_seurat, features = all.genes)
merged_seurat <- RunPCA(merged_seurat, features = VariableFeatures(object = merged_seurat))
merged_seurat <- RunHarmony(merged_seurat, group.by = "orig.ident", verbose = FALSE)
merged_seurat <- RunUMAP(merged_seurat, reduction = "harmony", dims = 1:30, verbose = FALSE)
merged_seurat <- FindNeighbors(merged_seurat, reduction = "harmony", dims = 1:30, verbose = FALSE)
merged_seurat <- FindClusters(merged_seurat, resolution = 0.8)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 39687

## Number of edges: 1408381

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.9038

## Number of communities: 22

## Elapsed time: 7 seconds

7 细胞类型标记、cluster的映射即细胞类型注释（R）

ENDO_marker <– c("ENG","PECAM1","VWF")
FIB_marker <– c("COL1A1","DCN")
SMC_marker <– c("ACTA2","TAGLN")
IMM_marker <– c("PTPRC","LYZ")
LYME_marker <– c("CCL21")
SGC_marker <– c("SCGB1B2P")
NEU_marker <– c("RBFOX3","NRXN1")
KRT_marker <– c("KRT14","KRT1","KRT5")
MELA_marker <– c("TYRP1","PMEL")
genes_to_check <– unique(c(ENDO_marker,FIB_marker,SMC_marker,IMM_marker,LYME_marker,SGC_marker,NEU_marker,KRT_marker,MELA_marker))
DotPlot(merged_seurat,group.by = 'seurat_clusters', features = genes_to_check, cluster.idents = TRUE) + coord_flip()

# cluster → celltype 映射
cluster_mapping <– list(
ENDO = c(0,3,6), FIB = c(1,4,7,13,15,17), SMC = c(2,9),
IMM = c(11,16), LYME = c(10), SGC = c(20), NEU = c(18),
KRT = c(5,8,14,19), MELA = c(21), Unknow = c(12)
)
celltype_info <– do.call(rbind, lapply(names(cluster_mapping), function(ctype) {
clusters <– cluster_mapping[[ctype]]
if(length(clusters) > 1) {
data.frame(Cluster = clusters, Celltype = ctype, SubCelltype = paste0(ctype, "-", seq_along(clusters)))
} else {
data.frame(Cluster = clusters, Celltype = ctype, SubCelltype = ctype)
}
}))
merged_seurat$Celltype <– celltype_info$Celltype[match(merged_seurat$seurat_clusters, celltype_info$Cluster)]
merged_seurat$SubCelltype <– celltype_info$SubCelltype[match(merged_seurat$seurat_clusters, celltype_info$Cluster)]
# UMAP + 中心编号
centers <– as.data.frame(Embeddings(merged_seurat, "umap")) %>%
cbind(cluster = merged_seurat$seurat_clusters) %>%
group_by(cluster) %>%
summarize(UMAP_1 = median(UMAP_1), UMAP_2 = median(UMAP_2), .groups = "drop")
legend_labels <– sapply(0:21, function(i) {
idx <– which(celltype_info$Cluster == i)
if (length(idx) > 0) {
ctype <– celltype_info$Celltype[idx]; subctype <– celltype_info$SubCelltype[idx]
if (ctype == subctype) paste0("C", i, "-", ctype) else paste0("C", i, "-", ctype, "(", subctype,")")
} else paste0("C", i)
})
merged_seurat$cluster_label <– factor(merged_seurat$seurat_clusters, labels = legend_labels)
merged_seurat$cluster_label <– factor(merged_seurat$cluster_label, levels = unique(merged_seurat$cluster_label[order(merged_seurat$SubCelltype)]))
DimPlot(merged_seurat, reduction = "umap", group.by = "cluster_label") +
guides(color = guide_legend(ncol = 1)) +
geom_text(data = centers, aes(x = UMAP_1, y = UMAP_2, label = cluster), size = 3)

8 各样本细胞类型组成比例、top10基因热图与标记基因的FeaturePlot（R）

# 各样本 Celltype 组成（堆叠柱）
prop_data <– merged_seurat@meta.data %>%
group_by(orig.ident, Celltype) %>% summarise(count = n(), .groups = "drop_last") %>%
group_by(orig.ident) %>% mutate(prop = count / sum(count))
celltype_colors <– c(
"ENDO" = "#1F77B4", "FIB" = "#FF7F0E", "SMC" = "#2CA02C",
"IMM" = "#D62728", "LYME" = "#9467BD", "SGC" = "#8C564B",
"NEU" = "#E377C2", "KRT" = "#7F7F7F", "MELA" = "#BCBD22",
"Unknow" = "#17BECF"
)
ggplot(prop_data, aes(x = orig.ident, y = prop, fill = Celltype)) +
geom_col(width = 0.8) +
scale_fill_manual(values = celltype_colors) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "orig.ident", y = "Proportion") +
theme_minimal() +
theme(panel.grid = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1))

# 每个 Celltype 的 DE 基因热图（Top10）
seurat_subset <– subset(merged_seurat, subset = Celltype != "Unknow")
DefaultAssay(seurat_subset) <– "RNA"
Idents(seurat_subset) <– "Celltype"
# —— 2) 计算 marker 并准备 features（去重 + 仅保留存在的基因）——
markers <– FindAllMarkers(seurat_subset, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top10 <– dplyr::group_by(markers, cluster) |> dplyr::slice_max(avg_log2FC, n = 10)
feats <– unique(top10$gene)
feats <– feats[feats %in% rownames(seurat_subset)]
stopifnot(length(feats) > 0)
suppressWarnings(seurat_subset <– ScaleData(seurat_subset, features = feats, verbose = FALSE))
missing_ct <– setdiff(levels(seurat_subset$Celltype), names(celltype_colors))
if (length(missing_ct)) {
add_cols <– setNames(rep("
#999999
", length(missing_ct)), missing_ct)
celltype_colors <– c(celltype_colors, add_cols)
}
top3_labels <– dplyr::group_by(top10, cluster) |> dplyr::slice_head(n = 3)
y_labels <– ifelse(feats %in% top3_labels$gene, feats, "")
label_colors <– ifelse(y_labels == "", "white", "black")
p <– DoHeatmap(
seurat_subset,
features = feats,
group.colors = celltype_colors,
slot = "scale.data",
disp.min = –2, disp.max = 2, label = FALSE
) +
scale_fill_gradient2(low = "white", mid = "white", high = "red") +
scale_y_discrete(labels = y_labels) +
theme(axis.text.y = element_text(color = label_colors, size = 8))
print(p)

dir.create("output", showWarnings = FALSE)
file.copy("/gpfs/hpc/home/chenchao/ranm/GSE163973/Heatmap_top10_Celltype.pdf", "output/Heatmap_top10_Celltype.pdf", overwrite = TRUE)

## [1] TRUE

knitr::include_graphics("output/Heatmap_top10_Celltype.pdf")

#rm
(seurat_subset)
# FeaturePlot（带轴线与标题配色）
features <– c("KRT14", "COL1A1", "TYRP1", "PTPRC", "SCGB1B2P",
"KRT1", "NRXN1", "CCL21", "ENG", "ACTA2")
celltype_colors <– c(
"KRT" = "#7F7F7F", "FIB" = "#FF7F0E", "MELA" = "#BCBD22",
"IMM" = "#D62728", "SGC" = "#8C564B", "NEU" = "#E377C2",
"LYME" = "#9467BD", "ENDO" = "#1F77B4", "SMC" = "#2CA02C")
titles <– c(
"KRT:KRT14", "FIB:COL1A1", "MELA:TYRP1", "IMM:PTPRC", "SGC:SCGB1B2P",
"KRT:KRT1", "NEU:NRXN1", "LYME:CCL21", "ENDO:ENG", "SMC:ACTA2")
plots <– FeaturePlot(merged_seurat, features = features,
cols = c("lightgrey", "red"), combine = FALSE)
for(i in seq_along(plots)) {
if(i %in% c(1, 6)) {
plots[[i]] <– plots[[i]] +
labs(x = NULL, y = "UMAP_2") +
theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_line())
} else if(i == 8) {
plots[[i]] <– plots[[i]] +
labs(x = "UMAP_1", y = NULL) +
theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_line())
} else {
plots[[i]] <– plots[[i]] +
labs(x = NULL, y = NULL) +
theme(axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_line())
}
celltype <– strsplit(titles[i], ":")[[1]][1]
gene_name <– strsplit(titles[i], ":")[[1]][2]
plots[[i]] <– plots[[i]] +
ggtitle(paste0("<span style='color:", celltype_colors[celltype], "'>", celltype, "</span>:<span style='color:black'>", gene_name, "</span>")) +
theme(plot.title = element_markdown(face = "bold"))
}
wrap_plots(plots, nrow = 2) + plot_layout(guides = "collect") & theme(legend.position = "none")

9 提取FIB子集、UMAP图、FIB亚聚类和差异表达基因的分半小提琴（R）

fib_seurat <- subset(merged_seurat, subset = Celltype == "FIB") |>
NormalizeData() |> FindVariableFeatures() |> ScaleData() |> RunPCA()
fib_seurat <- RunHarmony(fib_seurat, group.by = "orig.ident")
fib_seurat <- FindNeighbors(fib_seurat, dims = 1:20) |> FindClusters(resolution = 0.5) |> RunUMAP(dims = 1:20)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

##

## Number of nodes: 12423

## Number of edges: 417458

##

## Running Louvain algorithm…

## Maximum modularity in 10 random starts: 0.8796

## Number of communities: 12

## Elapsed time: 1 seconds

fib_seurat$sub_cluster <– paste0("sC", as.numeric(as.character(fib_seurat$seurat_clusters)) + 1)
fib_seurat$sub_cluster <– factor(fib_seurat$sub_cluster, levels = paste0("sC", 1:length(unique(fib_seurat$sub_cluster))))
DimPlot(fib_seurat, group.by = "sub_cluster", label = FALSE) + ggtitle("") +
guides(color = guide_legend(override.aes = list(size = 3)))

# 各样本子群比例
prop_data <- fib_seurat@meta.data %>%
group_by(orig.ident, sub_cluster) %>% summarise(count = n(), .groups = "drop_last") %>%
group_by(orig.ident) %>% mutate(prop = count / sum(count))
ggplot(prop_data, aes(x = orig.ident, y = prop, fill = sub_cluster)) +
geom_col(width = 0.8) + scale_y_continuous(limits = c(0, 1)) +
labs(x = "orig.ident", y = "Proportion") +
theme_minimal() + theme(panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_line(),
legend.position = "none")

# 分半小提琴：KL 左半 / NS 右半
features <– c("ASPN","POSTN","SFRP4","MDK","COL11A1","FN1","ADAM12","BGN")
fib_seurat$group <– ifelse(grepl("KL", fib_seurat$orig.ident), "KL", "NS")
d <– FetchData(fib_seurat, c(features, "group"))
df <– data.frame(expr = as.vector(as.matrix(d[features])),
gene = factor(rep(features, each = nrow(d)), levels = features),
group = factor(rep(d$group, times = length(features)), levels = c("KL","NS")))
library(ggh4x)
ggplot(df, aes(x = 1, y = expr, fill = group)) +
geom_half_violin(data = subset(df, group == "KL"), side = "l", trim = TRUE, linewidth = 0.3, colour = "black") +
geom_half_violin(data = subset(df, group == "NS"), side = "r", trim = TRUE, linewidth = 0.3, colour = "black") +
facet_wrap2(~ gene, ncol = 4, scales = "free_y", axes = "x") +
scale_fill_manual(values = c(KL = "
#1f77b4
", NS = "
#69b3a2
")) +
scale_x_continuous(breaks = c(0.9, 1.1), labels = c("KL","NS"), limits = c(0.6, 1.4)) +
labs(x = NULL, y = "Expression Level") +
theme_classic(base_size = 11) +
theme(legend.position = "none",
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.title.y = element_text(vjust = 0.5, angle = 90, margin = margin(r = 6)),
axis.title.x = element_blank())

10 导出KL/NS作为CellPhoneDB的输入文件（R）

merged_seurat$new_celltype <– ifelse(merged_seurat$Celltype == "FIB", merged_seurat$SubCelltype, merged_seurat$Celltype)
merged_seurat$group <– ifelse(grepl("KL", merged_seurat$orig.ident), "KL", "NS")
seurat_subset <– subset(merged_seurat, subset = new_celltype != "Unknow")
kl_seurat <– subset(seurat_subset, group == "KL")
ns_seurat <– subset(seurat_subset, group == "NS")
#SaveH5Seurat
(kl_seurat, filename = "/gpfs/hpc/home/GSE163973/kl_seurat.h5seurat")
#Convert
("/gpfs/hpc/home/GSE163973/kl_seurat.h5seurat", dest = "h5ad")
#SaveH5Seurat
(ns_seurat, filename = "/gpfs/hpc/home/GSE163973/NS_seurat.h5seurat")
#Convert
("/gpfs/hpc/home/GSE163973/NS_seurat.h5seurat", dest = "h5ad")

11 运行 CellPhoneDB（Python 代码块）

注意：下面是 Python 代码，两类样本分别作为Cellphonedb的输入文件，将输出文件保存到本地用于后续分析和画图。

import pandas as pd
import os
from cellphonedb.utils import db_releases_utils
from cellphonedb.src.core.methods import cpdb_statistical_analysis_method
cpdb_KL = cpdb_statistical_analysis_method.call(
cpdb_file_path = '/gpfs/hpc/home/GSE163973/cellphonedb.zip',
meta_file_path ='/gpfs/hpc/home/GSE163973/KL_meta.txt',
counts_file_path = '/gpfs/hpc/home/GSE163973/kl_seurat.h5ad',
counts_data = 'hgnc_symbol',
score_interactions = False,
iterations = 1000,
threshold = 0.25,
threads = 8,
debug_seed = 42,
result_precision = 3,
pvalue = 0.05,
separator = '|',
debug = False,
output_path = '/gpfs/hpc/home/chenchao/ranm/GSE163973/KL'
)
cpdb_NS = cpdb_statistical_analysis_method.call(
cpdb_file_path = '/gpfs/hpc/home/chenchao/ranm/GSE163973/cellphonedb.zip',
meta_file_path ='/gpfs/hpc/home/chenchao/ranm/GSE163973/NS_meta.txt',
counts_file_path = '/gpfs/hpc/home/chenchao/ranm/GSE163973/NS_seurat.h5ad',
counts_data = 'hgnc_symbol',
score_interactions = False,
iterations = 1000,
threshold = 0.25,
threads = 8,
debug_seed = 42,
result_precision = 3,
pvalue = 0.05,
separator = '|',
debug = False,
output_path = '/gpfs/hpc/home/chenchao/ranm/GSE163973/NS'
)
print("CellPhoneDB 统计完成。")

12 统计不同细胞类型间的显著互作对数 + 左下三角热图（R）

# 指定 CPDB 输出目录（若在 Windows 本地跑过，可改为 D:/wsw/Bioma_analysis/KL 和 NS）
dir_NS <– "/gpfs/hpc/home/GSE163973/NS"
dir_KL <– "/gpfs/hpc/home/GSE163973/KL"
get_pair_counts <– function(path, alpha = 0.05, expr_thr = 0) {
fs <– list.files(path, full.names = TRUE)
latest <– function(pat, exclude = NULL) {
idx <– grepl(pat, fs)
if (!is.null(exclude)) idx <– idx & !grepl(exclude, fs)
f <– fs[idx]
if (length(f)) f[which.max(file.mtime(f))] else character(0)
}
f_sig <– latest("statistical_analysis_significant_means")
if (length(f_sig)) {
X <– read.delim(f_sig, check.names = FALSE)
cols <– grep("\\\\|", names(X), value = TRUE)
cnt <– colSums(X[cols] > expr_thr, na.rm = TRUE)
} else {
f_means <– latest("statistical_analysis_means", exclude = "significant")
f_pvals <– latest("statistical_analysis_pvalues")
if (!length(f_means) || !length(f_pvals))
stop("未在目录中找到 CPDB 的 significant_means 或 means/pvalues 文件：", path)
M <– read.delim(f_means, check.names = FALSE)
P <– read.delim(f_pvals, check.names = FALSE)
cols <– intersect(grep("\\\\|", names(M), value = TRUE), grep("\\\\|", names(P), value = TRUE))
cnt <– colSums((P[cols] < alpha) & (M[cols] > expr_thr), na.rm = TRUE)
}
s <– sub("\\\\|.*", "", names(cnt)); r <– sub(".*\\\\|", "", names(cnt))
cells <– sort(unique(c(s, r)))
mat <– matrix(0L, length(cells), length(cells), dimnames = list(cells, cells))
mat[cbind(match(s, cells), match(r, cells))] <– as.integer(cnt)
as.data.frame.matrix(mat)
}
pair_count_NS <– get_pair_counts(dir_NS)
pair_count_KL <– get_pair_counts(dir_KL)
make_undirected_vectorized <– function(adj_matrix) {
mat <– as.matrix(adj_matrix)
undirected_mat <– mat + t(mat)
diag(undirected_mat) <– diag(mat)
return(as.data.frame(undirected_mat, check.names = FALSE))
}
pair_count_NS_undirect <– make_undirected_vectorized(pair_count_NS)
pair_count_KL_undirect <– make_undirected_vectorized(pair_count_KL)
plot_lower_tri <– function(mat, title = "normal scar") {
stopifnot(is.matrix(mat) || is.data.frame(mat))
mat <– as.data.frame(mat, check.names = FALSE)
stopifnot(!is.null(rownames(mat)), !is.null(colnames(mat)))
stopifnot(all(rownames(mat) %in% colnames(mat)), all(colnames(mat) %in% rownames(mat)))
mat <– mat[colnames(mat), colnames(mat), drop = FALSE] # 按列名顺序对齐
df <– reshape2::melt(cbind(row = rownames(mat), mat), id.vars = "row",
variable.name = "col", value.name = "n")
cells <– colnames(mat); n <– length(cells)
df$col <– factor(df$col, levels = cells) # X 轴：ENDO → SMC
df$row <– factor(df$row, levels = rev(cells)) # Y 轴：SMC → ENDO（下到上）
i <– as.integer(df$col); j <– as.integer(df$row)
df <– df[(i + j) <= (n + 1), ] # 仅保留下三角（含对角）
ggplot2::ggplot(df, ggplot2::aes(col, row, fill = n)) +
ggplot2::geom_tile(color = "grey60", linewidth = 0.4) +
ggplot2::geom_text(ggplot2::aes(label = n), size = 3) +
ggplot2::scale_fill_gradient(low = "
#2b6eb3
", high = "
#d7301f
") +
ggplot2::scale_x_discrete(drop = FALSE) +
ggplot2::scale_y_discrete(drop = FALSE) +
ggplot2::coord_fixed() +
ggplot2::labs(title = title, x = NULL, y = NULL) +
ggplot2::theme_minimal(base_size = 12) +
ggplot2::theme(legend.position = "none",
panel.grid = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5, face = "bold"),
axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, vjust = 1))
}
p1 <– plot_lower_tri(pair_count_NS_undirect, title = "normal scar")
p2 <– plot_lower_tri(pair_count_KL_undirect, title = "keloid")
p1+p2

欢迎致谢

如果以上内容对你有帮助，欢迎在文章的Acknowledgement中加上这一段，联系客服微信可以发放奖励：

Since Biomamba and his wechat public account team produce bioinformatics tutorials and share code with annotation, we thank Biomamba
for
their guidance
in
bioinformatics and data analysis
for
the current study.

欢迎在发文/毕业时向我们分享你的喜悦~