EDTAの出力データで色々やりたい
EDTAでトランスポゾンを検出したあとに色々やりたいが、それには細々した工夫が必要だったので忘備録としてまとめます。間違いなどあればご指摘いただけると幸いです。
※事前に解析に使うツールはcondaでインストールしてあるものとします。
※パソコンのスペックや目的によって最適なパラメータが違うのでオプションは各自いじっていただければと思います。
系統樹を描く
今回は系統樹を描いてみます。複数種のTElibをまとめて系統樹を書いたときに、A種とB種では共通のファミリーがゲノム中に存在するけどC種では全然違うファミリーがある、とか分かれば楽しいなと思っている。
あってるかわからないけど今回は{hoge}TElib.faを使って系統樹を描いてみます。
0. TElib.faを丸ごと使って下流の解析を行なってもいいと思うが、今回はseqkitを用いてGypsyだけ抽出してみる
seqkit faidx hoge.TElib.fa -r 'Gypsy' >hoge_Gypsy.TElib.fa
TElib.faはEDTAで使用したリファレンス配列が含まれるので、系統樹を描いたりする解析には使えなかった。てっきり新しく構築されたコンセンサスが出力されているのだと思っていたが違ったようだ。
1. TEsorterでdomain探索
TEsorter hoge_Gypsy.TElib.fa \
-db rexdb-plant \
-p 24
2. RH RT INTのdomainでアライメント+concatenate
python ~/miniconda3/envs/repeatmasker/pkgs/tesorter-1.3.0-py_0/site-packages/TEsorter/modules/concatenate_domains.py \
hoge_Gypsy.TElib.fa.rexdb-plant.cls.pep RH RT INT \
> hoge_Gypsy.TElib.fa.rexdb-plant.cls.pep.RH-RT-INT.aln
注意！
condaでTEsorterを入れた場合、TEsorterコマンドは正常に動くものの、concatenate_domains.pyは大元のスクリプトがあるところまでフルパスを記載する必要があるようです。
参考:
https://github.com/zhangrengang/TEsorter/issues/34
どのドメインでアライメントを行うかはお作法に従うのだと思われます。今回はLTR型なので以下の論文を参考にRH, RT, INTの3つのドメインを解析に使用しました。
参考:
https://mobilednajournal.biomedcentral.com/articles/10.1186/s13100-018-0144-1
3. iqtreeで系統樹推定
iqtree -s hoge_Gypsy.TElib.fa.rexdb-plant.cls.pep.RH-RT-INT.aln \
-bb 1000 -nt AUTO -T 20
出力されるtreefileで系統樹描く