2022/06/30 Seqkit のソート項目を追記
lastzについての日本語ドキュメントが少ないので備忘録として。
lastzはWhole genomeレベルでのアライメントができるツールである。
結構古くからあり有名なのでこの辺の説明は割愛。今も開発が続いている。
ダウンロードとインストール
GitHubが公開されているが、リリース版はこちらからダウンロードできる。
https://github.com/lastz/lastz/releases
今回はUbuntu 20.04での実行となる。現時点(2021/10/21)での最新版をダウンロードする。
$ wget https://github.com/lastz/lastz/archive/refs/tags/1.04.15.tar.gz $ tar -zxvf 1.04.15.tar.gz $ cd ./lastz-1.04.15/src $ make $ make install
これで完了。
ホームディレクトリ下にlastz-distrib/binができているので、そこに入っている。
また、今回は全ゲノムレベルの大きいファイルを必要とするため、追加でlastz_32もインストールしてやる。
$ make lastz_32 $ make install_32
あとはそのディレクトリにPATHを通しておく。
全ゲノムアライメントのための下準備
全ゲノムアライメントのための下準備。必須ではないが、やっておくと楽。今回はとある染色体レベルで繋がっているゲノムを使うのだが、染色体にマップされなかった領域のシーケンスも含まれている。
そこで染色体のみを切り出してやる必要がある。
そのためにはseqkitを使ってやる。seqkitはGitHubからバイナリをダウンロードできる他、biocondaでもインストール可能。インストール方法はリンク先を参照のこと。
18本染色体がある生物なので、頭からその分を切り出してやる(必ずしも染色体アセンブリが頭からあるわけではないと思うのでよく見ること)。
$ seqkit head -n 18 genome.fa > Only_chr_genome.fa
これで切り出せているかを確認する。
もし染色体が先に来ていないようであれば、最小の染色体の長さだけ控えて以下のようなコマンドで切り出す。例えば最小染色体の長さが100000000の場合は以下の通り。
$ seqkit seq -m 100000000 genome.fa > Only_chr_genome.fa
もしくは染色体レベルでつながっていない場合、同じくseqkitのSortでソートしてやってから先頭から適当なところを切り出す。
$ seqkit sort -lr2 [input.fa] > [out.fa]
helpには書かれていなかったが、標準では短い順に出力されるため、-r オプションを忘れずに。
また、sedコマンドを利用することで、わかりやすいヘッダに変えておくと今後楽である。
$ sed 's/^>.*chromosome: />Hge_chr/' chr_Hge13.fa > rename_Hge_chr.fa
同様のことを、もう一つのアライメントするゲノムについても行っておく。
次に、クエリとなるFASTAの分割を行う。Multiple vs Multiple FASTA はできない。
$ seqkit split -i chr_genome.fa
LASTZ_32によるアライメント
以下のコマンドで実行する。マルチコア、マルチスレッド非対応のため時間がかかる……。
$ lastz_32 genome1_chr1.fa genome2_chr.fa[multiple] --notransition --step=20 --nogapped --progress=1 --format=segments > genome1_vs_genome2.txt
これを染色体の数だけ繰り返す。ヤツメウナギとかだと大変そう……
オプションについては公式ドキュメントを参照のこと。上のやつはざっくりしたアライメントにしようすることができる。
また、例えばもっとかけ離れた染色体同士でギャップを許容してできる限り長くアライメントを取りたいという場合は以下のようなコマンドとなる。
$ lastz_32 genome1_chr1.fa genome2_chr.fa[multiple] --chain --gfextend --gapped --queryhspbest=500 --step=20
--queryhspbestはアライメントスコアの上位500を取ってくる……というものの様子。なぜか500も入っていないことが多い気がする。
Output
アウトプットの形式は色々選べる。今回はcircosでプロットすることを目的としているので、以下のようなコマンドで出力させた。
$ lastz_32 genome1_chr1.fa genome2_chr.fa[multiple] --chain --gfextend --gapped --queryhspbest=500 --step=20 --format=general:name1,start1,end1,name2,start2,end2 > out.txt
これでタブ区切りのcircosインプット互換の形式が得られる。
