转座子与重复序列注释（RepeatMasker / RepeatModeler / EDTA）¶

# 构建数据库
BuildDatabase -name genome_db -engine ncbi genome.fa

# 运行RepeatModeler2（包含LTRStruct）
RepeatModeler -database genome_db \
  -pa 16 \
  -LTRStruct

# 输出：
# genome_db-families.fa      - 重复序列共识库
# genome_db-families.stk     - Stockholm格式的种子比对
# RM_*/                      - 中间文件目录

# 查看结果统计
grep -c "^>" genome_db-families.fa
grep "^>" genome_db-families.fa | grep -oP "#\S+" | sort | uniq -c | sort -rn

# Step 1: LTRharvest候选鉴定
gt suffixerator -db genome.fa -indexname genome_idx -tis -suf -lcp -des -ssp -sds -dna
gt ltrharvest -index genome_idx -minlenltr 100 -maxlenltr 7000 \
  -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1 \
  -similar 85 -vic 10 -seed 20 -seqids yes \
  > ltrharvest_raw.gff3
gt ltrharvest -index genome_idx -minlenltr 100 -maxlenltr 7000 \
  -mintsd 4 -maxtsd 6 \
  -similar 85 -vic 10 -seed 20 -seqids yes \
  > ltrharvest_nomotif.gff3

# Step 2: LTR_FINDER候选鉴定（补充）
ltr_finder genome.fa -w 2 -C -D 15000 -d 1000 -L 7000 -l 100 -p 20 -M 0.85 > ltr_finder.scn

# Step 3: LTR_retriever精确过滤和分类
cat ltrharvest_raw.gff3 ltrharvest_nomotif.gff3 > ltrharvest_all.gff3

LTR_retriever -genome genome.fa \
  -inharvest ltrharvest_all.gff3 \
  -infinder ltr_finder.scn \
  -threads 16

# 输出：
# genome.fa.mod.LTRlib.fa      - 非冗余LTR库
# genome.fa.mod.out            - RepeatMasker格式注释
# genome.fa.mod.pass.list      - 完整LTR列表+分类+插入时间

# TIR-Learner
TIR-Learner.sh -g genome.fa -s species_name -t 16 -l 500

# MITE-Hunter
MITE_Hunter_manager.pl \
  -i genome.fa \
  -g species_name \
  -n 16 \
  -S 12345678

# 或使用MITE-Tracker
python3 -m MITETracker -g genome.fa -j mite_output -w 16

# HelitronScanner
java -jar HelitronScanner.jar scanHead \
  -genome genome.fa \
  -buffer_size 0 \
  -output head.txt

java -jar HelitronScanner.jar scanTail \
  -genome genome.fa \
  -buffer_size 0 \
  -output tail.txt

java -jar HelitronScanner.jar pairEnds \
  -head head.txt \
  -tail tail.txt \
  -genome genome.fa \
  -output helitrons.fa

# EDTA基本运行
EDTA.pl \
  --genome genome.fa \
  --threads 16 \
  --overwrite 1 \
  --sensitive 1

# 提供CDS信息以排除基因假阳性
EDTA.pl \
  --genome genome.fa \
  --cds gene_cds.fa \
  --threads 16 \
  --overwrite 1 \
  --sensitive 1 \
  --anno 1

# 如果已有curated TE库，可以合并
EDTA.pl \
  --genome genome.fa \
  --curatedlib curated_te.fa \
  --threads 16 \
  --anno 1

# 输出文件：
# genome.fa.mod.EDTA.TElib.fa      - 非冗余TE库
# genome.fa.mod.EDTA.TEanno.gff3   - 全基因组TE注释
# genome.fa.mod.EDTA.TEanno.sum    - 注释统计摘要
# genome.fa.mod.out                - RepeatMasker格式输出

# 使用EDTA生成的TE库进行屏蔽
RepeatMasker -pa 16 \
  -lib genome.fa.mod.EDTA.TElib.fa \
  -xsmall \
  -gff \
  -dir repeat_masker_output \
  genome.fa

# 使用Repbase+自定义库
RepeatMasker -pa 16 \
  -species "zea mays" \
  -lib custom_repeats.fa \
  -xsmall \
  -gff \
  -a \
  -dir repeat_masker_output \
  genome.fa

# -a: 生成比对文件(.align)，用于后续Kimura距离计算

# 查看RepeatMasker统计摘要
cat repeat_masker_output/genome.fa.tbl

# 或解析EDTA输出
cat genome.fa.mod.EDTA.TEanno.sum

# R中详细统计和可视化
library(tidyverse)

# 解析RepeatMasker .out文件
rm_out <- read.table("repeat_masker_output/genome.fa.out",
                     skip=3, fill=TRUE,
                     col.names=c("score","div","del","ins","query","begin","end",
                                 "left","strand","repeat_name","class","rbegin","rend","rleft","ID"))

# 提取大类信息
rm_out <- rm_out %>%
  mutate(
    length = end - begin + 1,
    superfamily = case_when(
      grepl("LTR/Gypsy", class) ~ "LTR/Gypsy",
      grepl("LTR/Copia", class) ~ "LTR/Copia",
      grepl("LTR", class) ~ "LTR/Other",
      grepl("LINE", class) ~ "LINE",
      grepl("SINE", class) ~ "SINE",
      grepl("DNA/TcMar", class) ~ "DNA/TcMar",
      grepl("DNA/hAT", class) ~ "DNA/hAT",
      grepl("DNA/MULE", class) ~ "DNA/MULE",
      grepl("DNA/CACTA", class) ~ "DNA/CACTA",
      grepl("DNA", class) ~ "DNA/Other",
      grepl("Helitron", class) ~ "Helitron",
      grepl("Simple_repeat", class) ~ "Simple_repeat",
      TRUE ~ "Other"
    )
  )

# 统计各超家族
genome_size <- 500000000  # 基因组大小(bp)

summary_table <- rm_out %>%
  group_by(superfamily) %>%
  summarise(
    count = n(),
    total_bp = sum(length),
    pct_genome = total_bp / genome_size * 100
  ) %>%
  arrange(desc(total_bp))

print(summary_table)

# 饼图
ggplot(summary_table, aes(x="", y=pct_genome, fill=superfamily)) +
  geom_col(width=1) +
  coord_polar(theta="y") +
  labs(title="基因组重复序列组成") +
  theme_void()

# LTR_retriever已经计算了插入时间
# 查看结果
head genome.fa.mod.pass.list
# 列：chr start end internal_id LTR_id TSD motif strand superfamily insertion_time

# 自定义计算（如需不同突变率）
# 提取LTR对，做配对比对
# 使用Kimura 2-parameter距离

# 插入时间分布可视化
ltr_data <- read.table("genome.fa.mod.pass.list", header=FALSE, sep="\t")
colnames(ltr_data) <- c("chr","start","end","id","ltr_id","tsd","motif",
                         "strand","superfamily","insertion_time")

# 按超家族绘制插入时间分布
ggplot(ltr_data, aes(x=insertion_time, fill=superfamily)) +
  geom_histogram(bins=50, alpha=0.7, position="stack") +
  labs(x="插入时间 (MYA)", y="LTR-RT数量",
       title="LTR逆转录转座子插入时间分布") +
  scale_fill_brewer(palette="Set2") +
  xlim(0, 10) +
  theme_minimal(base_size=14)

# 1. 提取RT结构域序列
# 使用TEsorter分类和提取
TEsorter genome.fa.mod.EDTA.TElib.fa -db rexdb -p 16

# 2. 提取Gypsy家族的RT序列
grep -A1 "Gypsy.*RT" TEsorter_output.fa > gypsy_RT.fa

# 3. 多序列比对
mafft --auto --thread 16 gypsy_RT.fa > gypsy_RT_aligned.fa

# 4. 构建进化树
FastTree -gamma gypsy_RT_aligned.fa > gypsy_RT.tree

# 或使用IQ-TREE（更准确）
iqtree -s gypsy_RT_aligned.fa -m MFP -bb 1000 -nt 16

# 1. 运行RepeatMasker时使用-a选项生成align文件
RepeatMasker -pa 16 -lib TE_lib.fa -xsmall -a -gff genome.fa

# 2. 使用RepeatMasker自带脚本计算Kimura距离
calcDivergenceFromAlign.pl \
  -s genome.fa.align.divsum \
  genome.fa.align

# 3. 生成Landscape数据
createRepeatLandscape.pl \
  -div genome.fa.align.divsum \
  -g 500000000 \
  > repeat_landscape.html

# R自定义Repeat Landscape图
library(ggplot2)
library(tidyverse)

# 解析divsum文件
parse_divsum <- function(divsum_file, genome_size) {
  lines <- readLines(divsum_file)
  # 找到coverage数据部分
  data_start <- grep("^Div", lines) + 1
  # 解析数据...
  # (具体格式根据文件版本调整)
}

# 或从RepeatMasker .out文件直接计算
landscape_data <- rm_out %>%
  mutate(kimura_bin = round(div)) %>%
  group_by(kimura_bin, superfamily) %>%
  summarise(total_bp = sum(length), .groups="drop") %>%
  mutate(pct_genome = total_bp / genome_size * 100)

# Repeat Landscape堆叠面积图
ggplot(landscape_data, aes(x=kimura_bin, y=pct_genome, fill=superfamily)) +
  geom_area(position="stack") +
  scale_fill_manual(values=c(
    "LTR/Gypsy"="darkblue", "LTR/Copia"="lightblue",
    "DNA/hAT"="orange", "DNA/CACTA"="red",
    "LINE"="green", "SINE"="yellow",
    "Helitron"="purple", "Other"="grey"
  )) +
  labs(x="Kimura substitution level (%)", y="基因组占比 (%)",
       title="Repeat Landscape") +
  xlim(0, 50) +
  theme_minimal(base_size=14) +
  theme(legend.position="right")

ggsave("repeat_landscape.pdf", width=12, height=6)

#!/bin/bash
# ====================================================
# 转座子与重复序列注释完整流程
# ====================================================

set -euo pipefail
THREADS=16
GENOME="genome.fa"
CDS="gene_cds.fa"  # 可选：已知CDS用于排除假阳性
OUTDIR="repeat_annotation"
GENOME_SIZE=$(grep -v "^>" ${GENOME} | tr -d '\n' | wc -c)

mkdir -p ${OUTDIR}/{edta,repeatmasker,analysis,figures}

echo "基因组大小: ${GENOME_SIZE} bp"

# ==== 方法1: EDTA一站式（推荐） ====
echo "=== EDTA注释 ==="
EDTA.pl \
  --genome ${GENOME} \
  --cds ${CDS} \
  --threads ${THREADS} \
  --overwrite 1 \
  --sensitive 1 \
  --anno 1

# 移动结果
cp ${GENOME}.mod.EDTA.TElib.fa ${OUTDIR}/edta/
cp ${GENOME}.mod.EDTA.TEanno.gff3 ${OUTDIR}/edta/
cp ${GENOME}.mod.EDTA.TEanno.sum ${OUTDIR}/edta/

echo "--- EDTA注释统计 ---"
cat ${GENOME}.mod.EDTA.TEanno.sum

# ==== 方法2: RepeatModeler + RepeatMasker（经典方法） ====
echo ""
echo "=== RepeatModeler + RepeatMasker ==="
# RepeatModeler
BuildDatabase -name ${OUTDIR}/repeatmasker/genome_db ${GENOME}
RepeatModeler -database ${OUTDIR}/repeatmasker/genome_db -pa ${THREADS} -LTRStruct

# RepeatMasker
RepeatMasker -pa ${THREADS} \
  -lib ${OUTDIR}/repeatmasker/genome_db-families.fa \
  -xsmall -gff -a \
  -dir ${OUTDIR}/repeatmasker \
  ${GENOME}

echo "--- RepeatMasker统计 ---"
cat ${OUTDIR}/repeatmasker/$(basename ${GENOME}).tbl

# ==== Repeat Landscape ====
echo ""
echo "=== Repeat Landscape ==="
calcDivergenceFromAlign.pl \
  -s ${OUTDIR}/analysis/divsum.txt \
  ${OUTDIR}/repeatmasker/$(basename ${GENOME}).align

createRepeatLandscape.pl \
  -div ${OUTDIR}/analysis/divsum.txt \
  -g ${GENOME_SIZE} \
  > ${OUTDIR}/figures/repeat_landscape.html

# ==== LTR插入时间 ====
echo ""
echo "=== LTR插入时间分析 ==="
if [ -f "${GENOME}.mod.pass.list" ]; then
  cp ${GENOME}.mod.pass.list ${OUTDIR}/analysis/
  echo "LTR-RT完整元素数: $(wc -l < ${OUTDIR}/analysis/$(basename ${GENOME}).mod.pass.list)"
fi

echo ""
echo "=== 注释完成 ==="
echo "TE库: ${OUTDIR}/edta/$(basename ${GENOME}).mod.EDTA.TElib.fa"
echo "注释GFF: ${OUTDIR}/edta/$(basename ${GENOME}).mod.EDTA.TEanno.gff3"
echo "Landscape: ${OUTDIR}/figures/repeat_landscape.html"