0%

手动安装RepeatMasker和RepeatModeler有点手疼,还好Dfam提供了基于docker的TE Tools。

TE Tools的官网:https://github.com/Dfam-consortium/TETools

使用docker,需要系统管理员将用户添加到docker用户组。

构建镜像

1
2
3
4
5
wget https://github.com/Dfam-consortium/TETools/archive/refs/tags/1.4.tar.gz
tar zxf 1.4.tar.gz
cd TETools-1.4/
./getsrc.sh
docker build -t dfam/tetools:1.4 .

构建成功后,可以使用docker image ls查看镜像。

构建镜像考验网速,可以下载我已经构建好的镜像。

1
2
3
wget https://www.biochen.org/public/software/tetools_1.4.tar.gz
gunzip tetools_1.4.tar.gz
docker load -i tetools_1.4.tar

在docker中使用RepeatMasker/RepeatModeler

1
docker run -it --rm --user $(id -u ${USER}):$(id -g ${USER}) -v /home/chenwen/data:/data dfam/tetools:1.4

参数说明

1
2
3
4
-it 为容器重新分配一个伪输入终端,以交互模式运行容器
--rm 容器退出后,自动清理容器内部的文件系统
--user $(id -u ${USER}):$(id -g ${USER}) 以当前用户和用户组运行docker
-v /home/chenwen/data:/data 挂载系统的/home/chenwen/data目录到容器的/data目录

在启动的新终端中,就可以运行RepeatMasker或者RepeatModeler,数据存放在宿主机的/home/chenwen/data目录中,对应到容器的/data目录。运行完之后,在终端输入exit退出容器。

RepeatMasker的用法参考之前的博文使用RepeatMasker屏蔽基因组重复序列

RepeatModeler的用法参考之前的博文使用RepeatModeler从头预测基因组重复序列

简介

官网:http://www.repeatmasker.org/RepeatModeler/

Repeatmasker 基于与已知的重复序列数据库比对来寻找重复序列,Repeatmodeler是通过重续序列的结构特征来进行从头注释,因此可以寻找一些物种特有的重复序列。

RepeatMasker依赖RepBase数据库和Dfam数据库来屏蔽重复序列,对于非模式生物来讲,这两个数据库覆盖有限。通常的做法就是RepeatModeler从头预测重复序列,然后将结果作为RepeatMasker的输入。

安装依赖

RepeatModeler依赖的RepeatMasker、TRF、RMBlast请参考我的上一篇博文使用RepeatMasker屏蔽基因组重复序列,接下来我们只安装剩下的依赖。

RECON - De Novo Repeat Finder

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz
tar zxf RECON-1.08.tar.gz
cd RECON-1.08/src/
make && make install
rm ~/software/RECON-1.08.tar.gz

RepeatScout - De Novo Repeat Finder

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatScout-1.0.6.tar.gz
tar zxf RepeatScout-1.0.6.tar.gz
cd RepeatScout-1.0.6/
make
rm ~/software/RepeatScout-1.0.6.tar.gz

运行LTR结构搜索,还需要安装如下软件。

LtrHarvest - The LtrHarvest program is part of the GenomeTools suite

1
2
3
4
cd ~/software/
wget http://genometools.org/pub/binary_distributions/gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz
tar zxf gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz
rm ~/software/gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz

CD-HIT - A sequence clustering package

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/weizhongli/cdhit/releases/download/V4.8.1/cd-hit-v4.8.1-2019-0228.tar.gz
tar zxf cd-hit-v4.8.1-2019-0228.tar.gz
cd cd-hit-v4.8.1-2019-0228/
make
cd cd-hit-auxtools/
make
rm ~/software/cd-hit-v4.8.1-2019-0228.tar.gz

Ltr_retriever - A LTR discovery post-processing and filtering tool

1
2
3
4
cd ~/software/
wget https://github.com/oushujun/LTR_retriever/archive/v2.9.0.tar.gz
tar zxf v2.9.0.tar.gz
rm ~/software/v2.9.0.tar.gz

LTR_retriever依赖TRF、BLAST+、CD-HIT、HMMER和RepeatMasker。CD-HIT的安装见本文,TRF、HMMER、RepeatMasker的安装见上一篇博文使用RepeatMasker屏蔽基因组重复序列,BLAST+的安装参考另外一篇博文本地BLAST

编辑~/software/LTR_retriever-2.9.0/目录下的paths文件,写入如下内容,指定这些依赖的安装路径。

1
2
3
4
BLAST+= /home/chenwen/software/ncbi-blast-2.12.0+/bin # a path that contains makeblastdb, blastn, blastx
RepeatMasker= /home/chenwen/software/RepeatMasker # a path that contains RepeatMasker
HMMER= /home/chenwen/software/hmmer-3.2.1/bin # a path that contains hmmsearch
CDHIT= /home/chenwen/software/cd-hit-v4.8.1-2019-0228 # a path that contains cd-hit-est (preferred). CDHIT and BLAST are replaceable

MAFFT - A multiple sequence alignment program.

1
2
3
cd ~/software/
wget https://mafft.cbrc.jp/alignment/software/mafft-7.487-with-extensions-src.tgz
tar zxf mafft-7.487-with-extensions-src.tgz

编辑~/software/mafft-7.487-with-extensions/core/目录下的Makefile文件,将PREFIX = /usr/local修改为PREFIX = /home/chenwen/software/mafft-7.487,然后继续执行如下命令。

1
2
3
4
5
cd ~/software/mafft-7.487-with-extensions/core/
make clean
make
make install
rm -rf ~/software/mafft-7.487-with-extensions*

Ninja - A tool for large-scale neighbor-joining phylogeny inference and clustering

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/TravisWheelerLab/NINJA/archive/0.95-cluster_only.tar.gz
tar zxf 0.95-cluster_only.tar.gz
cd NINJA-0.95-cluster_only/NINJA/
make all
cp Ninja ~/software/
rm -rf ~/software/NINJA-0.95-cluster_only
rm ~/software/0.95-cluster_only.tar.gz

UCSC TwoBit Tools

1
2
3
4
5
6
mkdir -p ~/software/kent/
cd ~/software/kent/
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitInfo
chmod +x *

安装和配置RepeatModeler

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatModeler/RepeatModeler-2.0.2a.tar.gz
tar zxf RepeatModeler-2.0.2a.tar.gz
rm ~/software/RepeatModeler-2.0.2a.tar.gz
cd RepeatModeler-2.0.2a/
perl ./configure

根据提示配置如下路径

1
2
3
4
5
6
7
8
9
10
11
12
REPEATMASKER_DIR /home/chenwen/software/RepeatMasker
RECON_DIR /home/chenwen/software/RECON-1.08/bin
RSCOUT_DIR /home/chenwen/software/RepeatScout-1.0.6
TRF_PRGM /home/chenwen/software/trf
CDHIT_DIR /home/chenwen/software/cd-hit-v4.8.1-2019-0228
UCSCTOOLS_DIR /home/chenwen/software/kent
RMBLAST_DIR /home/chenwen/software/rmblast-2.11.0/bin
ABBLAST_DIR /home/chenwen/software/ab-blast-20200317-linux-x64
GENOMETOOLS_DIR /home/chenwen/software/gt-1.6.2-Linux_x86_64-64bit-complete/bin
LTR_RETRIEVER_DIR /home/chenwen/software/LTR_retriever-2.9.0
MAFFT_DIR /home/chenwen/software/mafft-7.487/bin
NINJA_DIR /home/chenwen/software

如果配置成功,将会出现如下提示信息。

1
Congratulations!  RepeatModeler is now ready to use.

将RepeatModeler添加到环境变量

1
2
echo "export PATH=$HOME/software/RepeatModeler-2.0.2a:\$PATH" >> ~/.bashrc
source ~/.bashrc

运行RepeatModeler

这里我们以从Ensembl下载的鲤鱼基因组为例。

1
2
3
4
5
mkdir ~/test_RepeatModeler
cd ~/test_RepeatModeler/
wget http://ftp.ensembl.org/pub/release-104/fasta/cyprinus_carpio/dna/Cyprinus_carpio.common_carp_genome.dna.toplevel.fa.gz
gunzip Cyprinus_carpio.common_carp_genome.dna.toplevel.fa.gz Cyprinus_carpio.fa.gz
mv Cyprinus_carpio.common_carp_genome.dna.toplevel.fa Cyprinus_carpio.fa

建库

1
2
3
cd ~/test_RepeatModeler/
mkdir db
BuildDatabase -name db/Cyprinus_carpio Cyprinus_carpio.fa

参数说明

1
-name 库的名字

运行RepeatModeler

1
2
cd ~/test_RepeatModeler/
nohup RepeatModeler -database db/Cyprinus_carpio -pa 12 -LTRStruct >& run.out &

参数说明

1
2
3
-database 库的名字,与上一步一致
-pa 线程数
-LTRStruct 开启LTR结构搜索

RepeatModeler的搜索引擎已经只支持RMBlast了,如果使用-engine abblast,RepeatModeler会报出警告WARNING: "-engine abblast" is deprecated, this verison of RepeatModeler uses rmblast only.,然后继续使用RMBlast。

结果解读

我们查看一下日志文件run.out的末尾。

鲤鱼基因组大小约1.7G,RepeatModeler在我的i5 10400(6核12线程)电脑上面跑了约78小时。

RepeatModeler会生成一个名字像这样的文件夹RM_1208021.MonAug20735032021,这是工作目录,保存了中间结果,可以删除。主要结果有两个,db/Cyprinus_carpio-families.fa可以作为library用于RepeatMasker进行重复序列的屏蔽。db/Cyprinus_carpio-families.stk是Dfam兼容的Stockholm格式,可以上传到Dfam数据库。

运行RepeatMasker

我们还需要利用RepeatModeler生成的Cyprinus_carpio-families.fa,使用RepeatMasker屏蔽基因组重复序列。

1
2
cd ~/test_RepeatModeler/
RepeatMasker Cyprinus_carpio.fa -lib db/Cyprinus_carpio-families.fa -e rmblast -xsmall -s -gff -pa 12

使用了-lib参数,会使-species参数失效,不需要同时指定-species参数了。使用-lib db/Cyprinus_carpio-families.fa屏蔽了40.02%的重复序列,而使用-species 'Cyprinus carpio'也能屏蔽39.62%的重复序列。二者相差不大,可能是Dfam数据库的重复序列已经覆盖到了鲤鱼。有些非模式生物,单用RepeatMasker只能屏蔽约10%的重复序列,这时使用RepeatModeler就非常有必要了。

简介

复杂的基因组大部分是可转座元件和其他重复序列,在全基因组比对或者基因注释之前需要屏蔽重复序列,否则就会浪费计算资源,也会得到许多无意义的结果。

RepeatMasker(官网:http://www.repeatmasker.org/RMDownload.html)是一款识别基因组中重复序列和低复杂度序列的软件。RepeatMasker依赖序列搜索引擎(HMMER、Cross_Match、ABBlast或者RMBlast)和重复序列数据库(Dfam和Repbase)来屏蔽重复序列。

对于非模式生物,这两个数据库覆盖有限,因此需要先使用RepeatModeler从头预测重复序列,然后将RepeatModeler预测的序列作为RepeatMasker的输入,来完成最终的重复序列屏蔽。更多的信息请访问RepeatMasker的官网http://www.repeatmasker.org/RepeatModeler/

安装依赖

本文使用RepeatMasker-4.1.2-p1,它需要以下依赖:

  1. perl 5.8.0或者更高的版本
  2. Python 3和h5py库
  3. 序列搜索引擎
  4. TRF - Tandem Repeat Finder
  5. 重复序列数据库

我的HOME路径是/home/chenwen/,软件安装都安装在/home/chenwen/software/,请记得修改为适合自己的路径。

RepeatMasker支持序列4种搜索引擎,我们没有必要全都装上,只安装RMBlast也可以。

安装Cross_Match
按照官网http://www.phrap.org/说明,使用学术机构的邮箱申请软件。申请通过后,软件通过邮件发送给我们。这里提供我的副本phrap_cross_match_swat_1.090518.zip。

1
2
3
4
5
6
cd ~/software/
wget https://www.biochen.org/public/software/phrap_cross_match_swat_1.090518.zip
unzip phrap_cross_match_swat_1.090518.zip
cd phrap_cross_match_swat_1.090518/
make
rm ~/software/phrap_cross_match_swat_1.090518.zip

安装RMBlast
官网http://www.repeatmasker.org/RMBlast.html

1
2
3
4
cd ~/software/
wget http://www.repeatmasker.org/rmblast-2.11.0+-x64-linux.tar.gz
tar zxf rmblast-2.11.0+-x64-linux.tar.gz
rm rmblast-2.11.0+-x64-linux.tar.gz

安装HMMER
官网http://hmmer.org/,RepeatMasker官方推荐使用HMMER v3.2.1,而不是其他版本。

1
2
3
4
5
6
7
cd ~/
wget http://eddylab.org/software/hmmer/hmmer-3.2.1.tar.gz
tar zxf hmmer-3.2.1.tar.gz
cd hmmer-3.2.1/
./configure --prefix=$HOME/software/hmmer-3.2.1
make && make install
rm ~/hmmer-3.2.1* -rf

安装ABBlast
http://blast.advbiocomp.com/licensing/填写表格后,软件和使用许可通过邮件发送给我们。ABBlast解压即可用,使用许可license.xml需要放置到~/.config/ab-blast目录下。

1
2
3
4
5
6
7
8
cd ~/software/
wget https://www.biochen.org/public/software/ab-blast/ab-blast-20200317-linux-x64.tar.gz
tar zxf ab-blast-20200317-linux-x64.tar.gz
mkdir ~/.config/ab-blast
cd ~/.config/ab-blast/
wget https://www.biochen.org/public/software/ab-blast/license.xml
chmod 600 ~/.config/ab-blast/license.xml
rm ~/software/ab-blast-20200317-linux-x64.tar.gz

安装trf
官网https://github.com/Benson-Genomics-Lab/TRF

1
2
3
4
cd ~/software/
wget https://github.com/Benson-Genomics-Lab/TRF/releases/download/v4.09.1/trf409.linux64
mv trf409.linux64 trf
chmod 755 trf

安装和配置RepeatMasker

下载RepeatMasker

1
2
3
4
cd ~/software/
wget https://www.repeatmasker.org/RepeatMasker/RepeatMasker-4.1.2-p1.tar.gz
tar zxf RepeatMasker-4.1.2-p1.tar.gz
rm ~/software/RepeatMasker-4.1.2-p1.tar.gz

下载RepBase数据库

1
2
3
4
cd ~/software/RepeatMasker/
wget https://www.biochen.org/public/software/RepBaseRepeatMaskerEdition-20181026.tar.gz
tar zxf RepBaseRepeatMaskerEdition-20181026.tar.gz
rm ~/software/RepeatMasker/RepBaseRepeatMaskerEdition-20181026.tar.gz

更新Dfam数据库

RepeatMasker带一个小的Dfam数据库,我们需要重新下载完整版的数据库。解压后的Dfam.h5文件83G,非常大。

1
2
3
cd ~/software/RepeatMasker/Libraries/
wget https://www.dfam.org/releases/Dfam_3.3/families/Dfam.h5.gz
gunzip Dfam.h5.gz

安装h5py

1
pip install h5py

配置RepeatMasker

1
2
cd ~/software/RepeatMasker/
perl ./configure

配置过程中,必须配置trf路径,4种搜索引擎的路径可以选择配置。

1
2
3
4
5
TRF_PRGM: /home/chenwen/software/trf
CROSSMATCH_DIR: /home/chenwen/software/phrap_cross_match_swat_1.090518
RMBLAST_DIR: /home/chenwen/software/rmblast-2.11.0/bin
HMMER_DIR: /home/chenwen/software/hmmer-3.2.1/bin
ABBLAST_DIR: /home/chenwen/software/ab-blast-20200317-linux-x64

如果配置成功,将有如下提示。

配置成功之后,删掉巨大的Dfam.h5文件。

1
rm ~/software/RepeatMasker/Libraries/Dfam.h5

将RepeatMasker添加到环境变量

1
2
echo "export PATH=$HOME/software/RepeatMasker:\$PATH" >> ~/.bashrc
source ~/.bashrc

运行RepeatMasker

这里我们以从Ensembl下载的斑马鱼基因组为例。

1
2
3
4
5
mkdir ~/test_RepeatMasker
cd ~/test_RepeatMasker/
wget http://ftp.ensembl.org/pub/release-104/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
mv Danio_rerio.GRCz11.dna.primary_assembly.fa Danio_rerio.fa
1
RepeatMasker Danio_rerio.fa -species "Danio rerio" -e rmblast -xsmall -s -gff -pa 12

参数说明

1
2
3
4
5
6
7
8
-e rmblast 指定搜索引擎为rmblast,还可以选择crossmatch、abblast或者hmmer
-species "Danio rerio" 指定物种为斑马鱼,物种名必须在NCBI物种分类数据库(Taxonomy)中,推荐使用拉丁学名
-pa 12 使用12线程运行RepeatMasker
-s 慢速搜索,灵敏度比高0~5%,速度慢2-3倍
-xsmall 软屏蔽,将重复序列变成小写字母,而不是N或者X。
-gff 输出GFF文件
-nolow 忽略低重复的序列
-norna 忽略Small RNA

结果解读

斑马鱼基因组大小约1.4G,上述命令在我的i5 10400(6核12线程)电脑上面跑了约10.5小时。结果文件主要有:Danio_rerio.fa.tblDanio_rerio.fa.maskedDanio_rerio.fa.outDanio_rerio.fa.out.gffDanio_rerio.fa.cat.gz.cat.gz是重复序列与基因组序列的比对文件,其他文件我们打开看一下。

.tbl文件是重复序列的统计信息,我们可以看到斑马鱼有60%以上的重复序列。

.masked文件是屏蔽重复序列后的基因组文件,这里是软屏蔽,重复序列用小写字母表示。

.out文件是RepeatMasker默认的输出结果,.out.gff文件是相应的GFF文件。

可能让大家失望了,我的主力电脑不是服务器平台,也不是HEDT(high end desktop)平台。得益于技术进步,最普通的家用电脑,也能很好的完成我的生物信息分析。

配置单

1
2
3
4
5
6
CPU:英特尔 i5 10400
主板:华擎B460M Pro4
内存:光威 DDR4 2666 32G * 4 = 128G
硬盘:西部数据 SN550 1TB,西部数据紫盘 4T
电源:金河田 金牌 500W
机箱:先马 小黑洞

这个配置的价格约6000元,性价比非常高。

使用体验

盒装CPU自带的散热器,全铝鳍片,没有导热铜管,下压式而不是侧吹式,我担心因为散热的问题不能发挥出CPU的全部性能。实际测试,盒装自带的散热器完全够用,i5-10400能跑满全核心睿频4.0GHz,温度在62度左右。

双通道内存有15%左右的性能提升,多线程有12%左右的性能提升。

查看CPU运行频率

1
2
sudo apt install cpufrequtils
cpufreq-info |grep 'current CPU frequency'

查看CPU温度

1
2
sudo apt install lm-sensors
sensors

消费级CPU主频高、单核性能好,再配上NVMe固态硬盘,大部分时候的使用体验优于组里的戴尔T640服务器。目前消费级平台最大只支持128G内存,HEDT平台一般可以支持到256G内存,需要更高的内存的话,就需要上服务器平台了。128G内存除了跑大转录组和大基因组的从头组装不行外,其他任务都能胜任。一块4T硬盘不够用,之后又加了一块4T硬盘,以后还需要再加。所有的数据都不能删,硬盘存储是最大的压力。

0-based和1-based

生物信息文件格式中有很多格式是基于基因组坐标的,比如常见的BED格式或者GTF格式。然而对于对标系的定义,这两者有着截然的区别。BED格式第一个位置的下标是0,区间前开后闭;而GTF格式第一个位置的下标是1,区间都是闭的。不妨我们称前者为0-based,后者为1-based。0-based的优点是长度的计算很简单,直接相减就可以得到序列的长度;而1-based的优点是比较直观。

除了BED格式和GTF格式,下表列举了其他格式的情况。

长度计算

Length(0-based) = End(0-based) - Start(0-based)
Length(1-based) = End(1-based) - Start(1-based) + 1

坐标转换

0-based转1-based
Start(1-based) = Start(0-based) + 1
End(1-based) = End(0-based)

1-based转0-based
Start(0-based) = Start(1-based) - 1
End(0-based) = End(1-based)