RNA-seq 比对与表达定量
RNA-seq 上游分析把 FASTQ 读段变成可供统计的表达矩阵，包含读段质控、比对或伪比对、计数三个环节。路线一是用 STAR 或 HISAT2 比对到参考基因组，再用 featureCounts 按 GTF 计数，产生的 BAM 还能用于可变剪接、融合检测和 IGV 查看；路线二是用 Salmon 或 kallisto 对转录组做伪比对定量，速度快、资源占用低，再用 tximport 汇总到基因。最终产出原始计数矩阵和 TPM 矩阵，前者用于差异分析，后者用于展示。
什么时候用
- 拿到测序公司返回的 FASTQ，需要自己生成计数矩阵
- 从 SRA 下载公共原始数据，用统一流程重新定量以便与自有数据合并
- 需要 BAM 文件做可变剪接、融合基因检测或 IGV 查看
- 需要转录本水平的定量，用于异构体使用分析
- 排查测序公司结果异常，如比对率低、链方向设置错误
不适合的情况
- 没有参考基因组的非模式物种：应先从头组装转录组，再在组装结果上定量
- 长读长转录组（Iso-Seq、Nanopore）：应使用长读长比对和专用定量工具
- 小 RNA 测序：读段短、需去接头后按 miRNA 注释计数，流程不同
需要准备的数据
- FASTQ 文件
- 单端或双端 .fastq.gz，双端 R1、R2 成对且读段顺序一致；记录读长和测序平台
- 参考基因组 FASTA
- 与注释同一来源和版本，如 Ensembl 或 GENCODE 的 GRCh38，染色体命名须与 GTF 一致
- GTF 注释文件
- 决定计数单位和基因 ID 体系，版本号写入方法部分，后续富集和注释都依赖它
- 文库信息
- 是否链特异、链方向（如 dUTP 法）、polyA 富集还是去 rRNA，直接决定计数参数
分析步骤
- 1
读段质控与修剪
用 FastQC 和 MultiQC 检查碱基质量、接头含量和 GC 分布，接头明显时用 fastp 或 Trim Galore 修剪。RNA-seq 重复率偏高、前十几个碱基组成波动属于正常现象，细节见测序数据质控页。
- 2
构建索引
STAR 建索引时同时提供 GTF，sjdbOverhang 设为读长减 1，人类基因组建索引通常需要约 32 GB 内存；HISAT2 内存需求低得多。Salmon 建议把基因组序列作为 decoy 加入索引，减少错误分配。
- 3
比对或伪比对
STAR 加 --quantMode GeneCounts 可顺带输出基因计数；Salmon 用 -l A 自动推断文库类型。STAR 的 Log.final.out 中，人和小鼠常规文库唯一比对率经验上多在 70% 以上。
- 4
判断链特异性
用 RSeQC 的 infer_experiment.py，或比较 STAR ReadsPerGene.out.tab 第 3、4 列总数：一列远大于另一列为链特异文库，两列接近为非链特异。
- 5
基因计数
featureCounts 用 -s 指定链方向（0 不区分、1 正向、2 反向），双端数据加 -p。.summary 中 Assigned 比例低而 Unassigned_NoFeatures 高，多半是链方向或注释不匹配。
- 6
汇总与归一化
合并为基因 × 样本整数计数矩阵；Salmon 结果用 tximport 按转录本到基因的对应表汇总。另计算 TPM 用于作图。用 MultiQC 汇总各步指标，作为方法描述和补充材料。
关键参数与决策
链方向参数
featureCounts 的 -s、HTSeq-count 的 --stranded、Salmon 的 libType 须与文库一致。dUTP 法文库通常对应 -s 2，设错会使计数损失大半。
多重比对读段
featureCounts 默认不计多重比对读段，加 -M 才计入（可配合 --fraction 平均分配）。基因水平差异分析一般保持默认；研究重复元件或高度同源的基因家族时需要专门方法。
sjdbOverhang
STAR 建索引时理想值为读长减 1，读长不一时取最大读长减 1。默认值 100 对多数读长都能正常工作，但读长与之相差很大时，按实际读长重建索引更稳妥。
计数特征与 ID
featureCounts 默认 -t exon、-g gene_id，按外显子计数后汇总到基因。GENCODE 的 gene_id 带版本号后缀（如 .12），与其他数据合并前通常要去掉。
Salmon 偏好校正
建议加 --gcBias 和 --seqBias，分别校正片段 GC 含量偏好和读段起始序列偏好，双端数据中 GC 校正更有效。新版本 Salmon 默认使用选择性比对模式。
常见错误
- 基因组与 GTF 来源不一致：UCSC 写作 chr1，Ensembl 写作 1，混用时 featureCounts 几乎全判为 Unassigned。
- HTSeq-count 不显式设置 --stranded：它的默认值按链特异处理，非链特异文库会丢失约一半读段，而且程序不会报错。
- 按 DNA-seq 的习惯去 PCR 重复：高表达基因天然产生大量起点相同的读段，没有 UMI 时去重会压低高表达基因的计数。
- 双端数据按单端计数：featureCounts 不加 -p 时一个片段被计为两条读段，计数虚增，跨两个基因的片段还会被重复分配。
- 直接合并不同研究的 FPKM：注释版本和计数规则不同，数值不可比；合并公共数据应从 FASTQ 用同一流程重新定量。
- 不记录参考版本：GRCh37 与 GRCh38 坐标不同，后续与变异、峰或剪接位点比较时会出现系统性错位，难以追查。
结果怎么看
MultiQC 汇总报告
集中看各样本总读段数、唯一比对率、Assigned 比例和重复率。某样本在多项指标上同时离群，是后续剔除或在差异分析中标注的依据，单项略低通常不需处理。
STAR Log.final.out
关注唯一比对率、多重比对比例和 unmapped: too short 比例。too short 比例高常见于接头未去除、读段被过度修剪或物种不符；多重比对比例高提示 rRNA 残留较多。
featureCounts summary
NoFeatures 类高提示链方向或注释问题；Ambiguity 类高提示重叠基因多或链方向设为不区分；MultiMapping 类高提示 rRNA 或重复序列多。
基因体覆盖度
用 RSeQC 的 geneBody_coverage 看读段沿基因 5′ 到 3′ 的分布。理想情况较为均匀；明显偏向 3′ 端提示 RNA 降解，在 FFPE 样本中常见，差异分析时需考虑这一偏差。
计数矩阵与 TPM 矩阵
计数矩阵为整数，供 DESeq2、edgeR 使用；TPM 每个样本总和为一百万，适合作图和样本内基因比较。两种矩阵的基因 ID 体系和顺序应一致，并注明注释版本。
和相近方法怎么选
STAR 与 HISAT2
两者都是能识别剪接的基因组比对器。STAR 速度快、剪接位点检出灵敏，可直接输出基因计数和嵌合比对（融合检测依赖它），但内存需求高；HISAT2 内存占用小，适合普通工作站。需要做融合检测时优先 STAR。
基因组比对与 Salmon、kallisto
Salmon、kallisto 不产生基因组 BAM，速度快很多，并能按概率把多重比对读段分配到转录本，只需表达定量时很合适。需要看覆盖度、发现新剪接位点、检测融合或 RNA 变异时，仍需基因组比对。
TPM、FPKM 与原始计数
FPKM 先按测序深度再按长度归一化，各样本总和不固定；TPM 先按长度再按深度，各样本总和相同，样本间比较更合理。两者都不能作为差异分析输入，差异分析用原始计数，由工具内部做标准化。
featureCounts 与 HTSeq-count
两者都按 GTF 计数。featureCounts 速度快很多，支持多线程，一次处理多个 BAM 并输出统计摘要；HTSeq-count 规则直观但较慢。两者的默认链方向设置不同，切换工具时要重新核对。
常见问题
TPM 和 FPKM 有什么区别，用哪个？
做样本间展示优先 TPM，因为每个样本的 TPM 总和相同，FPKM 总和会随样本组成变化。但 TPM 只校正了深度和长度，不能消除组成偏差，差异分析仍要用原始计数。
RNA-seq 比对率多少算正常？
人和小鼠常规 polyA 文库在 STAR 中的唯一比对率经验上多在 70% 以上。偏低时依次检查：物种和基因组是否正确、接头是否去除、rRNA 比例、是否有细菌或支原体污染，可用 FastQ Screen 抽查污染来源。
featureCounts 的 -s 应该设 0、1 还是 2？
最直接的办法是实测：用 -s 1 和 -s 2 分别跑一个样本，Assigned 比例明显更高的就是正确方向；如果两者接近、都约为 -s 0 的一半，则是非链特异文库，应设 0。
没有服务器能做 RNA-seq 上游分析吗？
可以。Salmon 和 kallisto 对内存要求低，普通电脑通常即可完成人类样本定量；HISAT2 可下载官方预建索引后比对。STAR 建人类索引需要约 32 GB 内存，样本多时更适合放在服务器上批量运行。
Salmon 的结果怎么用于 DESeq2？
用 tximport 读取各样本的 quant.sf，并提供转录本到基因的对应表（tx2gene），汇总为基因水平后用 DESeqDataSetFromTximport 导入。DESeq2 会利用平均转录本长度作偏移，校正样本间的长度差异。
在 Codar 里怎么做
Codar 先读取 FASTQ 做质量评估；与你确认参考基因组和注释版本、文库链方向和定量路线后，执行读段比对与表达定量，生成基因计数矩阵和 TPM 矩阵。交付质控汇总、比对与计数统计、表达矩阵和可复现脚本。人类样本建索引内存需求高，专业版可通过远程 SSH 在服务器上运行。
我有 12 个人类双端链特异 RNA-seq 的 FASTQ，请先做质控，再用 GRCh38 比对，输出 featureCounts 计数矩阵和 TPM