植物基因组

TEsorter:TE的分类和可视化

2023-07-27  本文已影响0人  wo_monic

来源于:https://github.com/zhangrengang/TEsorter
TE的鉴定可以使用LTRreviewer 或/和 repeatMasker,此处是对鉴定的TE进行分类和可视化。
实质是对TE的进一步分类,因为目前的TE分类比较粗浅Copia和Gypsy,实质上可以分的更加详细。
TEsorter只是分类器,其实是基于dante开发的。分类的依据数据库可以是REXdb或GyDB,REXdb的分类更加详细,所以TEsorter作者推荐使用这个库。

dante本身就是用于基于domain的转座子注释工具。dante github

不同数据库对Copia的分类,来源于REXdb论文

image.png
不同数据库对Gypsy的分类
image.png

Copia包括这些小类Ale , AlesiaAlesia, AngelaAngela, BiancaBianca, BrycoBryco, LycoLyco, Gymco I–IV, IkerosIkeros, IvanaIvana, OsserOsser, SIRESIRE, TARTAR, and Tork clades and the
Gypsy包括CRM, ChlamyvirChlamyvir, GaladrielGaladriel, Tcn1Tcn1, ReinaReina, TekayTekay, AthilaAthila, Tat I–III, OgreOgre, RetandRetand, PhygyPhygy, and Selgy clades.

TEsorter把Copia和Gypsy分类为如下:

提取的时候注意使用的数据库不同关键词略有区别rexdb和gydb的域名有些不同:PROT(rexdb)=AP(gydb),RH(rexdb)=RNaseH(gydb)

安装TEsorter

git clone https://github.com/zhangrengang/TEsorter
cd TEsorter
python setup.py install

直接使用conda安装也可以conda install -c bioconda tesorter

测试TEsorter

TEsorter TEsorter/test/rice6.9.5.liban

使用参数讲解

从基因组中提取 TE 序列用于 TEsorter

当您只有基因组序列时,以下是从广泛使用的软件的输出中提取 TE 序列的示例。

  1. 从RepeatMasker输出中提取所有 TE 序列:
#基因组RepeatMasker的下游分析,TE的分类

#指定RepeatMasker的输出文件所在的路径里的基因组.out文件
genome_out="/share/home/repeatsequence/Repeat_all_result/C.s.genome.fa.out"
#基因组位置
genome="/share/home/repeatsequence/Repeat_all_result/C.s.genome.fa"
#指定输出的前缀
prefix="C.s"
pwd="/share/home/TEsort"

cd ${pwd}
if ! [ -e Genome_TEsort ];then
    mkdir Genome_TEsort
fi
cd Genome_TEsort

# run RepeatMasker, which will generate a *.out file.
#RepeatMasker [options] genome.fa

# extract sequences
RepeatMasker.py out2seqs $genome_out $genome > ${prefix}.whole_genome_te.fa

#获取分类信息(计算非常快)
TEsorter ${prefix}.whole_genome_te.fa -db rexdb-plant -p 20

#提取domains,此处是提取所有的。 也可以只提取一种,例如GAG 或PROT,此处只提取保守的三种RH,RT,INT
concatenate_domains.py ${prefix}.rexdb-plant.cls.pep INT > ${prefix}.rexdb-plant.cls.pep.INT.aln
concatenate_domains.py ${prefix}.rexdb-plant.cls.pep TPase > ${prefix}.rexdb-plant.cls.pep.TPase.aln
cat ${prefix}.rexdb-plant.cls.pep.INT.aln ${prefix}.rexdb-plant.cls.pep.TPase.aln > ${prefix}.rexdb-plant.cls.pep.INT_TPase.faa
mafft --auto ${prefix}.rexdb-plant.cls.pep.INT_TPase.faa > ${prefix}.rexdb-plant.cls.pep.INT_TPase.aln

#构建进化树
iqtree2 -s ${prefix}.rexdb-plant.cls.pep.INT_TPase.aln -B 1000 -T 20 -m MFP --bnni
#进化树的可视化(可视化可能会失败)
LTR_tree.R ${prefix}.rexdb-plant.cls.pep.INT_TPase.aln.treefile ${prefix}.rexdb-plant.cls.tsv ${prefix}.rexdb-plant.cls.pep.INT_TPase.aln.pdf
  1. 从LTR_retriever输出中提取所有完整的 LTR-RTs 序列:
run LTR_retriever, which generate two *.pass.list files.

LTR_retriever -genome genome.fa [options]
自动化脚本修改前3行的参数即可运行

#用于LTR-RTs的下游分析,TE的分类
LTR-RTs_genome="/share/home/LTR/C.s.genome.chr.fa" #指定LTRreviewer的输出文件所在的路径里的基因组文件
prefix="C.s" #指定输出的前缀

if ! [ -e TEsort ];then
    mkdir TEsort
fi
cd TEsort

LTR_retriever.py get_full_seqs ${LTR-RTs_genome} >${prefix}.ltr.fa
#获取LTR-RTs的分类信息(计算非常快)
TEsorter ${LTR-RTs_genome}.ltr.fa -db rexdb-plant -p 20

#提取domains,此处是提取所有的。 也可以只提取一种,例如GAG 或PROT或者RH
concatenate_domains.py ${prefix}.ltr.fa.rexdb-plant.cls.pep RH RT INT > ${prefix}.ltr.fa.rexdb-plant.cls.pep.full.aln
#构建进化树
iqtree2 -s ${prefix}.ltr.fa.rexdb-plant.cls.pep.full.aln -B 1000 -T 20 -m MFP --bnni

#进化树的可视化(可视化可能会失败)
LTR_tree.R ${prefix}.ltr.fa.rexdb-plant.cls.pep.full.aln.treefile ${prefix}.ltr.fa.rexdb-plant.cls.tsv ${prefix}.ltr.fa.rexdb-plant.cls.pep.full.aln.treefile.pdf

domains的不同种类的保守型不一致,RH RT INT的保守性比较好,所以一般用来构建进化树。GAG和PROT的保守性就没有前3者的稳定。来源于TEsorter和REXdb论文.


支持对 TE 多蛋白序列(示例)或基因蛋白序列进行分类:

TEsorter RepeatPeps.lib -st prot -p 20

从1.4版本可以直接使用基因组来鉴定

TEsorter genome.fasta -genome -p 20

输出结果

rice6.9.5.liban.rexdb.domtbl        HMMScan raw output
rice6.9.5.liban.rexdb.dom.faa       protein sequences of domain, which can be used for phylogenetic analysis.
rice6.9.5.liban.rexdb.dom.tsv       inner domains of TEs/LTR-RTs, which might be used to filter domains based on their scores and coverages.
rice6.9.5.liban.rexdb.dom.gff3      domain annotations in `gff3` format
rice6.9.5.liban.rexdb.cls.tsv       TEs/LTR-RTs classifications
    Column 1: raw id
    Column 2: Order, e.g. LTR
    Column 3: Superfamily, e.g. Copia
    Column 4: Clade, e.g. SIRE
    Column 5: Complete, "yes" means one LTR Copia/Gypsy element with full GAG-POL domains.
    Column 6: Strand, + or - or ?
    Column 7: Domains, e.g. GAG|SIRE PROT|SIRE INT|SIRE RT|SIRE RH|SIRE; `none` for pass-2 classifications
rice6.9.5.liban.rexdb.cls.lib      fasta library for RepeatMasker
rice6.9.5.liban.rexdb.cls.pep    the same sequences as `rice6.9.5.liban.rexdb.dom.faa`, but id is changed with classifications.

根据需要提取对应的domain来进行比对后构建进化树

提取的时候注意使用的数据库不同关键词略有区别rexdb和gydb的域名有些不同:PROT(rexdb)=AP(gydb),RH(rexdb)=RNaseH(gydb)
我此处使用的是rexdb,所以提取所有domain的时候使用的是PROT和RH

上一篇 下一篇

猜你喜欢

热点阅读