基因组注释01# | 重复序列注释（RepeatModler+RepeatMasker）
写在前面
进行这部分重复序列提取和注释，主要是为了屏蔽基因组的重复序列，以保证后面基因结构注释的准确性。
大致的流程如下：
1.首先使用RepeatModler对组装好的基因组进行重复序列预测，从头构建重复序列本地库(de novo repeat library)。
2.随后用RepeatMasker对重复序列进行比对注释。注释的流程又分两步：
- 基于dfam和repbase重复序列数据库对基因组进行近缘物种（species）的同源注释，拿到第一轮注释结果以及用于下一步的对重复序列屏蔽过的基因组
.fasta.masked; - 基于RepeatModler拿到的de novo重复序列本地库进行注释，但这次输入的基因组是上面的经过重复序列屏蔽过的
.fasta.masked，这样我们又得到第二轮注释结果。
3.后面就是合并两部分注释结果，生成我们需要一些文件用于后续分析，如.out, .align, .gff, .masked.fasta等。
数据准备
只需要组装好的基因组文件genome.fasta
##选择自己需要的基因组下载或上传服务器，我这里从NGDC数据库下载基因组序列文件，并解压出来
mkdir RepeatAnnoWorkdir
cd RepeatAnnoWorkdir
wget https://download.cncb.ac.cn/gwh/Animals/Silurus_lanzhouensis_Sl_YY_GWHERQF00000000/GWHERQF00000000.genome.fasta.gz
guzip GWHERQF00000000.genome.fasta.gz
软件
TETools主页：https://github.com/Dfam-consortium/TETools
这里用一个官方推荐的容器TETools，里面包含了RepeatModler、RepeatMasker等软件以及他们的诸多依赖。唯一比较麻烦的是这个容器里只有Dfam数据库而不包含RepBase重复序列数据库，因为RepBase是收费的。因此需要重新自定义RepeatMasker的数据库famDB,把这两部分数据库都添加进去。该容器的gihub网页上有相应的步骤。
GitHub上用的Docker，我这里用的singularity，根据其语法改了下，具体如下：
##1.下载容器，可修改容器名字，用singularity pull 去拉取Docker Hub里的容器
singularity pull dfam-tetools-latest .sif docker://dfam/tetools:latest
##2. 从容器中复制 RepeatMasker/Libraries/ 目录到宿主机文件系统
singularity shell --bind $(pwd):/work dfam-tetools-latest.sif # 启动交互式Singularity容器并挂载当前目录
ls /work # 确认挂载的当前目录可见
cp -r /opt/RepeatMasker/Libraries/ /work/ # 复制容器内的Libraries/到宿主机的当前目录
exit # 退出Singularity容器的交互模式
chown -R $USER ./Libraries/ # 修改宿主机上的Libraries/ 权限
##3. 添加Dfam数据库到宿主机的Libraries/famdb/里
Dfam数据库有一个必需的主分区partition 0，还有其他可选分区partition 1-16。里面包含生物界各个物种分类的重复序列数据库(https://www.dfam.org/releases/current/families/FamDB/README.txt)，根据自己研究物种进行选择，我这里除了partition 0，还选了partition 12 （脊椎动物Vertebrata）。
wget https://www.dfam.org/releases/current/families/FamDB/dfam39_full.0.h5.gz # 下载主分区
wget https://www.dfam.org/releases/current/families/FamDB/dfam39_full.12.h5.gz # 下载可选分区
gunzip -c dfam39_full.0.h5.gz > Libraries/famdb/dfam39_full.0.h5 # 解压到Libraries/famdb/目录里
gunzip -c dfam39_full.12.h5.gz > Libraries/famdb/dfam39_full.12.h5
singularity exec --bind ./Libraries:/opt/RepeatMasker/Libraries dfam-tetools-latest.sif /
bash -c "cd /opt/RepeatMasker && ./tetoolsDfamUpdate.pl" # 运行库更新脚本，让RepeatMasker识别这些文件并重新配置rmlib.config文件
##4. 整合RepBase数据库到宿主机的Libraries/famdb/里
RepBase数据库现在是收费的，但我在GitHub上找到一个别人上传的版本RepBaseRepeatMaskerEdition-20181026.tar.gz（https://github.com/yjx1217/RMRB)，将其下载上传到服务器上。
将RepBase数据库解压到Libraries/目录时，注意解压位置，解压出来的两个文件会被自动解压到Libraries/目录，这点一样要注意，一定得在Libraries/前一级解压
./Libraries/ # 主库目录
├───── famdb/ # Dfam库的h5数据
│ ├── dfam39_full.0.h5
│ ├── dfam39_full.12.h5
│ └── rmlib.config
├── RMRBSeqs.embl # RepBase库的EMBL格式数据
├── README.RMRBSeq
└── RepeatMasker.lib # 主库文件
tar -xvzf RepBaseRepeatMaskerEdition-20181026.tar.gz # 解压RepBase库
singularity exec --bind $(pwd)/Libraries:/opt/RepeatMasker/Libraries dfam-tetools-latest.sif /
bash -c "cd /opt/RepeatMasker && echo '/opt/RepeatMasker/Libraries' | ./tetoolsDfamUpdate.pl" #重新更新库
##5. 设置环境变量
# 将以下内容添加到 ~/.bashrc 或 ~/.bash_profile
export LIBDIR="/home/ug2117/wangtao_study/annotation/GenomeAnn_genek/container1/Libraries"
export SINGULARITY_BIND="$LIBDIR:/opt/RepeatMasker/Libraries"
由于以上将repbase库整合到TETools容器的内容是我后面弄的，在这之前我一直用的genek课程提供的容器TETools202309.sif。所以下面的操作都是用的这个容器。
RepeatModler预测重复序列，构建de novo repeat library
RepeatModler2的工作流程如下图，先通过RepeatScout/RECON进行1-6轮的denovo检测，再利用ltr_retriever进行LTR结构预测。前面6轮不一定全跑完，视物种的大小和复杂度而定。
# 1.建库BuildDatabase
singularity exec ../container/TETools202309.sif BuildDatabase /
-name SlYY / # 数据库名称
GWHERQF00000000.genome.fasta # 需要预测的基因组
# 2.运行RepeatModeler
singularity exec ../container/TETools202309.sif RepeatModeler /
-database SlYY /
-threads 50 / # 线程数
-LTRStruct / # 运行LTR鉴定
1>repeatmodeler.log 2>repeatmodeler.err
结果如下：
.fa文件就包含了预测到的所有重复序列，即de novo repeat library，另外我们这里RepeatModeler只跑了5轮。
RepeatMasker注释重复序列
1.基于dfam和repbase数据库进行第一轮同源注释
singularity exec ../container/TETools202309.sif RepeatMasker /
-e rmblast / # 比对软件
-pa 10 / # 并行任务数，实际占用线程为pa*4
-gff / # 输出gff
-s / # 低速模式，0-5%敏感度，比默认模式慢2-3倍
-a / # 输出align文件
-species Teleostei / # 物种范围，我这个物种是鱼，所以设置了硬骨鱼类
-dir ./repeatmasker_out_species / # 第一轮注释结果目录
GWHERQF00000000.genome.fasta / # 需要注释的基因组
1>repeatmasker_spe.log 2>&repeatmasker_spe.err
repeatmasker_out_species目录下的输出文件及.tbl文件的统计结果如下：
2.基于de novo repeat library进行第二轮从头注释
singularity exec ../container/TETools202309.sif RepeatMasker /
-e rmblast /
-pa 10 /
-gff -s -a /
-lib SlYY-families.fa / # 基于自己构建的de novo repeat library（SlYY-families.fa）进行注释
-dir ./repeatmasker_out_denovo / # 第二轮注释结果目录
repeatmasker_out_species/GWHERQF00000000.genome.fasta.masked / # 输入为上一轮注释输出的masked fasta
1>repeatmasker_denovo.log 2>repeatmasker_denovo.err
repeatmasker_out_denovo目录下的输出文件及.tbl文件的统计结果如下：
3.合并两轮注释结果
combineRMFiles.pl脚本只能合并两轮注释结果中的.out和.align文件，其他文件分别用其他几个脚本合并。
# 构建合并目录
mkdir repeatmaskerout_combine
# 合并
singularity exec ../container/TETools202309.sif /
perl /opt/RepeatMasker/util/combineRMFiles.pl / #combineRMFiles.pl脚本合并，注意使用容器内路径
./repeatmasker_out_species/GWHERQF00000000.genome.fasta / # 第一轮结果前缀
./repeatmasker_out_denovo/GWHERQF00000000.genome.fasta.masked / # 第二轮结果前缀
repeatmaskerout_combine/SlYYgenome # 合并结果附加前缀
合并输出.out和.align文件：
4.重新生成.gff重复序列注释文件
rmOutToGFF3.pl脚本将合并目录中的.out文件转换成.gff注释文件。
singularity exec ../container/TETools202309.sif rmOutToGFF3.pl /
repeatmaskerout_combine/SlYYgenome.out / # 合并结果中的.out文件作为输入
> repeatmaskerout_combine/SlYYgenome.out.gff # 输出到合并目录，命名.out.gff以区分前面两轮.gff文件
输出.out.gff文件：
5.重新生成屏蔽过重复序列的基因组注释文件.masked.fasta
屏蔽基因组重复序列有两种方式：
- hard mask：就是把基因组中的重复序列都转换成字母N；
- soft mask：就是把基因组中的重复序列转换成小写字母（一般基因组文件中的碱基序列都是大写字母）；
两种屏蔽方式都是用maskFile.pl脚本来转换，以原始的基因组文件为基础，利用合并目录中的.out文件里的重复序列信息将其转换成.softmask.fasta或.hardmask.fasta。两者的区别就是加不加-softmask参数，默认是hard mask，屏蔽过的基因组序列用于后续基因结构注释。
## soft mask
singularity exec ../container/TETools202309.sif maskFile.pl /
-fasta GWHERQF00000000.genome.fasta / # 原始的基因组文件作为基础
-annotations repeatmaskerout_combine/SlYYgenome.out / # 合并结果中的.out文件作为输入
-softmask # 该参数用于指定屏蔽方式
mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.softmask.fasta #默认会在当前文件夹中生成屏蔽过的序列文件.fasta.masked，将其挪到合并目录里并修改序列名
## hard mask
singularity exec ../container/TETools202309.sif maskFile.pl /
-fasta GWHERQF00000000.genome.fasta / # 原始的基因组文件作为基础
-annotations repeatmaskerout_combine/SlYYgenome.out / # 合并结果中的.out文件作为输入
mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.hardmask.fasta #默认会在当前文件夹中生成屏蔽过的序列文件.fasta.masked，将其挪到合并目录里并修改序列名
输出两种.masked.fasta文件：
注意
至此，我们拿到了用于后续基因结构注释所需的
.masked.fasta，但是这个文件是将所有的重复序列都屏蔽了，包括其中的Low_complexity和Simple_repeat重复元件，这些元件在基因的UTR，Intro，甚至Exon区域都有可能出现，尽管概率不高，但是屏蔽这些序列势必影响我们的基因结构注释时的比对识别。因此，我们在屏蔽重复序列的时候，选择性地不屏蔽Low_complexity和Simple_repeat元件。通过
grep命令滤掉合并目录中的.out文件里的Low_complexity和Simple_repeat行，然后再进行上述的mask步骤，那样输出的.masked.fasta文件就保留了这两种重复元件：## 1. 重新生成.out文件 grep -v "Low_complexity" ./repeatmaskerout_combine/SlYYgenome.out | / # grep -v过滤Low_complexity行 grep -v "Simple_repeat" / # grep -v进一步过滤Simple_repeat行 >./repeatmaskerout_combine/SlYYgenome.filter.out # 重新输出不带Low_complexity和Simple_repeat行的.filter.out文件 ## 2. 重新mask基因组 singularity exec ../container/TETools202309.sif maskFile.pl / -fasta GWHERQF00000000.genome.fasta / -annotations repeatmaskerout_combine/SlYYgenome.filter.out / # 重新生成的.filter.out文件作为输入 mv ./GWHERQF00000000.genome.fasta.masked repeatmaskerout_combine/SlYYgenome.filter_hardmask.fasta最后生成的
.filter_hardmask.fasta基因组序列文件即可用于后续的基因结构注释。
6.最后基于木村距离绘制重复序列景观图
首先，使用calcDivergenceFromAlign.pl脚本将合并目录中的.align文件里每个重复元件的木村距离信息Kimura (with divCpGMod)调出来，形成.div文件。
随后，使用createRepeatLandscape.pl脚本将进行图片绘制重复序列统计图。
## 1. 调取木村距离
singularity exec ../container/TETools202309.sif calcDivergenceFromAlign.pl /
-s SlYY-repeat.div / # 输出div文件
repeatmaskerout_combine/SlYYgenome.align # 输入合并目录里的.align文件
## 2. 作图
singularity exec ../container/TETools202309.sif createRepeatLandscape.pl /
-div SlYY-repeat.div / # 输入div文件
-g 776800000 / # 基因组大小
> SlYY-repeat.div.html # 输出html结果
html 结果：
写在后面
最后看一下重复序列注释输出的其他文件格式吧.align、 .out、 .gff ：
共有 0 条评论