基因组注释(2) |
您所在的位置:网站首页 › 基因序列注释怎么写的 › 基因组注释(2) |
前面我们注释了串联重复序列(Tandem repeat,TR),接下来是对散在重复序列(也称转座子,transposable element,TE)进行注释。注释之后我们对所有重复序列在基因组上进行屏蔽,就可以进行后面的结构基因预测和注释了。 1. 散在重复序列散在重复序列可以分为反转录转座子(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的时候要改下下参数设置。 2. RepeatModeler和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安装: 1conda install -c bioconda perl-devel-size最后需要修改一下环境变量(不修改运行的时候找不到pm文件),将RepeatMasker 安装路径添加到PERL5LIB环境变量中: 12# 打开 ~/.bashrcexport 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的参数是冲突的,需要自建数据库),因此注释所用的数据库将由以下三种数据库组成: LTR_retriever整合的LTR预测数据库(见这篇博客) 同源的(指该类群祖先和衍生节点)重复序列数据库 使用RepeatModeler从头预测序列,训练该物种的重复序列模型,构建预测的重复序列数据库需要注意这三种数据库都需要fasta格式,将三种数据库合并之后,使用RepeatMasker -lib指定自建数据库,预测TE序列。 4. 注释流程4.1 导出同源物种重复序列库前面2.1步骤将Repbase和Dfam数据库整合之后,RepeatMasker/Libraries目录下RepeatMaskerLib.h5这个文件为整合后构建的数据库文件,我们要在这个文件中导出同源物种的重复序列。 在RepeatMasker目录下提供了famdb.py这个程序查询目标近缘物种。如果你不知道自己的物种在什么分支上,我这里推荐一个查找已发表的植物基因组的网站Published Plant Genomes (plabipd.de),可以一级一级查看哪些近缘物种有人做过了。用以下命令查看物种重复序列否收录到库中: 1python famdb.py -i Libraries/RepeatMaskerLib.h5 lineage -ad lamiids # lamiids是我能查找到的最近的分支找到最近的分支后,导出最近分支的祖先节点和衍生节点物种的重复序列库,使用内置的perl软件转换成fasta格式: 1234python famdb.py -i Libraries/RepeatMaskerLib.h5 families -f embl -a -d lamiids > lamiids.embl # 查找近缘物种及其上祖先节点,其下所有类群repeat famlies,输出格式embl。 -a ancestor,-d descendentbuildRMLibFromEMBL.pl lamiids.embl > lamiids.fasta # 转换格式为fasta,方便后续合并 4.2 RepeatModeler从头预测新建一个目录,用于存放RepeatModeler的预测结果,写一个repeatmodeler.slurm脚本: 1234567#!/bin/bash#SBATCH -n 100#SBATCH -t 7200BuildDatabase -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说明一下,程序结束之后会给出如下四个文件: sample-families.fa de novo预测重复序列家族文件,也就是预测的重复序列库 sample-familes.stk Seed alignments RM_123456.XXXXXXXXX 中间文件(记录每一轮训练的流程和结果,仅用于中间程序崩了以后可以识别并继续跑流程) 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格式): 1cat 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脚本: 12345#!/bin/bash#SBATCH -n 100#SBATCH -t 7200RepeatMasker -nolow -no_is -pa 25 -lib final_luobuma_repeat.fasta -engine ncbi -gff -norna -dir luobuma /public/home/wlxie/NextPolish/luobuma_rundir/luobuma.fastaRepeatMasker的参数非常多,介绍一下这里用到的: -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分析,做一个基因组重复序列注释汇总表,我这里就不再演示了。 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |