简介本资源为基于单细胞测序数据的转座元件TEs表达量化源码包面向生物信息学研究者与单细胞数据分析人员旨在解决单细胞水平上TEs表达难以精确量化的问题。包内共42个文件以21个Python脚本和6个Shell脚本承担数据分析与自动化流程辅以2个Jupyter Notebook交互式分析示例、2个Markdown文档说明并包含gtf、bed、BAM、idx等关键生物信息学格式文件压缩包约34.63MB。源码借鉴自JiekaiLab的scTE项目提供从数据准备、表达量化到结果展示的完整流程其中scTE可执行文件可直接处理BAM文件并自动调用相关脚本example目录则演示了具体分析用法便于初学者快速上手。目前已有327人学习下载适合希望深入探究TEs在细胞异质性与基因表达调控中作用的研究者参考使用。1. 单细胞测序数据里转座元件表达量化到底难在哪单细胞测序分析流程走到表达定量这一步绝大多数人只盯着基因把转座元件TEs当成基因组里的“暗物质”直接跳过。但如果你做过肿瘤微环境、早期胚胎发育或者神经退行性疾病的单细胞项目就会发现一个反直觉的现象同一群细胞在基因层面看起来高度一致TEs 表达谱却能把它们切成好几个亚群。这正是转座元件表达量化在单细胞测序里的价值——它提供了一层基因表达之外的调控信息。问题在于TEs 的量化跟普通基因完全不是一回事。TEs 在基因组里以多拷贝家族形式存在同一个家族的不同拷贝序列高度相似短读长测序的 reads 很容易在多拷贝之间发生错误分配。再加上单细胞数据本身 dropout 率高、每个细胞只有几千条 UMITEs 的表达信号就更稀疏。所以“基于单细胞测序数据的转座元件表达量化设计源码”这个标题本质上要解决的是怎么在单细胞分辨率下把 TEs 的多拷贝分配、UMI 去重和表达定量串成一条可复现的流程并且把这条流程写成能跑起来的源码。这套方案适合两类人一是正在做单细胞转录组分析、想从基因层面往下挖一层调控机制的从业者二是需要把 TEs 量化能力集成到自己分析管线里的生信工程师。下面我从数据准备、多拷贝分配、UMI 定量到源码结构把这条路径拆开讲清楚。2. 转座元件量化的数据准备与参考构建2.1 为什么不能直接用基因注释的 GTF 跑 TEs常规单细胞流程里比对完直接拿基因注释 GTF 做 featureCounts 或 STARsolo这套逻辑对 TEs 完全不适用。原因有三个第一TEs 注释不在标准基因 GTF 里你需要额外引入 RepeatMasker 或 Dfam 的注释文件第二TEs 的注释是区间形式很多拷贝互相重叠featureCounts 默认的“一个 read 只分配给一个 feature”策略会直接丢掉大量 ambiguous reads第三单细胞数据需要 UMI 级别的定量而 TEs 的多拷贝特性让 UMI 去重变得复杂——同一个 UMI 可能来自不同拷贝但序列一样你没法简单去重。常见做法是先构建一个包含基因和 TEs 的联合参考然后在比对阶段就允许 reads 多位置比对把多拷贝分配留到定量阶段用 EM 算法或唯一分子标签来解决。我一般会推荐用 STAR 建索引时把 TEs 区间也加进去或者用 Cell Ranger 的 custom reference 流程把 TEs 注释转成 GTF 格式后合并。2.2 构建 TEs 注释 GTF 的具体步骤假设你手头有 RepeatMasker 的.out文件或者 Dfam 的.hmm结果第一步是转成标准 GTF。下面这段 Python 脚本把 RepeatMasker 输出转成 GTF并给每个 TEs 拷贝分配唯一 IDimport pandas as pd # 读取 RepeatMasker .out 文件跳过前两行注释 # 列顺序SW score, perc div, perc del, perc ins, query seq, begin, end, ... rm pd.read_csv(repeatmasker.out, sepr\s, skiprows2, headerNone) rm.columns [ sw_score, perc_div, perc_del, perc_ins, query_seq, q_begin, q_end, q_left, strand, repeat_name, repeat_class, repeat_family, r_begin, r_end, r_left, id ] # 只保留有明确家族分类的记录过滤低复杂度 rm rm[rm[repeat_class].notna()] # 构建 GTF 行seqname, source, feature, start, end, score, strand, frame, attribute gtf_lines [] for idx, row in rm.iterrows(): # 每个拷贝用 repeat_name 坐标做唯一 ID避免多拷贝混淆 te_id f{row[repeat_name]}_{row[query_seq]}_{row[q_begin]}_{row[q_end]} # GTF 坐标从 1 开始且 end 是闭区间 start int(row[q_begin]) end int(row[q_end]) strand row[strand] if row[strand] in [, -] else attr fgene_id {te_id}; transcript_id {te_id}; family {row[repeat_family]}; class {row[repeat_class]}; gtf_lines.append( f{row[query_seq]}\tRepeatMasker\texon\t{start}\t{end}\t.\t{strand}\t.\t{attr} ) with open(tes_annotation.gtf, w) as f: f.write(\n.join(gtf_lines))这段代码的逻辑是把每个 TEs 拷贝当成一个独立的“基因”来处理用repeat_name 染色体 起止坐标拼成唯一 ID。这样在后续定量时每个拷贝有独立的表达值你可以再按 family 或 class 聚合。参数上要注意q_begin和q_end是相对于 query 序列的坐标如果 RepeatMasker 跑的是全基因组query_seq 就是染色体名坐标可以直接用。如果跑的是 contig需要先做 liftOver 或者用原始基因组坐标。2.3 合并基因和 TEs 注释时的冲突处理把 TEs GTF 和基因 GTF 合并时最大的坑是区间重叠。很多 TEs 就插在基因内含子或 UTR 里如果直接合并STAR 或 Cell Ranger 建索引时会报重复 feature 或者把 reads 错误分配。我一般会做两步过滤第一去掉与基因外显子重叠超过 50% 的 TEs 拷贝第二对于基因内 TEs单独标记为intronic_TE在定量时单独统计。下面这段代码用 pybedtools 做重叠过滤import pybedtools genes pybedtools.BedTool(genes.bed) tes pybedtools.BedTool(tes_annotation.bed) # 找与基因外显子重叠的 TEs overlaps tes.intersect(genes, waTrue, wbTrue) # 计算重叠比例过滤掉重叠超过 50% 的 filtered_tes [] for interval in tes: te_len interval.end - interval.start ov interval.intersect(genes, waTrue, wbTrue) if ov.count() 0: filtered_tes.append(interval) else: max_ov max([min(interval.end, o.end) - max(interval.start, o.start) for o in ov]) if max_ov / te_len 0.5: filtered_tes.append(interval) pybedtools.BedTool(filtered_tes).saveas(tes_filtered.bed)参数说明intersect的wa和wb分别保留 A 和 B 的原始记录方便后续计算重叠长度。阈值 0.5 是个经验值如果你做的是胚胎发育早期TEs 表达活跃可以放宽到 0.7如果是体细胞组织建议收紧到 0.3。过滤后的 TEs 再转回 GTF和基因 GTF 合并用cat直接拼接即可注意去重。3. 多拷贝分配与 UMI 定量的源码实现3.1 多拷贝分配的 EM 算法核心逻辑TEs 量化的核心难点在于一条 read 比对到多个拷贝你怎么知道它来自哪一个常见做法是用 EM 算法迭代估计每个拷贝的表达比例。假设有 N 个 TEs 拷贝M 条 read每条 read 可能比对到多个拷贝EM 的步骤是先给每个拷贝一个初始表达值然后计算每条 read 来自每个拷贝的概率再根据概率重新估计表达值迭代到收敛。下面是一个简化版的 EM 实现输入是 read-拷贝比对矩阵和每个拷贝的有效长度import numpy as np def em_allocate(reads_by_te, te_lengths, max_iter100, tol1e-6): reads_by_te: list of lists, 每条 read 对应的候选 TE 索引 te_lengths: 每个 TE 的有效长度用于归一化 n_te len(te_lengths) # 初始化表达值用均匀分布 theta np.ones(n_te) / n_te # 有效长度归一化 eff_len np.array(te_lengths, dtypefloat) eff_len[eff_len 0] 1.0 for it in range(max_iter): # E 步计算每条 read 来自每个候选 TE 的概率 weights [] for candidates in reads_by_te: if len(candidates) 0: continue # 概率正比于 theta / eff_len probs np.array([theta[c] / eff_len[c] for c in candidates]) probs probs / probs.sum() weights.append((candidates, probs)) # M 步重新估计 theta new_theta np.zeros(n_te) for candidates, probs in weights: for c, p in zip(candidates, probs): new_theta[c] p new_theta new_theta / new_theta.sum() # 收敛判断 if np.abs(new_theta - theta).sum() tol: theta new_theta break theta new_theta return theta这段代码的关键参数是te_lengths它必须是有效长度而不是注释长度。有效长度等于注释长度减去 read 长度加 1如果小于 0 就设为 1。max_iter一般 100 次足够tol设 1e-6 在单细胞数据上比较稳。实际跑的时候reads_by_te会非常大建议用稀疏矩阵存储或者按细胞分批处理。3.2 UMI 去重与 TEs 定量的结合单细胞数据里每个细胞有独立的 UMI 标签。TEs 定量的正确做法是先按细胞和 UMI 分组对每组内的 reads 做多拷贝分配然后把分配结果累加到对应拷贝的表达值上。这里有个容易翻车的地方同一个 UMI 可能比对到多个拷贝如果你先做 UMI 去重再分配会丢失多拷贝信息如果先分配再去重又可能把不同拷贝的 reads 合并。我一般会先按 UMI 分组组内做 EM 分配然后把每个拷贝的期望计数作为该 UMI 对该拷贝的贡献。下面这段代码展示按细胞-UMI 分组后的定量流程from collections import defaultdict def quantify_tes_by_umi(bam_records, te_intervals): bam_records: 解析后的比对记录包含 cell_barcode, umi, ref_id te_intervals: TE 注释区间用于判断比对位置 # 按 cell umi 分组 umi_groups defaultdict(list) for rec in bam_records: key (rec[cell_barcode], rec[umi]) umi_groups[key].append(rec) # 每个细胞的 TE 表达向量 cell_te_counts defaultdict(lambda: defaultdict(float)) for (cell, umi), recs in umi_groups.items(): # 收集该 UMI 所有 reads 的候选 TE candidates set() for rec in recs: for te_id in rec[candidate_tes]: candidates.add(te_id) candidates list(candidates) if len(candidates) 0: continue if len(candidates) 1: # 唯一比对直接计数 cell_te_counts[cell][candidates[0]] 1.0 else: # 多拷贝用 EM 分配 reads_by_te [rec[candidate_tes] for rec in recs] te_lengths [te_intervals[c][length] for c in candidates] theta em_allocate(reads_by_te, te_lengths) for c, t in zip(candidates, theta): cell_te_counts[cell][c] t return cell_te_counts逻辑说明每个 UMI 组内如果只有唯一候选直接加 1如果有多个候选用 EM 算出的期望值累加。这样既保留了 UMI 去重的效果又不会丢失多拷贝信息。参数上要注意candidate_tes的生成——通常用比对工具输出的 NH 标签或者 MAPQ 过滤后的多位置比对记录。MAPQ 阈值建议设 10太低会引入大量错误比对太高会丢掉真实的多拷贝信号。3.3 源码目录结构与关键模块划分一套可复现的 TEs 量化源码目录结构应该清晰到别人能直接照着跑。我一般会这样组织te_quant/ ├── config/ │ └── params.yaml # 阈值、路径、线程数 ├── reference/ │ ├── build_te_gtf.py # RepeatMasker 转 GTF │ └── merge_annotation.py # 合并基因和 TEs 注释 ├── align/ │ └── run_star.sh # STAR 建索引和比对 ├── quantify/ │ ├── em_allocate.py # EM 多拷贝分配 │ ├── umi_count.py # UMI 分组定量 │ └── aggregate.py # 按 family/class 聚合 ├── utils/ │ ├── bam_parser.py # 解析 BAM提取 cell/UMI/候选 TE │ └── io.py # 读写矩阵 └── main.py # 主流程入口关键模块是bam_parser.py和em_allocate.py。bam_parser负责从 BAM 里提取每个 read 的 cell barcode、UMI 和候选 TEs候选 TEs 的获取依赖比对时的多位置输出。em_allocate是纯计算模块不依赖外部文件方便单独测试。main.py把流程串起来读params.yaml里的参数按细胞分批跑最后输出稀疏矩阵。参数文件params.yaml里我一般会放这些star: threads: 16 mapq_threshold: 10 multimap_max: 50 quantify: em_max_iter: 100 em_tol: 1e-6 min_umi_per_cell: 200 batch_size: 500 output: format: mtx aggregate_by: [family, class]mapq_threshold和multimap_max是最需要根据数据调整的。10x 数据一般 MAPQ 10 够用如果 TEs 表达特别低可以降到 5。multimap_max设 50 是防止某些 read 比对到太多拷贝导致 EM 太慢超过 50 的直接丢掉。4. 避坑与常见问题排查4.1 现象TEs 表达矩阵里大量细胞全为零原因单细胞数据本身稀疏TEs 又多拷贝如果 MAPQ 阈值设太高或者multimap_max设太小大量 reads 被过滤掉导致矩阵稀疏到没法用。解决先把 MAPQ 降到 5multimap_max提到 100看矩阵稀疏度是否改善。如果还是不行检查 UMI 去重逻辑——有些流程会把同一 UMI 的多条 reads 只保留一条这会直接丢掉多拷贝信号。正确做法是保留所有 reads在 EM 阶段处理。4.2 现象EM 迭代不收敛表达值震荡原因初始值给得太极端或者某些 TEs 的有效长度计算错误。解决初始值用均匀分布不要用基因表达值去初始化。有效长度检查一遍确保没有负数或零。如果还是震荡把tol放宽到 1e-4或者限制最大迭代次数到 50取最后一次结果。实际数据里震荡通常出现在 TEs 拷贝数特别多的家族可以考虑先按 family 聚合再跑 EM。4.3 现象同一细胞在不同批次间 TEs 表达差异巨大原因批次效应在 TEs 层面比基因层面更明显因为 TEs 表达受转座子调控网络影响不同建库批次的条件可能激活不同家族。解决在定量后做批次校正可以用 Harmony 或者 Seurat 的IntegrateLayers但注意 TEs 矩阵要单独做不要和基因矩阵混在一起。校正前先做 PCA看批次是否在 PC1/PC2 上分离。4.4 现象TEs 注释 GTF 和基因 GTF 合并后 STAR 建索引报错原因GTF 格式不标准比如属性字段缺少gene_id或transcript_id或者坐标超出染色体长度。解决用gtf_validator检查一遍确保每行都有gene_id和transcript_id。坐标超出的一般是 RepeatMasker 输出里的q_end大于染色体长度需要先和参考基因组的.fai文件比对截断或丢弃越界记录。4.5 现象跑完流程后 TEs 表达值全是整数没有期望计数原因EM 分配那一步被跳过了所有多拷贝 reads 都被当成唯一比对直接计数。解决检查bam_parser是否输出了candidate_tes列表如果列表长度总是 1说明比对时没有输出多位置信息。STAR 需要加--outMultimapperOrder Random和--outSAMattributes NH HI AS nMCell Ranger 需要改--nosecondary或者用 custom reference 时保留多位置比对。5. 从源码到产出验证 TEs 量化结果是否可信5.1 用已知 TEs 家族做阳性对照跑完流程后第一件事是拿几个已知在特定细胞类型里高表达的 TEs 家族做验证。比如在胚胎干细胞里LTR 类 TEs 通常活跃在神经元里LINE-1 家族有报道。你可以把量化结果按 family 聚合看这些家族的表达是否富集在预期细胞类型里。如果完全没信号说明流程有问题如果信号弥散在所有细胞里说明多拷贝分配太粗糙需要调 MAPQ 或 EM 参数。5.2 和 bulk RNA-seq 的 TEs 结果做一致性比较如果你有同一组织的 bulk RNA-seq 数据可以拿它跑一遍 TEs 量化然后和单细胞聚合后的伪 bulk 做相关性。相关系数在 0.6 以上算合理低于 0.4 说明单细胞流程里多拷贝分配或 UMI 去重有问题。注意 bulk 和单细胞的 TEs 注释要一致否则比较没意义。5.3 检查 TEs 表达和邻近基因表达的相关性TEs 插入位置往往影响邻近基因表达。你可以算每个 TEs 拷贝和它上下游 10kb 内基因的表达相关性如果大量 TEs 和邻近基因显著相关说明定量结果有生物学意义。如果全是随机相关可能是多拷贝分配把不同拷贝的信号混在一起了。这个验证不需要额外数据用你手头的单细胞矩阵就能做。5.4 一个具体技巧用唯一比对 reads 做校准多拷贝分配最怕的是把唯一比对 reads 也错误分配。我一般会留一部分唯一比对 reads 不参与 EM单独统计每个 TEs 拷贝的唯一 reads 数然后和 EM 分配后的总计数做比较。如果某个拷贝的唯一 reads 占比特别低说明它的表达几乎全靠多拷贝分配撑起来可信度低可以在下游分析里标记为低置信度。这个校准步骤不复杂但能帮你过滤掉一批不可信的 TEs 定量结果。我自己跑这类流程的血泪经验是参数不要一次调到位先拿一个细胞数少的小样本跑通确认 EM 收敛、UMI 分组正确、矩阵输出格式没问题再上全量数据。TEs 量化不像基因定量那样有成熟的金标准很多时候你得靠阳性对照和一致性检查来判断结果能不能用。希望帮到你。本文还有配套的精品资源点击获取 单细胞测序中转座元件表达量化：多拷贝分配与UMI定量源码实现
简介本资源为基于单细胞测序数据的转座元件TEs表达量化源码包面向生物信息学研究者与单细胞数据分析人员旨在解决单细胞水平上TEs表达难以精确量化的问题。包内共42个文件以21个Python脚本和6个Shell脚本承担数据分析与自动化流程辅以2个Jupyter Notebook交互式分析示例、2个Markdown文档说明并包含gtf、bed、BAM、idx等关键生物信息学格式文件压缩包约34.63MB。源码借鉴自JiekaiLab的scTE项目提供从数据准备、表达量化到结果展示的完整流程其中scTE可执行文件可直接处理BAM文件并自动调用相关脚本example目录则演示了具体分析用法便于初学者快速上手。目前已有327人学习下载适合希望深入探究TEs在细胞异质性与基因表达调控中作用的研究者参考使用。1. 单细胞测序数据里转座元件表达量化到底难在哪单细胞测序分析流程走到表达定量这一步绝大多数人只盯着基因把转座元件TEs当成基因组里的“暗物质”直接跳过。但如果你做过肿瘤微环境、早期胚胎发育或者神经退行性疾病的单细胞项目就会发现一个反直觉的现象同一群细胞在基因层面看起来高度一致TEs 表达谱却能把它们切成好几个亚群。这正是转座元件表达量化在单细胞测序里的价值——它提供了一层基因表达之外的调控信息。问题在于TEs 的量化跟普通基因完全不是一回事。TEs 在基因组里以多拷贝家族形式存在同一个家族的不同拷贝序列高度相似短读长测序的 reads 很容易在多拷贝之间发生错误分配。再加上单细胞数据本身 dropout 率高、每个细胞只有几千条 UMITEs 的表达信号就更稀疏。所以“基于单细胞测序数据的转座元件表达量化设计源码”这个标题本质上要解决的是怎么在单细胞分辨率下把 TEs 的多拷贝分配、UMI 去重和表达定量串成一条可复现的流程并且把这条流程写成能跑起来的源码。这套方案适合两类人一是正在做单细胞转录组分析、想从基因层面往下挖一层调控机制的从业者二是需要把 TEs 量化能力集成到自己分析管线里的生信工程师。下面我从数据准备、多拷贝分配、UMI 定量到源码结构把这条路径拆开讲清楚。2. 转座元件量化的数据准备与参考构建2.1 为什么不能直接用基因注释的 GTF 跑 TEs常规单细胞流程里比对完直接拿基因注释 GTF 做 featureCounts 或 STARsolo这套逻辑对 TEs 完全不适用。原因有三个第一TEs 注释不在标准基因 GTF 里你需要额外引入 RepeatMasker 或 Dfam 的注释文件第二TEs 的注释是区间形式很多拷贝互相重叠featureCounts 默认的“一个 read 只分配给一个 feature”策略会直接丢掉大量 ambiguous reads第三单细胞数据需要 UMI 级别的定量而 TEs 的多拷贝特性让 UMI 去重变得复杂——同一个 UMI 可能来自不同拷贝但序列一样你没法简单去重。常见做法是先构建一个包含基因和 TEs 的联合参考然后在比对阶段就允许 reads 多位置比对把多拷贝分配留到定量阶段用 EM 算法或唯一分子标签来解决。我一般会推荐用 STAR 建索引时把 TEs 区间也加进去或者用 Cell Ranger 的 custom reference 流程把 TEs 注释转成 GTF 格式后合并。2.2 构建 TEs 注释 GTF 的具体步骤假设你手头有 RepeatMasker 的.out文件或者 Dfam 的.hmm结果第一步是转成标准 GTF。下面这段 Python 脚本把 RepeatMasker 输出转成 GTF并给每个 TEs 拷贝分配唯一 IDimport pandas as pd # 读取 RepeatMasker .out 文件跳过前两行注释 # 列顺序SW score, perc div, perc del, perc ins, query seq, begin, end, ... rm pd.read_csv(repeatmasker.out, sepr\s, skiprows2, headerNone) rm.columns [ sw_score, perc_div, perc_del, perc_ins, query_seq, q_begin, q_end, q_left, strand, repeat_name, repeat_class, repeat_family, r_begin, r_end, r_left, id ] # 只保留有明确家族分类的记录过滤低复杂度 rm rm[rm[repeat_class].notna()] # 构建 GTF 行seqname, source, feature, start, end, score, strand, frame, attribute gtf_lines [] for idx, row in rm.iterrows(): # 每个拷贝用 repeat_name 坐标做唯一 ID避免多拷贝混淆 te_id f{row[repeat_name]}_{row[query_seq]}_{row[q_begin]}_{row[q_end]} # GTF 坐标从 1 开始且 end 是闭区间 start int(row[q_begin]) end int(row[q_end]) strand row[strand] if row[strand] in [, -] else attr fgene_id {te_id}; transcript_id {te_id}; family {row[repeat_family]}; class {row[repeat_class]}; gtf_lines.append( f{row[query_seq]}\tRepeatMasker\texon\t{start}\t{end}\t.\t{strand}\t.\t{attr} ) with open(tes_annotation.gtf, w) as f: f.write(\n.join(gtf_lines))这段代码的逻辑是把每个 TEs 拷贝当成一个独立的“基因”来处理用repeat_name 染色体 起止坐标拼成唯一 ID。这样在后续定量时每个拷贝有独立的表达值你可以再按 family 或 class 聚合。参数上要注意q_begin和q_end是相对于 query 序列的坐标如果 RepeatMasker 跑的是全基因组query_seq 就是染色体名坐标可以直接用。如果跑的是 contig需要先做 liftOver 或者用原始基因组坐标。2.3 合并基因和 TEs 注释时的冲突处理把 TEs GTF 和基因 GTF 合并时最大的坑是区间重叠。很多 TEs 就插在基因内含子或 UTR 里如果直接合并STAR 或 Cell Ranger 建索引时会报重复 feature 或者把 reads 错误分配。我一般会做两步过滤第一去掉与基因外显子重叠超过 50% 的 TEs 拷贝第二对于基因内 TEs单独标记为intronic_TE在定量时单独统计。下面这段代码用 pybedtools 做重叠过滤import pybedtools genes pybedtools.BedTool(genes.bed) tes pybedtools.BedTool(tes_annotation.bed) # 找与基因外显子重叠的 TEs overlaps tes.intersect(genes, waTrue, wbTrue) # 计算重叠比例过滤掉重叠超过 50% 的 filtered_tes [] for interval in tes: te_len interval.end - interval.start ov interval.intersect(genes, waTrue, wbTrue) if ov.count() 0: filtered_tes.append(interval) else: max_ov max([min(interval.end, o.end) - max(interval.start, o.start) for o in ov]) if max_ov / te_len 0.5: filtered_tes.append(interval) pybedtools.BedTool(filtered_tes).saveas(tes_filtered.bed)参数说明intersect的wa和wb分别保留 A 和 B 的原始记录方便后续计算重叠长度。阈值 0.5 是个经验值如果你做的是胚胎发育早期TEs 表达活跃可以放宽到 0.7如果是体细胞组织建议收紧到 0.3。过滤后的 TEs 再转回 GTF和基因 GTF 合并用cat直接拼接即可注意去重。3. 多拷贝分配与 UMI 定量的源码实现3.1 多拷贝分配的 EM 算法核心逻辑TEs 量化的核心难点在于一条 read 比对到多个拷贝你怎么知道它来自哪一个常见做法是用 EM 算法迭代估计每个拷贝的表达比例。假设有 N 个 TEs 拷贝M 条 read每条 read 可能比对到多个拷贝EM 的步骤是先给每个拷贝一个初始表达值然后计算每条 read 来自每个拷贝的概率再根据概率重新估计表达值迭代到收敛。下面是一个简化版的 EM 实现输入是 read-拷贝比对矩阵和每个拷贝的有效长度import numpy as np def em_allocate(reads_by_te, te_lengths, max_iter100, tol1e-6): reads_by_te: list of lists, 每条 read 对应的候选 TE 索引 te_lengths: 每个 TE 的有效长度用于归一化 n_te len(te_lengths) # 初始化表达值用均匀分布 theta np.ones(n_te) / n_te # 有效长度归一化 eff_len np.array(te_lengths, dtypefloat) eff_len[eff_len 0] 1.0 for it in range(max_iter): # E 步计算每条 read 来自每个候选 TE 的概率 weights [] for candidates in reads_by_te: if len(candidates) 0: continue # 概率正比于 theta / eff_len probs np.array([theta[c] / eff_len[c] for c in candidates]) probs probs / probs.sum() weights.append((candidates, probs)) # M 步重新估计 theta new_theta np.zeros(n_te) for candidates, probs in weights: for c, p in zip(candidates, probs): new_theta[c] p new_theta new_theta / new_theta.sum() # 收敛判断 if np.abs(new_theta - theta).sum() tol: theta new_theta break theta new_theta return theta这段代码的关键参数是te_lengths它必须是有效长度而不是注释长度。有效长度等于注释长度减去 read 长度加 1如果小于 0 就设为 1。max_iter一般 100 次足够tol设 1e-6 在单细胞数据上比较稳。实际跑的时候reads_by_te会非常大建议用稀疏矩阵存储或者按细胞分批处理。3.2 UMI 去重与 TEs 定量的结合单细胞数据里每个细胞有独立的 UMI 标签。TEs 定量的正确做法是先按细胞和 UMI 分组对每组内的 reads 做多拷贝分配然后把分配结果累加到对应拷贝的表达值上。这里有个容易翻车的地方同一个 UMI 可能比对到多个拷贝如果你先做 UMI 去重再分配会丢失多拷贝信息如果先分配再去重又可能把不同拷贝的 reads 合并。我一般会先按 UMI 分组组内做 EM 分配然后把每个拷贝的期望计数作为该 UMI 对该拷贝的贡献。下面这段代码展示按细胞-UMI 分组后的定量流程from collections import defaultdict def quantify_tes_by_umi(bam_records, te_intervals): bam_records: 解析后的比对记录包含 cell_barcode, umi, ref_id te_intervals: TE 注释区间用于判断比对位置 # 按 cell umi 分组 umi_groups defaultdict(list) for rec in bam_records: key (rec[cell_barcode], rec[umi]) umi_groups[key].append(rec) # 每个细胞的 TE 表达向量 cell_te_counts defaultdict(lambda: defaultdict(float)) for (cell, umi), recs in umi_groups.items(): # 收集该 UMI 所有 reads 的候选 TE candidates set() for rec in recs: for te_id in rec[candidate_tes]: candidates.add(te_id) candidates list(candidates) if len(candidates) 0: continue if len(candidates) 1: # 唯一比对直接计数 cell_te_counts[cell][candidates[0]] 1.0 else: # 多拷贝用 EM 分配 reads_by_te [rec[candidate_tes] for rec in recs] te_lengths [te_intervals[c][length] for c in candidates] theta em_allocate(reads_by_te, te_lengths) for c, t in zip(candidates, theta): cell_te_counts[cell][c] t return cell_te_counts逻辑说明每个 UMI 组内如果只有唯一候选直接加 1如果有多个候选用 EM 算出的期望值累加。这样既保留了 UMI 去重的效果又不会丢失多拷贝信息。参数上要注意candidate_tes的生成——通常用比对工具输出的 NH 标签或者 MAPQ 过滤后的多位置比对记录。MAPQ 阈值建议设 10太低会引入大量错误比对太高会丢掉真实的多拷贝信号。3.3 源码目录结构与关键模块划分一套可复现的 TEs 量化源码目录结构应该清晰到别人能直接照着跑。我一般会这样组织te_quant/ ├── config/ │ └── params.yaml # 阈值、路径、线程数 ├── reference/ │ ├── build_te_gtf.py # RepeatMasker 转 GTF │ └── merge_annotation.py # 合并基因和 TEs 注释 ├── align/ │ └── run_star.sh # STAR 建索引和比对 ├── quantify/ │ ├── em_allocate.py # EM 多拷贝分配 │ ├── umi_count.py # UMI 分组定量 │ └── aggregate.py # 按 family/class 聚合 ├── utils/ │ ├── bam_parser.py # 解析 BAM提取 cell/UMI/候选 TE │ └── io.py # 读写矩阵 └── main.py # 主流程入口关键模块是bam_parser.py和em_allocate.py。bam_parser负责从 BAM 里提取每个 read 的 cell barcode、UMI 和候选 TEs候选 TEs 的获取依赖比对时的多位置输出。em_allocate是纯计算模块不依赖外部文件方便单独测试。main.py把流程串起来读params.yaml里的参数按细胞分批跑最后输出稀疏矩阵。参数文件params.yaml里我一般会放这些star: threads: 16 mapq_threshold: 10 multimap_max: 50 quantify: em_max_iter: 100 em_tol: 1e-6 min_umi_per_cell: 200 batch_size: 500 output: format: mtx aggregate_by: [family, class]mapq_threshold和multimap_max是最需要根据数据调整的。10x 数据一般 MAPQ 10 够用如果 TEs 表达特别低可以降到 5。multimap_max设 50 是防止某些 read 比对到太多拷贝导致 EM 太慢超过 50 的直接丢掉。4. 避坑与常见问题排查4.1 现象TEs 表达矩阵里大量细胞全为零原因单细胞数据本身稀疏TEs 又多拷贝如果 MAPQ 阈值设太高或者multimap_max设太小大量 reads 被过滤掉导致矩阵稀疏到没法用。解决先把 MAPQ 降到 5multimap_max提到 100看矩阵稀疏度是否改善。如果还是不行检查 UMI 去重逻辑——有些流程会把同一 UMI 的多条 reads 只保留一条这会直接丢掉多拷贝信号。正确做法是保留所有 reads在 EM 阶段处理。4.2 现象EM 迭代不收敛表达值震荡原因初始值给得太极端或者某些 TEs 的有效长度计算错误。解决初始值用均匀分布不要用基因表达值去初始化。有效长度检查一遍确保没有负数或零。如果还是震荡把tol放宽到 1e-4或者限制最大迭代次数到 50取最后一次结果。实际数据里震荡通常出现在 TEs 拷贝数特别多的家族可以考虑先按 family 聚合再跑 EM。4.3 现象同一细胞在不同批次间 TEs 表达差异巨大原因批次效应在 TEs 层面比基因层面更明显因为 TEs 表达受转座子调控网络影响不同建库批次的条件可能激活不同家族。解决在定量后做批次校正可以用 Harmony 或者 Seurat 的IntegrateLayers但注意 TEs 矩阵要单独做不要和基因矩阵混在一起。校正前先做 PCA看批次是否在 PC1/PC2 上分离。4.4 现象TEs 注释 GTF 和基因 GTF 合并后 STAR 建索引报错原因GTF 格式不标准比如属性字段缺少gene_id或transcript_id或者坐标超出染色体长度。解决用gtf_validator检查一遍确保每行都有gene_id和transcript_id。坐标超出的一般是 RepeatMasker 输出里的q_end大于染色体长度需要先和参考基因组的.fai文件比对截断或丢弃越界记录。4.5 现象跑完流程后 TEs 表达值全是整数没有期望计数原因EM 分配那一步被跳过了所有多拷贝 reads 都被当成唯一比对直接计数。解决检查bam_parser是否输出了candidate_tes列表如果列表长度总是 1说明比对时没有输出多位置信息。STAR 需要加--outMultimapperOrder Random和--outSAMattributes NH HI AS nMCell Ranger 需要改--nosecondary或者用 custom reference 时保留多位置比对。5. 从源码到产出验证 TEs 量化结果是否可信5.1 用已知 TEs 家族做阳性对照跑完流程后第一件事是拿几个已知在特定细胞类型里高表达的 TEs 家族做验证。比如在胚胎干细胞里LTR 类 TEs 通常活跃在神经元里LINE-1 家族有报道。你可以把量化结果按 family 聚合看这些家族的表达是否富集在预期细胞类型里。如果完全没信号说明流程有问题如果信号弥散在所有细胞里说明多拷贝分配太粗糙需要调 MAPQ 或 EM 参数。5.2 和 bulk RNA-seq 的 TEs 结果做一致性比较如果你有同一组织的 bulk RNA-seq 数据可以拿它跑一遍 TEs 量化然后和单细胞聚合后的伪 bulk 做相关性。相关系数在 0.6 以上算合理低于 0.4 说明单细胞流程里多拷贝分配或 UMI 去重有问题。注意 bulk 和单细胞的 TEs 注释要一致否则比较没意义。5.3 检查 TEs 表达和邻近基因表达的相关性TEs 插入位置往往影响邻近基因表达。你可以算每个 TEs 拷贝和它上下游 10kb 内基因的表达相关性如果大量 TEs 和邻近基因显著相关说明定量结果有生物学意义。如果全是随机相关可能是多拷贝分配把不同拷贝的信号混在一起了。这个验证不需要额外数据用你手头的单细胞矩阵就能做。5.4 一个具体技巧用唯一比对 reads 做校准多拷贝分配最怕的是把唯一比对 reads 也错误分配。我一般会留一部分唯一比对 reads 不参与 EM单独统计每个 TEs 拷贝的唯一 reads 数然后和 EM 分配后的总计数做比较。如果某个拷贝的唯一 reads 占比特别低说明它的表达几乎全靠多拷贝分配撑起来可信度低可以在下游分析里标记为低置信度。这个校准步骤不复杂但能帮你过滤掉一批不可信的 TEs 定量结果。我自己跑这类流程的血泪经验是参数不要一次调到位先拿一个细胞数少的小样本跑通确认 EM 收敛、UMI 分组正确、矩阵输出格式没问题再上全量数据。TEs 量化不像基因定量那样有成熟的金标准很多时候你得靠阳性对照和一致性检查来判断结果能不能用。希望帮到你。本文还有配套的精品资源点击获取