当前位置：首页 > news >正文

29-植物单细胞RNA-seq分析教程3-2025年版

news 2026/5/12 19:23:33

本教程基于银白杨顶芽单细胞测序数据，聚焦 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 &# 3. 查看输出结果（关键文件：filtered_feature_bc_matrix）
ls /home/workingdirectory/bud1/outs/filtered_feature_bc_matrix/
# 包含3个文件：barcodes.tsv.gz（细胞条形码）、features.tsv.gz（基因信息）、matrix.mtx.gz（表达矩阵）

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()