基因组注释从0开始（1）：重复序列注释

20 Dec 2024 14578字 49分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

接手一个已经T2T组装完成的香椿基因组，接下来就是对基因组进行注释了，我本以为是比较简单的，然后查看了一些资料，发现并没有想象中的那么容易，又加上我是小白，所以感觉会比较复杂吧，那本系列笔记就作为我如何对组装的植物基因组进行注释吧。

基因组注释可以分为两部分：基因组的结构注释（重复序列识别、非编码基因预测、编码基因预测）和基因功能注释，结构注释是功能注释的基础。

这里先从结构注释中的重复序列注释开始。我们知道植物基因组多倍化频繁，且基因组中存在大量的重复序列（有的植物基因组中重复序列甚至能达到80%），这些重复序列控制植物表型调控中有非常重要的作用。基因组中的重复序列可以分为以下几种：

1

1. 重复序列注释

串联重复（Tandem Repeat, TR）指DNA中的一个或多个核苷酸前后相连接的重复。串联重复又分为卫星DNA(Satellite DNA)、小卫星（Minisatellite）、微卫星（Microsatellite）。微卫星在植物中一般称为SSR（Simple Sequence Repeats）SSR在植物基因组常被用做遗传标记使用。重复序列占基因组非常高的比例，对重复序列的注释一般是做基因组注释的第一步。常用的基因组重复序列注释软件有RepeatMasker, RepeatModeler, EDTA。

散在重复序列（也称转座子，transposable element，TE）是可以在基因组内改变位置的一段DNA序列，通常由DNA复制造成，TE是基因组的重要组成部分。

2. RepeatMasker和 RepeatModeler

2.1. RepeatMasker安装

RepeatMasker是基因组重复序列检测的常用工具。一般依赖于已有的重复序列参考库Repbase作同源预测。对于绝大部分目标真核物种，都收录在Repbase中。

RepeatMasker安装推荐手动安装（推荐），进入官网复制安装包下载连接，可以用wget下载后再解压到自己的软件目录下。

wget https://www.repeatmasker.org/RepeatMasker/RepeatMasker-4.1.7-p1.tar.gz 
tar -xvzf RepeatMasker-4.1.7-p1.tar.gz ./###

也可以conda安装：

conda install -c bioconda repeatmasker

下载Dfam3.8数据库

FamDB 现在支持通过分类群（taxonomic groups）来对数据库进行分区。这意味着用户可以选择只下载数据库的一部分，以便进行特定物种或分类群的工作，而不必下载整个数据库。尽管如此，root partition（根分区）是必须下载的，其包含所有其他分区的一些基本信息。

每个分区包含特定分类群的信息。以下是分区的具体信息：

​ • Partition 0 [dfam38-1_full.0.h5]: 包含根分区以及多个广泛的分类群，如哺乳动物（Mammalia）、细菌（Bacteria）、真菌（Fungi）、病毒（Viruses）等。

​ • Partition 1 [dfam38-1_full.1.h5]: 针对特定的分类群 Obtectomera。

​ • Partition 2 [dfam38-1_full.2.h5]: 针对 Euteleosteomorpha（硬骨鱼类）。

​ • Partition 3 [dfam38-1_full.3.h5]: 针对 Sarcopterygii（肉鳍鱼类）以及爬行动物、两栖动物等。

​ • Partition 4 [dfam38-1_full.4.h5]: 针对 Diptera（双翅目昆虫）。

​ • Partition 5 [dfam38-1_full.5.h5]: 针对 Viridiplantae（绿色植物）。

​ • Partition 6 [dfam38-1_full.6.h5]: 针对 Deuterostomia（后口动物），包括软骨鱼类、海鞘等。

​ • Partition 7 [dfam38-1_full.7.h5]: 针对 Hymenoptera（膜翅目昆虫）。

​ • Partition 8 [dfam38-1_full.8.h5]: 针对 Ecdysozoa（蜕皮动物），包括线虫、昆虫等。

比如我的是植物，需下载Partition 0、Partition 5

​ • Dfam 3.8 数据库是通过 HDF5 格式 分区的，每个分区包含不同的分类群（例如植物、哺乳动物、昆虫等）。

​ • Partition 0（dfam38-1_full.0.h5）是根分区，包含许多不同的生物类别，如哺乳动物、细菌、真菌、病毒等。这个分区是 必需的，如果你只下载了其他分区而没有 Partition 0，RepeatMasker 就会报错。

​ • 其他分区，如 Partition 5（dfam38-1_full.5.h5），仅包含与 Viridiplantae（绿色植物）相关的数据。你下载了这个文件，所以它适用于植物物种，但它 必须和 Partition 0 一起使用。

nohup wget -c -q https://www.dfam.org/releases/Dfam_3.8/families/FamDB/dfam38-1_full.0.h5.gz > dfam38-1_full.0.download.blog &
nohup wget -c https://www.dfam.org/releases/Dfam_3.8/families/FamDB/dfam38-1_full.5.h5.gz > dfam38-1_full.5.download.blog &

这边文件比较大，还是建议用电脑下载，然后再上传到服务器，服务器直接下载太慢了。

加载Repbase数据库：

RepeatMasker自带的重复序列数据库是Dfam数据库，这是一个转座子（TE）序列数据库，收录的物种比较少。Repbase是重复序列参考数据库，其中收录了大部分真核物种，适用于重复序列的同源预测。然而Repbase不是RepeatMasker自带的，需要额外下载，这里网上找到了一个20181026版本的Repbase下载地址：点击这里

下载Repbase数据库后用tar -xvf解压，将RMRBSeqs.embl和README.RMRBSeqs两个数据库文件放在RepeatMasker安装目录的Libraries目录下，注意不要修改后缀名。

tar -xvf RepBaseRepeatMaskerEdition-20181026.tar

在RepeatMasker安装目录下运行perl ./configure，一路回车确定路径，如果有缺失的依赖就用conda下载，一直到最后选择序列搜索比对的软件，我这里输入3回车，之后的界面再输入5回车确认，当看到提示信息Dfam和RBRM（也就是RepBase数据库）两个数据库版本的时候，就说明加载Repbase数据库成功了。

2

在安装目录下用RepeatMasker -h查看是否可以正常运行，如果正常输出说明已经安装成功。最后需要修改一下环境变量（不修改运行的时候找不到pm文件），将RepeatMasker 安装路径添加到PERL5LIB环境变量中：

  1. 确认 RepeatMasker 可执行文件路径

  2. 修改 ~/.bashrc文件**,在最后一行添加 RepeatMasker 路径：

    nano ~/.bashrc
export PERL5LIB=/**/**/SoftWare/RepeatMasker/lib:$PERL5LIB

  3. 保存更改并退出编辑器。如果你使用 nano 编辑器，按 Ctrl + O 保存，按 Ctrl + X 退出。

  4. 使修改生效，并验证配置，如果正确配置，应该显示 RepeatMasker 可执行文件的路径:

    source ~/.bashrc
which RepeatMasker

安装TRF

wget -c https://github.com/Benson-Genomics-Lab/TRF/releases/tag/v4.09.1
tar -zxvf TRF-4.09.1.tar.gz
cd TRF-4.09.1
./configure --prefix=to/you/path
make
make install

安装rmblast

wget -c https://www.repeatmasker.org/rmblast/rmblast-2.14.1+-x64-linux.tar.gz
tar rmblast-2.14.1+-x64-linux.tar.gz

安装hmmer

wget -c http://eddylab.org/software/hmmer/hmmer-3.2.1.tar.gz
tar -zxvf hmmer-3.2.1.tar.gz 
./configure --prefix=to/you/path
make
make install

2.2. RepeatMasker使用

2.2.1. 查看数据库中的物种

  1. Repbase数据库配置好以后，用RepeatMasker/util/queryRepeatDatabse.pl -tree >species.txt可以在species.txt中查看数据库中的物种列表。但RepeatMasker-4.1.1版本及以后不再提供queryRepeatDatabase.pl脚本，而是由famdb.py代替。比如查询Sapindales可以使用python famdb.py names Meliaceae来获取。
  2. 查看RepeatMasker/Libraries/taxonomy.dat文件，也有所有已收录物种的名称。
  • 找到数据库中已有训练的物种或者近缘种后，在使用RepeatMasker时用参数-species指定库中已有物种作为参考。

2.3. RepeatModeler及其依赖安装

2.3.1. RepeatModeler自我训练

  • 大部分情况下数据库中没有待注释物种的近缘种，此时最好用RepeatModeler等软件进行自我训练，用denovo预测来构建重复序列数据库。
  • 即使数据库中已有近缘种，有时候Repbase注释重复区的效果不是很好，也推荐进行自我训练。

2.3.2. RepeatModeler安装

  1. RepeatModeler依赖：
  • 配置好的RepeatMasker
  • RECON：Denovo Repeat Finder
  • RepeatScout：Denovo Repeat Finder
  • NSEG

  1. RepeatModeler支持conda安装：conda install -c bioconda repeatmodeler

  2. 或者下载源码编译，再一步步指定依赖软件路径，类似RepeatMasker的安装过程。这边由于https://repeatmasker.org/RepeatModeler/访问错误，选择github源码下载：

    git clone https://github.com/Dfam-consortium/RepeatModeler.git

  3. 进入RepeatModeler文件，perl ./configure，这边提示缺少一些依赖，可以用conda安装：

    3

    conda install -c conda-forge perl-file-which
conda install -c conda-forge perl-devel-size
conda install -c conda-forge perl-lwp-useragent

    有时候安装不了，要看看当前的perl是用的哪一个，一般用conda里面的，然后就是可以cpan安装，但是我碰到了make失败，比较麻烦，就可以通过cpan得到软件包，直接解压后手动安装就行了。

    tar -xvzf JSON-XS-4.03.tar.gz
cd JSON-XS-4.03
perl Makefile.PL
make
make install

  4. 重新运行配置perl ./configure

    安装recon

    踩坑，需要输入recon位置信息，结果使用conda install recon安装后，运行不了，后面我找到了recon官网，下载安装包重新安装了：

    5

​ 安装解压recon后，bin目录是空的，还需进行一下几个步骤：

  1. 运行 make 编译 RECON，进入 src 目录后，运行make编译一下，再运行make install然后安装到bin目录：

    cd recon/src
make
make install
cd ..

  2. 安装完成后，你需要配置 recon.pl 脚本，以便它能找到 bin 目录中的可执行文件。

    进入 scripts 目录，编辑 recon.pl 脚本，修改第三行，添加 bin 目录的路径

    cd scripts
vim recon.pl #修改bin目录路径，vscode修改方便些

    6

  3. 接下来就可以把bin目录的路径添加到RepeatModeler的./configure配置中了。

​ 安装RepeatScout

​ 在官网获取RepeatScout的安装包，同样的需要进入文件夹make编译一下，然后将路径添加到RepeatModeler的./configure配置中。

​ 安装TRF

​ 直接git,然后根据下面操作，一步步编译，最后将路径添加进去。

git clone https://github.com/Benson-Genomics-Lab/TRF.git#克隆代码仓库并进入项目目录
cd TRF
mkdir build#创建构建目录
cd build
../configure#运行配置脚本
make#编译源码
cp src/trf ../#手动复制可执行文件

​ 安装CD-Hit

​ 同样直接git，然后编译，添加路径

cd cdhitgit clone https://github.com/weizhongli/cdhit.git
cd cdhit
make

​ 安装UCSC，官网

​ 下载安装包，添加路径

#http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64
mkdir UCSCTOOLS
cd UCSCTOOLS
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitInfo
chmod 755 *
#将三个文件放在新建的UCSCTOOLS文件夹中

​ 安装RMBLAST，官网，解压完成后进入bin文件，添加bin文件路径

wget -c https://www.repeatmasker.org/rmblast/rmblast-2.14.1+-x64-linux.tar.gz#如果失效，返回上一目录寻找最新的
tar -xzvf rmblast-2.14.1+-x64-linux.tar.gz#解压

​ 这边安装的预编译版本有些问题，我这边下载源码重新安装一下：

wget -c https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.14.1/ncbi-blast-2.14.1+-src.tar.gz
wget -c https://www.repeatmasker.org/rmblast/isb-2.14.1+-rmblast.patch.gz #补丁
tar -zxvf ncbi-blast-2.14.1+-src.tar.gz
gunzip isb-2.14.1+-rmblast.patch.gz
cd ncbi-blast-2.14.1+-src
patch -p1 < ../isb-2.14.1+-rmblast.patch
./configure --with-mt --without-debug --without-krb5 --without-openssl --with-projects=scripts/projects/rmblastn/project.lst --prefix=/usr/local/rmblast

make -j4#启用4个核心，加速编译
make install

--with-mt：启用多线程支持。

--without-debug：禁用调试信息（优化构建）。

--without-krb5 和 --without-openssl：禁用 Kerberos 和 OpenSSL 支持。

--with-projects：指定 RMBlast 项目路径。

--prefix：设置安装路径。

安装genometools，官网

wget -c https://genometools.org/pub/genometools-1.6.2.tar.gz
tar -xzvf genometools-1.6.2.tar.gz
cd genometools-1.6.2
make#编译

​ 上面这个安装要安装一大堆依赖文件，我这边因为缺的比较多，所以选择了conda直接一步到位：

conda install broad-viral::genometools-genometools

​ Ltr_retriever安装，需要指定的v2.9.0版，手动安装：

wget -c https://github.com/oushujun/LTR_retriever/archive/refs/tags/v2.9.0.tar.gz
tar -zxvf v2.9.0.tar.gz

​ 这边也直接用conda安装：

conda install -c bioconda -c conda-forge ltr_retriever=2.9.0

​ 安装后把路径添加进去就行了。

​ MAFFT安装

wget -c https://mafft.cbrc.jp/alignment/software/mafft-7.487-with-extensions-src.tgz
tar -zvxf 
cd mafft-7.487-with-extensions/core/
vi Makefile (any other text editor is ok.)

#第一行
From:
PREFIX = /usr/local
To:
PREFIX = /home/your_home/somewhere (must be absolute path)
#第三行
From:
BINDIR = $(PREFIX)/bin
To:
BINDIR = /home/your_home/bin (or elsewhere in your command-search path)

#有sudo权限可以直接sudo安装
sudo apt install mafft

这边官网安装不上，我就直接conda安装了：

conda install -c bioconda mafft
conda list mafft

​ Ninja安装，官网

wget -c https://github.com/TravisWheelerLab/NINJA/archive/refs/tags/0.95-cluster_only.tar.gz
tzr -zxvf 0.95-cluster_only.tar.gz
cd NINJA
make all

7

全部配置好了后，就开始配置使用了。

3. 注释流程

3.1. 导出同源物种重复序列库

前面2.1步骤将Repbase和Dfam数据库整合之后，RepeatMasker/Libraries目录下RepeatMaskerLib.h5这个文件为整合后构建的数据库文件，我们要在这个文件中导出同源物种的重复序列。

在RepeatMasker目录下提供了famdb.py这个程序查询目标近缘物种。如果你不知道自己的物种在什么分支上，我这里推荐一个查找已发表的植物基因组的网站Published Plant Genomes (plabipd.de)，可以一级一级查看哪些近缘物种有人做过了。用以下命令查看物种重复序列否收录到库中：

8

然后运行famdb.py，查找：

python famdb.py -i /data/wangxingbin/Software/RepeatMasker/Libraries/famdb lineage -ad Meliaceae

这边数据需要两个数据库融合，Dfam数据库文件比较大，可以点击下载,解压后，改名，替换~/RepeatMasker/Libraries中的Dfam.h5。再然后结合2.1下载的Repbase数据库配置一下。

找到最近的分支后，导出最近分支的祖先节点和衍生节点物种的重复序列库，使用内置的perl软件转换成fasta格式：

python famdb.py -i /data/wangxingbin/Software/RepeatMasker/Libraries/famdb families -f embl -a -d Meliaceae > lamiids.embl	
# 查找近缘物种及其上祖先节点，其下所有类群repeat famlies，输出格式embl。 -a ancestor，-d descendent

buildRMLibFromEMBL.pl Meliaceae.embl > Meliaceae.fasta	# 转换格式为fasta，方便后续合并

也可以用下面的bash脚本：

#!/bin/bash
# run_repeat_family.sh

# 设定变量
FAMDB_DIR="/data/wangxingbin/Software/RepeatMasker/Libraries/famdb"
SPECIES="Meliaceae"
EMBL_FILE="${SPECIES}.embl"
FASTA_FILE="${SPECIES}.fasta"
WORK_DIR="01.repeat_family"
FAMDB_PY="/data/wangxingbin/Software/RepeatMasker/famdb.py"  # famdb.py 文件路径
BUILD_RMLIB_SCRIPT="/data/wangxingbin/Software/RepeatMasker/util/buildRMLibFromEMBL.pl"  # buildRMLibFromEMBL.pl 脚本路径

# 检查文件和目录是否存在
if [[ ! -f "$FAMDB_PY" ]]; then
    echo "错误: famdb.py 文件不存在，请检查路径：$FAMDB_PY"
    exit 1
fi

if [[ ! -f "$BUILD_RMLIB_SCRIPT" ]]; then
    echo "错误: buildRMLibFromEMBL.pl 脚本不存在，请检查路径：$BUILD_RMLIB_SCRIPT"
    exit 1
fi

if [[ ! -d "$FAMDB_DIR" ]]; then
    echo "错误: RepeatMasker famdb 目录不存在：$FAMDB_DIR"
    exit 1
fi

# 创建工作目录并进入
mkdir -p "$WORK_DIR"
cd "$WORK_DIR" || { echo "无法进入目录 $WORK_DIR"; exit 1; }

# 提取 Repeat family 信息
echo "[$(date)] 正在提取 Repeat family 数据..."
python "$FAMDB_PY" -i "$FAMDB_DIR" families -f embl -a -d "$SPECIES" > "$EMBL_FILE"

# 检查 EMBL 文件是否生成
if [[ ! -s "$EMBL_FILE" ]]; then
    echo "错误: 生成的 EMBL 文件为空或不存在！生成文件：$EMBL_FILE"
    exit 1
fi

# 转换 EMBL 到 FASTA 格式
echo "[$(date)] 正在转换 EMBL 文件为 FASTA 格式..."
perl "$BUILD_RMLIB_SCRIPT" "$EMBL_FILE" > "$FASTA_FILE"

# 检查 FASTA 文件是否生成
if [[ ! -s "$FASTA_FILE" ]]; then
    echo "错误: 生成的 FASTA 文件为空或不存在！生成文件：$FASTA_FILE"
    exit 1
fi

echo "[$(date)] 任务完成，输出文件：$FASTA_FILE"

3.2. RepeatModeler从头开始预测

3.2.1. 建库

mkdir 1.repeatmasker#创建一个文件夹来保存
cd 1.repeatmasker
mv /path/to/your/genome.fa ./sample.fa
mkdir db && cd db
BuildDatabase -name sample ../sample.fa

此步骤会生成以下文件：

  • sample.nhr
  • sample.nin
  • sample.nnd
  • sample.nni
  • sample.nog
  • sample.nsq
  • sample.translation

这些文件组成了一个 BLAST 格式的数据库，供 RepeatModeler 使用。

注意事项：

BuildDatabase 未添加到环境变量 即使安装了 RepeatModeler，若工具路径未添加到 $PATH，也会导致找不到命令。

使用以下命令搜索 BuildDatabase 的路径：

find /data/wangxingbin/SoftWare/ -name BuildDatabase

添加 RepeatModeler 到环境变量

将 RepeatModeler 目录添加到 $PATH：

echo 'export PATH=/you/path/SoftWare/RepeatModeler:$PATH' >> ~/.bashrc
source ~/.bashrc

3.2.2. 运行 RepeatModeler

在 db 文件夹的上级目录运行 RepeatModeler：

cd ..
nohup RepeatModeler -threads 60 -database db/sample -engine ncbi &

参数说明：

  • -threads 60: 指定使用 60个线程（根据服务器配置选择）。
  • -database db/sample: 指定上一步构建的数据库路径。
  • -engine ncbi: 使用 NCBI BLAST 引擎（默认）。

注意：

运行过程可能会耗时数天，根据你的基因组大小和服务器性能决定。

3.2.3. 结果文件解析

运行完成后，会生成以下重要文件：

  • sample-families.fa：包含预测到的重复序列库，可供 RepeatMasker 使用。
  • sample-families.stk：与 Dfam 数据库交互的文件。
  • consensi.fa.classified：包含所有重复序列及其分类状态。

目录内还会生成类似 RM_21945.SatAug11650032020 的子文件夹，其中包含中间文件和结果文件。

3.2.4. 提取可信的重复序列

RepeatModeler 的输出中，有些重复序列可能无法被 RepeatMasker 识别，通常标记为 Unknown。 为了提高结果的可信度，需将这些 Unknown 序列剔除。

提取可以识别的序列

seqkit grep -r -i -p "Unknown" -v ./RM_21945.SatAug11650032020/consensi.fa.classified > ModelerID.lib

-v：提取不包含 Unknown 的序列。

ModelerID.lib 是最终可信的重复序列数据库，可供其他软件（如 RepeatMasker）使用。

提取无法识别的序列

seqkit grep -r -i -p "Unknown" ./RM_21945.SatAug11650032020/consensi.fa.classified > Modelerunknown.lib

提取所有标记为 Unknown 的序列，供进一步分析或验证。

3.3. RepeatMasker重复序列屏蔽

鉴定基因组重复区域的方法有两种：一种基于文库(library)的同源(homology)方法，该文库收集了其他物种的某一种重复的一致性序列，通过相似性来鉴定重复；另一种是从头预测(de novo)，将序列和自己比较或者是高频K-mer来鉴定重复。

这边可以结合上面从头预测的结果作为参考数据库运行：

mkdir 2.repeatMasker
cd 2.repeatMasker
nohup RepeatMasker ../../sample.fa -species "Meliaceae" -pa 20 -lib consensi.fa.classified -poly -html -gff -dir repeatmasker &

因为从头预测太过于耗时，也可以根据同源序列直接进行屏蔽：

mkdir 00.repeatMasker
cd 00.repeatMasker
nohup RepeatMasker ../../sample.fa -species "Meliaceae" -pa 30 -poly -html -gff -dir repeatmasker_Meliaceae > nohup_Meliaceae.out &

3.3.1. 合并数据库

还可以合并从头预测以及famdb的数据库，然后来跑，这样完整些。

#!/bin/bash
# run_repeatmasker.sh

# 设置参数
THREADS=25  # 线程数x4
LIBRARY="/data/wangxingbin/00.Project/01.Toona_Annotation/hap2/01.repeatMasker/02.TE/01.repeatModeler/db/hap2-families.fa"
LIBRARY2="/data/wangxingbin/00.Project/01.Toona_Annotation/hap2/01.repeatMasker/02.TE/02.repeatMasker/repeat_family/Meliaceae.fasta"
GENOME="/data/wangxingbin/00.Project/01.Toona_Annotation/hap2_T2T.fasta"
OUTPUT_DIR="./repeatmasker_output"
LOG_FILE="./repeatmasker.log"
FINAL_LIBRARY="final_repeat_lib.fasta"

# 检查输入文件是否存在
if [[ ! -f "$LIBRARY" ]]; then
    echo "错误：RepeatModeler 生成的库文件不存在：$LIBRARY" | tee -a "$LOG_FILE"
    exit 1
fi

if [[ ! -f "$LIBRARY2" ]]; then
    echo "错误：第二个库文件不存在：$LIBRARY2" | tee -a "$LOG_FILE"
    exit 1
fi

if [[ ! -f "$GENOME" ]]; then
    echo "错误：输入基因组文件不存在：$GENOME" | tee -a "$LOG_FILE"
    exit 1
fi

# 合并数据库
echo "[$(date)] 合并同源数据库和RepeatModeler训练结果..." | tee -a "$LOG_FILE"
cat "$LIBRARY" "$LIBRARY2" > "$FINAL_LIBRARY"
if [[ ! -f "$FINAL_LIBRARY" ]]; then
    echo "错误：合并数据库失败，文件不存在：$FINAL_LIBRARY" | tee -a "$LOG_FILE"
    exit 1
fi

# 获取时间戳
TIMESTAMP=$(date +"%Y%m%d_%H%M%S")
# 将时间戳写入日志文件的第一行
echo "RepeatModeler 任务启动时间: $TIMESTAMP" > "$LOG_FILE"

# 创建输出目录
mkdir -p "$OUTPUT_DIR"
echo "[$(date)] 输出目录已创建：$OUTPUT_DIR" | tee -a "$LOG_FILE"

# 运行 RepeatMasker
echo "[$(date)] 开始运行 RepeatMasker..." | tee -a "$LOG_FILE"

nohup RepeatMasker \
    -nolow -no_is -pa "$THREADS" \
    -lib "$FINAL_LIBRARY" -engine ncbi -gff -norna \
    -dir "$OUTPUT_DIR" "$GENOME" > "$LOG_FILE" 2>&1 &

# 任务启动完成
echo "[$(date)] RepeatMasker 正在后台运行，日志输出到：$LOG_FILE" | tee -a "$LOG_FILE"

-nolow Does not mask low_complexity DNA or simple repeats 不屏蔽低复杂度DNA或简单重复序列（有的学者认为simple repeat不算严格意义上的重复序列类型）

-norna Does not mask small RNA (pseudo) genes 不屏蔽sRNA

-no_is Skips bacterial insertion element check 跳过细菌插入元件检查

-pa 和RepeatModeler一样，1 pa是4个线程

-lib 指定自建数据库（与-species冲突）

-gff 生成gff文件

-dir 指定输出目录

3.3.1. 参数说明：

sample.fa

  • 输入基因组序列文件，需要对其进行重复序列注释。
  • 示例中为 sample.fa。

-species "Sapindales"

  • 指定物种分类参考库，用于提高注释精度。
  • 示例中选择了被子植物中的无患子目（Sapindales）。

-pa 20

  • 指定并行线程数以加速运行。
  • 在使用 RMBlast 引擎时，每线程需要 4 核，因此 -pa 20 表示使用 80 核心线程（12 × 4 = 80）。

-lib consensi.fa.classified

  • 指定用户自定义的重复序列库文件路径，一般为前面的RepeatModeler自我训练的库。
  • 此参数优先级高于 -species，即 RepeatMasker 会优先使用该文件进行注释。

-poly

  • 将多态性的简单重复序列（如微卫星重复序列）单独提取到 .polyout 文件中。

  • 微卫星重复序列通常标注为 Simple_repeat 类型。

  • -html输出结果文件的 HTML 格式，便于使用网页浏览器查看注释结果。

  • -gff生成 GFF 格式的注释文件，方便与其他基因组注释工具（如 BEDtools、GENEMARK）联用。

  • -dir repeatmasker指定输出文件保存的目录，结果文件将存储在 repeatmasker 目录下。

3.3.2. 输出结果：

  • sample.fa.tbl重复序列统计报告，包含基因组长度、GC 含量、重复序列长度及类型占比等。

  • sample.fa.out重复序列详细注释文件，显示重复序列类型、起始/结束位置、碱基替换率、插入/缺失率等信息。

  • sample.fa.out.htmlsample.fa.out 的 HTML 格式，可通过网页浏览器查看。

  • sample.fa.out.gffGFF 格式的重复序列注释文件，便于与其他分析工具联用。

  • sample.fa.polyout微卫星重复序列的注释结果，提取自 sample.fa.out。

  • sample.fa.masked屏蔽重复序列后的基因组文件，重复序列替换为 N。

其中，sample.fa.out每一列含义：

  • 第一列：比对分值，SW score
  • 第二列：替代率 perc div.
  • 第三列：碱基缺失百分率
  • 第四列：在重复序列中碱基缺失百分率
  • 第五列：query sequence
  • 第六列：查询序列起始位置
  • 第七列：查询序列终止位置
  • 第八列：查询区域中超出比对区域碱基的数- 目，也就是没有比对上的碱基数
  • 第九列：+/-(C)
  • 第十列：比上的重复序列名称，类型命名
  • 第十一列：比上重复序列的分类，和- repeatmolder 中*.classed 是一样的
  • 第十二列：比上的在数据库中的起始位置
  • 第十三列：比上的在数据库中的终止位置
  • 第十四列：在第十列上超出比对区域碱基的数目，也就是没有比对上的碱基数
  • 第十五列：比对区域的ID，随机给的