UMI计数流程：从原始数据到计数矩阵
好，咱们正式进入UMI数据分析的核心。这一章我要把整个流程拆开，讲清楚数据是怎么从原始FASTQ文件，一步步变成最终那个计数矩阵的。
你想想看，如果手里拿着一堆测序数据，里面混着UMI序列、细胞条形码、还有实际的转录本序列，怎么把它们理清楚？这就是我们接下来要解决的问题。
原始数据长什么样？
先看看手里有什么。典型的单细胞UMI数据，一条read的结构大概是这样：
@read_id
[细胞条形码(16bp)] [UMI(10bp)] [polyT] [转录本序列...]
+
[质量分数...]嗯，这里要注意，不同平台和试剂盒的UMI长度不一样。10x Genomics常用10bp，BD Rhapsody用的是8bp，还有些定制方案会用12bp。我习惯在拿到数据的第一时间，先看一眼README或者试剂说明，确认UMI到底占多少碱基。
第一步：拆分与预处理
流程的第一步，就是把这些混在一起的信息拆开。说白了，就是把reads按细胞条形码分组，同时把UMI序列单独拎出来。
具体步骤：
- 提取细胞条形码：从read头部固定位置截取
- 提取UMI序列：紧跟着条形码的固定长度
- 过滤低质量reads：条形码或UMI中如果有太多N碱基，直接扔掉
- 去除接头和polyT：避免污染转录本比对
我在项目中遇到过一批数据，UMI序列里莫名其妙出现了连续三个N。最开始没在意，后来发现这些reads的计数异常偏高。排查了半天才发现是测序仪在UMI区域出了错。从那以后，我每次都会加一个UMI质量过滤的步骤：
# 过滤UMI中包含超过1个N的reads
if umi_sequence.count('N') > 1:
discard_read()第二步：比对与基因分配
拆出转录本序列后，下一步就是比对到参考基因组。这一步其实和普通RNA-seq没有本质区别，但有个坑——UMI数据往往包含大量PCR重复，比对软件如果处理不当，会浪费很多计算资源。
我常用的是STAR或HISAT2，不过针对UMI数据，我更推荐使用专门优化的比对器：
STAR --genomeDir /ref/hg38 \
--readFilesIn sample_R1.fastq.gz \
--outSAMtype BAM SortedByCoordinate \
--outFilterMultimapNmax 1 \
--outSAMattributes NH HI nM AS比对完成后，需要把每个read关联到具体的基因。这一步通常基于参考基因组的GTF注释文件。如果比对位置落在多个基因的重叠区域，我个人的处理原则是：只保留唯一匹配。虽然会损失一些数据，但能大幅降低假阳性。
第三步：UMI校正
终于到了UMI分析最核心的部分。为什么需要校正？我当年刚开始做这个的时候也纳闷：UMI不就是为了区分真实分子和PCR重复吗？怎么还需要校正？
原因其实很简单：测序会出错。UMI序列本身也会被测错，导致同一个原始分子被识别成不同的UMI。反过来，不同分子也可能因为测序错误变成同一个UMI序列。
常见的校正方法有下面几种：
| 方法 | 原理 | 适用场景 |
|---|---|---|
| 完全相同 | 只有序列完全一致的UMI才合并 | 测序质量高，UMI长度短（<8bp） |
| 汉明距离 | 最多允许1-2个碱基差异 | 最常用，平衡灵敏度与特异性 |
| 编辑距离 | 允许插入/删除，不仅限于替换 | 测序质量一般，有indel风险 |
| 邻近图聚类 | 基于序列相似性构建网络，自动聚类 | 高深度测序，适合自动优化 |
说说我的习惯。对于10x Genomics的10bp UMI，我一般先用汉明距离=1。如果发现UMI多样性异常低（也就是校正后分子数远小于UMI数），我会降到距离=0，看看是不是过度合并了。
一个典型的校正代码片段：
def correct_umis(umi_list, max_distance=1):
"""基于汉明距离合并相似的UMI"""
# 按丰度排序，丰度高的优先保留
sorted_umis = Counter(umi_list).most_common()
corrected = {}
for umi, count in sorted_umis:
# 如果这个UMI已经被合并，跳过
found = False
for ref_umi in corrected:
if hamming_distance(umi, ref_umi) <= max_distance:
corrected[ref_umi] += count
found = True
break
if not found:
corrected[umi] = count
return corrected第四步：生成计数矩阵
校正完之后，就可以统计每个细胞每个基因的UMI计数了。这就是最后要的东西——一个细胞×基因的矩阵。
格式大概长这样：
gene cell_1 cell_2 cell_3 cell_4
GAPDH 125 89 203 45
ACTB 78 156 34 212
TP53 12 3 8 0
... ... ... ... ...生成这个矩阵的工具很多。我最熟悉的是用Python的pandas来处理：
import pandas as pd
# 假设我们已经有了一个字典：{(cell_barcode, gene): corrected_umi_count}
count_data = []
for (cell, gene), count in corrected_counts.items():
count_data.append({'cell': cell, 'gene': gene, 'count': count})
umi_matrix = pd.DataFrame(count_data).pivot_table(
index='gene', columns='cell', values='count', fill_value=0
)但这里有一个隐藏问题。并不是所有细胞条形码都对应真正的细胞。有些是空的液滴，有些是死细胞。这就涉及到下一章要讲的内容——细胞过滤。嗯，你先记住这个节点就好。
常见错误类型
做了这么多年UMI分析，我踩过的坑真不少。整理一下最常见的几种：
- UMI长度识别错误：不同试剂盒的UMI位置和长度不同，搞错了直接导致全盘崩溃。我建议建立自己的备忘清单，每次开工前先核对。
- 校正阈值过松：我曾经把汉明距离设成2，结果把大量不同真实分子的UMI合到一起，导致基因表达量被严重低估。后来改成距离=1才正常。
- 忽略PCR偏差：如果文库扩增轮次太多，UMI多样性会被耗尽，这时候校正效果会大打折扣。需要结合扩增深度来调整策略。
- 细胞条形码污染：相邻液滴的细胞条形码有可能互换。我见过一个项目，5%的UMI分配到了错误的细胞里，直接影响了后续的聚类分析。