629_转座子插入分析¶

一句话概述： 转座子（TE）是基因组中能"跳跃"的DNA片段，占人类基因组的45%以上。TEtranscripts是分析TE转录活性的核心工具（能正确处理multi-mapped reads），RepeatMasker是TE注释的标准工具——两者结合可以从注释到表达定量完成TE分析全流程。

核心知识点速查表¶

概念	白话解释
Transposable Element (TE)	转座子——能在基因组中"跳跃"复制的DNA序列
Retrotransposon	反转录转座子——通过"复制粘贴"机制移动（如LINE、SINE、LTR）
DNA transposon	DNA转座子——通过"剪切粘贴"机制移动
LINE	长散在核元件——人类基因组中最多的TE（占17%）
SINE	短散在核元件——如Alu元件（占11%）
LTR	长末端重复——内源性逆转录病毒（如HERV）
RepeatMasker	TE注释的"金标准"工具
TEtranscripts	能正确处理多比对reads的TE表达定量工具
Multi-mapped reads	多比对reads——比对到基因组多个位置的reads
Repbase/Dfam	TE序列参考数据库

一、安装¶

# === RepeatMasker安装 ===
conda install -c bioconda repeatmasker  # conda安装
RepeatMasker -v  # 查看版本

# === TEtranscripts安装 ===
pip install TEtranscripts  # pip安装
TEcount --version  # 查看版本

# === 辅助工具 ===
conda install -c bioconda star        # STAR比对器（TEtranscripts推荐）
conda install -c bioconda samtools    # BAM处理
conda install -c bioconda bedtools    # BED操作

# === RepeatModeler（从头发现TE，可选）===
conda install -c bioconda repeatmodeler  # 从头TE建模

二、RepeatMasker TE注释¶

2.1 基本注释¶

# === RepeatMasker注释基因组中的TE ===
RepeatMasker \
    -species human \               # 物种（human/mouse/arabidopsis等）
    -pa 16 \                       # 16线程
    -gff \                         # 输出GFF格式注释
    -xsmall \                      # 用小写字母标记重复（而非N）
    -dir repeatmasker_output/ \    # 输出目录
    genome.fasta                   # 输入基因组

# === 输出文件说明 ===
# genome.fasta.out      — 详细TE注释（每个TE的位置、类型、得分）
# genome.fasta.masked   — 重复序列被mask后的基因组
# genome.fasta.tbl      — 统计表（各类TE的数量和占比）
# genome.fasta.out.gff  — GFF格式注释

# 查看TE统计
cat repeatmasker_output/genome.fasta.tbl

2.2 指定TE库¶

# === 使用自定义TE库 ===
# 从Dfam数据库下载（2026年最新）
# https://www.dfam.org/

# 使用自定义TE库注释
RepeatMasker \
    -lib custom_te_library.fa \    # 自定义TE序列库
    -pa 16 \
    -gff \
    -xsmall \
    -dir rm_custom/ \
    genome.fasta

# === 使用RepeatModeler从头发现TE（新物种）===
# 第1步：建数据库
BuildDatabase -name genome_db -engine ncbi genome.fasta

# 第2步：从头建模
RepeatModeler \
    -database genome_db \          # 基因组数据库
    -pa 16 \                       # 16线程
    -LTRStruct                     # 启用LTR结构搜索

# 输出：genome_db-families.fa（发现的TE一致性序列）
# 然后用这个库做RepeatMasker注释

2.3 解析RepeatMasker结果¶

# === 解析RepeatMasker输出 ===
import pandas as pd  # 数据处理
import matplotlib.pyplot as plt  # 画图

# 读取RepeatMasker .out文件
def parse_rm_out(filepath):
    """解析RepeatMasker .out文件"""
    records = []  # 存储记录
    with open(filepath) as f:
        for i, line in enumerate(f):
            if i < 3:  # 跳过头3行
                continue
            parts = line.split()  # 分割字段
            if len(parts) >= 15:  # 有效行
                records.append({
                    'score': int(parts[0]),        # 得分
                    'chr': parts[4],               # 染色体
                    'start': int(parts[5]),         # 起始位置
                    'end': int(parts[6]),           # 结束位置
                    'repeat_name': parts[9],        # TE名称
                    'repeat_class': parts[10],      # TE类/家族
                    'length': int(parts[6]) - int(parts[5])  # 长度
                })
    return pd.DataFrame(records)  # 返回数据框

# 解析
rm_data = parse_rm_out("repeatmasker_output/genome.fasta.out")

# 统计各类TE
class_stats = rm_data.groupby('repeat_class').agg(
    count=('repeat_name', 'count'),       # 数量
    total_bp=('length', 'sum')            # 总碱基数
).sort_values('total_bp', ascending=False)

print("=== TE类别统计 ===")
print(class_stats.head(20))

# 画饼图
top_classes = class_stats.head(8)  # 取前8类
plt.figure(figsize=(10, 8))
plt.pie(top_classes['total_bp'], 
        labels=top_classes.index,
        autopct='%1.1f%%',                # 显示百分比
        startangle=90)
plt.title("Transposable Element Composition")
plt.savefig("te_composition.pdf")
plt.close()

三、TEtranscripts TE表达定量¶

3.1 STAR比对（保留multi-mapped reads）¶

# === 用STAR比对，必须保留multi-mapped reads ===
# TEtranscripts需要multi-mapped reads信息来正确分配TE表达量

# 第1步：建STAR索引
STAR --runMode genomeGenerate \
    --genomeDir star_index/ \          # 索引输出目录
    --genomeFastaFiles genome.fasta \  # 参考基因组
    --sjdbGTFfile genes.gtf \          # 基因注释
    --runThreadN 16                    # 16线程

# 第2步：STAR比对（关键参数！）
STAR \
    --genomeDir star_index/ \          # 索引目录
    --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \  # 输入reads
    --readFilesCommand zcat \          # 解压命令
    --outSAMtype BAM SortedByCoordinate \  # 输出排序BAM
    --outMultimapperOrder Random \     # 多比对随机排序
    --winAnchorMultimapNmax 100 \      # 多比对锚点最大数（重要！）
    --outFilterMultimapNmax 100 \      # 允许最多100个比对位置（重要！）
    --outSAMmultNmax 1 \              # 每个位置只输出1个比对
    --runThreadN 16 \                  # 16线程
    --outFileNamePrefix sample_        # 输出前缀

# 索引BAM
samtools index sample_Aligned.sortedByCoord.out.bam

# 关键：--outFilterMultimapNmax 100 而非默认的10
# 因为TE序列在基因组中有很多拷贝，需要保留多比对信息

3.2 TEcount定量¶

# === TEcount：TE和基因的联合定量 ===
TEcount \
    -b sample_Aligned.sortedByCoord.out.bam \  # STAR输出的BAM
    --GTF genes.gtf \              # 基因注释GTF
    --TE te_annotation.gtf \       # TE注释GTF（从RepeatMasker生成）
    --mode multi \                 # 多比对reads处理模式（推荐multi）
    --project sample_te            # 输出前缀

# === 输出文件 ===
# sample_te.cntTable  — 基因+TE的计数表
# 格式：feature_name\tcount
# 前面是基因计数，后面是TE家族计数

# mode选项：
# "multi"  — 将multi-mapped reads按比例分配给各TE位点（推荐）
# "uniq"   — 只用uniquely-mapped reads（低估TE表达）

3.3 TEtranscripts差异表达分析¶

# === TEtranscripts差异表达（多样本比较）===
TEtranscripts \
    -t treatment1.bam treatment2.bam treatment3.bam \  # 处理组BAM
    -c control1.bam control2.bam control3.bam \        # 对照组BAM
    --GTF genes.gtf \              # 基因注释
    --TE te_annotation.gtf \       # TE注释
    --mode multi \                 # 多比对模式
    --minread 10 \                 # 最少reads数
    --project te_diff              # 输出前缀

# 输出：
# te_diff_gene_TE_analysis.txt — DESeq2差异分析结果
# 包含基因和TE家族的差异表达结果

3.4 准备TE注释GTF¶

# === 从RepeatMasker输出生成TEtranscripts需要的GTF ===
# TEtranscripts提供了预制的GTF文件
# 下载地址：http://hammelllab.labsites.cshl.edu/software/

# 或者从RepeatMasker .out文件自己生成
python3 -c "
import sys

# 读取RepeatMasker .out文件，转为GTF格式
with open('genome.fasta.out') as fin, open('te_annotation.gtf', 'w') as fout:
    for i, line in enumerate(fin):
        if i < 3:  # 跳过头3行
            continue
        parts = line.split()
        if len(parts) < 15:
            continue

        chrom = parts[4]           # 染色体
        start = parts[5]           # 起始
        end = parts[6]             # 结束
        strand = '+' if parts[8] == '+' else '-'  # 方向
        te_name = parts[9]         # TE名称
        te_class = parts[10]       # TE类别
        te_family = te_class.split('/')[0] if '/' in te_class else te_class

        # 写GTF行
        attrs = f'gene_id \"{te_name}\"; transcript_id \"{te_name}\"; family_id \"{te_family}\"; class_id \"{te_class}\";'
        fout.write(f'{chrom}\tRepeatMasker\texon\t{start}\t{end}\t.\t{strand}\t.\t{attrs}\n')

print('TE GTF注释生成完成')
"

四、TE表达结果分析¶

# === R语言分析TE差异表达 ===
library(ggplot2)  # 加载ggplot2
library(dplyr)    # 数据处理

# 读取TEtranscripts差异分析结果
te_results <- read.table(
    "te_diff_gene_TE_analysis.txt",
    header = TRUE, sep = "\t"
)

# 分离基因和TE结果
te_only <- te_results %>%
    dplyr::filter(grepl(":", gene)) %>%  # TE名称含冒号（如LINE:L1）
    dplyr::mutate(
        te_class = sub(":.*", "", gene),  # 提取TE类别
        te_family = sub(".*:", "", gene)  # 提取TE家族
    )

gene_only <- te_results %>%
    dplyr::filter(!grepl(":", gene))  # 基因名不含冒号

# 火山图：TE差异表达
ggplot(te_only, aes(x = log2FoldChange, y = -log10(padj))) +
    geom_point(aes(color = te_class), alpha = 0.6, size = 2) +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10(adjusted p-value)",
         title = "TE Differential Expression") +
    theme_minimal()

ggsave("te_volcano.pdf", width = 10, height = 8)

# TE类别水平的汇总
te_summary <- te_only %>%
    dplyr::group_by(te_class) %>%
    dplyr::summarise(
        total = n(),
        up = sum(log2FoldChange > 1 & padj < 0.05, na.rm = TRUE),
        down = sum(log2FoldChange < -1 & padj < 0.05, na.rm = TRUE)
    ) %>%
    dplyr::arrange(desc(up + down))

print(te_summary)

五、常见报错与解决¶

报错信息	原因	解决方案
`No TE annotation file`	TE GTF路径错误	下载预制GTF或从RepeatMasker生成
`Too few multi-mapped reads`	STAR过滤了多比对reads	设置`--outFilterMultimapNmax 100`
`RepeatMasker: no species`	物种不在Dfam/Repbase	用`-lib`指定自定义TE库
`Memory error in STAR`	基因组太大	增加`--limitGenomeGenerateRAM`
`TEcount very slow`	BAM文件太大	按染色体分别处理
`No results for TE families`	TE GTF格式不匹配	确认gene_id和transcript_id属性名

六、面试高频题¶

Q1：为什么TE表达分析需要特殊处理multi-mapped reads？¶

答： TE序列在基因组中有成千上万个拷贝，所以来自TE的reads往往能比对到多个位置——这就是multi-mapped reads。传统RNA-seq分析（如featureCounts）默认丢弃这些reads或只保留一个最佳比对，会严重低估TE的转录活性。TEtranscripts用统计方法（EM算法）将multi-mapped reads按比例分配给各个TE位点，同时联合考虑基因和TE的计数，避免基因和TE之间的reads错误分配。这是TEtranscripts的核心创新。

Q2：RepeatMasker的工作原理是什么？¶

答： RepeatMasker用序列相似性搜索来识别基因组中的重复序列。步骤：(1) 将基因组序列与已知TE序列数据库（Dfam/Repbase）比对；(2) 用交叉匹配（cross_match）、RMBlast或hmmer作为搜索引擎；(3) 根据相似度和覆盖度判断每个匹配是否是真正的TE；(4) 对重叠的匹配进行优先级排序和去冗余。输出包括每个TE的位置、类型、方向、与一致性序列的divergence等信息。对于新物种，需要先用RepeatModeler从头发现TE，然后用发现的TE库做注释。

Q3：TE在基因组中有什么生物学功能？¶

答： TE不是"垃圾DNA"：(1) 基因调控——TE序列中包含增强子、启动子等调控元件，可以影响附近基因的表达；(2) 基因组创新——TE插入可以产生新基因、新外显子、新可变剪接形式；(3) 基因组结构——着丝粒和端粒区域富含TE；(4) 疾病相关——TE异常激活与肿瘤发生、衰老、自身免疫有关（如L1在肿瘤中去甲基化后重新活跃）；(5) 进化驱动力——TE插入是基因组变异的重要来源。近年研究发现TE在干细胞分化、胚胎发育、免疫应答中都有重要角色。

Q4：TEtranscripts和SalmonTE有什么区别？¶

答： TEtranscripts基于基因组比对（需要先用STAR比对），用EM算法处理multi-mapped reads，输出每个TE家族的计数并做差异分析。优点是精确度高，缺点是需要比对步骤，速度较慢。SalmonTE基于伪比对（不需要基因组比对），直接将reads与TE一致性序列进行快速量化，速度快很多但精度稍低。选择建议：追求精确度用TEtranscripts，追求速度用SalmonTE，大规模筛选用SalmonTE然后对感兴趣的用TEtranscripts验证。

参考资料：TEtranscripts: Jin et al., Bioinformatics 2015 | RepeatMasker: Smit et al., repeatmasker.org | Dfam: Storer et al., Mobile DNA 2021 | TE biology综述: Chuong et al., Nat Rev Genet 2017 | TE Workshop 2024: Cold Spring Harbor