基因组注释（2）——散在重复序列注释

散在重复序列可以分为反转录转座子（class-I TEs）和DNA转座子（class-II TEs）

反转录转座子：通过RNA介导的copy and paste机制进行转座，主要由LTR（long terminal repeat）构成，而non-LTR根据长度又分为LINEs（long interspersed nuclear elements）和SINEs（short interspersed elements）。

DNA转座子：通过DNA介导的cut and paste机制进行转座。

这里我们用RepeatModeler和RepeatMasker两个软件跑一遍基因组散在重复序列注释的流程，需要注意下因为前面做了TRF注释串联重复序列，我们运行RepeatMasker的时候要改下下参数设置。

不建议用conda安装两款软件的本体（但是可以安装其他依赖）

RepeatMasker配置成功过是RepeatModeler配置的前提条件，且两者之间有版本关联（比如最新的RepeatModeler版本为2.0.4，需要最新的RepeatMasker版本4.1.4安装为前提），conda直接安装RepeatMasker会导致RepeatModeler无法找到RepeatMasker的路径，且输入正确路径也会提示找不到（不知道是不是我的原因）。

下载源码包编译，可以看官网RepeatMasker Home Page。

本篇博客所使用RepeatMasker版本为4.1.2，RepeatModeler版本为2.0.3

2.1 RepeatMasker安装

本体安装过程不多说，主要说一下加载Repbase数据库：

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

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

在RepeatMasker安装目录下运行perl ./configure，一路回车确定路径，如果有缺失的依赖就用conda下载，一直到最后选择序列搜索比对的软件，我这里输入3回车，之后的界面再输入5回车确认：

当看到提示信息Dfam和RBRM（也就是RepBase数据库）两个数据库版本的时候，就说明加载Repbase数据库成功了。

用RepeatMasker -h查看是否可以正常运行，如果提示Devel::Size这个perl模块缺失，可以用conda安装：

1

conda install -c bioconda perl-devel-size

最后需要修改一下环境变量（不修改运行的时候找不到pm文件），将RepeatMasker 安装路径添加到PERL5LIB环境变量中：

1
2

# 打开 ~/.bashrc
export PERL5LIB="/public/home/wlxie/miniconda3/envs/biosoft/share/RepeatMasker:$PERL5LIB"

2.2 RepeatModeler安装

安装过程与RepeatMasker差不多，有一个比较坑的地方是官方可选的一部分软件（比如CD-HIT）在configure过程中是必须指定的，所以还是按照github上的说明将所有依赖都用conda安装好。Dfam-consortium/RepeatModeler: De-Novo Repeat Discovery Tool (github.com)

接下来在RepeatModeler安装的目录下运行perl ./configure，同样是一路回车到底确定路径，最后会询问是否需要预测LTR结构，因为我在之前的求LAI指数的博客中已经做过LTR预测，因此这一步选择n跳过，后续我会说明如何利用LTR预测数据：

因为我要注释的生物是非模式生物，在Dfam库和Repbase库中均没有该物种信息（无法在RepeatMasker软件中指定特定的物种，-species 和 -lib的参数是冲突的，需要自建数据库），因此注释所用的数据库将由以下三种数据库组成：

1. LTR_retriever整合的LTR预测数据库（见这篇博客）
2. 同源的（指该类群祖先和衍生节点）重复序列数据库
3. 使用RepeatModeler从头预测序列，训练该物种的重复序列模型，构建预测的重复序列数据库

需要注意这三种数据库都需要fasta格式，将三种数据库合并之后，使用RepeatMasker -lib指定自建数据库，预测TE序列。

4.1 导出同源物种重复序列库

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

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

1

python famdb.py -i Libraries/RepeatMaskerLib.h5 lineage -ad lamiids	# lamiids是我能查找到的最近的分支

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

1
2
3
4

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

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

4.2 RepeatModeler从头预测

新建一个目录，用于存放RepeatModeler的预测结果，写一个repeatmodeler.slurm脚本：

1
2
3
4
5
6
7

#!/bin/bash
#SBATCH -n 100
#SBATCH -t 7200

BuildDatabase -name luobuma -engine ncbi /public/home/wlxie/NextPolish/luobuma_rundir/genome.nextpolish.fasta	# 用基因组组装结果构建数据库

RepeatModeler -pa 25 -database luobuma -engine ncbi	# 自训练

RepeatModeler以自身基因组数据做训练集，用三种重复序列分析软件（ RECON, RepeatScout 和 LtrHarvest/Ltr_retriever）进行预测，最后给出de novo预测结果。需要i说明一下，程序结束之后会给出如下四个文件：

1. sample-families.fa de novo预测重复序列家族文件，也就是预测的重复序列库
2. sample-familes.stk Seed alignments
3. RM_123456.XXXXXXXXX 中间文件（记录每一轮训练的流程和结果，仅用于中间程序崩了以后可以识别并继续跑流程）
4. sample-rmod.log log文件

最终得到的luobuma-families.fa文件是我们需要的，里面记录了各种de novo预测的重复序列家族。中间文件具体有什么可以参考官方的github文档，这里仅仅是起到Recover from a failure的作用，中间程序没有崩就不用管它。

注意下RepeatModeler -pa参数，1 pa可以运行4个线程，我申请了100个核，这里就是25 pa可以用完所有资源。

这一步运行时间最久，100个核对200Mbp大小的植物基因组进行de novo预测重复序列，跑了17个小时。

4.3 整合数据库

将4.1、4.2步骤的结果，以及前面做的LTR预测结果进行整合（都是fasta格式）：

1

cat lamiids.fasta luobuma-families.fa luobuma.fasta.mod.LTRlib.fa > final_luobuma_repeat.fasta	# 合并同源数据库、RepeatModeler训练结果和LTR预测结果

此时得到的final_luobuma_repeat.fasta就是后一步运行RepeatMarsker需要指定的自建数据库。

4.4 RepeatMasker搜索重复序列

根据需求确定参数，写一个repeatmasker.slurm脚本：

1
2
3
4
5

#!/bin/bash
#SBATCH -n 100
#SBATCH -t 7200

RepeatMasker -nolow -no_is -pa 25 -lib final_luobuma_repeat.fasta -engine ncbi -gff -norna -dir luobuma /public/home/wlxie/NextPolish/luobuma_rundir/luobuma.fasta 

RepeatMasker的参数非常多，介绍一下这里用到的：

-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 指定输出目录

在输出目录下可以找到以下几种格式的文件:

sample.fasta.cat.gz 基因组和数据库中参考重复序列比对详情，i代表碱基转换，v代表碱基颠换

sample.fasta.masked 重复序列屏蔽我iN后的序列

sample.fasta.out 预测的重复序列详细信息，Smith-Waterman 算法得分等等

sample.fasta.out.gff 上一个文件的gff格式

sample.fasta.tbl RepeatMasker的结果报告

主要看一下结果报告：

TE的预测结果被分为逆转录转座子、DNA转座子和Unclassified三类，总的转座子序列数量和在基因组的占比见Total interspersed repeats统计结果。做到这里可以再结合前面做的TR分析，做一个基因组重复序列注释汇总表，我这里就不再演示了。