转录组学

> 本章学习目标：掌握转录组学的核心概念、RNA-seq技术的完整分析流程、非编码RNA的分类与功能、单细胞转录组学的原理与方法，能够独立设计和执行转录组学研究项目。

中心法则与转录调控：
DNA →（转录）→ RNA →（翻译）→ 蛋白质
转录是基因表达的第一步，也是最关键的调控节点之一。生物体通过精确调控转录过程来实现：
- 时空特异性表达：不同组织、不同发育阶段表达不同的基因
- 环境响应：对外界刺激（温度、光照、营养、病原体等）做出快速反应
- 稳态维持：通过反馈调控维持细胞内环境稳定
转录调控的多层次机制：
1. 转录水平：转录因子结合、染色质可及性、DNA甲基化等
2. 转录后水平：mRNA剪接、编辑、稳定性调控
3. 翻译水平：翻译效率、miRNA调控
4. 翻译后水平：蛋白质修饰、定位、降解
转录组学主要聚焦于转录水平和转录后水平的研究。

第一代转录组技术：基因芯片（Microarray）
- 原理：将已知基因的寡核苷酸探针固定在芯片上，与标记的cDNA/cRNA杂交，通过荧光信号强度定量基因表达
- 局限：
- 只能检测已知序列（探针设计依赖已有基因组注释）
- 动态范围有限（通常2-3个数量级）
- 存在交叉杂交问题
- 难以检测新转录本和可变剪接
第二代转录组技术：RNA-seq
- 2008年首次报道（Mortazavi et al., Nature Methods）
- 原理：将RNA逆转录为cDNA，构建测序文库，进行高通量测序
- 优势：
- 不需要预先知道基因组序列
- 动态范围广（>5个数量级）
- 可同时检测已知和新转录本
- 可检测可变剪接、融合基因等
- 可检测SNP和等位基因特异性表达
第三代转录组技术：长读长RNA-seq
- PacBio Iso-Seq：产生全长转录本序列，无需组装即可确定完整异构体
- Nanopore direct RNA-seq：直接测序RNA分子，保留碱基修饰信息

（1）RNA提取与质控
- 使用Trizol或柱式提取法从样本中提取总RNA
- RNA完整性评估：Agilent Bioanalyzer检测RIN值（RNA Integrity Number）
- RIN > 7通常被认为是良好质量的RNA；RIN < 5可能导致3'端偏倚
（2）RNA富集或rRNA去除
- mRNA富集：使用oligo(dT)磁珠捕获带poly(A)尾的mRNA
- 优点：简单高效
- 缺点：会丢失不带poly(A)尾的RNA（如部分lncRNA、原核生物mRNA）
- rRNA去除：使用rRNA探针杂交去除rRNA（Ribo-Zero, Ribo-off等）
- 优点：保留所有非rRNA（包括mRNA、lncRNA、circRNA等）
- 缺点：操作稍复杂，成本略高
（3）RNA片段化与cDNA合成
- 片段化：使用二价阳离子（如Mg²⁺）在加热条件下将RNA随机打断成200-300 nt片段
- cDNA合成：
- 第一条链：使用随机六聚体或oligo(dT)引物，逆转录酶合成
- 第二条链：使用DNA聚合酶合成互补链（通常掺入dUTP以标记第二条链，用于链特异性文库）
（4）文库构建
- 末端修复：将cDNA片段末端补平
- A尾添加：在3'端添加单个A碱基
- 接头连接：连接带有barcode的测序接头
- PCR扩增：扩增文库至足够浓度
- 文库质控：检测文库大小分布（Agilent Bioanalyzer）和浓度（Qubit/qPCR）
（5）测序策略
- 单端测序（Single-end）：只测cDNA片段的一端，成本较低
- 双端测序（Paired-end）：测序两端，可更好地定位reads和识别可变剪接
- 链特异性测序（Strand-specific）：保留RNA链方向信息，可鉴定反义转录本

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, MultiQC)
    ↓
去除低质量和接头序列 (Trimmomatic, cutadapt)
    ↓
比对到参考基因组 (STAR, HISAT2)
    ↓
比对后处理 (排序, 索引)
    ↓
表达量定量 (featureCounts, HTSeq, Salmon)
    ↓
原始count矩阵
    ↓
数据探索 (PCA, 样本相关性热图)
    ↓
差异表达分析 (DESeq2, edgeR, limma)
    ↓
差异表达基因列表
    ↓
功能富集分析 (clusterProfiler, DAVID, GSEA)
    ↓
生物学解释和可视化

RNA-seq reads可能跨越一个或多个剪接位点（splice junction），这要求比对软件能够处理"分裂比对"（split alignment）。
STAR比对算法：
1. 种子生成：将read切分为多个短片段（seeds），在基因组上搜索精确匹配
2. 聚类：将基因组上接近的种子聚类为候选比对区域
3. 拼接：对于跨越内含子的reads，在候选区域之间构建剪接位点
4. 评分：对完整read进行比对评分，选择最优比对位置
STAR的特点：
- 速度快（比TopHat2快约50倍）
- 准确性高
- 内存需求大（需要约30-40 GB RAM用于人类基因组）
HISAT2比对算法：
- 使用基于图的比对方法（graph-based alignment）
- 利用FM-index和局部哈希实现快速搜索
- 对剪接变异的检测更敏感
- 内存效率高（<10 GB for human genome）

（1）基于count的方法
统计比对到每个基因区域的reads数。
- featureCounts：快速、准确，广泛使用的计数工具
- HTSeq-count：早期标准工具，灵活性高
（2）基于比对的定量（RPKM/FPKM/TPM）
原始count需要标准化才能进行样本内和样本间比较：
- RPKM（Reads Per Kilobase of exon per Million mapped reads）：
$$RPKM = \frac{ reads \times 10^9 }{ gene_length \times total_mapped_reads }$$
- FPKM（Fragments Per Kilobase of exon per Million mapped fragments）：
与RPKM类似，但使用fragment（pair-end测序中的一对reads）代替reads
- TPM（Transcripts Per Million）：
$$TPM_i = \frac{ RPKM_i }{ \sum_j RPKM_j } \times 10^6$$
TPM vs FPKM的关键区别：
- FPKM先标准化到总reads数，再除以基因长度
- TPM先除以基因长度得到RPKM-like值，再标准化到总和为10⁶
- TPM的优势：所有样本的TPM总和相同（10⁶），更适合样本间比较
- FPKM的问题：样本间总FPKM不同，受高表达基因影响大
（3）基于转录本估计的定量（Salmon, Kallisto, RSEM）
不依赖基因组比对，直接从reads推断转录本丰度：
- 使用伪比对（pseudoalignment）或EM算法
- 速度极快（Kallisto比传统方法快约100倍）
- 可直接输出TPM值

（1）为什么需要专门的统计方法
RNA-seq count数据具有以下统计特征：
- 非负整数（计数数据）
- 均值与方差相关（高表达基因的方差也大）
- 存在过度离散（overdispersion，方差>均值）
- 很多基因在部分样本中count为0（零膨胀）
因此不能使用普通t检验或ANOVA，需要专门的方法。
（2）DESeq2
基于负二项分布（Negative Binomial）模型：
1. 估计大小因子（Size Factor）：校正测序深度差异
$$size_factor_i = median_j \left( \frac{K_{ij}}{\left( \prod_{v=1}^{m} K_{vj} \right)^{1/m}} \right)$$
2. 估计离散度（Dispersion）：每个基因的变异程度
- 使用经验贝叶斯收缩（shrinkage）稳定低表达基因的离散度估计
3. 拟合广义线性模型（GLM）：
$$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$$
$$\log_2(\mu_{ij}) = \beta_0 + \beta_1 \cdot condition_j$$
4. Wald检验：检验系数β₁是否显著不为0
5. 多重检验校正：使用Benjamini-Hochberg方法控制FDR
（3）edgeR
同样基于负二项分布，但使用不同的离散度估计策略：
- 使用经验贝叶斯（empirical Bayes）从所有基因中借用信息
- 对小样本量（n<5）有较好的表现
差异表达基因的筛选标准：
- |log₂FC| ≥ 1（表达量变化≥2倍）
- padj < 0.05（校正后p值）

（1）火山图（Volcano Plot）
- X轴：log₂ Fold Change
- Y轴：-log₁₀(p-value)
- 每个点代表一个基因
- 同时展示效应大小（FC）和统计显著性（p-value）
- 右上和左上的点分别代表显著上调和下调的基因
（2）MA图（MA Plot）
- X轴：log₂(mean expression) —— 平均表达量
- Y轴：log₂ Fold Change
- 用于检查数据标准化效果：好的标准化应该使低表达和高表达基因的FC分布对称
（3）热图（Heatmap）
- 展示差异表达基因在不同样本中的表达模式
- 通常结合层次聚类，识别共表达基因模块
（4）PCA图
- 展示样本间的整体相似性
- 用于检测异常样本和批次效应

（1）ORA（Over-Representation Analysis）
- 输入：差异表达基因列表
- 方法：Fisher精确检验或超几何检验
- 输出：显著富集的GO term或KEGG pathway
- 工具：clusterProfiler, DAVID, g:Profiler
（2）GSEA（Gene Set Enrichment Analysis）
- 输入：所有基因的排序列表（如按log₂FC排序）
- 特点：不需要预设阈值，考虑所有基因
- 原理：检验预定义基因集（如某个通路）中的基因是否集中在排序列表的顶端或底端
- 输出：ES（Enrichment Score）和NES（Normalized Enrichment Score）
（3）GO语义相似性
- GO term之间存在层次关系
- 两个基因的功能相似性可以用其GO term在DAG中的距离来衡量

（1）动物miRNA的生物合成路径
1. 初级转录物（pri-miRNA）：由RNA聚合酶II转录，带有帽子结构和poly(A)尾，长度可达数千nt
2. 核内加工：Drosha（RNase III型核酸酶）在DGCR8辅助下，将pri-miRNA切割为~60-70 nt的pre-miRNA（发夹结构）
3. 核质运输：Exportin-5介导pre-miRNA从细胞核运输到细胞质
4. 成熟加工：Dicer酶（RNase III型）进一步切割pre-miRNA，形成~22 bp的miRNA:miRNA双链
5. RISC装载：AGO（Argonaute）蛋白结合双链，保留成熟miRNA链（通常5'端热力学稳定性较低的一条），降解另一条（miRNA）
（2）miRNA的调控机制
miRNA通过碱基配对识别靶mRNA：
- 种子区域（seed region）：miRNA的5'端第2-8位核苷酸，是靶识别的关键区域
- 完全互补（植物中常见）：导致靶mRNA的切割降解
- 不完全互补（动物中常见）：抑制翻译或促进mRNA降解（去腺苷酸化）
一个miRNA可以调控数百个靶基因，因此miRNA构成了复杂的基因调控网络。
（3）miRNA的计算预测
- 靶基因预测工具：TargetScan（基于种子互补性和保守性）、miRanda（基于自由能）、PITA（考虑靶位点可结合性）
- 新miRNA预测：基于发夹结构特征（MFEI > 0.85）和DCL-1加工依赖性

lncRNA的主要特征：
- 组织特异性和发育阶段特异性表达
- 表达量通常较低
- 物种间保守性相对较差（比蛋白质编码基因低）
- 可在表观遗传、转录及转录后等多层面调控基因表达
lncRNA的功能机制：
1. 顺式调控邻近基因：如XIST介导X染色体失活
2. 反式调控远端基因：如HOTAIR招募PRC2复合物到多个基因组位点
3. 作为miRNA海绵：lncRNA含有多个miRNA结合位点，竞争性吸附miRNA
4. 与蛋白质互作：作为支架RNA招募蛋白质复合物
5. 编码微肽：部分lncRNA含有短ORF，可翻译功能性微肽
lncRNA的鉴定：
1. 链特异性RNA-seq获取转录组数据
2. 转录本拼接（Cufflinks/StringTie）
3. 过滤已知注释基因和短转录本（<200 nt）
4. 编码潜能筛选（CPC2, CPAT, PLEK）
5. 排除具有显著蛋白质编码潜能的转录本

circRNA的生物合成：
1. 剪接体依赖的环化：下游外显子的5'端供体位点与上游的3'端受体位点结合
2. 内含子配对驱动：反向互补基序（如Alu元件）促进环化
3. 套索驱动：外显子跳过剪接产生的套索状副产物进一步环化
circRNA的独特性质：
- 高稳定性：无5'/3'端，对核酸外切酶免疫，半衰期可达数天（mRNA通常数小时）
- 组织特异性：具有较强的时间和组织表达特性
- 进化保守性：不同物种间具有一定保守性
circRNA的功能：
1. miRNA海绵：CDR1as含有63个miR-7结合位点
2. 调控亲本基因：与RNA结合蛋白互作，调节亲本基因表达
3. 翻译蛋白质：部分circRNA具有IRES和ORF，可翻译功能性蛋白
circRNA的计算鉴定：
- 主流方法：筛选不能线性比对的reads，检测反向剪接位点
- 工具：find_circ, CIRCexplorer, CIRI

早期方法（2011-2015）：
- Smart-seq/Smart-seq2：使用oligo(dT)引物和模板切换（template switching）逆转录，产生全长cDNA
- 优点：覆盖全转录本，可检测可变剪接
- 缺点：通量低（每次实验数十到数百个细胞），成本高
基于液滴的高通量方法（2015至今）：
- 10x Genomics Chromium：利用微流控技术将单个细胞与带有条形码的微珠包裹在油包水液滴中
- 优点：高通量（每次实验可测数千至数万个细胞），成本大幅下降
- 缺点：仅检测3'端（或5'端）序列，无法覆盖全转录本
- Drop-seq/InDrop：类似10x的原理，但由学术实验室开发

（1）单细胞悬液制备
- 组织解离为单细胞（机械+酶解法）
- 细胞质控：活性>90%，结团率<5%
- 细胞浓度调整至约1000 cells/μL
（2）GEM生成（Gel Bead-in-Emulsion）
使用微流控芯片同时导入三个流：
1. 单细胞悬液
2. Gel Bead（带有条形码寡核苷酸的凝胶珠）
3. 油相
在微流控通道中形成油包水液滴（GEM）：
- 每个液滴约含1个细胞 + 1个Gel Bead（泊松分布，实际捕获率约50-60%）
- Gel Bead上的寡核苷酸结构：
5'- [PCR引物] - [16bp细胞条形码] - [12bp UMI] - [30bp oligo(dT)] - 3'
（3）细胞内逆转录
- 液滴内，细胞裂解释放mRNA
- mRNA与Gel Bead上的oligo(dT)杂交
- 逆转录酶合成cDNA第一链，同时掺入细胞条形码和UMI
（4）破乳和cDNA扩增
- 破乳回收cDNA
- PCR扩增cDNA
（5）文库构建和测序
- 片段化、加接头、PCR扩增
- Illumina测序（通常2×150 bp）
（6）数据分析
- 根据reads中的细胞条形码将reads分配到不同细胞
- 根据UMI去重（同一个mRNA分子的多个PCR拷贝只计数一次）
- 生成基因×细胞的表达矩阵

Seurat是scRNA-seq分析最广泛使用的R包，由Satija Lab开发。
标准分析流程：

原始数据 (Cell Ranger输出)
    ↓
创建Seurat对象
    ↓
质控过滤 (nFeature, nCount, percent.mt)
    ↓
标准化 (LogNormalize, SCTransform)
    ↓
高变基因选择 (FindVariableFeatures)
    ↓
缩放 (ScaleData)
    ↓
线性降维 (PCA)
    ↓
非线性降维 (UMAP/t-SNE)
    ↓
聚类 (FindNeighbors + FindClusters, Louvain算法)
    ↓
细胞类型注释 (marker基因/自动注释)
    ↓
差异表达分析 (FindAllMarkers)
    ↓
可视化 (DimPlot, FeaturePlot, VlnPlot)

（1）质控
- nFeature_RNA：每个细胞检测到的基因数。过低（<200）可能是死细胞/空液滴；过高（>5000-8000）可能是双细胞
- nCount_RNA：每个细胞的总UMI数。与nFeature相关
- percent.mt：线粒体基因比例。过高（>10-20%）提示细胞膜破损，细胞质RNA流失
（2）标准化
- LogNormalize：对每个细胞的总UMI数进行标准化，然后取log1p
$$norm_count = \log\left( \frac{count}{total_UMI} \times 10000 + 1 \right)$$
- SCTransform：使用正则化负二项回归进行标准化，同时去除技术变异（如测序深度影响），效果通常优于简单标准化
（3）降维与聚类
- PCA：线性降维，保留主要变异方向
- UMAP（Uniform Manifold Approximation and Projection）：非线性降维，保留局部结构的同时拉开不同聚类
- Louvain聚类：基于共享最近邻（SNN）图的社区检测算法
（4）细胞类型注释
- 根据已知marker基因的表达进行手动注释
- 或使用自动注释工具（SingleR, scCATCH, celldex）
- 常用marker：
- T细胞：CD3D, CD3E
- B细胞：CD79A, MS4A1
- 单核/巨噬细胞：CD14, LYZ
- NK细胞：NKG7, GNLY

当多个样本或批次的数据需要合并分析时，批次效应是一个严重问题。
整合方法：
- CCA（Canonical Correlation Analysis）：Seurat v3的整合方法，找到跨批次共享的相关基因模块
- Harmony：快速、灵活的整合方法，适用于大规模数据
- scVI：基于变分自编码器的深度学习方法
- ComBat/sva：统计学校正方法

拟时序分析推断细胞在生物学过程中的发展轨迹：
- Monocle：使用反向图嵌入（RGE）构建发育轨迹
- Slingshot：基于降维后的聚类结果推断分支轨迹
- RNA velocity：利用未剪接和已剪接mRNA的比例推断细胞状态转换的方向

转录组学作为连接静态基因组和动态功能世界的桥梁，在过去十五年中经历了从基因芯片到RNA-seq再到单细胞测序的技术跃迁。每一次技术革新都带来了全新的生物学发现维度，深刻改变了我们对基因表达调控、细胞异质性和生命过程动态性的理解。

本章核心要点回顾

RNA-seq基础分析：从原始FASTQ数据到生物学洞见的完整流程包括质控、比对、定量、差异分析和功能富集。剪接比对（STAR/HISAT2）是处理真核生物数据的关键，TPM是推荐的标准化表达量指标，DESeq2基于负二项分布模型进行差异分析，而GSEA等富集分析则将基因列表转化为可解释的生物学发现。
非编码RNA：从曾被视为"转录噪音"到如今被确认为基因调控网络的核心参与者，非编码RNA的研究彻底改变了我们对基因组功能的认识。miRNA通过种子区域调控数百个靶基因，lncRNA通过多种机制在多个层次参与基因调控，circRNA则以其独特的稳定性开辟了新的调控维度。
单细胞转录组学：scRNA-seq技术使我们首次能够在单细胞分辨率上解析转录组异质性。10x Genomics的液滴微流控平台实现了高通量单细胞捕获，Seurat分析流程提供了从质控到注释的系统化方法。拟时序分析、RNA velocity等计算工具进一步使我们能够推断细胞状态转换的动态过程。

未来展望

空间转录组学：将转录组信息与组织空间位置结合，揭示基因表达的空间模式。Visium、Stereo-seq等技术正在将转录组学从"单细胞"推向"空间单细胞"时代。
多模态单细胞技术：同时测量单个细胞的转录组、表观基因组（scATAC-seq）、蛋白质组（CITE-seq）和免疫组库（scTCR/BCR-seq），实现多维度细胞状态的系统解析。
长读长RNA-seq：PacBio Iso-Seq和Nanopore direct RNA-seq能够产生全长转录本序列，彻底解决了短读长无法精确确定异构体的问题，将推动可变剪接和融合基因的精准研究。
单细胞多组学整合分析：结合基因组、转录组、表观遗传组和空间信息的整合分析方法，将使我们能够构建从基因型到表型的完整调控图谱。
临床转化：转录组生物标志物正在从研究走向临床。基于血液转录组的疾病诊断、基于单细胞图谱的肿瘤免疫分型、基于lncRNA/circRNA的新型标志物，都有望在精准医学中发挥重要作用。
转录组学的发展轨迹清晰地展示了技术创新如何驱动科学发现。从测量细胞群体的平均表达，到解析单细胞的异质性，再到定位空间位置信息，我们对生命活动动态过程的理解正在不断深化。作为生命科学的" read-out "技术，转录组学将继续在基础研究和临床转化中发挥核心作用。

本章思考题

1. 比较bulk RNA-seq和scRNA-seq在数据特征和分析方法上的主要差异。为什么scRNA-seq需要专门的分析工具？
2. 如果你要研究一种新型lncRNA在癌症中的功能，请设计一个结合RNA-seq、RIP-seq和CRISPR的实验方案。
3. 讨论单细胞转录组学面临的主要技术限制（如dropout、双细胞、批次效应），以及当前的研究如何尝试克服这些限制。
4. 空间转录组学如何改变了我们对组织生物学的理解？与传统scRNA-seq相比，空间信息带来了哪些额外的生物学洞见？
5. 展望10年后的转录组学，你认为哪些技术或方法会成为主流？转录组学将如何融入精准医学的临床实践？

推荐阅读
1. Mortazavi A, Williams B A, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq [J]. Nature methods, 2008, 5: 621-628. （RNA-seq奠基之作）
2. Trapnell C, Roberts A, Goff L, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks [J]. Nature protocols, 2012, 7: 562-578.
3. Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [J]. Genome biology, 2014, 15: 550.
4. Satija R, Farrell J A, Gennert D, et al. Spatial reconstruction of single-cell gene expression data [J]. Nature biotechnology, 2015, 33: 495-502. （Seurat经典论文）
5. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity [J]. Nature reviews immunology, 2018, 18: 35-45.

> 本章学习目标：掌握转录组学的核心概念、RNA-seq技术的完整分析流程、非编码RNA的分类与功能、单细胞转录组学的原理与方法，能够独立设计和执行转录组学研究项目。

中心法则与转录调控：
DNA →（转录）→ RNA →（翻译）→ 蛋白质
转录是基因表达的第一步，也是最关键的调控节点之一。生物体通过精确调控转录过程来实现：
- 时空特异性表达：不同组织、不同发育阶段表达不同的基因
- 环境响应：对外界刺激（温度、光照、营养、病原体等）做出快速反应
- 稳态维持：通过反馈调控维持细胞内环境稳定
转录调控的多层次机制：
1. 转录水平：转录因子结合、染色质可及性、DNA甲基化等
2. 转录后水平：mRNA剪接、编辑、稳定性调控
3. 翻译水平：翻译效率、miRNA调控
4. 翻译后水平：蛋白质修饰、定位、降解
转录组学主要聚焦于转录水平和转录后水平的研究。

第一代转录组技术：基因芯片（Microarray）
- 原理：将已知基因的寡核苷酸探针固定在芯片上，与标记的cDNA/cRNA杂交，通过荧光信号强度定量基因表达
- 局限：
- 只能检测已知序列（探针设计依赖已有基因组注释）
- 动态范围有限（通常2-3个数量级）
- 存在交叉杂交问题
- 难以检测新转录本和可变剪接
第二代转录组技术：RNA-seq
- 2008年首次报道（Mortazavi et al., Nature Methods）
- 原理：将RNA逆转录为cDNA，构建测序文库，进行高通量测序
- 优势：
- 不需要预先知道基因组序列
- 动态范围广（>5个数量级）
- 可同时检测已知和新转录本
- 可检测可变剪接、融合基因等
- 可检测SNP和等位基因特异性表达
第三代转录组技术：长读长RNA-seq
- PacBio Iso-Seq：产生全长转录本序列，无需组装即可确定完整异构体
- Nanopore direct RNA-seq：直接测序RNA分子，保留碱基修饰信息

（1）RNA提取与质控
- 使用Trizol或柱式提取法从样本中提取总RNA
- RNA完整性评估：Agilent Bioanalyzer检测RIN值（RNA Integrity Number）
- RIN > 7通常被认为是良好质量的RNA；RIN < 5可能导致3'端偏倚
（2）RNA富集或rRNA去除
- mRNA富集：使用oligo(dT)磁珠捕获带poly(A)尾的mRNA
- 优点：简单高效
- 缺点：会丢失不带poly(A)尾的RNA（如部分lncRNA、原核生物mRNA）
- rRNA去除：使用rRNA探针杂交去除rRNA（Ribo-Zero, Ribo-off等）
- 优点：保留所有非rRNA（包括mRNA、lncRNA、circRNA等）
- 缺点：操作稍复杂，成本略高
（3）RNA片段化与cDNA合成
- 片段化：使用二价阳离子（如Mg²⁺）在加热条件下将RNA随机打断成200-300 nt片段
- cDNA合成：
- 第一条链：使用随机六聚体或oligo(dT)引物，逆转录酶合成
- 第二条链：使用DNA聚合酶合成互补链（通常掺入dUTP以标记第二条链，用于链特异性文库）
（4）文库构建
- 末端修复：将cDNA片段末端补平
- A尾添加：在3'端添加单个A碱基
- 接头连接：连接带有barcode的测序接头
- PCR扩增：扩增文库至足够浓度
- 文库质控：检测文库大小分布（Agilent Bioanalyzer）和浓度（Qubit/qPCR）
（5）测序策略
- 单端测序（Single-end）：只测cDNA片段的一端，成本较低
- 双端测序（Paired-end）：测序两端，可更好地定位reads和识别可变剪接
- 链特异性测序（Strand-specific）：保留RNA链方向信息，可鉴定反义转录本

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, MultiQC)
    ↓
去除低质量和接头序列 (Trimmomatic, cutadapt)
    ↓
比对到参考基因组 (STAR, HISAT2)
    ↓
比对后处理 (排序, 索引)
    ↓
表达量定量 (featureCounts, HTSeq, Salmon)
    ↓
原始count矩阵
    ↓
数据探索 (PCA, 样本相关性热图)
    ↓
差异表达分析 (DESeq2, edgeR, limma)
    ↓
差异表达基因列表
    ↓
功能富集分析 (clusterProfiler, DAVID, GSEA)
    ↓
生物学解释和可视化

RNA-seq reads可能跨越一个或多个剪接位点（splice junction），这要求比对软件能够处理"分裂比对"（split alignment）。
STAR比对算法：
1. 种子生成：将read切分为多个短片段（seeds），在基因组上搜索精确匹配
2. 聚类：将基因组上接近的种子聚类为候选比对区域
3. 拼接：对于跨越内含子的reads，在候选区域之间构建剪接位点
4. 评分：对完整read进行比对评分，选择最优比对位置
STAR的特点：
- 速度快（比TopHat2快约50倍）
- 准确性高
- 内存需求大（需要约30-40 GB RAM用于人类基因组）
HISAT2比对算法：
- 使用基于图的比对方法（graph-based alignment）
- 利用FM-index和局部哈希实现快速搜索
- 对剪接变异的检测更敏感
- 内存效率高（<10 GB for human genome）

（1）基于count的方法
统计比对到每个基因区域的reads数。
- featureCounts：快速、准确，广泛使用的计数工具
- HTSeq-count：早期标准工具，灵活性高
（2）基于比对的定量（RPKM/FPKM/TPM）
原始count需要标准化才能进行样本内和样本间比较：
- RPKM（Reads Per Kilobase of exon per Million mapped reads）：
$$RPKM = \frac{ reads \times 10^9 }{ gene_length \times total_mapped_reads }$$
- FPKM（Fragments Per Kilobase of exon per Million mapped fragments）：
与RPKM类似，但使用fragment（pair-end测序中的一对reads）代替reads
- TPM（Transcripts Per Million）：
$$TPM_i = \frac{ RPKM_i }{ \sum_j RPKM_j } \times 10^6$$
TPM vs FPKM的关键区别：
- FPKM先标准化到总reads数，再除以基因长度
- TPM先除以基因长度得到RPKM-like值，再标准化到总和为10⁶
- TPM的优势：所有样本的TPM总和相同（10⁶），更适合样本间比较
- FPKM的问题：样本间总FPKM不同，受高表达基因影响大
（3）基于转录本估计的定量（Salmon, Kallisto, RSEM）
不依赖基因组比对，直接从reads推断转录本丰度：
- 使用伪比对（pseudoalignment）或EM算法
- 速度极快（Kallisto比传统方法快约100倍）
- 可直接输出TPM值

（1）为什么需要专门的统计方法
RNA-seq count数据具有以下统计特征：
- 非负整数（计数数据）
- 均值与方差相关（高表达基因的方差也大）
- 存在过度离散（overdispersion，方差>均值）
- 很多基因在部分样本中count为0（零膨胀）
因此不能使用普通t检验或ANOVA，需要专门的方法。
（2）DESeq2
基于负二项分布（Negative Binomial）模型：
1. 估计大小因子（Size Factor）：校正测序深度差异
$$size_factor_i = median_j \left( \frac{K_{ij}}{\left( \prod_{v=1}^{m} K_{vj} \right)^{1/m}} \right)$$
2. 估计离散度（Dispersion）：每个基因的变异程度
- 使用经验贝叶斯收缩（shrinkage）稳定低表达基因的离散度估计
3. 拟合广义线性模型（GLM）：
$$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$$
$$\log_2(\mu_{ij}) = \beta_0 + \beta_1 \cdot condition_j$$
4. Wald检验：检验系数β₁是否显著不为0
5. 多重检验校正：使用Benjamini-Hochberg方法控制FDR
（3）edgeR
同样基于负二项分布，但使用不同的离散度估计策略：
- 使用经验贝叶斯（empirical Bayes）从所有基因中借用信息
- 对小样本量（n<5）有较好的表现
差异表达基因的筛选标准：
- |log₂FC| ≥ 1（表达量变化≥2倍）
- padj < 0.05（校正后p值）

（1）火山图（Volcano Plot）
- X轴：log₂ Fold Change
- Y轴：-log₁₀(p-value)
- 每个点代表一个基因
- 同时展示效应大小（FC）和统计显著性（p-value）
- 右上和左上的点分别代表显著上调和下调的基因
（2）MA图（MA Plot）
- X轴：log₂(mean expression) —— 平均表达量
- Y轴：log₂ Fold Change
- 用于检查数据标准化效果：好的标准化应该使低表达和高表达基因的FC分布对称
（3）热图（Heatmap）
- 展示差异表达基因在不同样本中的表达模式
- 通常结合层次聚类，识别共表达基因模块
（4）PCA图
- 展示样本间的整体相似性
- 用于检测异常样本和批次效应

（1）ORA（Over-Representation Analysis）
- 输入：差异表达基因列表
- 方法：Fisher精确检验或超几何检验
- 输出：显著富集的GO term或KEGG pathway
- 工具：clusterProfiler, DAVID, g:Profiler
（2）GSEA（Gene Set Enrichment Analysis）
- 输入：所有基因的排序列表（如按log₂FC排序）
- 特点：不需要预设阈值，考虑所有基因
- 原理：检验预定义基因集（如某个通路）中的基因是否集中在排序列表的顶端或底端
- 输出：ES（Enrichment Score）和NES（Normalized Enrichment Score）
（3）GO语义相似性
- GO term之间存在层次关系
- 两个基因的功能相似性可以用其GO term在DAG中的距离来衡量

（1）动物miRNA的生物合成路径
1. 初级转录物（pri-miRNA）：由RNA聚合酶II转录，带有帽子结构和poly(A)尾，长度可达数千nt
2. 核内加工：Drosha（RNase III型核酸酶）在DGCR8辅助下，将pri-miRNA切割为~60-70 nt的pre-miRNA（发夹结构）
3. 核质运输：Exportin-5介导pre-miRNA从细胞核运输到细胞质
4. 成熟加工：Dicer酶（RNase III型）进一步切割pre-miRNA，形成~22 bp的miRNA:miRNA双链
5. RISC装载：AGO（Argonaute）蛋白结合双链，保留成熟miRNA链（通常5'端热力学稳定性较低的一条），降解另一条（miRNA）
（2）miRNA的调控机制
miRNA通过碱基配对识别靶mRNA：
- 种子区域（seed region）：miRNA的5'端第2-8位核苷酸，是靶识别的关键区域
- 完全互补（植物中常见）：导致靶mRNA的切割降解
- 不完全互补（动物中常见）：抑制翻译或促进mRNA降解（去腺苷酸化）
一个miRNA可以调控数百个靶基因，因此miRNA构成了复杂的基因调控网络。
（3）miRNA的计算预测
- 靶基因预测工具：TargetScan（基于种子互补性和保守性）、miRanda（基于自由能）、PITA（考虑靶位点可结合性）
- 新miRNA预测：基于发夹结构特征（MFEI > 0.85）和DCL-1加工依赖性

lncRNA的主要特征：
- 组织特异性和发育阶段特异性表达
- 表达量通常较低
- 物种间保守性相对较差（比蛋白质编码基因低）
- 可在表观遗传、转录及转录后等多层面调控基因表达
lncRNA的功能机制：
1. 顺式调控邻近基因：如XIST介导X染色体失活
2. 反式调控远端基因：如HOTAIR招募PRC2复合物到多个基因组位点
3. 作为miRNA海绵：lncRNA含有多个miRNA结合位点，竞争性吸附miRNA
4. 与蛋白质互作：作为支架RNA招募蛋白质复合物
5. 编码微肽：部分lncRNA含有短ORF，可翻译功能性微肽
lncRNA的鉴定：
1. 链特异性RNA-seq获取转录组数据
2. 转录本拼接（Cufflinks/StringTie）
3. 过滤已知注释基因和短转录本（<200 nt）
4. 编码潜能筛选（CPC2, CPAT, PLEK）
5. 排除具有显著蛋白质编码潜能的转录本

circRNA的生物合成：
1. 剪接体依赖的环化：下游外显子的5'端供体位点与上游的3'端受体位点结合
2. 内含子配对驱动：反向互补基序（如Alu元件）促进环化
3. 套索驱动：外显子跳过剪接产生的套索状副产物进一步环化
circRNA的独特性质：
- 高稳定性：无5'/3'端，对核酸外切酶免疫，半衰期可达数天（mRNA通常数小时）
- 组织特异性：具有较强的时间和组织表达特性
- 进化保守性：不同物种间具有一定保守性
circRNA的功能：
1. miRNA海绵：CDR1as含有63个miR-7结合位点
2. 调控亲本基因：与RNA结合蛋白互作，调节亲本基因表达
3. 翻译蛋白质：部分circRNA具有IRES和ORF，可翻译功能性蛋白
circRNA的计算鉴定：
- 主流方法：筛选不能线性比对的reads，检测反向剪接位点
- 工具：find_circ, CIRCexplorer, CIRI

早期方法（2011-2015）：
- Smart-seq/Smart-seq2：使用oligo(dT)引物和模板切换（template switching）逆转录，产生全长cDNA
- 优点：覆盖全转录本，可检测可变剪接
- 缺点：通量低（每次实验数十到数百个细胞），成本高
基于液滴的高通量方法（2015至今）：
- 10x Genomics Chromium：利用微流控技术将单个细胞与带有条形码的微珠包裹在油包水液滴中
- 优点：高通量（每次实验可测数千至数万个细胞），成本大幅下降
- 缺点：仅检测3'端（或5'端）序列，无法覆盖全转录本
- Drop-seq/InDrop：类似10x的原理，但由学术实验室开发

（1）单细胞悬液制备
- 组织解离为单细胞（机械+酶解法）
- 细胞质控：活性>90%，结团率<5%
- 细胞浓度调整至约1000 cells/μL
（2）GEM生成（Gel Bead-in-Emulsion）
使用微流控芯片同时导入三个流：
1. 单细胞悬液
2. Gel Bead（带有条形码寡核苷酸的凝胶珠）
3. 油相
在微流控通道中形成油包水液滴（GEM）：
- 每个液滴约含1个细胞 + 1个Gel Bead（泊松分布，实际捕获率约50-60%）
- Gel Bead上的寡核苷酸结构：
5'- [PCR引物] - [16bp细胞条形码] - [12bp UMI] - [30bp oligo(dT)] - 3'
（3）细胞内逆转录
- 液滴内，细胞裂解释放mRNA
- mRNA与Gel Bead上的oligo(dT)杂交
- 逆转录酶合成cDNA第一链，同时掺入细胞条形码和UMI
（4）破乳和cDNA扩增
- 破乳回收cDNA
- PCR扩增cDNA
（5）文库构建和测序
- 片段化、加接头、PCR扩增
- Illumina测序（通常2×150 bp）
（6）数据分析
- 根据reads中的细胞条形码将reads分配到不同细胞
- 根据UMI去重（同一个mRNA分子的多个PCR拷贝只计数一次）
- 生成基因×细胞的表达矩阵

Seurat是scRNA-seq分析最广泛使用的R包，由Satija Lab开发。
标准分析流程：

原始数据 (Cell Ranger输出)
    ↓
创建Seurat对象
    ↓
质控过滤 (nFeature, nCount, percent.mt)
    ↓
标准化 (LogNormalize, SCTransform)
    ↓
高变基因选择 (FindVariableFeatures)
    ↓
缩放 (ScaleData)
    ↓
线性降维 (PCA)
    ↓
非线性降维 (UMAP/t-SNE)
    ↓
聚类 (FindNeighbors + FindClusters, Louvain算法)
    ↓
细胞类型注释 (marker基因/自动注释)
    ↓
差异表达分析 (FindAllMarkers)
    ↓
可视化 (DimPlot, FeaturePlot, VlnPlot)

（1）质控
- nFeature_RNA：每个细胞检测到的基因数。过低（<200）可能是死细胞/空液滴；过高（>5000-8000）可能是双细胞
- nCount_RNA：每个细胞的总UMI数。与nFeature相关
- percent.mt：线粒体基因比例。过高（>10-20%）提示细胞膜破损，细胞质RNA流失
（2）标准化
- LogNormalize：对每个细胞的总UMI数进行标准化，然后取log1p
$$norm_count = \log\left( \frac{count}{total_UMI} \times 10000 + 1 \right)$$
- SCTransform：使用正则化负二项回归进行标准化，同时去除技术变异（如测序深度影响），效果通常优于简单标准化
（3）降维与聚类
- PCA：线性降维，保留主要变异方向
- UMAP（Uniform Manifold Approximation and Projection）：非线性降维，保留局部结构的同时拉开不同聚类
- Louvain聚类：基于共享最近邻（SNN）图的社区检测算法
（4）细胞类型注释
- 根据已知marker基因的表达进行手动注释
- 或使用自动注释工具（SingleR, scCATCH, celldex）
- 常用marker：
- T细胞：CD3D, CD3E
- B细胞：CD79A, MS4A1
- 单核/巨噬细胞：CD14, LYZ
- NK细胞：NKG7, GNLY

当多个样本或批次的数据需要合并分析时，批次效应是一个严重问题。
整合方法：
- CCA（Canonical Correlation Analysis）：Seurat v3的整合方法，找到跨批次共享的相关基因模块
- Harmony：快速、灵活的整合方法，适用于大规模数据
- scVI：基于变分自编码器的深度学习方法
- ComBat/sva：统计学校正方法

拟时序分析推断细胞在生物学过程中的发展轨迹：
- Monocle：使用反向图嵌入（RGE）构建发育轨迹
- Slingshot：基于降维后的聚类结果推断分支轨迹
- RNA velocity：利用未剪接和已剪接mRNA的比例推断细胞状态转换的方向

转录组学作为连接静态基因组和动态功能世界的桥梁，在过去十五年中经历了从基因芯片到RNA-seq再到单细胞测序的技术跃迁。每一次技术革新都带来了全新的生物学发现维度，深刻改变了我们对基因表达调控、细胞异质性和生命过程动态性的理解。

本章核心要点回顾

RNA-seq基础分析：从原始FASTQ数据到生物学洞见的完整流程包括质控、比对、定量、差异分析和功能富集。剪接比对（STAR/HISAT2）是处理真核生物数据的关键，TPM是推荐的标准化表达量指标，DESeq2基于负二项分布模型进行差异分析，而GSEA等富集分析则将基因列表转化为可解释的生物学发现。
非编码RNA：从曾被视为"转录噪音"到如今被确认为基因调控网络的核心参与者，非编码RNA的研究彻底改变了我们对基因组功能的认识。miRNA通过种子区域调控数百个靶基因，lncRNA通过多种机制在多个层次参与基因调控，circRNA则以其独特的稳定性开辟了新的调控维度。
单细胞转录组学：scRNA-seq技术使我们首次能够在单细胞分辨率上解析转录组异质性。10x Genomics的液滴微流控平台实现了高通量单细胞捕获，Seurat分析流程提供了从质控到注释的系统化方法。拟时序分析、RNA velocity等计算工具进一步使我们能够推断细胞状态转换的动态过程。

未来展望

空间转录组学：将转录组信息与组织空间位置结合，揭示基因表达的空间模式。Visium、Stereo-seq等技术正在将转录组学从"单细胞"推向"空间单细胞"时代。
多模态单细胞技术：同时测量单个细胞的转录组、表观基因组（scATAC-seq）、蛋白质组（CITE-seq）和免疫组库（scTCR/BCR-seq），实现多维度细胞状态的系统解析。
长读长RNA-seq：PacBio Iso-Seq和Nanopore direct RNA-seq能够产生全长转录本序列，彻底解决了短读长无法精确确定异构体的问题，将推动可变剪接和融合基因的精准研究。
单细胞多组学整合分析：结合基因组、转录组、表观遗传组和空间信息的整合分析方法，将使我们能够构建从基因型到表型的完整调控图谱。
临床转化：转录组生物标志物正在从研究走向临床。基于血液转录组的疾病诊断、基于单细胞图谱的肿瘤免疫分型、基于lncRNA/circRNA的新型标志物，都有望在精准医学中发挥重要作用。
转录组学的发展轨迹清晰地展示了技术创新如何驱动科学发现。从测量细胞群体的平均表达，到解析单细胞的异质性，再到定位空间位置信息，我们对生命活动动态过程的理解正在不断深化。作为生命科学的" read-out "技术，转录组学将继续在基础研究和临床转化中发挥核心作用。

本章思考题

1. 比较bulk RNA-seq和scRNA-seq在数据特征和分析方法上的主要差异。为什么scRNA-seq需要专门的分析工具？
2. 如果你要研究一种新型lncRNA在癌症中的功能，请设计一个结合RNA-seq、RIP-seq和CRISPR的实验方案。
3. 讨论单细胞转录组学面临的主要技术限制（如dropout、双细胞、批次效应），以及当前的研究如何尝试克服这些限制。
4. 空间转录组学如何改变了我们对组织生物学的理解？与传统scRNA-seq相比，空间信息带来了哪些额外的生物学洞见？
5. 展望10年后的转录组学，你认为哪些技术或方法会成为主流？转录组学将如何融入精准医学的临床实践？

推荐阅读
1. Mortazavi A, Williams B A, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq [J]. Nature methods, 2008, 5: 621-628. （RNA-seq奠基之作）
2. Trapnell C, Roberts A, Goff L, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks [J]. Nature protocols, 2012, 7: 562-578.
3. Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [J]. Genome biology, 2014, 15: 550.
4. Satija R, Farrell J A, Gennert D, et al. Spatial reconstruction of single-cell gene expression data [J]. Nature biotechnology, 2015, 33: 495-502. （Seurat经典论文）
5. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity [J]. Nature reviews immunology, 2018, 18: 35-45.

> 本章学习目标：掌握转录组学的核心概念、RNA-seq技术的完整分析流程、非编码RNA的分类与功能、单细胞转录组学的原理与方法，能够独立设计和执行转录组学研究项目。

中心法则与转录调控：
DNA →（转录）→ RNA →（翻译）→ 蛋白质
转录是基因表达的第一步，也是最关键的调控节点之一。生物体通过精确调控转录过程来实现：
- 时空特异性表达：不同组织、不同发育阶段表达不同的基因
- 环境响应：对外界刺激（温度、光照、营养、病原体等）做出快速反应
- 稳态维持：通过反馈调控维持细胞内环境稳定
转录调控的多层次机制：
1. 转录水平：转录因子结合、染色质可及性、DNA甲基化等
2. 转录后水平：mRNA剪接、编辑、稳定性调控
3. 翻译水平：翻译效率、miRNA调控
4. 翻译后水平：蛋白质修饰、定位、降解
转录组学主要聚焦于转录水平和转录后水平的研究。

第一代转录组技术：基因芯片（Microarray）
- 原理：将已知基因的寡核苷酸探针固定在芯片上，与标记的cDNA/cRNA杂交，通过荧光信号强度定量基因表达
- 局限：
- 只能检测已知序列（探针设计依赖已有基因组注释）
- 动态范围有限（通常2-3个数量级）
- 存在交叉杂交问题
- 难以检测新转录本和可变剪接
第二代转录组技术：RNA-seq
- 2008年首次报道（Mortazavi et al., Nature Methods）
- 原理：将RNA逆转录为cDNA，构建测序文库，进行高通量测序
- 优势：
- 不需要预先知道基因组序列
- 动态范围广（>5个数量级）
- 可同时检测已知和新转录本
- 可检测可变剪接、融合基因等
- 可检测SNP和等位基因特异性表达
第三代转录组技术：长读长RNA-seq
- PacBio Iso-Seq：产生全长转录本序列，无需组装即可确定完整异构体
- Nanopore direct RNA-seq：直接测序RNA分子，保留碱基修饰信息

（1）RNA提取与质控
- 使用Trizol或柱式提取法从样本中提取总RNA
- RNA完整性评估：Agilent Bioanalyzer检测RIN值（RNA Integrity Number）
- RIN > 7通常被认为是良好质量的RNA；RIN < 5可能导致3'端偏倚
（2）RNA富集或rRNA去除
- mRNA富集：使用oligo(dT)磁珠捕获带poly(A)尾的mRNA
- 优点：简单高效
- 缺点：会丢失不带poly(A)尾的RNA（如部分lncRNA、原核生物mRNA）
- rRNA去除：使用rRNA探针杂交去除rRNA（Ribo-Zero, Ribo-off等）
- 优点：保留所有非rRNA（包括mRNA、lncRNA、circRNA等）
- 缺点：操作稍复杂，成本略高
（3）RNA片段化与cDNA合成
- 片段化：使用二价阳离子（如Mg²⁺）在加热条件下将RNA随机打断成200-300 nt片段
- cDNA合成：
- 第一条链：使用随机六聚体或oligo(dT)引物，逆转录酶合成
- 第二条链：使用DNA聚合酶合成互补链（通常掺入dUTP以标记第二条链，用于链特异性文库）
（4）文库构建
- 末端修复：将cDNA片段末端补平
- A尾添加：在3'端添加单个A碱基
- 接头连接：连接带有barcode的测序接头
- PCR扩增：扩增文库至足够浓度
- 文库质控：检测文库大小分布（Agilent Bioanalyzer）和浓度（Qubit/qPCR）
（5）测序策略
- 单端测序（Single-end）：只测cDNA片段的一端，成本较低
- 双端测序（Paired-end）：测序两端，可更好地定位reads和识别可变剪接
- 链特异性测序（Strand-specific）：保留RNA链方向信息，可鉴定反义转录本

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, MultiQC)
    ↓
去除低质量和接头序列 (Trimmomatic, cutadapt)
    ↓
比对到参考基因组 (STAR, HISAT2)
    ↓
比对后处理 (排序, 索引)
    ↓
表达量定量 (featureCounts, HTSeq, Salmon)
    ↓
原始count矩阵
    ↓
数据探索 (PCA, 样本相关性热图)
    ↓
差异表达分析 (DESeq2, edgeR, limma)
    ↓
差异表达基因列表
    ↓
功能富集分析 (clusterProfiler, DAVID, GSEA)
    ↓
生物学解释和可视化

RNA-seq reads可能跨越一个或多个剪接位点（splice junction），这要求比对软件能够处理"分裂比对"（split alignment）。
STAR比对算法：
1. 种子生成：将read切分为多个短片段（seeds），在基因组上搜索精确匹配
2. 聚类：将基因组上接近的种子聚类为候选比对区域
3. 拼接：对于跨越内含子的reads，在候选区域之间构建剪接位点
4. 评分：对完整read进行比对评分，选择最优比对位置
STAR的特点：
- 速度快（比TopHat2快约50倍）
- 准确性高
- 内存需求大（需要约30-40 GB RAM用于人类基因组）
HISAT2比对算法：
- 使用基于图的比对方法（graph-based alignment）
- 利用FM-index和局部哈希实现快速搜索
- 对剪接变异的检测更敏感
- 内存效率高（<10 GB for human genome）

（1）基于count的方法
统计比对到每个基因区域的reads数。
- featureCounts：快速、准确，广泛使用的计数工具
- HTSeq-count：早期标准工具，灵活性高
（2）基于比对的定量（RPKM/FPKM/TPM）
原始count需要标准化才能进行样本内和样本间比较：
- RPKM（Reads Per Kilobase of exon per Million mapped reads）：
$$RPKM = \frac{ reads \times 10^9 }{ gene_length \times total_mapped_reads }$$
- FPKM（Fragments Per Kilobase of exon per Million mapped fragments）：
与RPKM类似，但使用fragment（pair-end测序中的一对reads）代替reads
- TPM（Transcripts Per Million）：
$$TPM_i = \frac{ RPKM_i }{ \sum_j RPKM_j } \times 10^6$$
TPM vs FPKM的关键区别：
- FPKM先标准化到总reads数，再除以基因长度
- TPM先除以基因长度得到RPKM-like值，再标准化到总和为10⁶
- TPM的优势：所有样本的TPM总和相同（10⁶），更适合样本间比较
- FPKM的问题：样本间总FPKM不同，受高表达基因影响大
（3）基于转录本估计的定量（Salmon, Kallisto, RSEM）
不依赖基因组比对，直接从reads推断转录本丰度：
- 使用伪比对（pseudoalignment）或EM算法
- 速度极快（Kallisto比传统方法快约100倍）
- 可直接输出TPM值

（1）为什么需要专门的统计方法
RNA-seq count数据具有以下统计特征：
- 非负整数（计数数据）
- 均值与方差相关（高表达基因的方差也大）
- 存在过度离散（overdispersion，方差>均值）
- 很多基因在部分样本中count为0（零膨胀）
因此不能使用普通t检验或ANOVA，需要专门的方法。
（2）DESeq2
基于负二项分布（Negative Binomial）模型：
1. 估计大小因子（Size Factor）：校正测序深度差异
$$size_factor_i = median_j \left( \frac{K_{ij}}{\left( \prod_{v=1}^{m} K_{vj} \right)^{1/m}} \right)$$
2. 估计离散度（Dispersion）：每个基因的变异程度
- 使用经验贝叶斯收缩（shrinkage）稳定低表达基因的离散度估计
3. 拟合广义线性模型（GLM）：
$$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$$
$$\log_2(\mu_{ij}) = \beta_0 + \beta_1 \cdot condition_j$$
4. Wald检验：检验系数β₁是否显著不为0
5. 多重检验校正：使用Benjamini-Hochberg方法控制FDR
（3）edgeR
同样基于负二项分布，但使用不同的离散度估计策略：
- 使用经验贝叶斯（empirical Bayes）从所有基因中借用信息
- 对小样本量（n<5）有较好的表现
差异表达基因的筛选标准：
- |log₂FC| ≥ 1（表达量变化≥2倍）
- padj < 0.05（校正后p值）

（1）火山图（Volcano Plot）
- X轴：log₂ Fold Change
- Y轴：-log₁₀(p-value)
- 每个点代表一个基因
- 同时展示效应大小（FC）和统计显著性（p-value）
- 右上和左上的点分别代表显著上调和下调的基因
（2）MA图（MA Plot）
- X轴：log₂(mean expression) —— 平均表达量
- Y轴：log₂ Fold Change
- 用于检查数据标准化效果：好的标准化应该使低表达和高表达基因的FC分布对称
（3）热图（Heatmap）
- 展示差异表达基因在不同样本中的表达模式
- 通常结合层次聚类，识别共表达基因模块
（4）PCA图
- 展示样本间的整体相似性
- 用于检测异常样本和批次效应

（1）ORA（Over-Representation Analysis）
- 输入：差异表达基因列表
- 方法：Fisher精确检验或超几何检验
- 输出：显著富集的GO term或KEGG pathway
- 工具：clusterProfiler, DAVID, g:Profiler
（2）GSEA（Gene Set Enrichment Analysis）
- 输入：所有基因的排序列表（如按log₂FC排序）
- 特点：不需要预设阈值，考虑所有基因
- 原理：检验预定义基因集（如某个通路）中的基因是否集中在排序列表的顶端或底端
- 输出：ES（Enrichment Score）和NES（Normalized Enrichment Score）
（3）GO语义相似性
- GO term之间存在层次关系
- 两个基因的功能相似性可以用其GO term在DAG中的距离来衡量

（1）动物miRNA的生物合成路径
1. 初级转录物（pri-miRNA）：由RNA聚合酶II转录，带有帽子结构和poly(A)尾，长度可达数千nt
2. 核内加工：Drosha（RNase III型核酸酶）在DGCR8辅助下，将pri-miRNA切割为~60-70 nt的pre-miRNA（发夹结构）
3. 核质运输：Exportin-5介导pre-miRNA从细胞核运输到细胞质
4. 成熟加工：Dicer酶（RNase III型）进一步切割pre-miRNA，形成~22 bp的miRNA:miRNA双链
5. RISC装载：AGO（Argonaute）蛋白结合双链，保留成熟miRNA链（通常5'端热力学稳定性较低的一条），降解另一条（miRNA）
（2）miRNA的调控机制
miRNA通过碱基配对识别靶mRNA：
- 种子区域（seed region）：miRNA的5'端第2-8位核苷酸，是靶识别的关键区域
- 完全互补（植物中常见）：导致靶mRNA的切割降解
- 不完全互补（动物中常见）：抑制翻译或促进mRNA降解（去腺苷酸化）
一个miRNA可以调控数百个靶基因，因此miRNA构成了复杂的基因调控网络。
（3）miRNA的计算预测
- 靶基因预测工具：TargetScan（基于种子互补性和保守性）、miRanda（基于自由能）、PITA（考虑靶位点可结合性）
- 新miRNA预测：基于发夹结构特征（MFEI > 0.85）和DCL-1加工依赖性

lncRNA的主要特征：
- 组织特异性和发育阶段特异性表达
- 表达量通常较低
- 物种间保守性相对较差（比蛋白质编码基因低）
- 可在表观遗传、转录及转录后等多层面调控基因表达
lncRNA的功能机制：
1. 顺式调控邻近基因：如XIST介导X染色体失活
2. 反式调控远端基因：如HOTAIR招募PRC2复合物到多个基因组位点
3. 作为miRNA海绵：lncRNA含有多个miRNA结合位点，竞争性吸附miRNA
4. 与蛋白质互作：作为支架RNA招募蛋白质复合物
5. 编码微肽：部分lncRNA含有短ORF，可翻译功能性微肽
lncRNA的鉴定：
1. 链特异性RNA-seq获取转录组数据
2. 转录本拼接（Cufflinks/StringTie）
3. 过滤已知注释基因和短转录本（<200 nt）
4. 编码潜能筛选（CPC2, CPAT, PLEK）
5. 排除具有显著蛋白质编码潜能的转录本

circRNA的生物合成：
1. 剪接体依赖的环化：下游外显子的5'端供体位点与上游的3'端受体位点结合
2. 内含子配对驱动：反向互补基序（如Alu元件）促进环化
3. 套索驱动：外显子跳过剪接产生的套索状副产物进一步环化
circRNA的独特性质：
- 高稳定性：无5'/3'端，对核酸外切酶免疫，半衰期可达数天（mRNA通常数小时）
- 组织特异性：具有较强的时间和组织表达特性
- 进化保守性：不同物种间具有一定保守性
circRNA的功能：
1. miRNA海绵：CDR1as含有63个miR-7结合位点
2. 调控亲本基因：与RNA结合蛋白互作，调节亲本基因表达
3. 翻译蛋白质：部分circRNA具有IRES和ORF，可翻译功能性蛋白
circRNA的计算鉴定：
- 主流方法：筛选不能线性比对的reads，检测反向剪接位点
- 工具：find_circ, CIRCexplorer, CIRI

早期方法（2011-2015）：
- Smart-seq/Smart-seq2：使用oligo(dT)引物和模板切换（template switching）逆转录，产生全长cDNA
- 优点：覆盖全转录本，可检测可变剪接
- 缺点：通量低（每次实验数十到数百个细胞），成本高
基于液滴的高通量方法（2015至今）：
- 10x Genomics Chromium：利用微流控技术将单个细胞与带有条形码的微珠包裹在油包水液滴中
- 优点：高通量（每次实验可测数千至数万个细胞），成本大幅下降
- 缺点：仅检测3'端（或5'端）序列，无法覆盖全转录本
- Drop-seq/InDrop：类似10x的原理，但由学术实验室开发

（1）单细胞悬液制备
- 组织解离为单细胞（机械+酶解法）
- 细胞质控：活性>90%，结团率<5%
- 细胞浓度调整至约1000 cells/μL
（2）GEM生成（Gel Bead-in-Emulsion）
使用微流控芯片同时导入三个流：
1. 单细胞悬液
2. Gel Bead（带有条形码寡核苷酸的凝胶珠）
3. 油相
在微流控通道中形成油包水液滴（GEM）：
- 每个液滴约含1个细胞 + 1个Gel Bead（泊松分布，实际捕获率约50-60%）
- Gel Bead上的寡核苷酸结构：
5'- [PCR引物] - [16bp细胞条形码] - [12bp UMI] - [30bp oligo(dT)] - 3'
（3）细胞内逆转录
- 液滴内，细胞裂解释放mRNA
- mRNA与Gel Bead上的oligo(dT)杂交
- 逆转录酶合成cDNA第一链，同时掺入细胞条形码和UMI
（4）破乳和cDNA扩增
- 破乳回收cDNA
- PCR扩增cDNA
（5）文库构建和测序
- 片段化、加接头、PCR扩增
- Illumina测序（通常2×150 bp）
（6）数据分析
- 根据reads中的细胞条形码将reads分配到不同细胞
- 根据UMI去重（同一个mRNA分子的多个PCR拷贝只计数一次）
- 生成基因×细胞的表达矩阵

Seurat是scRNA-seq分析最广泛使用的R包，由Satija Lab开发。
标准分析流程：

原始数据 (Cell Ranger输出)
    ↓
创建Seurat对象
    ↓
质控过滤 (nFeature, nCount, percent.mt)
    ↓
标准化 (LogNormalize, SCTransform)
    ↓
高变基因选择 (FindVariableFeatures)
    ↓
缩放 (ScaleData)
    ↓
线性降维 (PCA)
    ↓
非线性降维 (UMAP/t-SNE)
    ↓
聚类 (FindNeighbors + FindClusters, Louvain算法)
    ↓
细胞类型注释 (marker基因/自动注释)
    ↓
差异表达分析 (FindAllMarkers)
    ↓
可视化 (DimPlot, FeaturePlot, VlnPlot)

（1）质控
- nFeature_RNA：每个细胞检测到的基因数。过低（<200）可能是死细胞/空液滴；过高（>5000-8000）可能是双细胞
- nCount_RNA：每个细胞的总UMI数。与nFeature相关
- percent.mt：线粒体基因比例。过高（>10-20%）提示细胞膜破损，细胞质RNA流失
（2）标准化
- LogNormalize：对每个细胞的总UMI数进行标准化，然后取log1p
$$norm_count = \log\left( \frac{count}{total_UMI} \times 10000 + 1 \right)$$
- SCTransform：使用正则化负二项回归进行标准化，同时去除技术变异（如测序深度影响），效果通常优于简单标准化
（3）降维与聚类
- PCA：线性降维，保留主要变异方向
- UMAP（Uniform Manifold Approximation and Projection）：非线性降维，保留局部结构的同时拉开不同聚类
- Louvain聚类：基于共享最近邻（SNN）图的社区检测算法
（4）细胞类型注释
- 根据已知marker基因的表达进行手动注释
- 或使用自动注释工具（SingleR, scCATCH, celldex）
- 常用marker：
- T细胞：CD3D, CD3E
- B细胞：CD79A, MS4A1
- 单核/巨噬细胞：CD14, LYZ
- NK细胞：NKG7, GNLY

当多个样本或批次的数据需要合并分析时，批次效应是一个严重问题。
整合方法：
- CCA（Canonical Correlation Analysis）：Seurat v3的整合方法，找到跨批次共享的相关基因模块
- Harmony：快速、灵活的整合方法，适用于大规模数据
- scVI：基于变分自编码器的深度学习方法
- ComBat/sva：统计学校正方法

拟时序分析推断细胞在生物学过程中的发展轨迹：
- Monocle：使用反向图嵌入（RGE）构建发育轨迹
- Slingshot：基于降维后的聚类结果推断分支轨迹
- RNA velocity：利用未剪接和已剪接mRNA的比例推断细胞状态转换的方向

转录组学作为连接静态基因组和动态功能世界的桥梁，在过去十五年中经历了从基因芯片到RNA-seq再到单细胞测序的技术跃迁。每一次技术革新都带来了全新的生物学发现维度，深刻改变了我们对基因表达调控、细胞异质性和生命过程动态性的理解。

本章核心要点回顾

RNA-seq基础分析：从原始FASTQ数据到生物学洞见的完整流程包括质控、比对、定量、差异分析和功能富集。剪接比对（STAR/HISAT2）是处理真核生物数据的关键，TPM是推荐的标准化表达量指标，DESeq2基于负二项分布模型进行差异分析，而GSEA等富集分析则将基因列表转化为可解释的生物学发现。
非编码RNA：从曾被视为"转录噪音"到如今被确认为基因调控网络的核心参与者，非编码RNA的研究彻底改变了我们对基因组功能的认识。miRNA通过种子区域调控数百个靶基因，lncRNA通过多种机制在多个层次参与基因调控，circRNA则以其独特的稳定性开辟了新的调控维度。
单细胞转录组学：scRNA-seq技术使我们首次能够在单细胞分辨率上解析转录组异质性。10x Genomics的液滴微流控平台实现了高通量单细胞捕获，Seurat分析流程提供了从质控到注释的系统化方法。拟时序分析、RNA velocity等计算工具进一步使我们能够推断细胞状态转换的动态过程。

未来展望

空间转录组学：将转录组信息与组织空间位置结合，揭示基因表达的空间模式。Visium、Stereo-seq等技术正在将转录组学从"单细胞"推向"空间单细胞"时代。
多模态单细胞技术：同时测量单个细胞的转录组、表观基因组（scATAC-seq）、蛋白质组（CITE-seq）和免疫组库（scTCR/BCR-seq），实现多维度细胞状态的系统解析。
长读长RNA-seq：PacBio Iso-Seq和Nanopore direct RNA-seq能够产生全长转录本序列，彻底解决了短读长无法精确确定异构体的问题，将推动可变剪接和融合基因的精准研究。
单细胞多组学整合分析：结合基因组、转录组、表观遗传组和空间信息的整合分析方法，将使我们能够构建从基因型到表型的完整调控图谱。
临床转化：转录组生物标志物正在从研究走向临床。基于血液转录组的疾病诊断、基于单细胞图谱的肿瘤免疫分型、基于lncRNA/circRNA的新型标志物，都有望在精准医学中发挥重要作用。
转录组学的发展轨迹清晰地展示了技术创新如何驱动科学发现。从测量细胞群体的平均表达，到解析单细胞的异质性，再到定位空间位置信息，我们对生命活动动态过程的理解正在不断深化。作为生命科学的" read-out "技术，转录组学将继续在基础研究和临床转化中发挥核心作用。

本章思考题

1. 比较bulk RNA-seq和scRNA-seq在数据特征和分析方法上的主要差异。为什么scRNA-seq需要专门的分析工具？
2. 如果你要研究一种新型lncRNA在癌症中的功能，请设计一个结合RNA-seq、RIP-seq和CRISPR的实验方案。
3. 讨论单细胞转录组学面临的主要技术限制（如dropout、双细胞、批次效应），以及当前的研究如何尝试克服这些限制。
4. 空间转录组学如何改变了我们对组织生物学的理解？与传统scRNA-seq相比，空间信息带来了哪些额外的生物学洞见？
5. 展望10年后的转录组学，你认为哪些技术或方法会成为主流？转录组学将如何融入精准医学的临床实践？

推荐阅读
1. Mortazavi A, Williams B A, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq [J]. Nature methods, 2008, 5: 621-628. （RNA-seq奠基之作）
2. Trapnell C, Roberts A, Goff L, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks [J]. Nature protocols, 2012, 7: 562-578.
3. Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [J]. Genome biology, 2014, 15: 550.
4. Satija R, Farrell J A, Gennert D, et al. Spatial reconstruction of single-cell gene expression data [J]. Nature biotechnology, 2015, 33: 495-502. （Seurat经典论文）
5. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity [J]. Nature reviews immunology, 2018, 18: 35-45.

> 本章学习目标：掌握转录组学的核心概念、RNA-seq技术的完整分析流程、非编码RNA的分类与功能、单细胞转录组学的原理与方法，能够独立设计和执行转录组学研究项目。

中心法则与转录调控：
DNA →（转录）→ RNA →（翻译）→ 蛋白质
转录是基因表达的第一步，也是最关键的调控节点之一。生物体通过精确调控转录过程来实现：
- 时空特异性表达：不同组织、不同发育阶段表达不同的基因
- 环境响应：对外界刺激（温度、光照、营养、病原体等）做出快速反应
- 稳态维持：通过反馈调控维持细胞内环境稳定
转录调控的多层次机制：
1. 转录水平：转录因子结合、染色质可及性、DNA甲基化等
2. 转录后水平：mRNA剪接、编辑、稳定性调控
3. 翻译水平：翻译效率、miRNA调控
4. 翻译后水平：蛋白质修饰、定位、降解
转录组学主要聚焦于转录水平和转录后水平的研究。

第一代转录组技术：基因芯片（Microarray）
- 原理：将已知基因的寡核苷酸探针固定在芯片上，与标记的cDNA/cRNA杂交，通过荧光信号强度定量基因表达
- 局限：
- 只能检测已知序列（探针设计依赖已有基因组注释）
- 动态范围有限（通常2-3个数量级）
- 存在交叉杂交问题
- 难以检测新转录本和可变剪接
第二代转录组技术：RNA-seq
- 2008年首次报道（Mortazavi et al., Nature Methods）
- 原理：将RNA逆转录为cDNA，构建测序文库，进行高通量测序
- 优势：
- 不需要预先知道基因组序列
- 动态范围广（>5个数量级）
- 可同时检测已知和新转录本
- 可检测可变剪接、融合基因等
- 可检测SNP和等位基因特异性表达
第三代转录组技术：长读长RNA-seq
- PacBio Iso-Seq：产生全长转录本序列，无需组装即可确定完整异构体
- Nanopore direct RNA-seq：直接测序RNA分子，保留碱基修饰信息

（1）RNA提取与质控
- 使用Trizol或柱式提取法从样本中提取总RNA
- RNA完整性评估：Agilent Bioanalyzer检测RIN值（RNA Integrity Number）
- RIN > 7通常被认为是良好质量的RNA；RIN < 5可能导致3'端偏倚
（2）RNA富集或rRNA去除
- mRNA富集：使用oligo(dT)磁珠捕获带poly(A)尾的mRNA
- 优点：简单高效
- 缺点：会丢失不带poly(A)尾的RNA（如部分lncRNA、原核生物mRNA）
- rRNA去除：使用rRNA探针杂交去除rRNA（Ribo-Zero, Ribo-off等）
- 优点：保留所有非rRNA（包括mRNA、lncRNA、circRNA等）
- 缺点：操作稍复杂，成本略高
（3）RNA片段化与cDNA合成
- 片段化：使用二价阳离子（如Mg²⁺）在加热条件下将RNA随机打断成200-300 nt片段
- cDNA合成：
- 第一条链：使用随机六聚体或oligo(dT)引物，逆转录酶合成
- 第二条链：使用DNA聚合酶合成互补链（通常掺入dUTP以标记第二条链，用于链特异性文库）
（4）文库构建
- 末端修复：将cDNA片段末端补平
- A尾添加：在3'端添加单个A碱基
- 接头连接：连接带有barcode的测序接头
- PCR扩增：扩增文库至足够浓度
- 文库质控：检测文库大小分布（Agilent Bioanalyzer）和浓度（Qubit/qPCR）
（5）测序策略
- 单端测序（Single-end）：只测cDNA片段的一端，成本较低
- 双端测序（Paired-end）：测序两端，可更好地定位reads和识别可变剪接
- 链特异性测序（Strand-specific）：保留RNA链方向信息，可鉴定反义转录本

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, MultiQC)
    ↓
去除低质量和接头序列 (Trimmomatic, cutadapt)
    ↓
比对到参考基因组 (STAR, HISAT2)
    ↓
比对后处理 (排序, 索引)
    ↓
表达量定量 (featureCounts, HTSeq, Salmon)
    ↓
原始count矩阵
    ↓
数据探索 (PCA, 样本相关性热图)
    ↓
差异表达分析 (DESeq2, edgeR, limma)
    ↓
差异表达基因列表
    ↓
功能富集分析 (clusterProfiler, DAVID, GSEA)
    ↓
生物学解释和可视化

RNA-seq reads可能跨越一个或多个剪接位点（splice junction），这要求比对软件能够处理"分裂比对"（split alignment）。
STAR比对算法：
1. 种子生成：将read切分为多个短片段（seeds），在基因组上搜索精确匹配
2. 聚类：将基因组上接近的种子聚类为候选比对区域
3. 拼接：对于跨越内含子的reads，在候选区域之间构建剪接位点
4. 评分：对完整read进行比对评分，选择最优比对位置
STAR的特点：
- 速度快（比TopHat2快约50倍）
- 准确性高
- 内存需求大（需要约30-40 GB RAM用于人类基因组）
HISAT2比对算法：
- 使用基于图的比对方法（graph-based alignment）
- 利用FM-index和局部哈希实现快速搜索
- 对剪接变异的检测更敏感
- 内存效率高（<10 GB for human genome）

（1）基于count的方法
统计比对到每个基因区域的reads数。
- featureCounts：快速、准确，广泛使用的计数工具
- HTSeq-count：早期标准工具，灵活性高
（2）基于比对的定量（RPKM/FPKM/TPM）
原始count需要标准化才能进行样本内和样本间比较：
- RPKM（Reads Per Kilobase of exon per Million mapped reads）：
$$RPKM = \frac{ reads \times 10^9 }{ gene_length \times total_mapped_reads }$$
- FPKM（Fragments Per Kilobase of exon per Million mapped fragments）：
与RPKM类似，但使用fragment（pair-end测序中的一对reads）代替reads
- TPM（Transcripts Per Million）：
$$TPM_i = \frac{ RPKM_i }{ \sum_j RPKM_j } \times 10^6$$
TPM vs FPKM的关键区别：
- FPKM先标准化到总reads数，再除以基因长度
- TPM先除以基因长度得到RPKM-like值，再标准化到总和为10⁶
- TPM的优势：所有样本的TPM总和相同（10⁶），更适合样本间比较
- FPKM的问题：样本间总FPKM不同，受高表达基因影响大
（3）基于转录本估计的定量（Salmon, Kallisto, RSEM）
不依赖基因组比对，直接从reads推断转录本丰度：
- 使用伪比对（pseudoalignment）或EM算法
- 速度极快（Kallisto比传统方法快约100倍）
- 可直接输出TPM值

（1）为什么需要专门的统计方法
RNA-seq count数据具有以下统计特征：
- 非负整数（计数数据）
- 均值与方差相关（高表达基因的方差也大）
- 存在过度离散（overdispersion，方差>均值）
- 很多基因在部分样本中count为0（零膨胀）
因此不能使用普通t检验或ANOVA，需要专门的方法。
（2）DESeq2
基于负二项分布（Negative Binomial）模型：
1. 估计大小因子（Size Factor）：校正测序深度差异
$$size_factor_i = median_j \left( \frac{K_{ij}}{\left( \prod_{v=1}^{m} K_{vj} \right)^{1/m}} \right)$$
2. 估计离散度（Dispersion）：每个基因的变异程度
- 使用经验贝叶斯收缩（shrinkage）稳定低表达基因的离散度估计
3. 拟合广义线性模型（GLM）：
$$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$$
$$\log_2(\mu_{ij}) = \beta_0 + \beta_1 \cdot condition_j$$
4. Wald检验：检验系数β₁是否显著不为0
5. 多重检验校正：使用Benjamini-Hochberg方法控制FDR
（3）edgeR
同样基于负二项分布，但使用不同的离散度估计策略：
- 使用经验贝叶斯（empirical Bayes）从所有基因中借用信息
- 对小样本量（n<5）有较好的表现
差异表达基因的筛选标准：
- |log₂FC| ≥ 1（表达量变化≥2倍）
- padj < 0.05（校正后p值）

（1）火山图（Volcano Plot）
- X轴：log₂ Fold Change
- Y轴：-log₁₀(p-value)
- 每个点代表一个基因
- 同时展示效应大小（FC）和统计显著性（p-value）
- 右上和左上的点分别代表显著上调和下调的基因
（2）MA图（MA Plot）
- X轴：log₂(mean expression) —— 平均表达量
- Y轴：log₂ Fold Change
- 用于检查数据标准化效果：好的标准化应该使低表达和高表达基因的FC分布对称
（3）热图（Heatmap）
- 展示差异表达基因在不同样本中的表达模式
- 通常结合层次聚类，识别共表达基因模块
（4）PCA图
- 展示样本间的整体相似性
- 用于检测异常样本和批次效应

（1）ORA（Over-Representation Analysis）
- 输入：差异表达基因列表
- 方法：Fisher精确检验或超几何检验
- 输出：显著富集的GO term或KEGG pathway
- 工具：clusterProfiler, DAVID, g:Profiler
（2）GSEA（Gene Set Enrichment Analysis）
- 输入：所有基因的排序列表（如按log₂FC排序）
- 特点：不需要预设阈值，考虑所有基因
- 原理：检验预定义基因集（如某个通路）中的基因是否集中在排序列表的顶端或底端
- 输出：ES（Enrichment Score）和NES（Normalized Enrichment Score）
（3）GO语义相似性
- GO term之间存在层次关系
- 两个基因的功能相似性可以用其GO term在DAG中的距离来衡量

（1）动物miRNA的生物合成路径
1. 初级转录物（pri-miRNA）：由RNA聚合酶II转录，带有帽子结构和poly(A)尾，长度可达数千nt
2. 核内加工：Drosha（RNase III型核酸酶）在DGCR8辅助下，将pri-miRNA切割为~60-70 nt的pre-miRNA（发夹结构）
3. 核质运输：Exportin-5介导pre-miRNA从细胞核运输到细胞质
4. 成熟加工：Dicer酶（RNase III型）进一步切割pre-miRNA，形成~22 bp的miRNA:miRNA双链
5. RISC装载：AGO（Argonaute）蛋白结合双链，保留成熟miRNA链（通常5'端热力学稳定性较低的一条），降解另一条（miRNA）
（2）miRNA的调控机制
miRNA通过碱基配对识别靶mRNA：
- 种子区域（seed region）：miRNA的5'端第2-8位核苷酸，是靶识别的关键区域
- 完全互补（植物中常见）：导致靶mRNA的切割降解
- 不完全互补（动物中常见）：抑制翻译或促进mRNA降解（去腺苷酸化）
一个miRNA可以调控数百个靶基因，因此miRNA构成了复杂的基因调控网络。
（3）miRNA的计算预测
- 靶基因预测工具：TargetScan（基于种子互补性和保守性）、miRanda（基于自由能）、PITA（考虑靶位点可结合性）
- 新miRNA预测：基于发夹结构特征（MFEI > 0.85）和DCL-1加工依赖性

lncRNA的主要特征：
- 组织特异性和发育阶段特异性表达
- 表达量通常较低
- 物种间保守性相对较差（比蛋白质编码基因低）
- 可在表观遗传、转录及转录后等多层面调控基因表达
lncRNA的功能机制：
1. 顺式调控邻近基因：如XIST介导X染色体失活
2. 反式调控远端基因：如HOTAIR招募PRC2复合物到多个基因组位点
3. 作为miRNA海绵：lncRNA含有多个miRNA结合位点，竞争性吸附miRNA
4. 与蛋白质互作：作为支架RNA招募蛋白质复合物
5. 编码微肽：部分lncRNA含有短ORF，可翻译功能性微肽
lncRNA的鉴定：
1. 链特异性RNA-seq获取转录组数据
2. 转录本拼接（Cufflinks/StringTie）
3. 过滤已知注释基因和短转录本（<200 nt）
4. 编码潜能筛选（CPC2, CPAT, PLEK）
5. 排除具有显著蛋白质编码潜能的转录本

circRNA的生物合成：
1. 剪接体依赖的环化：下游外显子的5'端供体位点与上游的3'端受体位点结合
2. 内含子配对驱动：反向互补基序（如Alu元件）促进环化
3. 套索驱动：外显子跳过剪接产生的套索状副产物进一步环化
circRNA的独特性质：
- 高稳定性：无5'/3'端，对核酸外切酶免疫，半衰期可达数天（mRNA通常数小时）
- 组织特异性：具有较强的时间和组织表达特性
- 进化保守性：不同物种间具有一定保守性
circRNA的功能：
1. miRNA海绵：CDR1as含有63个miR-7结合位点
2. 调控亲本基因：与RNA结合蛋白互作，调节亲本基因表达
3. 翻译蛋白质：部分circRNA具有IRES和ORF，可翻译功能性蛋白
circRNA的计算鉴定：
- 主流方法：筛选不能线性比对的reads，检测反向剪接位点
- 工具：find_circ, CIRCexplorer, CIRI

早期方法（2011-2015）：
- Smart-seq/Smart-seq2：使用oligo(dT)引物和模板切换（template switching）逆转录，产生全长cDNA
- 优点：覆盖全转录本，可检测可变剪接
- 缺点：通量低（每次实验数十到数百个细胞），成本高
基于液滴的高通量方法（2015至今）：
- 10x Genomics Chromium：利用微流控技术将单个细胞与带有条形码的微珠包裹在油包水液滴中
- 优点：高通量（每次实验可测数千至数万个细胞），成本大幅下降
- 缺点：仅检测3'端（或5'端）序列，无法覆盖全转录本
- Drop-seq/InDrop：类似10x的原理，但由学术实验室开发

（1）单细胞悬液制备
- 组织解离为单细胞（机械+酶解法）
- 细胞质控：活性>90%，结团率<5%
- 细胞浓度调整至约1000 cells/μL
（2）GEM生成（Gel Bead-in-Emulsion）
使用微流控芯片同时导入三个流：
1. 单细胞悬液
2. Gel Bead（带有条形码寡核苷酸的凝胶珠）
3. 油相
在微流控通道中形成油包水液滴（GEM）：
- 每个液滴约含1个细胞 + 1个Gel Bead（泊松分布，实际捕获率约50-60%）
- Gel Bead上的寡核苷酸结构：
5'- [PCR引物] - [16bp细胞条形码] - [12bp UMI] - [30bp oligo(dT)] - 3'
（3）细胞内逆转录
- 液滴内，细胞裂解释放mRNA
- mRNA与Gel Bead上的oligo(dT)杂交
- 逆转录酶合成cDNA第一链，同时掺入细胞条形码和UMI
（4）破乳和cDNA扩增
- 破乳回收cDNA
- PCR扩增cDNA
（5）文库构建和测序
- 片段化、加接头、PCR扩增
- Illumina测序（通常2×150 bp）
（6）数据分析
- 根据reads中的细胞条形码将reads分配到不同细胞
- 根据UMI去重（同一个mRNA分子的多个PCR拷贝只计数一次）
- 生成基因×细胞的表达矩阵

Seurat是scRNA-seq分析最广泛使用的R包，由Satija Lab开发。
标准分析流程：

原始数据 (Cell Ranger输出)
    ↓
创建Seurat对象
    ↓
质控过滤 (nFeature, nCount, percent.mt)
    ↓
标准化 (LogNormalize, SCTransform)
    ↓
高变基因选择 (FindVariableFeatures)
    ↓
缩放 (ScaleData)
    ↓
线性降维 (PCA)
    ↓
非线性降维 (UMAP/t-SNE)
    ↓
聚类 (FindNeighbors + FindClusters, Louvain算法)
    ↓
细胞类型注释 (marker基因/自动注释)
    ↓
差异表达分析 (FindAllMarkers)
    ↓
可视化 (DimPlot, FeaturePlot, VlnPlot)

（1）质控
- nFeature_RNA：每个细胞检测到的基因数。过低（<200）可能是死细胞/空液滴；过高（>5000-8000）可能是双细胞
- nCount_RNA：每个细胞的总UMI数。与nFeature相关
- percent.mt：线粒体基因比例。过高（>10-20%）提示细胞膜破损，细胞质RNA流失
（2）标准化
- LogNormalize：对每个细胞的总UMI数进行标准化，然后取log1p
$$norm_count = \log\left( \frac{count}{total_UMI} \times 10000 + 1 \right)$$
- SCTransform：使用正则化负二项回归进行标准化，同时去除技术变异（如测序深度影响），效果通常优于简单标准化
（3）降维与聚类
- PCA：线性降维，保留主要变异方向
- UMAP（Uniform Manifold Approximation and Projection）：非线性降维，保留局部结构的同时拉开不同聚类
- Louvain聚类：基于共享最近邻（SNN）图的社区检测算法
（4）细胞类型注释
- 根据已知marker基因的表达进行手动注释
- 或使用自动注释工具（SingleR, scCATCH, celldex）
- 常用marker：
- T细胞：CD3D, CD3E
- B细胞：CD79A, MS4A1
- 单核/巨噬细胞：CD14, LYZ
- NK细胞：NKG7, GNLY

当多个样本或批次的数据需要合并分析时，批次效应是一个严重问题。
整合方法：
- CCA（Canonical Correlation Analysis）：Seurat v3的整合方法，找到跨批次共享的相关基因模块
- Harmony：快速、灵活的整合方法，适用于大规模数据
- scVI：基于变分自编码器的深度学习方法
- ComBat/sva：统计学校正方法

拟时序分析推断细胞在生物学过程中的发展轨迹：
- Monocle：使用反向图嵌入（RGE）构建发育轨迹
- Slingshot：基于降维后的聚类结果推断分支轨迹
- RNA velocity：利用未剪接和已剪接mRNA的比例推断细胞状态转换的方向

转录组学作为连接静态基因组和动态功能世界的桥梁，在过去十五年中经历了从基因芯片到RNA-seq再到单细胞测序的技术跃迁。每一次技术革新都带来了全新的生物学发现维度，深刻改变了我们对基因表达调控、细胞异质性和生命过程动态性的理解。

本章核心要点回顾

RNA-seq基础分析：从原始FASTQ数据到生物学洞见的完整流程包括质控、比对、定量、差异分析和功能富集。剪接比对（STAR/HISAT2）是处理真核生物数据的关键，TPM是推荐的标准化表达量指标，DESeq2基于负二项分布模型进行差异分析，而GSEA等富集分析则将基因列表转化为可解释的生物学发现。
非编码RNA：从曾被视为"转录噪音"到如今被确认为基因调控网络的核心参与者，非编码RNA的研究彻底改变了我们对基因组功能的认识。miRNA通过种子区域调控数百个靶基因，lncRNA通过多种机制在多个层次参与基因调控，circRNA则以其独特的稳定性开辟了新的调控维度。
单细胞转录组学：scRNA-seq技术使我们首次能够在单细胞分辨率上解析转录组异质性。10x Genomics的液滴微流控平台实现了高通量单细胞捕获，Seurat分析流程提供了从质控到注释的系统化方法。拟时序分析、RNA velocity等计算工具进一步使我们能够推断细胞状态转换的动态过程。

未来展望

空间转录组学：将转录组信息与组织空间位置结合，揭示基因表达的空间模式。Visium、Stereo-seq等技术正在将转录组学从"单细胞"推向"空间单细胞"时代。
多模态单细胞技术：同时测量单个细胞的转录组、表观基因组（scATAC-seq）、蛋白质组（CITE-seq）和免疫组库（scTCR/BCR-seq），实现多维度细胞状态的系统解析。
长读长RNA-seq：PacBio Iso-Seq和Nanopore direct RNA-seq能够产生全长转录本序列，彻底解决了短读长无法精确确定异构体的问题，将推动可变剪接和融合基因的精准研究。
单细胞多组学整合分析：结合基因组、转录组、表观遗传组和空间信息的整合分析方法，将使我们能够构建从基因型到表型的完整调控图谱。
临床转化：转录组生物标志物正在从研究走向临床。基于血液转录组的疾病诊断、基于单细胞图谱的肿瘤免疫分型、基于lncRNA/circRNA的新型标志物，都有望在精准医学中发挥重要作用。
转录组学的发展轨迹清晰地展示了技术创新如何驱动科学发现。从测量细胞群体的平均表达，到解析单细胞的异质性，再到定位空间位置信息，我们对生命活动动态过程的理解正在不断深化。作为生命科学的" read-out "技术，转录组学将继续在基础研究和临床转化中发挥核心作用。

本章思考题

1. 比较bulk RNA-seq和scRNA-seq在数据特征和分析方法上的主要差异。为什么scRNA-seq需要专门的分析工具？
2. 如果你要研究一种新型lncRNA在癌症中的功能，请设计一个结合RNA-seq、RIP-seq和CRISPR的实验方案。
3. 讨论单细胞转录组学面临的主要技术限制（如dropout、双细胞、批次效应），以及当前的研究如何尝试克服这些限制。
4. 空间转录组学如何改变了我们对组织生物学的理解？与传统scRNA-seq相比，空间信息带来了哪些额外的生物学洞见？
5. 展望10年后的转录组学，你认为哪些技术或方法会成为主流？转录组学将如何融入精准医学的临床实践？

推荐阅读
1. Mortazavi A, Williams B A, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq [J]. Nature methods, 2008, 5: 621-628. （RNA-seq奠基之作）
2. Trapnell C, Roberts A, Goff L, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks [J]. Nature protocols, 2012, 7: 562-578.
3. Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [J]. Genome biology, 2014, 15: 550.
4. Satija R, Farrell J A, Gennert D, et al. Spatial reconstruction of single-cell gene expression data [J]. Nature biotechnology, 2015, 33: 495-502. （Seurat经典论文）
5. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity [J]. Nature reviews immunology, 2018, 18: 35-45.

> 本章学习目标：掌握转录组学的核心概念、RNA-seq技术的完整分析流程、非编码RNA的分类与功能、单细胞转录组学的原理与方法，能够独立设计和执行转录组学研究项目。

中心法则与转录调控：
DNA →（转录）→ RNA →（翻译）→ 蛋白质
转录是基因表达的第一步，也是最关键的调控节点之一。生物体通过精确调控转录过程来实现：
- 时空特异性表达：不同组织、不同发育阶段表达不同的基因
- 环境响应：对外界刺激（温度、光照、营养、病原体等）做出快速反应
- 稳态维持：通过反馈调控维持细胞内环境稳定
转录调控的多层次机制：
1. 转录水平：转录因子结合、染色质可及性、DNA甲基化等
2. 转录后水平：mRNA剪接、编辑、稳定性调控
3. 翻译水平：翻译效率、miRNA调控
4. 翻译后水平：蛋白质修饰、定位、降解
转录组学主要聚焦于转录水平和转录后水平的研究。

第一代转录组技术：基因芯片（Microarray）
- 原理：将已知基因的寡核苷酸探针固定在芯片上，与标记的cDNA/cRNA杂交，通过荧光信号强度定量基因表达
- 局限：
- 只能检测已知序列（探针设计依赖已有基因组注释）
- 动态范围有限（通常2-3个数量级）
- 存在交叉杂交问题
- 难以检测新转录本和可变剪接
第二代转录组技术：RNA-seq
- 2008年首次报道（Mortazavi et al., Nature Methods）
- 原理：将RNA逆转录为cDNA，构建测序文库，进行高通量测序
- 优势：
- 不需要预先知道基因组序列
- 动态范围广（>5个数量级）
- 可同时检测已知和新转录本
- 可检测可变剪接、融合基因等
- 可检测SNP和等位基因特异性表达
第三代转录组技术：长读长RNA-seq
- PacBio Iso-Seq：产生全长转录本序列，无需组装即可确定完整异构体
- Nanopore direct RNA-seq：直接测序RNA分子，保留碱基修饰信息

（1）RNA提取与质控
- 使用Trizol或柱式提取法从样本中提取总RNA
- RNA完整性评估：Agilent Bioanalyzer检测RIN值（RNA Integrity Number）
- RIN > 7通常被认为是良好质量的RNA；RIN < 5可能导致3'端偏倚
（2）RNA富集或rRNA去除
- mRNA富集：使用oligo(dT)磁珠捕获带poly(A)尾的mRNA
- 优点：简单高效
- 缺点：会丢失不带poly(A)尾的RNA（如部分lncRNA、原核生物mRNA）
- rRNA去除：使用rRNA探针杂交去除rRNA（Ribo-Zero, Ribo-off等）
- 优点：保留所有非rRNA（包括mRNA、lncRNA、circRNA等）
- 缺点：操作稍复杂，成本略高
（3）RNA片段化与cDNA合成
- 片段化：使用二价阳离子（如Mg²⁺）在加热条件下将RNA随机打断成200-300 nt片段
- cDNA合成：
- 第一条链：使用随机六聚体或oligo(dT)引物，逆转录酶合成
- 第二条链：使用DNA聚合酶合成互补链（通常掺入dUTP以标记第二条链，用于链特异性文库）
（4）文库构建
- 末端修复：将cDNA片段末端补平
- A尾添加：在3'端添加单个A碱基
- 接头连接：连接带有barcode的测序接头
- PCR扩增：扩增文库至足够浓度
- 文库质控：检测文库大小分布（Agilent Bioanalyzer）和浓度（Qubit/qPCR）
（5）测序策略
- 单端测序（Single-end）：只测cDNA片段的一端，成本较低
- 双端测序（Paired-end）：测序两端，可更好地定位reads和识别可变剪接
- 链特异性测序（Strand-specific）：保留RNA链方向信息，可鉴定反义转录本

原始测序数据 (FASTQ)
    ↓
质量控制 (FastQC, MultiQC)
    ↓
去除低质量和接头序列 (Trimmomatic, cutadapt)
    ↓
比对到参考基因组 (STAR, HISAT2)
    ↓
比对后处理 (排序, 索引)
    ↓
表达量定量 (featureCounts, HTSeq, Salmon)
    ↓
原始count矩阵
    ↓
数据探索 (PCA, 样本相关性热图)
    ↓
差异表达分析 (DESeq2, edgeR, limma)
    ↓
差异表达基因列表
    ↓
功能富集分析 (clusterProfiler, DAVID, GSEA)
    ↓
生物学解释和可视化

RNA-seq reads可能跨越一个或多个剪接位点（splice junction），这要求比对软件能够处理"分裂比对"（split alignment）。
STAR比对算法：
1. 种子生成：将read切分为多个短片段（seeds），在基因组上搜索精确匹配
2. 聚类：将基因组上接近的种子聚类为候选比对区域
3. 拼接：对于跨越内含子的reads，在候选区域之间构建剪接位点
4. 评分：对完整read进行比对评分，选择最优比对位置
STAR的特点：
- 速度快（比TopHat2快约50倍）
- 准确性高
- 内存需求大（需要约30-40 GB RAM用于人类基因组）
HISAT2比对算法：
- 使用基于图的比对方法（graph-based alignment）
- 利用FM-index和局部哈希实现快速搜索
- 对剪接变异的检测更敏感
- 内存效率高（<10 GB for human genome）

（1）基于count的方法
统计比对到每个基因区域的reads数。
- featureCounts：快速、准确，广泛使用的计数工具
- HTSeq-count：早期标准工具，灵活性高
（2）基于比对的定量（RPKM/FPKM/TPM）
原始count需要标准化才能进行样本内和样本间比较：
- RPKM（Reads Per Kilobase of exon per Million mapped reads）：
$$RPKM = \frac{ reads \times 10^9 }{ gene_length \times total_mapped_reads }$$
- FPKM（Fragments Per Kilobase of exon per Million mapped fragments）：
与RPKM类似，但使用fragment（pair-end测序中的一对reads）代替reads
- TPM（Transcripts Per Million）：
$$TPM_i = \frac{ RPKM_i }{ \sum_j RPKM_j } \times 10^6$$
TPM vs FPKM的关键区别：
- FPKM先标准化到总reads数，再除以基因长度
- TPM先除以基因长度得到RPKM-like值，再标准化到总和为10⁶
- TPM的优势：所有样本的TPM总和相同（10⁶），更适合样本间比较
- FPKM的问题：样本间总FPKM不同，受高表达基因影响大
（3）基于转录本估计的定量（Salmon, Kallisto, RSEM）
不依赖基因组比对，直接从reads推断转录本丰度：
- 使用伪比对（pseudoalignment）或EM算法
- 速度极快（Kallisto比传统方法快约100倍）
- 可直接输出TPM值

（1）为什么需要专门的统计方法
RNA-seq count数据具有以下统计特征：
- 非负整数（计数数据）
- 均值与方差相关（高表达基因的方差也大）
- 存在过度离散（overdispersion，方差>均值）
- 很多基因在部分样本中count为0（零膨胀）
因此不能使用普通t检验或ANOVA，需要专门的方法。
（2）DESeq2
基于负二项分布（Negative Binomial）模型：
1. 估计大小因子（Size Factor）：校正测序深度差异
$$size_factor_i = median_j \left( \frac{K_{ij}}{\left( \prod_{v=1}^{m} K_{vj} \right)^{1/m}} \right)$$
2. 估计离散度（Dispersion）：每个基因的变异程度
- 使用经验贝叶斯收缩（shrinkage）稳定低表达基因的离散度估计
3. 拟合广义线性模型（GLM）：
$$K_{ij} \sim NB(\mu_{ij}, \alpha_i)$$
$$\log_2(\mu_{ij}) = \beta_0 + \beta_1 \cdot condition_j$$
4. Wald检验：检验系数β₁是否显著不为0
5. 多重检验校正：使用Benjamini-Hochberg方法控制FDR
（3）edgeR
同样基于负二项分布，但使用不同的离散度估计策略：
- 使用经验贝叶斯（empirical Bayes）从所有基因中借用信息
- 对小样本量（n<5）有较好的表现
差异表达基因的筛选标准：
- |log₂FC| ≥ 1（表达量变化≥2倍）
- padj < 0.05（校正后p值）

（1）火山图（Volcano Plot）
- X轴：log₂ Fold Change
- Y轴：-log₁₀(p-value)
- 每个点代表一个基因
- 同时展示效应大小（FC）和统计显著性（p-value）
- 右上和左上的点分别代表显著上调和下调的基因
（2）MA图（MA Plot）
- X轴：log₂(mean expression) —— 平均表达量
- Y轴：log₂ Fold Change
- 用于检查数据标准化效果：好的标准化应该使低表达和高表达基因的FC分布对称
（3）热图（Heatmap）
- 展示差异表达基因在不同样本中的表达模式
- 通常结合层次聚类，识别共表达基因模块
（4）PCA图
- 展示样本间的整体相似性
- 用于检测异常样本和批次效应

（1）ORA（Over-Representation Analysis）
- 输入：差异表达基因列表
- 方法：Fisher精确检验或超几何检验
- 输出：显著富集的GO term或KEGG pathway
- 工具：clusterProfiler, DAVID, g:Profiler
（2）GSEA（Gene Set Enrichment Analysis）
- 输入：所有基因的排序列表（如按log₂FC排序）
- 特点：不需要预设阈值，考虑所有基因
- 原理：检验预定义基因集（如某个通路）中的基因是否集中在排序列表的顶端或底端
- 输出：ES（Enrichment Score）和NES（Normalized Enrichment Score）
（3）GO语义相似性
- GO term之间存在层次关系
- 两个基因的功能相似性可以用其GO term在DAG中的距离来衡量

（1）动物miRNA的生物合成路径
1. 初级转录物（pri-miRNA）：由RNA聚合酶II转录，带有帽子结构和poly(A)尾，长度可达数千nt
2. 核内加工：Drosha（RNase III型核酸酶）在DGCR8辅助下，将pri-miRNA切割为~60-70 nt的pre-miRNA（发夹结构）
3. 核质运输：Exportin-5介导pre-miRNA从细胞核运输到细胞质
4. 成熟加工：Dicer酶（RNase III型）进一步切割pre-miRNA，形成~22 bp的miRNA:miRNA双链
5. RISC装载：AGO（Argonaute）蛋白结合双链，保留成熟miRNA链（通常5'端热力学稳定性较低的一条），降解另一条（miRNA）
（2）miRNA的调控机制
miRNA通过碱基配对识别靶mRNA：
- 种子区域（seed region）：miRNA的5'端第2-8位核苷酸，是靶识别的关键区域
- 完全互补（植物中常见）：导致靶mRNA的切割降解
- 不完全互补（动物中常见）：抑制翻译或促进mRNA降解（去腺苷酸化）
一个miRNA可以调控数百个靶基因，因此miRNA构成了复杂的基因调控网络。
（3）miRNA的计算预测
- 靶基因预测工具：TargetScan（基于种子互补性和保守性）、miRanda（基于自由能）、PITA（考虑靶位点可结合性）
- 新miRNA预测：基于发夹结构特征（MFEI > 0.85）和DCL-1加工依赖性

lncRNA的主要特征：
- 组织特异性和发育阶段特异性表达
- 表达量通常较低
- 物种间保守性相对较差（比蛋白质编码基因低）
- 可在表观遗传、转录及转录后等多层面调控基因表达
lncRNA的功能机制：
1. 顺式调控邻近基因：如XIST介导X染色体失活
2. 反式调控远端基因：如HOTAIR招募PRC2复合物到多个基因组位点
3. 作为miRNA海绵：lncRNA含有多个miRNA结合位点，竞争性吸附miRNA
4. 与蛋白质互作：作为支架RNA招募蛋白质复合物
5. 编码微肽：部分lncRNA含有短ORF，可翻译功能性微肽
lncRNA的鉴定：
1. 链特异性RNA-seq获取转录组数据
2. 转录本拼接（Cufflinks/StringTie）
3. 过滤已知注释基因和短转录本（<200 nt）
4. 编码潜能筛选（CPC2, CPAT, PLEK）
5. 排除具有显著蛋白质编码潜能的转录本

circRNA的生物合成：
1. 剪接体依赖的环化：下游外显子的5'端供体位点与上游的3'端受体位点结合
2. 内含子配对驱动：反向互补基序（如Alu元件）促进环化
3. 套索驱动：外显子跳过剪接产生的套索状副产物进一步环化
circRNA的独特性质：
- 高稳定性：无5'/3'端，对核酸外切酶免疫，半衰期可达数天（mRNA通常数小时）
- 组织特异性：具有较强的时间和组织表达特性
- 进化保守性：不同物种间具有一定保守性
circRNA的功能：
1. miRNA海绵：CDR1as含有63个miR-7结合位点
2. 调控亲本基因：与RNA结合蛋白互作，调节亲本基因表达
3. 翻译蛋白质：部分circRNA具有IRES和ORF，可翻译功能性蛋白
circRNA的计算鉴定：
- 主流方法：筛选不能线性比对的reads，检测反向剪接位点
- 工具：find_circ, CIRCexplorer, CIRI

早期方法（2011-2015）：
- Smart-seq/Smart-seq2：使用oligo(dT)引物和模板切换（template switching）逆转录，产生全长cDNA
- 优点：覆盖全转录本，可检测可变剪接
- 缺点：通量低（每次实验数十到数百个细胞），成本高
基于液滴的高通量方法（2015至今）：
- 10x Genomics Chromium：利用微流控技术将单个细胞与带有条形码的微珠包裹在油包水液滴中
- 优点：高通量（每次实验可测数千至数万个细胞），成本大幅下降
- 缺点：仅检测3'端（或5'端）序列，无法覆盖全转录本
- Drop-seq/InDrop：类似10x的原理，但由学术实验室开发

（1）单细胞悬液制备
- 组织解离为单细胞（机械+酶解法）
- 细胞质控：活性>90%，结团率<5%
- 细胞浓度调整至约1000 cells/μL
（2）GEM生成（Gel Bead-in-Emulsion）
使用微流控芯片同时导入三个流：
1. 单细胞悬液
2. Gel Bead（带有条形码寡核苷酸的凝胶珠）
3. 油相
在微流控通道中形成油包水液滴（GEM）：
- 每个液滴约含1个细胞 + 1个Gel Bead（泊松分布，实际捕获率约50-60%）
- Gel Bead上的寡核苷酸结构：
5'- [PCR引物] - [16bp细胞条形码] - [12bp UMI] - [30bp oligo(dT)] - 3'
（3）细胞内逆转录
- 液滴内，细胞裂解释放mRNA
- mRNA与Gel Bead上的oligo(dT)杂交
- 逆转录酶合成cDNA第一链，同时掺入细胞条形码和UMI
（4）破乳和cDNA扩增
- 破乳回收cDNA
- PCR扩增cDNA
（5）文库构建和测序
- 片段化、加接头、PCR扩增
- Illumina测序（通常2×150 bp）
（6）数据分析
- 根据reads中的细胞条形码将reads分配到不同细胞
- 根据UMI去重（同一个mRNA分子的多个PCR拷贝只计数一次）
- 生成基因×细胞的表达矩阵

Seurat是scRNA-seq分析最广泛使用的R包，由Satija Lab开发。
标准分析流程：

原始数据 (Cell Ranger输出)
    ↓
创建Seurat对象
    ↓
质控过滤 (nFeature, nCount, percent.mt)
    ↓
标准化 (LogNormalize, SCTransform)
    ↓
高变基因选择 (FindVariableFeatures)
    ↓
缩放 (ScaleData)
    ↓
线性降维 (PCA)
    ↓
非线性降维 (UMAP/t-SNE)
    ↓
聚类 (FindNeighbors + FindClusters, Louvain算法)
    ↓
细胞类型注释 (marker基因/自动注释)
    ↓
差异表达分析 (FindAllMarkers)
    ↓
可视化 (DimPlot, FeaturePlot, VlnPlot)

（1）质控
- nFeature_RNA：每个细胞检测到的基因数。过低（<200）可能是死细胞/空液滴；过高（>5000-8000）可能是双细胞
- nCount_RNA：每个细胞的总UMI数。与nFeature相关
- percent.mt：线粒体基因比例。过高（>10-20%）提示细胞膜破损，细胞质RNA流失
（2）标准化
- LogNormalize：对每个细胞的总UMI数进行标准化，然后取log1p
$$norm_count = \log\left( \frac{count}{total_UMI} \times 10000 + 1 \right)$$
- SCTransform：使用正则化负二项回归进行标准化，同时去除技术变异（如测序深度影响），效果通常优于简单标准化
（3）降维与聚类
- PCA：线性降维，保留主要变异方向
- UMAP（Uniform Manifold Approximation and Projection）：非线性降维，保留局部结构的同时拉开不同聚类
- Louvain聚类：基于共享最近邻（SNN）图的社区检测算法
（4）细胞类型注释
- 根据已知marker基因的表达进行手动注释
- 或使用自动注释工具（SingleR, scCATCH, celldex）
- 常用marker：
- T细胞：CD3D, CD3E
- B细胞：CD79A, MS4A1
- 单核/巨噬细胞：CD14, LYZ
- NK细胞：NKG7, GNLY

当多个样本或批次的数据需要合并分析时，批次效应是一个严重问题。
整合方法：
- CCA（Canonical Correlation Analysis）：Seurat v3的整合方法，找到跨批次共享的相关基因模块
- Harmony：快速、灵活的整合方法，适用于大规模数据
- scVI：基于变分自编码器的深度学习方法
- ComBat/sva：统计学校正方法

拟时序分析推断细胞在生物学过程中的发展轨迹：
- Monocle：使用反向图嵌入（RGE）构建发育轨迹
- Slingshot：基于降维后的聚类结果推断分支轨迹
- RNA velocity：利用未剪接和已剪接mRNA的比例推断细胞状态转换的方向

转录组学作为连接静态基因组和动态功能世界的桥梁，在过去十五年中经历了从基因芯片到RNA-seq再到单细胞测序的技术跃迁。每一次技术革新都带来了全新的生物学发现维度，深刻改变了我们对基因表达调控、细胞异质性和生命过程动态性的理解。

本章核心要点回顾

RNA-seq基础分析：从原始FASTQ数据到生物学洞见的完整流程包括质控、比对、定量、差异分析和功能富集。剪接比对（STAR/HISAT2）是处理真核生物数据的关键，TPM是推荐的标准化表达量指标，DESeq2基于负二项分布模型进行差异分析，而GSEA等富集分析则将基因列表转化为可解释的生物学发现。
非编码RNA：从曾被视为"转录噪音"到如今被确认为基因调控网络的核心参与者，非编码RNA的研究彻底改变了我们对基因组功能的认识。miRNA通过种子区域调控数百个靶基因，lncRNA通过多种机制在多个层次参与基因调控，circRNA则以其独特的稳定性开辟了新的调控维度。
单细胞转录组学：scRNA-seq技术使我们首次能够在单细胞分辨率上解析转录组异质性。10x Genomics的液滴微流控平台实现了高通量单细胞捕获，Seurat分析流程提供了从质控到注释的系统化方法。拟时序分析、RNA velocity等计算工具进一步使我们能够推断细胞状态转换的动态过程。

未来展望

空间转录组学：将转录组信息与组织空间位置结合，揭示基因表达的空间模式。Visium、Stereo-seq等技术正在将转录组学从"单细胞"推向"空间单细胞"时代。
多模态单细胞技术：同时测量单个细胞的转录组、表观基因组（scATAC-seq）、蛋白质组（CITE-seq）和免疫组库（scTCR/BCR-seq），实现多维度细胞状态的系统解析。
长读长RNA-seq：PacBio Iso-Seq和Nanopore direct RNA-seq能够产生全长转录本序列，彻底解决了短读长无法精确确定异构体的问题，将推动可变剪接和融合基因的精准研究。
单细胞多组学整合分析：结合基因组、转录组、表观遗传组和空间信息的整合分析方法，将使我们能够构建从基因型到表型的完整调控图谱。
临床转化：转录组生物标志物正在从研究走向临床。基于血液转录组的疾病诊断、基于单细胞图谱的肿瘤免疫分型、基于lncRNA/circRNA的新型标志物，都有望在精准医学中发挥重要作用。
转录组学的发展轨迹清晰地展示了技术创新如何驱动科学发现。从测量细胞群体的平均表达，到解析单细胞的异质性，再到定位空间位置信息，我们对生命活动动态过程的理解正在不断深化。作为生命科学的" read-out "技术，转录组学将继续在基础研究和临床转化中发挥核心作用。

本章思考题

1. 比较bulk RNA-seq和scRNA-seq在数据特征和分析方法上的主要差异。为什么scRNA-seq需要专门的分析工具？
2. 如果你要研究一种新型lncRNA在癌症中的功能，请设计一个结合RNA-seq、RIP-seq和CRISPR的实验方案。
3. 讨论单细胞转录组学面临的主要技术限制（如dropout、双细胞、批次效应），以及当前的研究如何尝试克服这些限制。
4. 空间转录组学如何改变了我们对组织生物学的理解？与传统scRNA-seq相比，空间信息带来了哪些额外的生物学洞见？
5. 展望10年后的转录组学，你认为哪些技术或方法会成为主流？转录组学将如何融入精准医学的临床实践？

推荐阅读
1. Mortazavi A, Williams B A, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq [J]. Nature methods, 2008, 5: 621-628. （RNA-seq奠基之作）
2. Trapnell C, Roberts A, Goff L, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks [J]. Nature protocols, 2012, 7: 562-578.
3. Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 [J]. Genome biology, 2014, 15: 550.
4. Satija R, Farrell J A, Gennert D, et al. Spatial reconstruction of single-cell gene expression data [J]. Nature biotechnology, 2015, 33: 495-502. （Seurat经典论文）
5. Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity [J]. Nature reviews immunology, 2018, 18: 35-45.