本教程基于银白杨顶芽单细胞测序数据，聚焦 CYC 和 CDK 基因家族参与叶片发育的分子机制，适配有少量生信基础（了解 Linux 基本操作、R 语言入门）的读者。教程共4节，遵循 “软件准备→基础分析→细胞注释→高级挖掘” 的分析流程，代码均来自实际研究，注释结合发表文章核心结果，确保 “代码可复现、结果可解读”。
第 2 节：单细胞数据预处理与基础分析
2.1 数据定量（cellranger count）
对两个生物学重复（bud1、bud2）进行单细胞定量，输出细胞 - 基因表达矩阵：
# 1. 数据整理（测序文件放在指定文件夹）
mkdir -p bud1_fastq bud2_fastq bud1 bud2
# 将bud1的测序文件（R1/R2/I1）放入bud1_fastq文件夹中，bud2同理
# 2. 定量分析（后台运行，输出日志）
nohup cellranger count --id=bud1 \
--transcriptome=Poalb.genome \
--fastqs=bud1_fastq/ \
--sample=bud1 > bud1.log 2>&1&
nohup cellranger count --id=bud2 \
--transcriptome=Poalb.genome \
--fastqs=bud2_fastq/ \
--sample=bud2 > bud2.log 2>&1&
2.2 Seurat 对象创建与 QC 筛选
# 1. 加载包与数据
library(Seurat)
library(dplyr)
library(patchwork)
setwd("/home/workingdirectory")
# 读取bud1数据（bud2流程完全一致）
bud1.data <- Read10X(data.dir ="bud1/outs/filtered_feature_bc_matrix")
# 创建Seurat对象（过滤低质量基因：至少在3个细胞中表达；过滤低质量细胞：至少含200个基因）
bud1 <- CreateSeuratObject(counts =bud1.data, project = "bud1", min.cells = 3, min.features = 200)
# 2. 计算QC指标（重点：线粒体污染比例）
# 银白杨线粒体基因以"gene-E"开头，需根据实际基因组注释调整
bud1[["percent.mt"]] <-PercentageFeatureSet(bud1, pattern = "^gene-E")
# 3. QC可视化（评估数据质量）
pdf("VlnPlot_QC.bud1.pdf", height= 5, width = 10)
VlnPlot(bud1, features =c("nFeature_RNA", "nCount_RNA", "percent.mt"),ncol = 3)
dev.off()
pdf("FeatureScatter.bud1.pdf",height = 5, width = 10)
plot1 <- FeatureScatter(bud1, feature1 ="nCount_RNA", feature2 = "percent.mt") #总表达量与线粒体污染的关系
plot2 <- FeatureScatter(bud1, feature1 ="nCount_RNA", feature2 = "nFeature_RNA") #总表达量与基因数的关系
plot1 + plot2
dev.off()
# 4. 筛选高质量细胞
# 筛选标准：基因数200-7500，线粒体比例<30%，总表达量<40000
bud1 <- subset(bud1, subset =nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 0.3& nCount_RNA < 40000)
bud1 #输出结果：27052
features across 16872 samples（仅过滤22个低质量细胞）
2.3 数据标准化与降维聚类
# 1. 标准化（LogNormalize：总表达量归一化后log转换）
bud1 <- NormalizeData(bud1,normalization.method = "LogNormalize", scale.factor = 10000)
# 2. 鉴定高变基因（用于后续PCA降维，默认2000个）
bud1 <- FindVariableFeatures(bud1,selection.method = "vst", nfeatures = 2000)
# 可视化高变基因（标记top10基因）
top10 <- head(VariableFeatures(bud1),10)
pdf("VariableFeaturePlot.bud1.pdf",height = 5, width = 10)
plot1 <- VariableFeaturePlot(bud1)
plot2 <- LabelPoints(plot = plot1,points = top10, repel = TRUE)
plot1 + plot2
dev.off()
# 3. 数据缩放（消除批次效应与高表达基因影响）
bud1 <- ScaleData(bud1)
# 4. PCA降维（仅用高变基因，选择前20个主成分）
bud1 <- RunPCA(bud1, features =VariableFeatures(object = bud1))
# 可视化PCA结果
pdf("DimPlot_PCA.bud1.pdf",height = 5, width = 5)
DimPlot(bud1, reduction = "pca")
dev.off()
# 5. 确定最佳聚类分辨率（用clustree可视化）
library(clustree)
bud1 <- FindNeighbors(bud1, dims = 1:20)
obj <- FindClusters(bud1, resolution =seq(0.5, 1.2, by = 0.1))
pdf("clustree.bud1.pdf", height =5, width = 7)
clustree(obj)
dev.off() #论文选择resolution=0.9（得到17个cluster）
2.4 降维聚类
# 6. 最终聚类与非线性降维（UMAP/tSNE）
bud1 <- FindClusters(bud1, resolution =0.9)
# UMAP可视化（细胞聚类结果）
bud1 <- RunUMAP(bud1, dims = 1:20)
pdf("DimPlot_UMAP.bud1.pdf",height = 5, width = 5)
DimPlot(bud1, reduction = "umap",label = TRUE)
dev.off()
# tSNE可视化（可选，用于验证聚类稳定性）
bud1 <- RunTSNE(bud1, dims = 1:20)
pdf("DimPlot_tSNE.bud1.pdf",height = 5, width = 5)
DimPlot(bud1, reduction = "tsne",label = TRUE)
dev.off()
#封面：憧憬成为魔法少女
