企业官网建设流程全解析

简介：本资源面向单细胞测序与生物信息学分析人员，提供一套基于单细胞测序数据的转座元件（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用于读 BAM，scipy用于稀疏矩阵运算，numpy和pandas做数据处理，scanpy用于下游的可视化。R 那边主要用Matrix和Seurat做二次验证。建议用 conda 建一个独立环境，避免和现有分析环境冲突。

# 创建独立环境，Python 版本建议 3.9 或 3.10 conda create -n te_quant python=3.10 conda activate te_quant # 安装核心依赖，pysam 对 htslib 版本敏感，用 conda 装更稳 conda install -c bioconda pysam=0.22 pip install numpy pandas scipy scanpy # R 侧依赖，如果只跑 Python 部分可以跳过 conda install -c conda-forge r-base=4.3 R -e "install.packages(c('Matrix', 'Seurat'), repos='https://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', names=cols, 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', header=False, index=False) 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_genes=50) sc.pp.filter_genes(adata, min_cells=10) # 归一化和降维 sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes=2000) sc.tl.pca(adata, n_comps=30) sc.pp.neighbors(adata, n_neighbors=15) sc.tl.leiden(adata, resolution=0.8) sc.tl.umap(adata) # 可视化 sc.pl.umap(adata, color=['leiden'], save='_te_clusters.png')

这里用.T转置是因为源码输出的矩阵是 TE×细胞，而scanpy默认要求细胞×基因。min_genes=50这个阈值比常规 RNA 分析低很多，因为 TEs 的表达丰度普遍低于蛋白编码基因，阈值设太高会把真实细胞过滤掉。n_top_genes=2000也是类似考虑，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 现象：量化过程内存溢出，被系统 kill

TEs 的多比对 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_key='cell_type') # 筛选在某个细胞类型里得分显著高的 TE # 比如看 Excitatory 神经元里富集的 TE ex_te = scores[scores['cell_type'] == 'Excitatory'].sort_values('score', ascending=False) 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 和注释文件，都强制走一遍坐标一致性检查，确认无误再开始量化。希望帮到你。

本文还有配套的精品资源，点击获取