20260912 093449 Seurat V5 数据可视化方法与 Sctransform 归一化教程

library(Seurat)
library(ggplot2)
library(patchwork)

# 加载阶段1保存的最终对象
pbmc <- readRDS("output/pbmc3k_final.rds")

# 添加随机分组变量（用于演示 split.by 功能）
set.seed(42)
pbmc$groups <- sample(c("group1", "group2"), size = ncol(pbmc), replace = TRUE)

# 演示用基因列表
features <- c("LYZ", "CCL5", "IL32", "PTPRCAP", "FCGR3A", "PF4")

p_ridge <- RidgePlot(pbmc, features = features, ncol = 2)

p_vln <- VlnPlot(pbmc, features = features)

p_feature <- FeaturePlot(pbmc, features = features)

p_dot <- DotPlot(pbmc, features = features) + RotatedAxis()

p_heat <- DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)

FeaturePlot(pbmc, features = "MS4A1")

FeaturePlot(pbmc, features = "MS4A1", min.cutoff = 1, max.cutoff = 3)

FeaturePlot(pbmc, features = c("MS4A1", "PTPRCAP"), min.cutoff = "q10", max.cutoff = "q90")

FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), blend = TRUE)

FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), split.by = "groups")

VlnPlot(pbmc, features = "percent.mt", split.by = "groups")

DotPlot(pbmc, features = features, split.by = "groups") + RotatedAxis()

DimPlot(pbmc)

DoHeatmap(pbmc, features = VariableFeatures(pbmc)[1:100], cells = 1:500, size = 4, angle = 90) + NoLegend()

baseplot <- DimPlot(pbmc, reduction = "umap")
baseplot + labs(title = "Clustering of 2,700 PBMCs")

baseplot + DarkTheme()

baseplot + DarkTheme() +
    theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18)) +
    guides(colour = guide_legend(override.aes = list(size = 10)))

plot1 <- DimPlot(pbmc)
plot2 <- FeatureScatter(pbmc, feature1 = "LYZ", feature2 = "CCL5")
plot1 + plot2

(plot1 + plot2) & NoLegend()

p_label <- DimPlot(pbmc, reduction = "pca") + NoLegend()
p_label <- LabelClusters(plot = p_label, id = "ident")

library(Seurat)
library(ggplot2)
library(sctransform)

# 从原始数据开始（独立于阶段1的标准流程）
pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data)

# 计算线粒体百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# QC 过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

# SCTransform 一步完成归一化 + 变量特征选择 + 标准化
# 同时回归线粒体百分比（去除技术变异）
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

> DefaultAssay(pbmc)
[1] "SCT"

> pbmc[["SCT"]]
SCTAssay data with 12519 features for 2638 cells, and 1 SCTModel(s)
Top 10 variable features:
 S100A9, GNLY, LYZ, S100A8, NKG7, FTL, GZMB, IGLL5, FTH1, CCL5

pbmc <- RunPCA(pbmc, verbose = FALSE)
DimPlot(pbmc, reduction = "pca")

pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)
DimPlot(pbmc, reduction = "umap")

ElbowPlot(pbmc, ndims = 50)

pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindClusters(pbmc, verbose = FALSE)
DimPlot(pbmc, reduction = "umap", label = TRUE)

VlnPlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7", "ISG15", "CD3D"),
         pt.size = 0.2, ncol = 4)

FeaturePlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7"),
            pt.size = 0.2, ncol = 3)

FeaturePlot(pbmc, features = c("CD3D", "ISG15", "TCL1A", "FCER2", "XCL1", "FCGR3A"),
            pt.size = 0.2, ncol = 3)

length(VariableFeatures(pbmc))  # 3000

DimPlot(pbmc, reduction = "umap", group.by = "seurat_clusters") +
    ggtitle("SCT: 30 PCs (推荐)")

library(dplyr)
pbmc <- CreateSeuratObject(pbmc_data) %>%
    PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
    SCTransform(vars.to.regress = "percent.mt") %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    RunUMAP(dims = 1:30) %>%
    FindClusters()

###############################################################################
# Seurat v5 阶段2-A：数据可视化方法
# 基于官方教程: https://satijalab.org/seurat/articles/visualization_vignette
# 使用阶段1保存的 pbmc3k_final.rds 对象
###############################################################################

library(Seurat)
library(ggplot2)
library(patchwork)

# === 路径配置 ===
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
viz_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output/visualization"
dir.create(viz_dir, showWarnings = FALSE, recursive = TRUE)

# === 加载阶段1的最终对象 ===
pbmc <- readRDS(file.path(output_dir, "pbmc3k_final.rds"))
cat("已加载 pbmc3k_final.rds\n")
print(pbmc)

# 添加分组变量（用于演示 split.by 功能）
set.seed(42)
pbmc$groups <- sample(c("group1", "group2"), size = ncol(pbmc), replace = TRUE)

# 定义演示用基因列表
features <- c("LYZ", "CCL5", "IL32", "PTPRCAP", "FCGR3A", "PF4")


# ============================================================================
# 1. 五种可视化方法展示 Marker 基因表达
# ============================================================================

cat("\n========== 1. 五种可视化方法 ==========\n")

# 1.1 Ridge Plot（山脊图）- 展示各聚类的基因表达分布
p_ridge <- RidgePlot(pbmc, features = features, ncol = 2)
ggsave(file.path(viz_dir, "01_ridge_plot.png"), plot = p_ridge, width = 12, height = 10, dpi = 150)
cat("已保存: 01_ridge_plot.png\n")

# 1.2 Violin Plot（小提琴图）- 展示各聚类的基因表达分布
p_vln <- VlnPlot(pbmc, features = features)
ggsave(file.path(viz_dir, "02_violin_plot.png"), plot = p_vln, width = 14, height = 8, dpi = 150)
cat("已保存: 02_violin_plot.png\n")

# 1.3 Feature Plot - 在低维空间展示基因表达
p_feature <- FeaturePlot(pbmc, features = features)
ggsave(file.path(viz_dir, "03_feature_plot.png"), plot = p_feature, width = 14, height = 10, dpi = 150)
cat("已保存: 03_feature_plot.png\n")

# 1.4 Dot Plot（气泡图）- 点大小=表达比例，颜色=平均表达水平
p_dot <- DotPlot(pbmc, features = features) + RotatedAxis()
ggsave(file.path(viz_dir, "04_dot_plot.png"), plot = p_dot, width = 10, height = 6, dpi = 150)
cat("已保存: 04_dot_plot.png\n")

# 1.5 Heatmap（热图）- 降采样后展示基因表达
p_heat <- DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)
ggsave(file.path(viz_dir, "05_heatmap.png"), plot = p_heat, width = 10, height = 6, dpi = 150)
cat("已保存: 05_heatmap.png\n")


# ============================================================================
# 2. 自定义 FeaturePlot
# ============================================================================

cat("\n========== 2. 自定义 FeaturePlot ==========\n")

# 2.1 带图例的 FeaturePlot
p_fp1 <- FeaturePlot(pbmc, features = "MS4A1")
ggsave(file.path(viz_dir, "06_featureplot_legend.png"), plot = p_fp1, width = 7, height = 6, dpi = 150)
cat("已保存: 06_featureplot_legend.png\n")

# 2.2 调整对比度（设定 min/max cutoff）
p_fp2 <- FeaturePlot(pbmc, features = "MS4A1", min.cutoff = 1, max.cutoff = 3)
ggsave(file.path(viz_dir, "07_featureplot_contrast.png"), plot = p_fp2, width = 7, height = 6, dpi = 150)
cat("已保存: 07_featureplot_contrast.png\n")

# 2.3 分位数 cutoff（适合多基因同时绘图）
p_fp3 <- FeaturePlot(pbmc, features = c("MS4A1", "PTPRCAP"), min.cutoff = "q10", max.cutoff = "q90")
ggsave(file.path(viz_dir, "08_featureplot_quantile.png"), plot = p_fp3, width = 12, height = 5, dpi = 150)
cat("已保存: 08_featureplot_quantile.png\n")

# 2.4 混合展示两个基因的共表达（blend 模式）
p_fp4 <- FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), blend = TRUE)
ggsave(file.path(viz_dir, "09_featureplot_blend.png"), plot = p_fp4, width = 14, height = 5, dpi = 150)
cat("已保存: 09_featureplot_blend.png\n")

# 2.5 按分组分割展示
p_fp5 <- FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), split.by = "groups")
ggsave(file.path(viz_dir, "10_featureplot_split.png"), plot = p_fp5, width = 14, height = 10, dpi = 150)
cat("已保存: 10_featureplot_split.png\n")


# ============================================================================
# 3. 其他可视化功能
# ============================================================================

cat("\n========== 3. 其他可视化功能 ==========\n")

# 3.1 按分组分割的小提琴图
p_vln_split <- VlnPlot(pbmc, features = "percent.mt", split.by = "groups")
ggsave(file.path(viz_dir, "11_violin_split.png"), plot = p_vln_split, width = 8, height = 5, dpi = 150)
cat("已保存: 11_violin_split.png\n")

# 3.2 按分组分割的气泡图
p_dot_split <- DotPlot(pbmc, features = features, split.by = "groups") + RotatedAxis()
ggsave(file.path(viz_dir, "12_dotplot_split.png"), plot = p_dot_split, width = 12, height = 6, dpi = 150)
cat("已保存: 12_dotplot_split.png\n")

# 3.3 DimPlot 默认使用 UMAP
p_dim1 <- DimPlot(pbmc)
ggsave(file.path(viz_dir, "13_dimplot_umap.png"), plot = p_dim1, width = 8, height = 7, dpi = 150)
cat("已保存: 13_dimplot_umap.png\n")

# 3.4 DoHeatmap 带分组条
p_heat2 <- DoHeatmap(pbmc, features = VariableFeatures(pbmc)[1:100], cells = 1:500, size = 4, angle = 90) + NoLegend()
ggsave(file.path(viz_dir, "14_heatmap_grouped.png"), plot = p_heat2, width = 10, height = 8, dpi = 150)
cat("已保存: 14_heatmap_grouped.png\n")


# ============================================================================
# 4. 应用主题（Theme）
# ============================================================================

cat("\n========== 4. 应用主题 ==========\n")

baseplot <- DimPlot(pbmc, reduction = "umap")

# 4.1 添加标题
p_theme1 <- baseplot + labs(title = "Clustering of 2,700 PBMCs")
ggsave(file.path(viz_dir, "15_theme_title.png"), plot = p_theme1, width = 8, height = 7, dpi = 150)
cat("已保存: 15_theme_title.png\n")

# 4.2 暗色主题
p_theme2 <- baseplot + DarkTheme()
ggsave(file.path(viz_dir, "16_theme_dark.png"), plot = p_theme2, width = 8, height = 7, dpi = 150)
cat("已保存: 16_theme_dark.png\n")

# 4.3 链式主题
p_theme3 <- baseplot + DarkTheme() +
    theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18)) +
    guides(colour = guide_legend(override.aes = list(size = 10)))
ggsave(file.path(viz_dir, "17_theme_chained.png"), plot = p_theme3, width = 10, height = 8, dpi = 150)
cat("已保存: 17_theme_chained.png\n")


# ============================================================================
# 5. 图组合并
# ============================================================================

cat("\n========== 5. 图组合并 ==========\n")

plot1 <- DimPlot(pbmc)
plot2 <- FeatureScatter(pbmc, feature1 = "LYZ", feature2 = "CCL5")

# 5.1 拼接两图
p_combine1 <- plot1 + plot2
ggsave(file.path(viz_dir, "18_combine_plots.png"), plot = p_combine1, width = 14, height = 6, dpi = 150)
cat("已保存: 18_combine_plots.png\n")

# 5.2 移除所有图例
p_combine2 <- (plot1 + plot2) & NoLegend()
ggsave(file.path(viz_dir, "19_combine_nolegend.png"), plot = p_combine2, width = 14, height = 6, dpi = 150)
cat("已保存: 19_combine_nolegend.png\n")


# ============================================================================
# 6. 聚类标签
# ============================================================================

cat("\n========== 6. 聚类标签 ==========\n")

# 6.1 标注聚类
p_label1 <- DimPlot(pbmc, reduction = "pca") + NoLegend()
p_label1 <- LabelClusters(plot = p_label1, id = "ident")
ggsave(file.path(viz_dir, "20_label_clusters.png"), plot = p_label1, width = 8, height = 7, dpi = 150)
cat("已保存: 20_label_clusters.png\n")


# ============================================================================
# 总结
# ============================================================================

cat("\n")
cat("=======================================================\n")
cat("  阶段2-A 可视化教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", viz_dir, "\n")
cat("  图片数量: 20 张\n")
cat("  覆盖内容:\n")
cat("    01-05: 五种核心可视化方法\n")
cat("    06-10: FeaturePlot 自定义技巧\n")
cat("    11-14: 分组分割与热图\n")
cat("    15-17: 主题定制\n")
cat("    18-19: 图组合并\n")
cat("    20:    聚类标签\n")
cat("=======================================================\n")

###############################################################################
# Seurat v5 阶段2-B：SCTransform 归一化分析
# 基于官方教程: https://satijalab.org/seurat/articles/sctransform_vignette
# SCTransform 替代 NormalizeData + ScaleData + FindVariableFeatures 三步
# 使用 PBMC 3k 数据从头创建新对象（独立于阶段1的标准流程）
###############################################################################

library(Seurat)
library(ggplot2)
library(sctransform)
library(patchwork)

# === 路径配置 ===
data_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/data/filtered_gene_bc_matrices/hg19"
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
sct_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output/sctransform"
dir.create(sct_dir, showWarnings = FALSE, recursive = TRUE)

# ============================================================================
# 1. 加载数据并创建 Seurat 对象
# ============================================================================
cat("\n========== 1. 加载数据 ==========\n")

pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data)
cat("已创建 Seurat 对象\n")
print(pbmc)

# ============================================================================
# 2. 计算线粒体百分比并应用 SCTransform
# ============================================================================
cat("\n========== 2. SCTransform 归一化 ==========\n")

# 计算线粒体基因百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
cat("已计算线粒体百分比\n")

# QC 过滤（与阶段1一致）
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("QC 过滤后细胞数:\n")
cat(ncol(pbmc), "\n")

# 应用 SCTransform（同时回归线粒体百分比）
# SCTransform 一步完成：归一化 + 变量特征选择 + 标准化
cat("正在运行 SCTransform（可能需要1-2分钟）...\n")
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = TRUE)
cat("SCTransform 完成!\n")
cat("当前默认 assay:\n")
print(DefaultAssay(pbmc))

# ============================================================================
# 3. 降维分析：PCA + UMAP
# ============================================================================
cat("\n========== 3. 降维分析 ==========\n")

# PCA（使用 SCT residuals）
pbmc <- RunPCA(pbmc, verbose = FALSE)
cat("PCA 完成\n")

# UMAP（使用前30个PC，SCTransform 推荐使用更多PC）
pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)
cat("UMAP 完成\n")

# 3.1 PCA 可视化
p_pca <- DimPlot(pbmc, reduction = "pca")
ggsave(file.path(sct_dir, "01_pca.png"), plot = p_pca, width = 8, height = 7, dpi = 150)
cat("已保存: 01_pca.png\n")

# 3.2 UMAP 可视化
p_umap <- DimPlot(pbmc, reduction = "umap")
ggsave(file.path(sct_dir, "02_umap.png"), plot = p_umap, width = 8, height = 7, dpi = 150)
cat("已保存: 02_umap.png\n")

# 3.3 PCA Elbow Plot（选择PC数量参考）
p_elbow <- ElbowPlot(pbmc, ndims = 50)
ggsave(file.path(sct_dir, "03_elbow_plot.png"), plot = p_elbow, width = 8, height = 5, dpi = 150)
cat("已保存: 03_elbow_plot.png\n")

# ============================================================================
# 4. 聚类分析
# ============================================================================
cat("\n========== 4. 聚类分析 ==========\n")

pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindClusters(pbmc, verbose = FALSE)
cat("聚类完成\n")

# 4.1 带标签的 UMAP
p_umap_label <- DimPlot(pbmc, reduction = "umap", label = TRUE)
ggsave(file.path(sct_dir, "04_umap_labeled.png"), plot = p_umap_label, width = 8, height = 7, dpi = 150)
cat("已保存: 04_umap_labeled.png\n")

# ============================================================================
# 5. SCTransform 优势展示：更精细的细胞亚群
# ============================================================================
cat("\n========== 5. SCTransform 优势展示 ==========\n")

# 5.1 T 细胞亚群标志基因小提琴图
# SCTransform 能更清晰地区分 CD8 T 细胞亚群（naive/memory/effector）
p_vln_t <- VlnPlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7", "ISG15", "CD3D"),
                   pt.size = 0.2, ncol = 4)
ggsave(file.path(sct_dir, "05_violin_tcell_markers.png"), plot = p_vln_t, width = 16, height = 10, dpi = 150)
cat("已保存: 05_violin_tcell_markers.png\n")

# 5.2 T 细胞标志基因 FeaturePlot
p_fp_t <- FeaturePlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7"),
                       pt.size = 0.2, ncol = 3)
ggsave(file.path(sct_dir, "06_featureplot_tcell.png"), plot = p_fp_t, width = 14, height = 10, dpi = 150)
cat("已保存: 06_featureplot_tcell.png\n")

# 5.3 B 细胞与 NK 细胞标志基因 FeaturePlot
p_fp_bnk <- FeaturePlot(pbmc, features = c("CD3D", "ISG15", "TCL1A", "FCER2", "XCL1", "FCGR3A"),
                         pt.size = 0.2, ncol = 3)
ggsave(file.path(sct_dir, "07_featureplot_bnk.png"), plot = p_fp_bnk, width = 14, height = 10, dpi = 150)
cat("已保存: 07_featureplot_bnk.png\n")

# ============================================================================
# 6. SCTransform 与标准归一化对比
# ============================================================================
cat("\n========== 6. 对比分析 ==========\n")

# 6.1 SCTransform 变量特征数量（默认3000个，比标准流程的2000更多）
sct_features <- VariableFeatures(pbmc)
cat("SCTransform 变量特征数量:", length(sct_features), "\n")

# 6.2 展示前20个变量特征
p_var <- LabelPoints(plot = VariableFeaturePlot(pbmc), points = sct_features[1:20], repel = TRUE)
ggsave(file.path(sct_dir, "08_variable_features.png"), plot = p_var, width = 10, height = 8, dpi = 150)
cat("已保存: 08_variable_features.png\n")

# 6.3 SCT assay 结构说明
cat("\nSCT Assay 层级结构:\n")
print(pbmc[["SCT"]])

# 6.4 标准流程对比图：使用前10个PC（标准流程推荐）vs 前30个PC（SCT推荐）
p_umap_10 <- DimPlot(pbmc, reduction = "umap", group.by = "seurat_clusters") +
    ggtitle("SCT: 30 PCs (推荐)")
ggsave(file.path(sct_dir, "09_umap_sct_30pc.png"), plot = p_umap_10, width = 8, height = 7, dpi = 150)
cat("已保存: 09_umap_sct_30pc.png\n")

# ============================================================================
# 7. 保存 SCTransform 结果
# ============================================================================
cat("\n========== 7. 保存结果 ==========\n")

saveRDS(pbmc, file.path(output_dir, "pbmc3k_sct.rds"))
cat("已保存 SCTransform 对象: pbmc3k_sct.rds\n")

# ============================================================================
# 总结
# ============================================================================
cat("\n")
cat("=======================================================\n")
cat("  阶段2-B SCTransform 教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", sct_dir, "\n")
cat("  图片数量: 9 张\n")
cat("  覆盖内容:\n")
cat("    01-03: PCA + UMAP + Elbow Plot\n")
cat("    04:    聚类标注\n")
cat("    05-07: T/B/NK 细胞亚群标志基因可视化\n")
cat("    08:    SCTransform 变量特征\n")
cat("    09:    SCT UMAP（30 PCs）\n")
cat("  SCT vs 标准流程:\n")
cat("    - 变量特征: 3000（SCT） vs 2000（标准）\n")
cat("    - 推荐PC数: 30（SCT） vs 10（标准）\n")
cat("    - 归一化方法: 负二项回归 vs LogNormalize\n")
cat("    - 能揭示更精细的细胞亚群（CD8 naive/memory/effector）\n")
cat("=======================================================\n")