单细胞高变基因

在单细胞分析里,基因表达矩阵通常包含几万个基因,但其中只有一小部分基因在不同细胞间表现出足够大的生物学差异,能够帮助我们区分细胞类型或状态。选取“高变基因”(Highly Variable Genes, HVGs)这一步,就是为了从海量基因中筛出那部分“最具信息量”的基因,避免噪声和低变基因拖累后续的降维、聚类和可视化。


高变基因选择的“做什么”——分步解析#

  1. 输入数据:归一化 & 对数化后的表达

    • 在此之前,已经做了 normalize_total(将每个细胞的总 UMI 统一到 1e4)和 log1p 转换,得到“log1p_norm_count”层。

    • 目的是让不同细胞的测序深度可比,同时压缩极端值、近似正态分布。

  2. 计算每个基因的平均表达(μ)和方差(σ²)

    • 对所有细胞,分别统计每个基因在多少细胞中被检测到(用于过滤),以及在这些细胞里的均值与方差。
  3. 拟合均值–方差关系曲线

    • 生物学上,表达量很低的基因方差往往也小;表达量很高的基因方差也大,这是一个“漂移”趋势。

    • 我们要找的,是“比同等级表达量下更高方差”的基因,因此先用曲线(如 loess 或 Seurat 的方差标准化方法)拟合全局的 μ–σ² 关系。

  4. 计算“标准化离散度”

    • 对每个基因,求出它的实际方差与拟合曲线上的预期方差之差,或更常见地,计算 “归一化后的离散度指标”(standardized dispersion/z-score)。

    • 这个值越大,说明该基因在同样表达水平下,离散得越厉害——更可能反映细胞间真实的生物学差异。

  5. (可选)按批次分组再汇总

    • 如果你指定了 batch_key(比如 BatchSample_Site),Scanpy 会先在每个组(批次)内独立做上面几步,得到每个基因在每个组内的标准化离散度;

    • 然后再取各组离散度的中位数或其它汇总指标,对基因做全局排名,以确保所选 HVG 既在各批次中都有表现,而不是被某一个批次“带跑偏”了。

  6. 按照“标准化离散度”从高到低取前 N 个基因

    • 最常见的做法是直接指定 n_top_genes=5000(或其它数字),把排名靠前的 5,000 个基因标记为 adata.var['highly_variable'] = True

    • 这些基因集合就是你后续做 PCA/UMAP/聚类时的 “特征基因集”。


为什么要做这一步?#

  • 去噪:排除那些在所有细胞里都表现“很平均”、或几乎不表达的基因,它们对分辨细胞类型毫无帮助。

  • 减维:基因维度从数万降到几千,不仅加快计算,也聚焦在最有信号的区域。

  • 提升可解释性:高变基因往往包含各细胞类型或状态的标志性基因,更容易在后续可视化图(如特征基因热图、DotPlot)中找到生物学意义。


小结#

sc.pp.highly_variable_genes 这一步,本质上是在问:

“在你给定的表达矩阵里,哪些基因的细胞间变异度,明显超出了它们自身平均表达量所能解释的那部分?”

通过拟合均值–方差关系、计算标准化离散度,再按批次(可选)汇总,Scanpy 就能帮你自动挑出那批“最会变”的基因,作为后续降维和聚类的高信息量特征。

ScripyTCR术语表

克隆 (Clone) 和克隆型 (Clonotype): 在免疫学中,特别是T细胞或B细胞的语境下,一个“克隆”指的是一群起源于同一个祖细胞的细胞。这些细胞因为共享相同的T细胞受体 (TCR) 或B细胞受体 (BCR) 而被定义为一个“克隆型”。这意味着它们能够识别和响应同一种特定的抗原。

扩增 (Expansion): 当一个T细胞(或B细胞)遇到它能够识别的特定抗原(比如来自病原体或肿瘤细胞的抗原)并被激活时,它会开始快速分裂和增殖,产生大量具有相同TCR(即属于同一克隆型)的子细胞。这个过程就叫做“克隆扩增”。

克隆型大小 (Clonotype Size): 就是每个克隆型中包含的细胞数量。一个克隆型的细胞数量越多,说明这个克隆型发生的“扩增”程度越高。

AIRR

Adaptive Immune Receptor Repertoire.
Within the Scirpy documentation, we simply speak of immune receptors (IR).
The AIRR community defines standards around AIRR data. Scirpy supports the AIRR Rearrangement schema and complies with the AIRR Software Guidelines.

Alellically included B-cells

A B cell with two pairs of IG chains.
See Dual IR.

awkward array

TCR在癌症免疫治疗中研究实例

我们先来看第一个方面:在癌症免疫治疗中的发现。

癌症免疫治疗,尤其是免疫检查点抑制剂(比如PD-1/PD-L1抑制剂)的出现,彻底改变了很多癌症的治疗格局。而单细胞TCR分析在其中扮演了“侦探”和“导航员”的关键角色,帮助科学家们理解治疗为什么有效、为什么会失败,以及如何改进治疗策略。

以下是一些应用实例:

  1. 精准识别并描绘“抗癌特种兵”——肿瘤浸润淋巴细胞(TILs)

    • 发现:肿瘤组织中其实混杂着很多T细胞,它们被称为肿瘤浸润淋巴细胞(TILs)。但并非所有TILs都在积极对抗癌症。单细胞TCR分析能够同时读取单个TIL的TCR序列和其完整的基因表达谱。
    • 应用
      • 找出“功勋卓著”的T细胞克隆:通过分析TCR序列,可以识别出在肿瘤中显著扩增的T细胞克隆,这些通常是针对肿瘤抗原的“主力军”。
      • 描绘这些“主力军”的“画像”:结合基因表达数据,科学家可以了解这些关键T细胞克隆处于什么样的状态。比如,它们是充满活力的效应T细胞(正在积极杀伤癌细胞),还是已经“筋疲力尽”的耗竭T细胞(Exhausted T cells,虽然识别肿瘤但功能受损),或者是能够长期提供保护的记忆T细胞。
      • 案例:研究发现,在对免疫治疗有反应的患者肿瘤中,往往能找到一群特定的、具有“耗竭前体细胞 (precursor exhausted)”表型的CD8+ T细胞克隆。这些细胞虽然也表达一些耗竭标记,但仍保留了一定的增殖和杀伤潜力,并且在PD-1抑制剂治疗后能够被“重新激活”并扩增,从而有效控制肿瘤。单细胞TCR分析使得我们能够精确地将这些特定TCR克隆与这种关键的细胞状态联系起来。
  2. 追踪免疫治疗过程中T细胞克隆的“人生轨迹”

    • 发现:免疫治疗的效果不是一蹴而就的,T细胞的反应是一个动态变化的过程。
    • 应用
      • 监测克隆演变:通过在治疗前、治疗中、治疗后对患者的血液或肿瘤样本进行单细胞TCR分析,科学家可以追踪特定T细胞克隆(由其TCR序列定义)的数量变化和状态转变。
      • 揭示响应机制:例如,研究显示,在对PD-1抑制剂有反应的黑色素瘤患者中,治疗后肿瘤内不仅有原有T细胞克隆的扩增,有时还会出现新的、之前未被检测到的T细胞克隆的涌现和浸润。单细胞TCR分析还能告诉我们这些扩增或新出现的克隆具体是什么类型的T细胞,它们表达了哪些关键的效应分子或检查点分子。
      • 理解耐药机制:对于治疗无效或产生耐药的患者,单细胞TCR分析可以帮助揭示原因,比如有效的抗肿瘤T细胞克隆未能成功扩增,或者它们虽然扩增但迅速进入了深度耗竭状态,无法发挥作用。
  3. 指导新型T细胞疗法的开发

    • 发现:如果我们能找到那些TCR识别肿瘤特别“精准狠”的T细胞,就可以把这些TCR利用起来。
    • 应用
      • 筛选“精英TCR”:通过单细胞TCR分析,可以从那些成功清除肿瘤或对免疫治疗反应良好的患者的T细胞中,筛选出具有高效肿瘤识别能力的TCR序列。
      • 赋能工程化T细胞:一旦找到了这些“精英TCR”,就可以通过基因工程技术将编码这些TCR的基因导入患者自身的T细胞(或者健康的供体T细胞)中,制造出大量能特异性攻击癌细胞的“定制T细胞”,这就是所谓的TCR-T细胞疗法。单细胞TCR分析不仅提供了TCR序列,其偶联的转录组数据还能帮助评估这个TCR所在原始细胞的功能特性,进一步指导筛选。
  4. 发现新的免疫治疗靶点和生物标志物

    • 发现:通过比较响应者和非响应者体内T细胞的TCR库特征及其基因表达谱,可以找到与治疗效果相关的潜在靶点或预测标志物。
    • 应用:例如,如果发现某一类具有特定基因表达特征(比如高表达某些细胞因子或趋化因子受体)的T细胞克隆在响应者中显著富集,那么这些基因或通路就可能成为未来增强免疫治疗效果的新靶点。同时,特定的TCR库特征(如克隆多样性、特定克隆的丰度)也可能作为预测患者是否会对免疫治疗产生反应的生物标志物。

UMI

UMI,全称是Unique Molecular Identifier,中文意思是“独特分子标识符”。它在单细胞测序(尤其是scRNA-seq)中扮演着纠正定量偏差、实现更精确分子计数的“幕后英雄”角色。

  1. UMI是什么?

    UMI通常是一段短的(比如6-12个核苷酸)、随机合成的DNA序列。在单细胞实验的早期步骤中,这些UMI会被整合到我们想要研究的分子(比如mRNA)的反转录产物(cDNA)上。

  2. UMI是如何被添加的?

    在单细胞RNA测序流程中,当细胞裂解后释放出mRNA,这些mRNA会被带有oligo-dT(用于捕获mRNA的polyA尾)、细胞条形码(标记细胞身份)以及UMI的引物捕获。在反转录酶的作用下,当mRNA被逆转录成cDNA时,每一条原始的mRNA分子(理论上)就会在其对应的cDNA拷贝上“随机分配”并连接上一个独特的UMI序列。

  3. UMI的核心作用:消除PCR扩增偏好,实现精确计数

    这是UMI最重要的功能。在后续的实验步骤中,为了获得足够的DNA量进行测序,我们需要对产生的cDNA进行PCR扩增。然而,PCR扩增过程并非完美均一:

    • 扩增偏差 (Amplification Bias):有些cDNA分子可能因为序列特性或其他随机因素,比其他cDNA分子更容易或更有效地被扩增。这意味着,如果一条原始的mRNA分子A被扩增了100次,而另一条原始的mRNA分子B只被扩增了10次,仅仅通过计算测序后A和B的读长数量,我们就会错误地认为细胞中A的原始数量是B的10倍,而实际上它们可能都只有1个原始分子。

    UMI如何解决这个问题?

    • 因为(理想情况下)每个原始的mRNA分子在反转录时都带上了一个独特的UMI,所以即使经过多轮PCR扩增,所有源自同一个原始mRNA分子的扩增产物(cDNA拷贝)都会共享相同的细胞条形码相同的UMI。
    • 在数据分析时,对于某个特定基因,在某个特定细胞(由细胞条形码确定)中,我们不再是简单地计算所有映射到这个基因的测序读长数量。相反,我们去统计有多少种不同的UMI序列与这个基因的读长相关联。
    • 每种不同的UMI就代表了细胞中一个原始的mRNA分子

    举个例子: 假设在细胞X中,对于基因Y,我们有:

    • 一个原始的mRNA分子,我们称之为 mRNA_Y1。在反转录时,它被打上了一个独特的分子标识符,我们称之为 UMI_A。经过PCR扩增后,可能产生了50个都带有这个UMI_A的cDNA拷贝。
    • 细胞X中,基因Y的另一个原始mRNA分子,我们称之为 mRNA_Y2。在反转录时,它被打上了另一个独特的分子标识符,我们称之为 UMI_B。经过PCR扩增后,可能产生了20个都带有这个UMI_B的cDNA拷贝。

    在测序和数据分析后:

    • 如果我们仅仅计算所有与基因Y相关的测序读长数量,那么我们会得到 50 (来自mRNA_Y1的扩增产物) + 20 (来自mRNA_Y2的扩增产物) = 70条读长。
    • 但是,因为我们使用了UMI,我们会去识别这些读长所携带的UMI种类。我们发现,所有这70条读长,实际上只对应着两种不同的UMI序列:UMI_AUMI_B

    因此,尽管总共检测到了70条基因Y的测序读长,但通过识别不同的UMI,我们就知道在细胞X中,基因Y最初实际上只有2个原始的mRNA分子(一个对应UMI_A,另一个对应UMI_B)。

  4. 更准确的基因表达定量:

    通过这种方式,UMI能够帮助我们更准确地估计每个细胞中每个基因的原始转录本数量,从而得到更可靠的基因表达谱。这对于比较不同细胞间或不同条件下基因表达的细微差异至关重要。

总结一下,UMI就像是给每个原始RNA分子在变成cDNA时贴上了一个独一无二的小标签。无论这个分子后来被复制了多少次,通过识别这些独特的小标签,我们就能追溯到最初有多少个原始分子,从而实现更精确的定量。

R-seurat

R-seurat#

1. Seurat对象的整体结构#

Seurat对象主要包含以下几大组件:

  • Assays:存储表达矩阵和相关的数据(如空间转录组数据,通常为"Spatial")。

  • Meta data:样本、细胞、或者空间spot级别的元数据。

  • Dimensional Reductions:PCA、UMAP、tSNE等降维分析的结果。

  • Images:空间转录组中特有的组件,存储组织图像信息和空间坐标。

2. 典型结构:10x Visium#

SeuratObject
├── assays
│   └── Spatial(默认)
│       ├── counts:原始计数矩阵
│       ├── data:标准化表达矩阵(log normalized)
│       └── scale.data:经过缩放的表达数据(比如scale之后)
├── meta.data
│   └── 空间spot相关信息(细胞类型、样本条件等)
├── reductions
│   ├── pca:主成分分析结果
│   └── umap/tSNE:降维可视化结果
└── images
    └── slice1(空间图像信息)
        ├── coordinates:空间坐标数据(x, y位置)
        ├── image:组织的原始切片图像
        ├── scale factors:图像与空间坐标之间的尺度因子
        └── spot radius:空间spots的半径

3. 访问Seurat对象内容的示例代码#

library(Seurat)

# 假设数据对象名为obj
# 查看Seurat对象结构
str(obj)

# 提取空间坐标
coords <- GetTissueCoordinates(obj)

# 提取表达矩阵
expr_matrix <- GetAssayData(obj, slot = "data", assay = "Spatial")

# 提取空间图片
image <- obj@images$slice1@image

Python环境(scanpy/Squidpy)#

1.AnnData对象的整体结构#

AnnData对象
├── X:主表达矩阵(通常为基因表达)
├── obs:观测单元元数据(细胞或spot的元信息)
├── var:基因层面的信息
├── uns:非结构化数据(通常存储如图像、空间坐标元数据)
├── obsm:观测单元的多维信息(如PCA、空间坐标)
├── layers:可选项,存储多个版本的表达数据(如raw counts、normalized数据)
└── obsp/varp:可选,存储pairwise关系数据(如空间邻域图)

2. 空间转录组的典型结构#

以10x Visium数据为例:

Cellchat

得到cellchat对象后的可视化分析#

需要的包#

library(CellChat)
library(zellkonverter)
library(reticulate)
library(SingleCellExperiment)
library(patchwork)
library(presto)
library(Matrix)
library(future)
library(NMF)

读取和合并数据#

Hot_cell = readRDS("/share/home/sunLab/yanghc/project/TNBC_scRNA/GSE246613/cellchat/Hot_2.rds")
nonHot_cell = readRDS("/share/home/sunLab/yanghc/project/TNBC_scRNA/GSE246613/cellchat/nonHot_2.rds")

# 1. 合并两个 CellChat 对象  
object.list <- list(Hot = Hot_cell, nonHot = nonHot_cell)  
cellchat <- mergeCellChat(object.list, add.names = names(object.list))  # 合并后自动在 cellchat@net 等 slot 中保存两个数据集的信息 :contentReference[oaicite:0]{index=0}

分组间的dotplot可视化#

这里注意,netVisual_bubble()函数中有两个参数,第一个是comparison = c(1, 2)它会决定数据集的先后位置,例如在上一步的合并中 list(Hot = Hot_cell, nonHot = nonHot_cell)这样代表,Hot在1这个位置,nonHot在2这个位置。此时使用第二个参数max.dataset = 1,就会显示在Hot组中更强的信号,反之=2的话则显示在nonHot组中更强的信号也就是上调信息。

对于不同的对比画图需求,只需要改动第二个参数max.dataset = 1即可。如果同时改动则等于没有改动。总结如下:

  • 当你写 comparison = c(1, 2) 时,第 1 个子数据集是 Hot,第 2 个是 nonHot;此时 max.dataset = 2 就会高亮 nonHot 组中更强的信号。

  • 如果你改成 comparison = c(2, 1),那在这一对比较里第 1 个就变成 nonHot,第 2 个变成 Hot;这时候再用 max.dataset = 2,就会高亮 Hot 组中更强的信号。