基因组注释-重复序列注释

1.介绍

Hite是一款基于python开发的专注于基因组中高精度检测和注释全长转座元件(full-length TEs)的生物信息学软件。相比于RepeatModeler、EDTA等软件，Hite基于动态边界调整识别更完整的TEs结构、更长的TEs。

2.安装

#将软件包下载到目录下/path
git clone https://github.com/CSU-KangHu/HiTE.git
#配置环境依赖
conda env create --name HiTE -f environment.yml
#激活环境
conda activate HiTE
#运行程序
python /path/HiTE/main.py --help

environment.yml

3.使用

3.1 流程简介

流程图

1.先将基因组分割为不同的chunk
2.TRF对每个chunk做Tandem repeats的屏蔽
3.将每个chunk分割成1M并使用blastn两两比对，根据FMEA算法找到最大的TEs的边界
4.根据HybridLTR鉴定LTR-RTs(也可以根据LTR_harvest and LTR_Finder鉴定，再用LTR_retriever整合，但是这个方法费时间，两种方法鉴定的LTR-RTs基本相近)
5.根据TE Finder的itrsearch寻找TIR元素
6.根据EAHelitron和HelitronScanner鉴定Helitron元素
7.使用HiTE-NonLTR算法利用RepeatMasker(v4.1.1)的Dfam库中已知的LINE和SINE建库，然后根据序列比对结果鉴定non-LTR元素

3.2 脚本

#source conda环境
source ~/conda_init.sh
#激活HiTE环境
conda activate HiTE
#导入依赖
export PATH=/path/HiTE/module:/path/HiTE/tools:$PATH

#输入参数
genome=$PWD/genome.fasta
out_dir=$PWD/Hite_out
thread=50
### 1: true, 0: false
isplant=0                      ### Is it a plant genome
isremove_nested=1              ### Whether to remove nested TE
isdomain=1                     ### Whether to obtain TE domains, HiTE uses RepeatPeps.lib from RepeatMasker to obtain TE domains
isrecover=1                    ### Whether to enable recovery mode to avoid starting from the beginning
isannotate=1                   ### Whether to annotate the genome using the TE library generated

python /path/HiTE/main.py \
       --genome ${genome} \
       --out_dir ${out_dir} \
       --thread ${thread} \
       --plant ${isplant} \
       --remove_nested ${isremove_nested} \
       --domain ${isdomain} \
       --recover ${isrecover} \
       --annotate ${isannotate}

参数设置：
1.chunk_size：默认是将基因组按照400Mb切分为一个chunk，后续的每个chunk都会执行Coarse-grained boundary mapping、determine fine-grained TIR、determine fine-grained Helitron、determine fine-grained Non-LTR四个步骤，会生成对应chunk名字开头的文件。如果你的基因组>3G，可以将该值设置小一点。
2.curated_lib:有信任的库可以添加，没有的话主要影响在于LINE的比例(参考issue49)
3.如果想对HiTE的注释结果和EDTA、repeatmodeler软件的结果做基准测试，可以添加参数--BM_RM2、--BM_EDTA、--BM_HiTE；
4.如果想加速比对，可以将--chrom_seg_length参数设置小一点，默认是1M；
Tips：导入依赖不建议写入.bashrc文件，随着生信分析的增多，.bashrc下面会有很多的软件依赖(比如perl、python等不同版本之间会存在冲突)。

4.输出

4.1 输出文件

1.chr_name.map：第一列为输出文件genome.rename.fa的名字，第二列输入基因组的序列ID；
2.confident_TE.cons.fa：HiTE鉴定的TE库
3.low_confident_TE.cons.fa：具有低全长拷贝数的未分类的TE候选序列
3.all_TE.fa：2和3的总和
4.HiTE.out, HiTE.gff, and HiTE.tbl：基因组的注释文件，HiTE.out类似repeatmasker的.out文件，HiTE.tbl是统计文件
5.intact_LTR.list：具有完整LTR-RTs结构的注释信息
6.confident_ltr.internal.fa：根据intact_LTR.list提取分不为Unknown的internal序列
7.confident_ltr.terminal.fa：根据intact_LTR.list提取分不为Unknown的terminal序列

4.2 细节展示

#输出每个亚家族的TEs和每条染色体被屏蔽的长度和比例
HiTE/tools/buildSummary.pl HiTE.out > HiTE.detail.tbl

4.3 成对LTR提取

写了一个脚本用于根据HiTE输出提取成对的LTR计算LTR插入时间或者计算成对的LTR之间的序列分歧度等，用法：python3 Hite_out_extract_intact_LTR.py -i intact_LTR.list -g genome.fasta -o output，如果将3'和5'LTR放入一个文件夹，则不设置该参数(--keep-split)

import argparse
from pyfaidx import Fasta
from pathlib import Path

def parse_arguments():
    parser = argparse.ArgumentParser(description='Extract and merge 5\'LTR and 3\'LTR sequences.')
    parser.add_argument('-i', '--input', required=True, help='Input annotation file')
    parser.add_argument('-g', '--genome', required=True, help='Genome FASTA file')
    parser.add_argument('-o', '--output', default='output', help='Output directory')
    parser.add_argument('--keep-split', action='store_true', help='Keep individual 5LTR/3LTR files')
    return parser.parse_args()

def extract_and_merge_ltr(annotation_file, genome_file, output_dir, keep_split=False):
    Path(output_dir).mkdir(parents=True, exist_ok=True)
    genome = Fasta(genome_file)  # 自动建立索引

    with open(annotation_file, 'r') as f:
        for line in f:
            if line.startswith('#') or not line.strip():
                continue

            parts = line.strip().split()
            if len(parts) < 7:
                continue

            # 解析染色体和坐标
            chr_coords = parts[0]
            chr_info, coords = chr_coords.split(':')
            start, end = map(int, coords.split('..'))

            # 解析 IN 坐标
            in_part = parts[6]
            if not in_part.startswith('IN:'):
                continue

            in_coords = in_part[3:]
            in_start, in_end = map(int, in_coords.split('..'))

            # 计算 5'LTR 和 3'LTR 区间
            ltr5_start, ltr5_end = start, in_start - 1
            ltr3_start, ltr3_end = in_end + 1, end

            # 提取序列
            ltr5_seq = genome[chr_info][ltr5_start-1:ltr5_end].seq  # pyfaidx 是 0-based
            ltr3_seq = genome[chr_info][ltr3_start-1:ltr3_end].seq

            # 生成文件名
            base_name = f"{chr_info}_{start}_{end}"
            merged_file = Path(output_dir) / f"{base_name}_LTR.fa"

            # 写入合并后的文件
            with open(merged_file, 'w') as f_out:
                f_out.write(f">{base_name}_5LTR\n{ltr5_seq}\n")
                f_out.write(f">{base_name}_3LTR\n{ltr3_seq}\n")

            # 可选：保留单独的 5LTR/3LTR 文件
            if keep_split:
                with open(Path(output_dir) / f"{base_name}_5LTR.fa", 'w') as f5:
                    f5.write(f">{base_name}_5LTR\n{ltr5_seq}\n")
                with open(Path(output_dir) / f"{base_name}_3LTR.fa", 'w') as f3:
                    f3.write(f">{base_name}_3LTR\n{ltr3_seq}\n")

def main():
    args = parse_arguments()
    extract_and_merge_ltr(args.input, args.genome, args.output, args.keep_split)
    print(f"Extraction and merging complete. Results are in {args.output} directory.")

if __name__ == '__main__':
    main()