正文

蛋白质转录组分析怎么做 单细胞测序技术解读基因表达差异实验流程与数据分析指南

/0 浏览量
0910

嘿，朋友。看到“单细胞测序”这几个字，你是不是脑子里已经浮现出那些密密麻麻的热图、复杂的降维聚类，还有满屏看不懂的代码？别紧张。其实，单细胞转录组（scRNA-seq）并没有那么玄乎。它本质上就是回答一个问题：在同一个组织里，为什么有的细胞是A，有的是B，而有的又是C？

今天，我们不整那些虚头巴脑的教科书定义，我就像个大伙儿身边的老学长，带着你一步步把这套流程摸透。从湿实验拿到那一管细胞悬液，到最后跑出漂亮的差异表达基因，咱们一个一个过。

一、 为什么我们非要“单细胞”？

在讲怎么做之前，你得先明白为什么要做。

以前的 Bulk RNA-seq（Bulk转录组），就像是你把一筐水果全都打成汁，然后尝尝这杯汁。你知道这筐水果平均是甜的还是酸的，但你永远不知道里面有没有坏掉的果子，也不知道哪颗苹果特别甜，哪颗柠檬特别酸。

而 scRNA-seq，就是给每一颗水果单独拍照。

你得到的不再是一个平均值，而是成千上万个独立的数据点。每一个点，就是一个细胞。你可以看到：

  • 这个组织里到底有哪些细胞类型？（发现新细胞）
  • 同一种细胞在不同状态下，基因表达有什么细微差别？（细胞亚群）
  • 细胞之间是怎么“聊天”的？（细胞通讯）

所以，单细胞技术的核心价值在于分辨率。它让我们看到了组织内部的异质性，这是 Bulk 测序永远做不到的。

二、 湿实验流程：从组织到条形码

数据分析的第一步，其实在实验室里就已经决定了。如果你拿到的原始数据质量不好，后面算法再神也救不回来。

1. 样本制备：最关键的一步

这一步的核心目标是：把组织变成单个细胞的悬液，并且尽量不让细胞“受伤”。

  • 新鲜组织 vs. 冷冻样本：
    • 理想情况是新鲜组织。冷冻和解冻过程会破坏细胞膜，导致大量细胞死亡，释放出 RNA，污染其他细胞，产生所谓的“环境 RNA”（ambient RNA）。
    • 如果必须用冷冻样本，要确保冷冻速度够快，解冻过程够温和。
  • 酶消化：
    • 不同的组织用不同的酶。肝脏、胰腺这种结缔组织多的，可能需要胶原蛋白酶、透明质酸酶混合消化。
    • 心脏、肌肉需要更强的机械分散配合酶消化。
    • 注意：消化时间不能太长，温度也不能太高（通常 37℃），否则细胞应激反应会改变基因表达谱。
  • 过滤与离心：
    • 用 40μm 或 70μm 的细胞过滤器去除团块和杂质。
    • 离心后，用冰冷的 PBS 或专用缓冲液重悬细胞。

2. 活率检测与清洗

这一步至关重要。你必须确保你的细胞是活的，而且没有碎片。

  • 使用台盼蓝或自动细胞计数仪（如 Countess 或 NucleoCounter）检测活率。
  • 目标：活率最好 >85%，甚至 >90%。死细胞会破裂，释放 RNA，被微流控芯片误捕获，导致数据污染。
  • 如果活率不够，可以用磁性分选珠（如 Dead Cell Removal Kit）去除死细胞。

3. 微流控芯片加载（以 10x Genomics 为例）

目前最主流的平台是 10x Genomics。它的原理很巧妙：

  • GEMs（微滴）：每个微滴里包含一个 Gel Bead（凝胶珠，上面连着独特的条形码 RNA）、一个油包水液滴、以及最多一个细胞。
  • 裂解与条形码：在微滴里，细胞被裂解，RNA 释放出来，与 Gel Bead 上的条形码结合。这样，每一条 RNA 分子都被打上了“我是从哪个微滴（哪个细胞）里来的”标签。
  • 逆转录：在微滴内完成逆转录，生成带条形码的 cDNA。

关键点：

  • 载荷效率：你需要调整细胞浓度，使得大部分微滴里只有一个细胞，而没有微滴（空滴）或双细胞微滴（两个细胞在一个微滴里）。通常目标捕获细胞数是预期细胞数的 1.5 倍左右，以减少空滴。
  • 双细胞比例：如果加载细胞太多，双细胞比例会升高，影响后续分析。一般建议双细胞比例 <5-10%。

4. cDNA 扩增与文库构建

  • 将微滴破乳，收集 cDNA。
  • 对 cDNA 进行 PCR 扩增。
  • 片段化、加接头、PCR 扩增，最终构建测序文库。
  • 质检：用 Bioanalyzer 或 TapeStation 检测文库片段大小分布和浓度。

5. 测序

  • 通常使用 Illumina 测序平台。
  • 测序深度：每个细胞建议测序深度在 20,000 - 50,000 reads 以上。如果研究的是低丰度转录本或想要更精细的亚群划分，可能需要更深（100,000+ reads）。
  • 读长：通常是 2x90bp 或 2x150bp。Read 1 包含细胞条形码和 UMI（唯一分子标识符），Read 2 包含 cDNA 序列。

三、 数据质控：垃圾进，垃圾出

拿到 FASTQ 文件后，别急着跑分析。先看看数据质量怎么样。这一步决定了你后面所有结果的可靠性。

1. 原始数据质控（FastQC）

使用 FastQC 工具检查每个样本的原始序列质量。

  • Per base sequence quality：每个位置的碱基质量是否达标（Q30 以上）。
  • Per sequence GC content：GC 含量分布是否异常。
  • Adapter contamination：是否有接头污染。

如果发现质量问题，可能需要用 Trimmomatic 或 Cutadapt 进行剪接和过滤。

2. 细胞过滤：什么是“好”细胞？

这是 scRNA-seq 分析中最重要的一步。我们需要定义什么是“好细胞”，什么是“噪音”。

常用的质控指标有三个：

  • nCount_RNA（检测到的基因数）：反映细胞捕获的 RNA 分子数量。太低可能是空滴或破损细胞；太高可能是双细胞或多细胞。
  • nFeature_RNA（检测到的特征基因数）：与 nCount 类似，但更关注基因种类。
  • percent.mt（线粒体基因百分比）：线粒体基因占比高，通常意味着细胞膜破损，细胞质 RNA 流失，只剩下线粒体 RNA。这是死细胞或应激细胞的标志。

典型的过滤策略（以 Seurat 为例）：

# 加载 Seurat 对象
library(Seurat)

# 查看每个细胞的三个指标分布
VlnPlot(object = scRNA_seq, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"), ncol = 3)

# 设定阈值进行过滤
# 假设我们观察到：
# nFeature_RNA 在 200-2500 之间比较合理
# percent.mt 小于 5-10% 比较合理（取决于组织类型，血液可能更低，心肌可能更高）
scRNA_seq <- subset(scRNA_seq, 
                    subset = nFeature_RNA > 200 & 
                             nFeature_RNA < 2500 & 
                             percent.mt < 10)

注意：不同的组织类型，质控阈值不同。

  • 血液细胞：线粒体基因占比通常很低（%），因为红细胞没有细胞核，白细胞代谢相对温和。
  • 心肌、肝脏细胞：线粒体占比可能较高（10-20%），因为这些细胞能量需求高。
  • 肿瘤细胞：可能因代谢重编程而线粒体基因占比异常。

3. 去除双细胞（Doublets）

即使你控制了加载浓度，仍会有一部分微滴包含两个细胞。这些“双细胞”会干扰聚类分析，因为它们可能表达两种不同细胞类型的标志基因，看起来像是一个“过渡态”或“杂交”细胞。

检测方法：

  • Simulation-based approaches：如 Scrublet 或 DoubletFinder。它们通过模拟人工双细胞，计算每个真实细胞与模拟双细胞的相似度，从而预测双细胞。
  • 基于表达的方法：如果一个细胞高表达两个完全不相关的细胞类型标志基因，很可能是双细胞。
# 使用 Scrublet 检测双细胞（Python 环境）
import scrublet as scr
# 假设 counts_matrix 是基因表达矩阵
doublet_scores, expected_doublet_rate = scrublet.scrub_doublets(counts_matrix)
# 设定阈值，标记高双细胞分数为双细胞

4. 标准化与归一化

不同细胞捕获的 RNA 分子数量不同（测序深度差异），直接比较基因表达量是没有意义的。我们需要进行归一化，消除技术噪音。

  • Log-normalization：最常用。先将计数除以总计数（得到比例），再加 1 后取对数。 $\( \text{normalized} = \log(1 + \text{count} / \text{total\_counts} \times 10000) \)$ 这相当于将每个细胞的总表达量缩放到相同的水平（通常是 10,000）。
# Seurat 默认流程
scRNA_seq <- NormalizeData(scRNA_seq, normalization.method = "LogNormalize", scale.factor = 10000)

四、 特征基因筛选与降维

1. 筛选高变基因（HVGs）

并不是所有基因都有分析价值。看那些在所有细胞中都稳定表达的“管家基因”（如 GAPDH, ACTB），它们不能帮助我们区分细胞类型。我们需要找的是高变基因（Highly Variable Genes, HVGs），即在不同细胞间表达差异很大的基因。

这些基因往往与细胞类型特异性、细胞状态有关。

# 在 Seurat 中识别 HVGs
scRNA_seq <- FindVariableFeatures(scRNA_seq, selection.method = "vst", nfeatures = 2000)
# 绘制 HVGs 图，看看哪些基因变异系数最大
Visualization::FeatureScatter(scRNA_seq, feature1 = "mean.expr", feature2 = "dispersion")

通常选择前 2000-3000 个 HVGs 进入后续分析。

2. 数据缩放（Scaling）

将基因表达数据标准化，使均值为 0，方差为 1。这一步是为了让后续的主成分分析（PCA）能够公平地对待每个基因，避免高表达基因主导结果。

scRNA_seq <- ScaleData(scRNA_seq, features = all.genes)

3. 主成分分析（PCA）

PCA 是一种线性降维技术。它将高维的基因表达数据（比如 20000 个基因）投影到低维空间（比如 50 个主成分），同时保留尽可能多的方差信息。

  • 前几个 PCs 通常捕获的是技术噪音（如测序深度、线粒体效应）。
  • 中间的 PCs 通常包含生物学信号（细胞类型差异）。
  • 后面的 PCs 可能又是噪音或细微的生物学差异。

如何确定用哪些 PCs？

  • 查看 Elbow Plot（手肘图）：寻找拐点。
  • 查看 PC Heatmap：观察哪些 PC 与生物学变量相关。
  • 目前常用的方法是结合 JackStraw 或 BIC 方法，或者简单地使用前 10-30 个 PCs。
# 计算 PCA
scRNA_seq <- RunPCA(scRNA_seq, features = variable.features)

# 绘制 Elbow Plot
ElbowPlot(object = scRNA_seq)

# 选择 PCs（例如选择前 20 个）
DimReduce(object = scRNA_seq, dims = 1:20)

4. 非线性降维：UMAP 和 t-SNE

PCA 是线性的，可能无法捕捉复杂的细胞结构。UMAP（Uniform Manifold Approximation and Projection）和 t-SNE 是两种流行的非线性降维方法，用于将数据可视化到 2D 或 3D 空间。

  • t-SNE：早期常用，能很好地保留局部结构，但计算慢，且全局结构有时难以解读。
  • UMAP：近年来更流行。它在保留局部结构的同时，也能更好地保持全局结构，且计算速度更快。

注意：UMAP/t-SNE 的结果高度依赖于前面的 PCA 结果。所以，先做好 PCA，再用 UMAP/t-SNE 可视化。

# 运行 UMAP
scRNA_seq <- RunUMAP(scRNA_seq, dims = 1:20)

# 绘制 UMAP 图
DimPlot(scRNA_seq, reduction = "umap")

五、 细胞聚类与注释

1. 聚类（Clustering）

聚类的目标是将相似的细胞分到同一组。常用的算法是 Louvain 或 Leiden 算法。

  • 分辨率（Resolution）：这是一个超参数。分辨率越高，聚类越细，细胞群分得越碎；分辨率越低，聚类越粗，细胞群合并得越多。
  • 通常从一个较低的分辨率（如 0.5）开始，逐步增加（0.8, 1.0, 1.2…），观察聚类结果的变化。
# 构建 KNN 图并聚类
scRNA_seq <- FindNeighbors(scRNA_seq, dims = 1:20)
scRNA_seq <- FindClusters(scRNA_seq, resolution = 0.5)

2. 细胞注释（Cell Type Annotation）

聚类完成后，你需要知道每个簇代表什么细胞类型。这通常通过标记基因（Marker Genes）来实现。

方法一：已知标记基因对照 查阅文献，找到各个细胞类型的经典标志基因，看它们在哪些簇中高度表达。

# 寻找每个簇的标记基因
cluster_markers <- FindAllMarkers(scRNA_seq, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

# 查看某个簇的特异性表达基因
VlnPlot(scRNA_seq, features = c("CD3D", "MS4A1", "PPBP", "FCGR3A"), group.by = "seurat_clusters")
  • CD3D, CD2, ITM2C：T 细胞
  • MS4A1 (CD20), CD79A：B 细胞
  • PPBP, CBX3：血小板/巨核细胞
  • FCGR3A, MS4A7：单核/巨噬细胞
  • EPCAM：上皮细胞
  • PECAM1, CLDN5：内皮细胞

方法二：自动注释工具 现在有很多自动注释工具，如 SingleR、CellTypist、scCATCH 等。它们利用已知的参考数据集（如 Human Cell Atlas），将你的细胞与参考数据比对，自动给出注释建议。

# 使用 SingleR 进行自动注释
library(SingleR)
# 参考数据集
ref <- HumanPrimaryCellAtlasData()
# 预测细胞类型
predictions <- SingleR(test = scRNA_seq@assays$RNA@counts, ref = ref, labels = ref$label.main)
# 将预测结果添加到 Seurat 对象
scRNA_seq$predicted.cell.type <- predictions$labels

方法三：人工审阅 自动注释可能会有错误。一定要结合 UMAP 图、VlnPlot、FeaturePlot 等多种可视化手段，人工审阅每个簇的标记基因，确认注释的合理性。

六、 差异表达分析：寻找关键基因

当你确定了不同的细胞类型或状态后，下一步就是找出它们之间的差异。

1. 组间差异分析

比较不同条件（如：处理组 vs. 对照组）下，同一细胞类型中基因表达的变化。

方法：

  • 提取特定细胞类型的所有细胞。
  • 使用 Wilcoxon 秩和检验、Mann-Whitney U 检验或 DESeq2/edgeR（如果考虑批量效应）进行差异分析。
  • Seurat 的 FindMarkers 函数默认使用 Wilcoxon 检验。
# 比较 T 细胞在对照组和处理组之间的差异
t_cells <- subset(scRNA_seq, idents = "T cells")
t_cells$condition <- metadata(t_cells)$condition # 假设 metadata 中有条件信息

diff_markers <- FindMarkers(t_cells, 
                            ident.1 = "Treatment", 
                            ident.2 = "Control",
                            min.pct = 0.1,
                            logfc.threshold = 0.25)

解读结果：

  • LogFC：对数折叠变化，表示基因表达上调或下调的程度。
  • p.val：p 值，表示差异的显著性。
  • avg_log2FC：平均 log2 折叠变化。
  • pct.1 和 pct.2：分别在两个组中表达该基因的细胞比例。

通常，我们会选择 LogFC > 0.25 且 p.adj < 0.05 的基因作为显著差异基因。

2. 差异基因可视化

  • 火山图：展示所有基因的 LogFC 和 -log10(p值) 的关系，一眼看出哪些基因是显著上调或下调的。
  • 热图：展示差异基因在不同样本中的表达模式。
  • 小提琴图：展示关键差异基因在不同组间的表达分布。

”`r

绘制火山图

EnhancedVolcano(diff_markers,

            lab = rownames(diff_markers),
            x = 'p_val
-- 展开阅读全文 --