第三阶段 3.2：RNA 分析方向

对应大纲第三阶段 3.2 RNA 分析方向 目标：理解 RNA-seq 和 scRNA-seq 的完整分析流程，能从 FASTQ 跑到差异表达结果和细胞分群图 预计完成时间：3-4 周

3.2-rna-analysis-tutorial.md 结构

批量 RNA-seq 完整流程：

• 从"为什么不能用 BWA"切入，用 mRNA 跨外显子读段解释剪接感知比对的必要性
• STAR 比对——建索引、运行命令、Log.final.out 关键指标解读
• featureCounts 定量——输出格式、为什么不能直接比原始计数
• DESeq2（R）——从创建对象到提取 DEG，每个步骤的参数逻辑，padj vs pvalue 的区别
• 火山图 + 热图可视化代码
• GO/KEGG 富集分析（clusterProfiler），结果表格逐字段解读
• Salmon 伪比对方案（快 10 倍，大样本推荐）

单细胞 RNA-seq：

• 批量 vs 单细胞的核心区别（平均 vs 单细胞分辨率）
• 液滴法（10x Genomics）原理和 UMI 去重逻辑
• Cell Ranger 流程和 metrics_summary 质控指标
• Scanpy 完整流程——读取数据→质控（线粒体比例）→标准化→高变基因→PCA→UMAP→Leiden 聚类→标志基因
• 细胞类型注释（手动标记 + CellTypist 自动注释）
• RNA velocity 和 LIANA 细胞通讯的代码入口
• Seurat（R）等价实现
• 批量 vs 单细胞对比总结表

---

先建立直觉：RNA-seq 在测什么？

DNA 测序回答的是"基因组里有什么变异"，是静态信息。

RNA-seq回答的是"此时此刻，哪些基因正在被激活，激活程度如何"，是动态信息。

状态 A（正常肝细胞）         状态 B（肿瘤肝细胞）
  ↓ 提取所有 mRNA              ↓ 提取所有 mRNA
  ↓ 反转录为 cDNA              ↓ 反转录为 cDNA
  ↓ 测序                       ↓ 测序
  基因A: 1000 reads             基因A: 50 reads   ← 下调 20 倍
  基因B:  200 reads             基因B: 4000 reads  ← 上调 20 倍
  基因C:  500 reads             基因C: 480 reads   ← 无变化

reads 数量代表该基因的 mRNA 丰度，丰度越高表达量越高。RNA-seq 的核心任务就是：

1. 数清楚每个基因有多少 reads（定量）
2. 比较两组样本之间哪些基因定量差异显著（差异表达）

---

3.2.1 RNA-seq：批量转录组测序

流程总览

FASTQ（原始测序读段）
  ↓ 质控：FastQC + Trimmomatic/fastp
FASTQ（清洁数据）
  ↓ 比对：STAR（识别可变剪接）
BAM（读段定位到基因组坐标）
  ↓ 定量：featureCounts / Salmon
Count Matrix（基因 × 样本的计数矩阵）
  ↓ 差异表达：DESeq2 / edgeR（R）
DEG 列表（差异表达基因 + 统计量）
  ↓ 功能解读：KEGG/GO 富集分析
Pathway 结果（哪条信号通路被激活/抑制）

为什么 RNA-seq 比对不能用 BWA？

这是 RNA 分析里最常见的新手错误。

DNA 是连续的，读段可以无缝比对。但 mRNA 经过了剪接（Splicing），内含子已经被去掉。如果一条 150bp 的 mRNA 读段恰好跨越两个外显子的边界，它在基因组上对应的位置之间有一段内含子间隔（可能几千 bp）：

mRNA 读段（150bp）：
ATCGATCG...（外显子1末端 75bp）TGCATGCA...（外显子2开头 75bp）

基因组实际结构：
外显子1...[ 内含子，1000bp ]...外显子2

BWA 看到：这条读段中间有 1000bp 的空缺，不知道怎么比对 → 丢弃
STAR 看到：识别剪接位点，正确对齐到两个外显子 → 保留，CIGAR 为 75M1000N75M

所以 RNA-seq 必须用剪接感知（splice-aware）的比对工具，STAR 是目前最主流的选择。

第一步：质控

RNA-seq 的质控和 DNA-seq 类似，但需要额外关注：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc/

# RNA-seq 特有的质控工具：RSeQC（检查读段在基因上的分布）
pip install RSeQC
infer_experiment.py -r genes.bed -i sample.bam    # 判断是否为链特异性测序
read_distribution.py -r genes.bed -i sample.bam  # 读段落在基因哪些区域

RNA-seq 质控的特殊关注点：

指标	正常范围	异常含义
rRNA 比例	<5%（polyA 选择建库）	>20% 说明 rRNA 去除不彻底
外显子覆盖均匀度	均匀	5' 或 3' 端偏斜说明 RNA 降解
重复读段比例	20-40%（正常）	>60% 说明文库复杂度不足
链特异性	按建库方式确定	判断错误会导致反义链基因被计入

第二步：STAR 比对

# 第一步：建立 STAR 基因组索引（只需一次，约 30-60 分钟）
STAR \
    --runMode genomeGenerate \
    --genomeDir ./hg38_star_index/ \
    --genomeFastaFiles hg38.fa \
    --sjdbGTFfile gencode.v44.annotation.gtf \  # 提供基因注释，提升剪接识别精度
    --sjdbOverhang 149 \                         # 读段长度 - 1（150bp 的读段填 149）
    --runThreadN 16

# 第二步：比对
STAR \
    --runMode alignReads \
    --genomeDir ./hg38_star_index/ \
    --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
    --readFilesCommand zcat \               # 输入是 gz 压缩的
    --outSAMtype BAM SortedByCoordinate \   # 直接输出排序后的 BAM
    --outSAMattributes NH HI AS NM \        # 常用标签
    --quantMode GeneCounts \               # 同时输出基因计数（简单情况可用）
    --runThreadN 16 \
    --outFileNamePrefix ./star_output/sample_

STAR 输出文件：

文件	内容
*Aligned.sortedByCoord.out.bam	主要输出：排序后的 BAM 文件
*Log.final.out	比对统计（最重要的质控文件）
*SJ.out.tab	检测到的剪接位点
*ReadsPerGene.out.tab	基因计数（--quantMode GeneCounts 时产生）

Log.final.out 关键指标：

Uniquely mapped reads %  |  95.23%      ← 核心指标，正常 >85%
Multi-mapped reads %     |  3.15%       ← 比对到多处，较低为好
Unmapped reads %         |  1.62%       ← 正常 <5%
% of reads unmapped: too short | 0.87%  ← 通常是 Trimming 过度
Splicing:
  Number of splices: Total | 52,341,021
  Number of annotated splices | 98.7%   ← 注释内的剪接比例，正常应 >90%

第三步：featureCounts 定量

featureCounts 根据 GTF 注释，统计每个基因被多少读段覆盖：

featureCounts \
    -T 8 \                           # 8 线程
    -p \                             # 双端测序（paired-end）
    -s 2 \                           # 链特异性（0=无，1=正链，2=反链，需根据建库方式确定）
    -a gencode.v44.annotation.gtf \ # 基因注释文件
    -o counts_matrix.txt \           # 输出文件
    sample1.bam sample2.bam sample3.bam sample4.bam  # 可以一次处理多个样本

# 提取计数矩阵（去掉前几行注释）
cut -f1,7- counts_matrix.txt | tail -n +2 > counts_clean.txt

输出格式：

Geneid           sample1  sample2  sample3  sample4
ENSG00000223972  0        0        1        0
ENSG00000227232  156      143      178      162
ENSG00000243485  8        12       7        11
ENSG00000237613  0        0        0        0
...

这个矩阵就是后续所有分析的输入：行是基因，列是样本，值是原始 read 计数。

重要：为什么不能直接比较原始计数？

基因 A 在 样本1 中有 1000 reads，在 样本2 中有 500 reads
→ 能说样本2 中基因 A 表达量低一半吗？
→ 不能，因为样本1 可能总共测了 5000 万 reads，样本2 只测了 2000 万

需要标准化处理，消除测序深度和基因长度的影响。

第四步：DESeq2 差异表达分析（R）

DESeq2 是 RNA-seq 差异表达分析的行业标准工具，使用负二项分布模型处理计数数据。

library(DESeq2)

# ── 1. 准备输入数据 ──
count_matrix <- read.table("counts_clean.txt", header=TRUE, row.names=1)
# 列名：sample1_tumor, sample2_tumor, sample3_normal, sample4_normal
count_matrix <- round(count_matrix)  # 必须是整数

# 样本分组信息（必须和矩阵列顺序一致）
sample_info <- data.frame(
    row.names = colnames(count_matrix),
    condition = c("tumor", "tumor", "normal", "normal")
)

# ── 2. 创建 DESeqDataSet 对象 ──
dds <- DESeqDataSetFromMatrix(
    countData = count_matrix,
    colData   = sample_info,
    design    = ~ condition           # 比较 condition 变量
)

# ── 3. 过滤低表达基因（至少 3 个样本中有 ≥ 10 reads） ──
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]
cat("过滤后保留基因数：", sum(keep), "\n")

# ── 4. 运行 DESeq2 ──
dds <- DESeq(dds)
# 内部步骤：估计大小因子（标准化）→ 估计离散度 → 负二项 GLM → Wald 检验

# ── 5. 提取结果 ──
res <- results(dds,
    contrast = c("condition", "tumor", "normal"),  # 比较方向：tumor vs normal
    alpha    = 0.05                                 # FDR 阈值
)
summary(res)
# 输出：
# LFC > 0 (up)   : 1234（肿瘤中上调）
# LFC < 0 (down) : 987（肿瘤中下调）
# outliers        : 23
# low counts      : 4521

# ── 6. 结果排序和过滤 ──
res_df <- as.data.frame(res)
res_df$gene_id <- rownames(res_df)

# 筛选显著差异基因（DEG）
deg <- res_df[
    !is.na(res_df$padj) &
    res_df$padj < 0.05 &
    abs(res_df$log2FoldChange) >= 1,   # |log2FC| ≥ 1 即变化倍数 ≥ 2
]
cat("显著差异基因数：", nrow(deg), "\n")
write.csv(deg, "DEG_results.csv", row.names=FALSE)

结果字段解读：

列名	含义	判断标准
baseMean	该基因在所有样本中的平均表达量（标准化后）	太低（<10）的结果不可靠
log2FoldChange	以 2 为底的表达倍数变化（tumor/normal）	>1 表示上调 2 倍，<-1 表示下调 2 倍
lfcSE	log2FC 的标准误差	越小越精确
stat	Wald 检验统计量
pvalue	原始 p 值	不用原始 p 值，用校正后的
padj	BH 多重检验校正后的 p 值（FDR）	< 0.05 为显著

为什么用 padj 而不是 pvalue？ RNA-seq 同时检验几万个基因，随机期望 5% 的基因（即约 1000 个）会得到 p < 0.05，即使它们实际无差异。FDR 校正控制的是"在所有显著结果中，假阳性的预期比例"。

第五步：可视化

library(ggplot2)
library(pheatmap)

# ── 火山图（Volcano Plot）──
res_df$significance <- "不显著"
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange > 1]  <- "上调"
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange < -1] <- "下调"

ggplot(res_df[!is.na(res_df$padj), ], aes(x=log2FoldChange, y=-log10(padj))) +
    geom_point(aes(color=significance), alpha=0.6, size=1) +
    scale_color_manual(values=c("上调"="red", "下调"="blue", "不显著"="grey")) +
    geom_vline(xintercept=c(-1, 1), linetype="dashed") +
    geom_hline(yintercept=-log10(0.05), linetype="dashed") +
    labs(title="Tumor vs Normal 差异表达火山图",
         x="log2(Fold Change)", y="-log10(adjusted p-value)") +
    theme_bw()

# ── 热图（Heatmap）——Top 50 差异基因 ──
top50 <- head(deg[order(deg$padj), ], 50)
# 提取标准化表达量（vst 转换，适合可视化）
vst_data <- assay(vst(dds))
mat <- vst_data[top50$gene_id, ]
mat <- t(scale(t(mat)))  # 每行（基因）Z-score 标准化

pheatmap(mat,
    annotation_col = sample_info,
    show_rownames  = FALSE,
    cluster_cols   = TRUE,
    cluster_rows   = TRUE,
    color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
    main = "Top 50 差异基因表达热图"
)

第六步：富集分析（KEGG/GO Pathway）

找到 DEG 之后，更重要的是理解这些基因共同代表什么生物学过程：

library(clusterProfiler)
library(org.Hs.eg.db)   # 人类基因注释数据库

# 把 Ensembl ID 转换为 Entrez ID（KEGG 需要）
gene_ids <- bitr(
    geneID   = deg$gene_id,
    fromType = "ENSEMBL",
    toType   = c("ENTREZID", "SYMBOL"),
    OrgDb    = org.Hs.eg.db
)

# ── GO 富集分析 ──
go_result <- enrichGO(
    gene          = gene_ids$ENTREZID,
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",          # BP=生物学过程, MF=分子功能, CC=细胞组分
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    qvalueCutoff  = 0.2,
    readable      = TRUE           # 把 Entrez ID 转回基因名
)

# ── KEGG 通路富集 ──
kegg_result <- enrichKEGG(
    gene          = gene_ids$ENTREZID,
    organism      = "hsa",         # hsa = Homo sapiens
    pvalueCutoff  = 0.05
)

# ── 可视化 ──
dotplot(go_result, showCategory=20, title="GO Biological Process")
barplot(kegg_result, showCategory=15, title="KEGG Pathway")

富集分析结果解读：

Description              GeneRatio  BgRatio   pvalue    p.adjust  Count
cell cycle               42/1234    124/18000  1.2e-15   3.4e-13   42
DNA repair               28/1234    67/18000   8.9e-13   1.3e-10   28
immune response          65/1234    312/18000  2.1e-10   2.0e-08   65

• GeneRatio：你的 DEG 中有多少比例属于这个通路（42/1234 = 3.4%）
• BgRatio：全基因组中这个通路的基因比例（124/18000 = 0.69%）
• 两者比值越大（富集倍数越高）、p.adjust 越小，说明该通路越显著

替代方案：Salmon 伪比对定量（更快）

对于大样本量项目，可以跳过 BAM 比对，直接用 Salmon 做定量：

# 建立转录组索引（使用转录本序列，而非基因组）
salmon index \
    -t gencode.v44.transcripts.fa.gz \
    -d decoys.txt \
    -i salmon_index/ \
    --gencode

# 定量（每个样本约 5 分钟，远快于 STAR）
salmon quant \
    -i salmon_index/ \
    -l A \                          # 自动检测文库类型
    -1 sample_R1.fastq.gz \
    -2 sample_R2.fastq.gz \
    --validateMappings \
    -o salmon_output/sample/ \
    -p 8

Salmon 输出 quant.sf 文件，用 tximport 导入 R 后接 DESeq2：

library(tximport)
# 读取多个样本的 Salmon 结果
files <- list.files("salmon_output/", pattern="quant.sf", recursive=TRUE, full.names=TRUE)
txi <- tximport(files, type="salmon", tx2gene=tx2gene_map)  # tx2gene_map: 转录本→基因的对应表
dds <- DESeqDataSetFromTximport(txi, colData=sample_info, design=~condition)

---

3.2.2 单细胞 RNA-seq（scRNA-seq）

批量测序 vs 单细胞测序

批量 RNA-seq（Bulk RNA-seq）测的是一管细胞的平均：

1000 个细胞一起测 → 基因A = (细胞1的A + 细胞2的A + ... + 细胞1000的A) / 1000

问题：如果你的样本里有 800 个正常细胞和 200 个肿瘤细胞，
     基因A 在肿瘤细胞里高表达，但在批量测序里被"稀释"了

单细胞 RNA-seq（scRNA-seq）对每个细胞单独测序：

细胞1: 基因A=高, 基因B=低  → 肿瘤细胞特征
细胞2: 基因A=低, 基因B=高  → 免疫细胞特征
细胞3: 基因A=低, 基因B=低  → 正常上皮细胞
...
10000 个细胞 → 10000 行，每行是一个细胞的表达谱

scRNA-seq 能做批量测序做不到的事：

• 发现样本中有哪些细胞类型（细胞分型）
• 找稀有细胞亚群（如肿瘤干细胞，可能只占 0.1%）
• 重建细胞从一种状态到另一种状态的发育轨迹（Trajectory Analysis）
• 分析细胞间的通讯网络（Cell-Cell Communication）

实验技术：液滴法（10x Genomics）

目前最主流的 scRNA-seq 平台是 10x Genomics Chromium，原理：

1. 将细胞和含有 barcode 的凝胶珠包裹在油滴（液滴）中
2. 每个液滴理想状态下包含一个细胞和一颗珠
3. 珠上的 barcode 序列会标记该液滴中的 mRNA
4. 裂解细胞，mRNA 与珠上的 oligo-dT 结合，逆转录
5. 最终每条 cDNA 序列携带：
   - Cell Barcode（16bp）：标识来自哪个细胞
   - UMI（12bp，Unique Molecular Identifier）：标识原始 mRNA 分子
   - cDNA 序列：来自哪个基因

UMI 的作用：去除 PCR 扩增偏差

原始 mRNA 分子（1个）
  ↓ PCR 扩增 10 倍
10 条 reads → 如果没有 UMI，会认为这个基因有 10 个 mRNA 分子

有 UMI：10 条 reads 的 UMI 完全相同 → 去重后只计为 1 个 mRNA 分子

Cell Ranger：原始数据处理

10x Genomics 提供官方处理工具 Cell Ranger：

# 比对并生成细胞×基因计数矩阵
cellranger count \
    --id=sample1_output \
    --transcriptome=refdata-gex-GRCh38-2020-A/ \  # 10x 提供的参考数据库
    --fastqs=./fastq_dir/ \
    --sample=sample1 \
    --localcores=16 \
    --localmem=64

# 多样本合并（不同患者/条件）
cellranger aggr \
    --id=combined_output \
    --csv=aggregation.csv \
    --normalize=mapped

Cell Ranger 输出（在 outs/ 目录下）：

filtered_feature_bc_matrix/   ← 主要输出：过滤后的细胞矩阵
    barcodes.tsv.gz            ← 细胞 barcode 列表（行）
    features.tsv.gz            ← 基因列表（列）
    matrix.mtx.gz              ← 稀疏矩阵（大多数值是 0）
raw_feature_bc_matrix/        ← 未过滤版本（包含空液滴）
metrics_summary.csv           ← 关键质控指标
web_summary.html              ← 可视化报告

metrics_summary.csv 关键指标：

指标	正常范围	含义
Estimated Number of Cells	1000-10000（取决于实验设计）	检测到的细胞数
Median Genes per Cell	1500-4000	每个细胞检测到的基因数，越高越好
Median UMI Counts per Cell	3000-15000	每个细胞的 UMI 总数
Reads Mapped to Genome	>90%	比对率
Sequencing Saturation	40-80%	测序深度是否饱和

Scanpy：Python 单细胞分析（推荐路径）

import scanpy as sc
import numpy as np
import pandas as pd

sc.settings.verbosity = 2
sc.settings.figdir   = './figures/'

# ── 1. 读取数据 ──
adata = sc.read_10x_mtx(
    "filtered_feature_bc_matrix/",
    var_names="gene_symbols",   # 用基因名而非 Ensembl ID
    cache=True
)
adata                           # AnnData object: 8000 cells × 33000 genes

# ── 2. 质控：过滤低质量细胞和基因 ──
# 标记线粒体基因（MT- 开头，线粒体 mRNA 比例过高说明细胞死亡/破损）
adata.var["mt"] = adata.var_names.str.startswith("MT-")

sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)

# 可视化质控指标分布，目测确定阈值
sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
             jitter=0.4, multi_panel=True)

# 过滤
adata = adata[adata.obs.n_genes_by_counts > 200, :]  # 至少检测到 200 个基因
adata = adata[adata.obs.n_genes_by_counts < 5000, :]  # 排除双联体（doublet）
adata = adata[adata.obs.pct_counts_mt < 20, :]        # 线粒体比例 <20%
sc.pp.filter_genes(adata, min_cells=3)                # 基因在至少 3 个细胞中有表达

print(f"过滤后：{adata.n_obs} 细胞，{adata.n_vars} 基因")

# ── 3. 标准化 ──
sc.pp.normalize_total(adata, target_sum=1e4)  # 每个细胞标准化到总计数 10000
sc.pp.log1p(adata)                             # log(x+1) 转换，压缩分布
adata.raw = adata                              # 保存标准化后数据，供后续差异分析用

# ── 4. 高变基因筛选 ──
# 不是所有 33000 个基因都有信息量，只保留在细胞间变异最大的基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)
adata = adata[:, adata.var.highly_variable]    # 通常筛选出约 2000-3000 个基因

# ── 5. PCA 降维 ──
sc.pp.scale(adata, max_value=10)  # 零均值，单位方差（Z-score）
sc.tl.pca(adata, svd_solver="arpack", n_comps=50)
sc.pl.pca_variance_ratio(adata, n_pcs=50)  # 看多少个 PC 解释了主要变异

# ── 6. 构建邻居图 ──
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)  # 用前 30 个 PC，15 近邻

# ── 7. 聚类（Leiden 算法）──
sc.tl.leiden(adata, resolution=0.5)  # resolution 越大，cluster 越多越细
# 结果存在 adata.obs["leiden"]

# ── 8. UMAP 可视化 ──
sc.tl.umap(adata)
sc.pl.umap(adata, color=["leiden", "n_genes_by_counts", "pct_counts_mt"],
           frameon=False, ncols=3)

# ── 9. 寻找每个 cluster 的标志基因 ──
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)

# 查看 cluster 3 的 Top10 标志基因
sc.get.rank_genes_groups_df(adata, group="3").head(10)

细胞类型注释

找到 cluster 后，需要根据标志基因判断每个 cluster 是什么细胞类型：

# 已知的细胞类型标志基因（来自文献/数据库）
marker_genes = {
    "T 细胞":      ["CD3D", "CD3E", "CD3G"],
    "CD8 T 细胞":  ["CD8A", "CD8B"],
    "CD4 T 细胞":  ["CD4", "IL7R", "CCR7"],
    "B 细胞":      ["CD19", "MS4A1", "CD79A"],
    "NK 细胞":     ["GNLY", "NKG7", "KLRD1"],
    "Monocyte":    ["CD14", "LYZ", "CST3"],
    "DC":          ["FCER1A", "CST3"],
    "巨噬细胞":    ["CD68", "CD163", "MRC1"],
    "上皮细胞":    ["EPCAM", "KRT8", "KRT18"],
    "内皮细胞":    ["PECAM1", "VWF", "CDH5"],
    "成纤维细胞":  ["COL1A1", "COL1A2", "DCN"],
}

# 气泡图：每个 cluster 中各标志基因的表达
sc.pl.dotplot(adata, marker_genes, groupby="leiden", dendrogram=True)

# 小提琴图：特定基因在各 cluster 的表达分布
sc.pl.stacked_violin(adata, marker_genes, groupby="leiden", rotation=90)

# 手动给 cluster 打上细胞类型标签
cluster_to_celltype = {
    "0": "CD4 T 细胞",
    "1": "CD8 T 细胞",
    "2": "上皮细胞",
    "3": "Monocyte",
    "4": "B 细胞",
    "5": "NK 细胞",
    "6": "成纤维细胞",
    "7": "内皮细胞",
}
adata.obs["cell_type"] = adata.obs["leiden"].map(cluster_to_celltype)
sc.pl.umap(adata, color="cell_type", frameon=False, legend_loc="on data")

自动注释工具（减少主观判断）：

# CellTypist：基于机器学习的自动注释
import celltypist
from celltypist import models

model = models.download_models(model="Immune_All_Low.pkl")
predictions = celltypist.annotate(adata, model=model, majority_voting=True)
adata = predictions.to_adata()
sc.pl.umap(adata, color="majority_voting")

进阶分析：轨迹分析和细胞间通讯

# 伪时序分析（Pseudotime）：重建细胞分化轨迹
import scvelo as scv  # RNA velocity（需要含 unspliced reads 的数据）

scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)
scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap")

# 细胞间通讯分析：哪些细胞类型在"说话"
import liana
liana.mt.rank_aggregate(adata, groupby="cell_type", resource_name="consensus")
liana.pl.dotplot(adata, colour="magnitude_rank", size="specificity_rank",
                 source_labels=["Cancer", "T cells"], target_labels=["Macrophage"])

Seurat（R）：同等功能的 R 实现

library(Seurat)

# 读取数据
data <- Read10X("filtered_feature_bc_matrix/")
seurat_obj <- CreateSeuratObject(counts=data, min.cells=3, min.features=200)

# 质控
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern="^MT-")
seurat_obj <- subset(seurat_obj,
    subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 20)

# 标准流程
seurat_obj <- NormalizeData(seurat_obj)
seurat_obj <- FindVariableFeatures(seurat_obj, selection.method="vst", nfeatures=2000)
seurat_obj <- ScaleData(seurat_obj)
seurat_obj <- RunPCA(seurat_obj, npcs=30)
seurat_obj <- FindNeighbors(seurat_obj, dims=1:30)
seurat_obj <- FindClusters(seurat_obj, resolution=0.5)
seurat_obj <- RunUMAP(seurat_obj, dims=1:30)

# 可视化
DimPlot(seurat_obj, reduction="umap", label=TRUE)
FeaturePlot(seurat_obj, features=c("CD3D", "CD8A", "CD4", "MS4A1"))
VlnPlot(seurat_obj, features=c("nFeature_RNA", "percent.mt"), ncol=2)

# 标志基因
markers <- FindAllMarkers(seurat_obj, only.pos=TRUE, min.pct=0.25, logfc.threshold=0.25)
top10   <- markers %>% group_by(cluster) %>% slice_max(avg_log2FC, n=10)
DoHeatmap(seurat_obj, features=top10$gene)

---

3.2.3 批量 vs 单细胞对比总结

维度	批量 RNA-seq	scRNA-seq
分辨率	样本平均	单细胞
每样本细胞数	百万级	通常 5000-20000
检测基因数/细胞	N/A（全部）	约 1000-5000（稀疏）
测序深度	高（每基因几百 reads）	低（每细胞总 UMI~5000）
主要输出	DEG 列表 + Pathway	细胞类型图谱 + cluster 标志基因
每样本成本	~500 元	~5000-10000 元
数据分析难度	中等	较高（稀疏性、批次效应）
主要用途	两组条件比较	细胞组成、发育轨迹、罕见细胞
推荐工具（Python）	DESeq2（R）	Scanpy
推荐工具（R）	DESeq2	Seurat

---

本阶段完成标志

• 能解释为什么 RNA-seq 必须用 STAR 而不是 BWA
• 知道 UMI 在 scRNA-seq 中的作用是什么
• 能从 DESeq2 结果中正确筛选差异表达基因（padj < 0.05 且 |log2FC| ≥ 1）
• 能解释 Leiden clustering 和 UMAP 各自的作用
• 能根据标志基因判断一个 scRNA-seq cluster 是什么细胞类型
• 用 Scanpy 跑过完整的 scRNA-seq 分析（读取数据 → 质控 → 聚类 → 注释）