拓冰网站定制 
简介本资源面向单细胞测序与生物信息学分析人员提供一套基于单细胞测序数据的转座元件TEs表达量化源码旨在解决单细胞水平上TEs表达难以精确量化的问题。资源包共42个文件约34.63MB包含21个Python脚本、6个Shell脚本、2个Jupyter Notebook、1个gtf基因组注释、1个bed文件、1个BAM文件及索引文件等覆盖数据准备、表达量化到结果展示的完整流程。其中可执行文件scTE可直接处理BAM文件并自动调用相关脚本example目录则演示了从聚类、标准化学习到差异表达与标记基因分析的具体用法便于初学者快速上手。目前已有327人学习下载适合具备一定单细胞分析基础、希望深入探究TEs在细胞异质性与基因表达调控中作用的研究者参考使用。1. 单细胞测序数据里被忽略的转座元件表达这套源码能解决什么做单细胞转录组分析的人十有八九把全部注意力放在蛋白编码基因上转座元件TEs的表达信号往往在比对阶段就被丢弃了。原因很直接TEs 在基因组里高度重复多比对 reads 一大堆常规流程要么直接扔掉要么用 featureCounts 的默认参数把它们过滤掉。但 TEs 恰恰是研究细胞异质性、早期胚胎发育、肿瘤微环境时越来越绕不开的一类调控元件。这套「基于单细胞测序数据的转座元件表达量化设计源码」解决的就是从原始比对结果出发把 TEs 的表达量在单细胞层面算出来这件事。它适合已经跑过 Cell Ranger 或 STARsolo、手里有 BAM 文件、想进一步挖 TEs 信号的生信从业者也适合做课程设计、需要一套完整可复现流程的学生。源码本身不依赖商业软件核心逻辑用 Python 和 R 串起来拿到就能改。2. 转座元件量化的技术底座为什么不能直接套用基因表达流程2.1 多比对 reads 的处理策略决定了 TEs 量化的成败普通基因表达量化只保留唯一比对unique mapping的 reads因为一个 read 只应该来自一个基因。但 TEs 不同同一个 read 可能同时匹配到多个 TE 拷贝如果直接丢弃TEs 的表达量会被严重低估。常见做法是保留多比对 reads然后用 EM 算法或最大似然估计把 reads 分配到各个 TE 拷贝上。这套源码里用的是「先保留多比对、再按拷贝数加权分配」的策略比直接丢弃多比对 reads 的召回率高出一大截。具体来说源码在 BAM 处理阶段设置了一个--keep-multi开关默认开启。开启后MAPQ 低于阈值的 reads 不会被直接扔掉而是进入一个待分配池。分配池里的 reads 会根据每个 TE 拷贝的有效长度和比对得分做加权最终把分数累加到对应的 TE 上。这个逻辑在te_quant/assign.py里实现核心是一个迭代加权的循环。注意如果你的数据来自 10x Genomics 的 3 或 5 建库UMI 信息必须保留否则重复 reads 会被重复计数TEs 表达量会虚高。2.2 从 BAM 到表达矩阵源码的模块划分与数据流源码整体分成三个模块预处理、量化、下游分析。预处理模块负责从 BAM 里提取比对信息、过滤低质量 reads、保留 UMI 和细胞条形码量化模块负责把 reads 分配到 TE 拷贝并汇总成细胞×TE 的表达矩阵下游分析模块提供了一些基础的降维聚类和差异表达函数方便快速验证结果。数据流的起点是 Cell Ranger 或 STARsolo 输出的possorted_genome_bam.bam终点是一个稀疏矩阵文件.mtx格式和一个对应的 TE 注释文件。中间过程全部用 Python 脚本串联不需要手动干预。如果你用的是其他比对工具只要 BAM 里包含标准的 CB细胞条形码和 UMI 标签也能接入。2.3 环境准备与依赖安装源码依赖的 Python 包不多但版本要卡准。pysam用于读 BAMscipy用于稀疏矩阵运算numpy和pandas做数据处理scanpy用于下游的可视化。R 那边主要用Matrix和Seurat做二次验证。建议用 conda 建一个独立环境避免和现有分析环境冲突。# 创建独立环境Python 版本建议 3.9 或 3.10 conda create -n te_quant python3.10 conda activate te_quant # 安装核心依赖pysam 对 htslib 版本敏感用 conda 装更稳 conda install -c bioconda pysam0.22 pip install numpy pandas scipy scanpy # R 侧依赖如果只跑 Python 部分可以跳过 conda install -c conda-forge r-base4.3 R -e install.packages(c(Matrix, Seurat), reposhttps://cloud.r-project.org)这里把pysam放在 conda 里装而不是 pip是因为 pip 版的pysam经常在编译时找不到htslib的头文件尤其是在没有 root 权限的服务器上。conda 版自带预编译的htslib省去很多麻烦。scanpy用 pip 装最新版即可它和pysam没有直接依赖冲突。3. 跑通量化流程从 BAM 到细胞×TE 表达矩阵3.1 准备 TE 注释文件格式要求与常见来源源码需要一个 BED 格式的 TE 注释文件每行至少包含染色体、起始位置、终止位置、TE 名称、拷贝编号。推荐用 RepeatMasker 的输出转成 BED或者直接下载 UCSC 的rmsk.txt自己转。转换脚本源码里带了一个utils/rmsk_to_bed.py可以直接用。# utils/rmsk_to_bed.py 的核心逻辑 import pandas as pd def rmsk_to_bed(rmsk_file, output_bed): # rmsk.txt 是制表符分隔列名固定 cols [bin, swScore, milliDiv, milliDel, milliIns, genoName, genoStart, genoEnd, genoLeft, strand, repName, repClass, repFamily, repStart, repEnd, repLeft, id] df pd.read_csv(rmsk_file, sep\t, namescols, comment#) # BED 需要 0-based 起始位置rmsk 是 1-based df[genoStart] df[genoStart] - 1 # 只保留有明确家族分类的 TE df df[df[repClass].notna()] # 输出六列标准 BED bed df[[genoName, genoStart, genoEnd, repName, repClass, strand]] bed.to_csv(output_bed, sep\t, headerFalse, indexFalse) if __name__ __main__: rmsk_to_bed(rmsk.txt, te_annotations.bed)这段代码做了三件事读入 RepeatMasker 的原始输出、把坐标从 1-based 转成 0-based、过滤掉没有家族分类的条目。repClass这一列很关键后面做 TE 家族层面的汇总时会用到。如果你的注释文件来自其他来源只要保证有染色体、起始、终止、名称、链方向这五列就能直接替换。3.2 运行主量化脚本参数含义与输出解读主脚本是te_quant/quantify.py调用方式如下python te_quant/quantify.py \ --bam possorted_genome_bam.bam \ --bed te_annotations.bed \ --output te_matrix \ --keep-multi \ --min-mapq 10 \ --threads 8参数逐个说清楚--bam指定输入 BAM必须是包含 CB 和 UMI 标签的--bed是上一步生成的 TE 注释--output是输出目录会生成matrix.mtx、barcodes.tsv、features.tsv三个文件--keep-multi决定是否保留多比对 reads默认开启如果你的数据多比对率特别高超过 50%可以关掉试试对比--min-mapq是比对质量阈值低于这个值的 reads 不进入分配池默认 10对于 TEs 来说可以适当放宽到 5--threads是并行线程数建议设为 CPU 核数的 70% 左右。跑完之后输出目录里会有三个文件。matrix.mtx是稀疏矩阵行是 TE列是细胞barcodes.tsv是细胞条形码列表features.tsv是 TE 名称和家族信息。用scanpy读进来就能直接做下游分析。3.3 用 scanpy 快速验证量化结果拿到矩阵后先做个基础的质控和聚类看看 TEs 的表达能不能区分出已知的细胞类型。这一步不是必须的但能帮你判断量化结果是否合理。import scanpy as sc # 读入源码输出的矩阵 adata sc.read_mtx(te_matrix/matrix.mtx).T adata.obs_names [l.strip() for l in open(te_matrix/barcodes.tsv)] adata.var_names [l.strip().split(\t)[0] for l in open(te_matrix/features.tsv)] # 基础质控过滤掉表达 TE 数太少的细胞 sc.pp.filter_cells(adata, min_genes50) sc.pp.filter_genes(adata, min_cells10) # 归一化和降维 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes2000) sc.tl.pca(adata, n_comps30) sc.pp.neighbors(adata, n_neighbors15) sc.tl.leiden(adata, resolution0.8) sc.tl.umap(adata) # 可视化 sc.pl.umap(adata, color[leiden], save_te_clusters.png)这里用.T转置是因为源码输出的矩阵是 TE×细胞而scanpy默认要求细胞×基因。min_genes50这个阈值比常规 RNA 分析低很多因为 TEs 的表达丰度普遍低于蛋白编码基因阈值设太高会把真实细胞过滤掉。n_top_genes2000也是类似考虑TEs 的高变基因数量通常比 mRNA 少。跑完 UMAP 后如果能看到清晰的聚类结构说明量化结果至少没有大的偏差。4. 避坑与排查TEs 量化中五个血泪教训4.1 现象表达矩阵里大量细胞全为零原因通常是细胞条形码不匹配。Cell Ranger 输出的 BAM 里CB 标签是CB:Z:xxxxx格式但有些版本的 STARsolo 用的是CR:Z:xxxxx。源码默认读CB如果标签不对所有 reads 都会被判为无效。解决办法是在quantify.py里加一个--cb-tag参数指定正确的标签名。跑之前先用samtools view看一眼 BAM 的标签格式确认无误再往下走。4.2 现象TEs 表达量整体偏高和已知生物学预期不符多比对 reads 的分配权重没设对。源码默认用的是「比对得分加权」但如果你的 BAM 里所有多比对 reads 的得分都是 0有些比对工具不输出 AS 标签权重就退化成均匀分配导致每个 TE 拷贝都分到相同的 reads表达量虚高。解决办法是改用「有效长度加权」在assign.py里把weight_method改成effective_length。这个参数在配置文件里就能改不用动代码。4.3 现象跑完量化后下游聚类完全分不开细胞类型先检查是不是把 TE 家族和 TE 拷贝搞混了。源码默认输出的是拷贝级别的矩阵同一个家族的不同拷贝是分开的。如果你的生物学问题关注的是家族层面的变化需要在features.tsv里按repClass汇总。源码提供了一个utils/collapse_to_family.py脚本跑一下就能把拷贝矩阵合并成家族矩阵。家族矩阵的稀疏性更低聚类效果通常更好。4.4 现象量化过程内存溢出被系统 killTEs 的多比对 reads 数量可能非常大尤其是植物或大型哺乳动物基因组。源码默认把所有待分配 reads 读进内存数据量大时确实会爆。解决办法是加--batch-size参数分批处理。每批处理 100 万条 reads处理完一批写一次中间结果最后再合并。这个参数在quantify.py里已经预留了默认是 0不分批手动设成1000000就行。4.5 现象和 bulk RNA-seq 的 TEs 结果对不上单细胞和 bulk 的 TEs 表达量本来就不应该完全一致但如果是趋势相反大概率是 UMI 去重没做好。单细胞数据里UMI 相同的 reads 应该只算一次但如果 BAM 里的 UMI 标签有多个比如UB和UB同时存在源码可能会重复计数。检查方法是随机抽一个细胞数一下它的总 UMI 数和总 reads 数如果比值明显偏低就是去重出了问题。在quantify.py里加--umi-tag UB明确指定标签即可。5. 进阶技巧用 TEs 表达矩阵做细胞类型注释的验证5.1 把 TEs 信号和已知 marker 基因做相关性分析TEs 本身不是经典的细胞类型 marker但某些 TE 家族在特定细胞类型里会有特异性表达。一个实用的技巧是先用 mRNA 数据做好细胞类型注释然后把注释标签映射到 TEs 矩阵上看哪些 TE 在特定细胞类型里显著富集。源码的下游模块里有一个te_marker_score函数输入是 TEs 矩阵和细胞类型标签输出是每个 TE 在每个细胞类型里的富集得分。from te_quant.downstream import te_marker_score # adata_te 是 TEs 矩阵adata_rna 是 mRNA 矩阵 # cell_type_key 是 mRNA 注释里的细胞类型列名 scores te_marker_score(adata_te, adata_rna, cell_type_keycell_type) # 筛选在某个细胞类型里得分显著高的 TE # 比如看 Excitatory 神经元里富集的 TE ex_te scores[scores[cell_type] Excitatory].sort_values(score, ascendingFalse) print(ex_te.head(20))这个函数的逻辑是对每个 TE计算它在目标细胞类型和其他细胞类型之间的表达差异然后用 Wilcoxon 秩和检验算显著性。score列是差异倍数和显著性的综合得分排在前面的就是候选 marker TE。我一般会取前 20 个然后手动检查一下这些 TE 的家族分类看看有没有已知的生物学关联。5.2 用 TEs 做批次效应评估单细胞数据整合时批次效应是个绕不开的问题。常规做法是看 mRNA 的整合效果但 mRNA 有时候会掩盖批次效应。TEs 的表达模式对批次更敏感可以用来做辅助评估。具体做法是在整合前后分别计算 TEs 矩阵的批次混合熵如果整合后熵值明显下降说明批次效应被校正了如果没降反升说明整合参数可能过拟合了。源码里提供了一个batch_entropy函数输入是 TEs 矩阵和批次标签输出是每个批次的混合熵。这个指标不是绝对的但可以作为 mRNA 评估的补充。我自己的习惯是mRNA 的整合指标和 TEs 的混合熵都看一遍两者趋势一致才认为整合是可靠的。5.3 一个容易忽略的细节TE 注释的版本一致性最后说一个我踩过的坑。TE 注释文件一定要和参考基因组的版本匹配。比如你用 hg38 的 BAM就必须用 hg38 的 RepeatMasker 注释。如果用 hg19 的注释坐标全错量化结果完全没有意义。更隐蔽的是有些注释文件虽然基因组版本对但 RepeatMasker 的版本不同TE 家族的命名会有差异。我现在的习惯是每次跑新数据之前先用bedtools intersect检查一下注释文件和 BAM 的染色体命名是否一致再随机抽几个 TE 看看坐标能不能对上。这个检查花不了两分钟但能省掉后面几小时的排查。从那以后我每次拿到新的 BAM 和注释文件都强制走一遍坐标一致性检查确认无误再开始量化。希望帮到你。本文还有配套的精品资源点击获取 单细胞测序中转座元件表达量化源码：从BAM到细胞×TE矩阵
简介本资源面向单细胞测序与生物信息学分析人员提供一套基于单细胞测序数据的转座元件TEs表达量化源码旨在解决单细胞水平上TEs表达难以精确量化的问题。资源包共42个文件约34.63MB包含21个Python脚本、6个Shell脚本、2个Jupyter Notebook、1个gtf基因组注释、1个bed文件、1个BAM文件及索引文件等覆盖数据准备、表达量化到结果展示的完整流程。其中可执行文件scTE可直接处理BAM文件并自动调用相关脚本example目录则演示了从聚类、标准化学习到差异表达与标记基因分析的具体用法便于初学者快速上手。目前已有327人学习下载适合具备一定单细胞分析基础、希望深入探究TEs在细胞异质性与基因表达调控中作用的研究者参考使用。1. 单细胞测序数据里被忽略的转座元件表达这套源码能解决什么做单细胞转录组分析的人十有八九把全部注意力放在蛋白编码基因上转座元件TEs的表达信号往往在比对阶段就被丢弃了。原因很直接TEs 在基因组里高度重复多比对 reads 一大堆常规流程要么直接扔掉要么用 featureCounts 的默认参数把它们过滤掉。但 TEs 恰恰是研究细胞异质性、早期胚胎发育、肿瘤微环境时越来越绕不开的一类调控元件。这套「基于单细胞测序数据的转座元件表达量化设计源码」解决的就是从原始比对结果出发把 TEs 的表达量在单细胞层面算出来这件事。它适合已经跑过 Cell Ranger 或 STARsolo、手里有 BAM 文件、想进一步挖 TEs 信号的生信从业者也适合做课程设计、需要一套完整可复现流程的学生。源码本身不依赖商业软件核心逻辑用 Python 和 R 串起来拿到就能改。2. 转座元件量化的技术底座为什么不能直接套用基因表达流程2.1 多比对 reads 的处理策略决定了 TEs 量化的成败普通基因表达量化只保留唯一比对unique mapping的 reads因为一个 read 只应该来自一个基因。但 TEs 不同同一个 read 可能同时匹配到多个 TE 拷贝如果直接丢弃TEs 的表达量会被严重低估。常见做法是保留多比对 reads然后用 EM 算法或最大似然估计把 reads 分配到各个 TE 拷贝上。这套源码里用的是「先保留多比对、再按拷贝数加权分配」的策略比直接丢弃多比对 reads 的召回率高出一大截。具体来说源码在 BAM 处理阶段设置了一个--keep-multi开关默认开启。开启后MAPQ 低于阈值的 reads 不会被直接扔掉而是进入一个待分配池。分配池里的 reads 会根据每个 TE 拷贝的有效长度和比对得分做加权最终把分数累加到对应的 TE 上。这个逻辑在te_quant/assign.py里实现核心是一个迭代加权的循环。注意如果你的数据来自 10x Genomics 的 3 或 5 建库UMI 信息必须保留否则重复 reads 会被重复计数TEs 表达量会虚高。2.2 从 BAM 到表达矩阵源码的模块划分与数据流源码整体分成三个模块预处理、量化、下游分析。预处理模块负责从 BAM 里提取比对信息、过滤低质量 reads、保留 UMI 和细胞条形码量化模块负责把 reads 分配到 TE 拷贝并汇总成细胞×TE 的表达矩阵下游分析模块提供了一些基础的降维聚类和差异表达函数方便快速验证结果。数据流的起点是 Cell Ranger 或 STARsolo 输出的possorted_genome_bam.bam终点是一个稀疏矩阵文件.mtx格式和一个对应的 TE 注释文件。中间过程全部用 Python 脚本串联不需要手动干预。如果你用的是其他比对工具只要 BAM 里包含标准的 CB细胞条形码和 UMI 标签也能接入。2.3 环境准备与依赖安装源码依赖的 Python 包不多但版本要卡准。pysam用于读 BAMscipy用于稀疏矩阵运算numpy和pandas做数据处理scanpy用于下游的可视化。R 那边主要用Matrix和Seurat做二次验证。建议用 conda 建一个独立环境避免和现有分析环境冲突。# 创建独立环境Python 版本建议 3.9 或 3.10 conda create -n te_quant python3.10 conda activate te_quant # 安装核心依赖pysam 对 htslib 版本敏感用 conda 装更稳 conda install -c bioconda pysam0.22 pip install numpy pandas scipy scanpy # R 侧依赖如果只跑 Python 部分可以跳过 conda install -c conda-forge r-base4.3 R -e install.packages(c(Matrix, Seurat), reposhttps://cloud.r-project.org)这里把pysam放在 conda 里装而不是 pip是因为 pip 版的pysam经常在编译时找不到htslib的头文件尤其是在没有 root 权限的服务器上。conda 版自带预编译的htslib省去很多麻烦。scanpy用 pip 装最新版即可它和pysam没有直接依赖冲突。3. 跑通量化流程从 BAM 到细胞×TE 表达矩阵3.1 准备 TE 注释文件格式要求与常见来源源码需要一个 BED 格式的 TE 注释文件每行至少包含染色体、起始位置、终止位置、TE 名称、拷贝编号。推荐用 RepeatMasker 的输出转成 BED或者直接下载 UCSC 的rmsk.txt自己转。转换脚本源码里带了一个utils/rmsk_to_bed.py可以直接用。# utils/rmsk_to_bed.py 的核心逻辑 import pandas as pd def rmsk_to_bed(rmsk_file, output_bed): # rmsk.txt 是制表符分隔列名固定 cols [bin, swScore, milliDiv, milliDel, milliIns, genoName, genoStart, genoEnd, genoLeft, strand, repName, repClass, repFamily, repStart, repEnd, repLeft, id] df pd.read_csv(rmsk_file, sep\t, namescols, comment#) # BED 需要 0-based 起始位置rmsk 是 1-based df[genoStart] df[genoStart] - 1 # 只保留有明确家族分类的 TE df df[df[repClass].notna()] # 输出六列标准 BED bed df[[genoName, genoStart, genoEnd, repName, repClass, strand]] bed.to_csv(output_bed, sep\t, headerFalse, indexFalse) if __name__ __main__: rmsk_to_bed(rmsk.txt, te_annotations.bed)这段代码做了三件事读入 RepeatMasker 的原始输出、把坐标从 1-based 转成 0-based、过滤掉没有家族分类的条目。repClass这一列很关键后面做 TE 家族层面的汇总时会用到。如果你的注释文件来自其他来源只要保证有染色体、起始、终止、名称、链方向这五列就能直接替换。3.2 运行主量化脚本参数含义与输出解读主脚本是te_quant/quantify.py调用方式如下python te_quant/quantify.py \ --bam possorted_genome_bam.bam \ --bed te_annotations.bed \ --output te_matrix \ --keep-multi \ --min-mapq 10 \ --threads 8参数逐个说清楚--bam指定输入 BAM必须是包含 CB 和 UMI 标签的--bed是上一步生成的 TE 注释--output是输出目录会生成matrix.mtx、barcodes.tsv、features.tsv三个文件--keep-multi决定是否保留多比对 reads默认开启如果你的数据多比对率特别高超过 50%可以关掉试试对比--min-mapq是比对质量阈值低于这个值的 reads 不进入分配池默认 10对于 TEs 来说可以适当放宽到 5--threads是并行线程数建议设为 CPU 核数的 70% 左右。跑完之后输出目录里会有三个文件。matrix.mtx是稀疏矩阵行是 TE列是细胞barcodes.tsv是细胞条形码列表features.tsv是 TE 名称和家族信息。用scanpy读进来就能直接做下游分析。3.3 用 scanpy 快速验证量化结果拿到矩阵后先做个基础的质控和聚类看看 TEs 的表达能不能区分出已知的细胞类型。这一步不是必须的但能帮你判断量化结果是否合理。import scanpy as sc # 读入源码输出的矩阵 adata sc.read_mtx(te_matrix/matrix.mtx).T adata.obs_names [l.strip() for l in open(te_matrix/barcodes.tsv)] adata.var_names [l.strip().split(\t)[0] for l in open(te_matrix/features.tsv)] # 基础质控过滤掉表达 TE 数太少的细胞 sc.pp.filter_cells(adata, min_genes50) sc.pp.filter_genes(adata, min_cells10) # 归一化和降维 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes2000) sc.tl.pca(adata, n_comps30) sc.pp.neighbors(adata, n_neighbors15) sc.tl.leiden(adata, resolution0.8) sc.tl.umap(adata) # 可视化 sc.pl.umap(adata, color[leiden], save_te_clusters.png)这里用.T转置是因为源码输出的矩阵是 TE×细胞而scanpy默认要求细胞×基因。min_genes50这个阈值比常规 RNA 分析低很多因为 TEs 的表达丰度普遍低于蛋白编码基因阈值设太高会把真实细胞过滤掉。n_top_genes2000也是类似考虑TEs 的高变基因数量通常比 mRNA 少。跑完 UMAP 后如果能看到清晰的聚类结构说明量化结果至少没有大的偏差。4. 避坑与排查TEs 量化中五个血泪教训4.1 现象表达矩阵里大量细胞全为零原因通常是细胞条形码不匹配。Cell Ranger 输出的 BAM 里CB 标签是CB:Z:xxxxx格式但有些版本的 STARsolo 用的是CR:Z:xxxxx。源码默认读CB如果标签不对所有 reads 都会被判为无效。解决办法是在quantify.py里加一个--cb-tag参数指定正确的标签名。跑之前先用samtools view看一眼 BAM 的标签格式确认无误再往下走。4.2 现象TEs 表达量整体偏高和已知生物学预期不符多比对 reads 的分配权重没设对。源码默认用的是「比对得分加权」但如果你的 BAM 里所有多比对 reads 的得分都是 0有些比对工具不输出 AS 标签权重就退化成均匀分配导致每个 TE 拷贝都分到相同的 reads表达量虚高。解决办法是改用「有效长度加权」在assign.py里把weight_method改成effective_length。这个参数在配置文件里就能改不用动代码。4.3 现象跑完量化后下游聚类完全分不开细胞类型先检查是不是把 TE 家族和 TE 拷贝搞混了。源码默认输出的是拷贝级别的矩阵同一个家族的不同拷贝是分开的。如果你的生物学问题关注的是家族层面的变化需要在features.tsv里按repClass汇总。源码提供了一个utils/collapse_to_family.py脚本跑一下就能把拷贝矩阵合并成家族矩阵。家族矩阵的稀疏性更低聚类效果通常更好。4.4 现象量化过程内存溢出被系统 killTEs 的多比对 reads 数量可能非常大尤其是植物或大型哺乳动物基因组。源码默认把所有待分配 reads 读进内存数据量大时确实会爆。解决办法是加--batch-size参数分批处理。每批处理 100 万条 reads处理完一批写一次中间结果最后再合并。这个参数在quantify.py里已经预留了默认是 0不分批手动设成1000000就行。4.5 现象和 bulk RNA-seq 的 TEs 结果对不上单细胞和 bulk 的 TEs 表达量本来就不应该完全一致但如果是趋势相反大概率是 UMI 去重没做好。单细胞数据里UMI 相同的 reads 应该只算一次但如果 BAM 里的 UMI 标签有多个比如UB和UB同时存在源码可能会重复计数。检查方法是随机抽一个细胞数一下它的总 UMI 数和总 reads 数如果比值明显偏低就是去重出了问题。在quantify.py里加--umi-tag UB明确指定标签即可。5. 进阶技巧用 TEs 表达矩阵做细胞类型注释的验证5.1 把 TEs 信号和已知 marker 基因做相关性分析TEs 本身不是经典的细胞类型 marker但某些 TE 家族在特定细胞类型里会有特异性表达。一个实用的技巧是先用 mRNA 数据做好细胞类型注释然后把注释标签映射到 TEs 矩阵上看哪些 TE 在特定细胞类型里显著富集。源码的下游模块里有一个te_marker_score函数输入是 TEs 矩阵和细胞类型标签输出是每个 TE 在每个细胞类型里的富集得分。from te_quant.downstream import te_marker_score # adata_te 是 TEs 矩阵adata_rna 是 mRNA 矩阵 # cell_type_key 是 mRNA 注释里的细胞类型列名 scores te_marker_score(adata_te, adata_rna, cell_type_keycell_type) # 筛选在某个细胞类型里得分显著高的 TE # 比如看 Excitatory 神经元里富集的 TE ex_te scores[scores[cell_type] Excitatory].sort_values(score, ascendingFalse) print(ex_te.head(20))这个函数的逻辑是对每个 TE计算它在目标细胞类型和其他细胞类型之间的表达差异然后用 Wilcoxon 秩和检验算显著性。score列是差异倍数和显著性的综合得分排在前面的就是候选 marker TE。我一般会取前 20 个然后手动检查一下这些 TE 的家族分类看看有没有已知的生物学关联。5.2 用 TEs 做批次效应评估单细胞数据整合时批次效应是个绕不开的问题。常规做法是看 mRNA 的整合效果但 mRNA 有时候会掩盖批次效应。TEs 的表达模式对批次更敏感可以用来做辅助评估。具体做法是在整合前后分别计算 TEs 矩阵的批次混合熵如果整合后熵值明显下降说明批次效应被校正了如果没降反升说明整合参数可能过拟合了。源码里提供了一个batch_entropy函数输入是 TEs 矩阵和批次标签输出是每个批次的混合熵。这个指标不是绝对的但可以作为 mRNA 评估的补充。我自己的习惯是mRNA 的整合指标和 TEs 的混合熵都看一遍两者趋势一致才认为整合是可靠的。5.3 一个容易忽略的细节TE 注释的版本一致性最后说一个我踩过的坑。TE 注释文件一定要和参考基因组的版本匹配。比如你用 hg38 的 BAM就必须用 hg38 的 RepeatMasker 注释。如果用 hg19 的注释坐标全错量化结果完全没有意义。更隐蔽的是有些注释文件虽然基因组版本对但 RepeatMasker 的版本不同TE 家族的命名会有差异。我现在的习惯是每次跑新数据之前先用bedtools intersect检查一下注释文件和 BAM 的染色体命名是否一致再随机抽几个 TE 看看坐标能不能对上。这个检查花不了两分钟但能省掉后面几小时的排查。从那以后我每次拿到新的 BAM 和注释文件都强制走一遍坐标一致性检查确认无误再开始量化。希望帮到你。本文还有配套的精品资源点击获取