macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

ゲノムの相同性の高い領域の網羅的な検索 LAST

2019 6/9 biocondaインストール追記

2021 2/19 曖昧な日本語を修正

2021 7/28 maf-convert追記

 

 ゲノム配列のアラインメントは多くの研究の基礎を成している。ゲノムアラインメントは、さまざまな日常的ではあるが重要な選択、たとえばリピートをマスクする方法や使用するスコアパラメータに依存する。驚くべきことに、実際のゲノムデータを用いたこれらの選択の大規模な評価はない。さらに、疑わしいアライメントの速度を制御するための厳密な手順は採用されていない。

 動物、植物、真菌ゲノムのアラインメントについて、495のスコアパラメータの組み合わせを評価した。我々の正確さのゴールドスタンダードとして、我々はタンパク質と構造RNAのマルチプルアラインメントによって暗示されるゲノムアラインメントを使用した。UCSCゲノムデータベースにおけるアラインメントの基礎をなすHOXDスコアリングスキームが最適からは程遠いことを見出し、そしてより良いパラメーターを示唆する。 X-dropパラメータの値が大きいほど、必ずしも良いとは限らない。 E値は疑わしいアライメントの割合を正確に示すが、これはタンデムリピートが非標準的な方法でマスクされている場合に限られる。最後に、γ-セントロイド(確率的)アラインメントがアラインメントされた塩基の信頼性の高いサブセットを見つけることができることを示す。

 

 

マニュアル

http://last.cbrc.jp/doc/last-tutorial.html

 

インストール

condaやbrewワンライナー導入できる。

 #bioconda (link)
cponda install -c bioconda -y last

#homebrew
brew install LAST

 

ラン 

相同性の高い領域を検索するには、はじめに比較するリファンレスゲノム(ref.fa)のインデックスを作る必要がある。

lastdb -cR01 db ref.fa

いくつか.dbファイルができる。

 -cR01でcacacacacacacacacacacacaなど単純リピートのアライメントが除かれる。

 

比較したいゲノム配列(compare.fa)を指定して以下のコマンドを打つ。

lastal db qurry.fa -P0 > myalns.maf

 比較結果がmyalns.mafに保存される。

 -P: CPU数。maxまで使うなら-P0

 

結果を表示するには-Sをつけてlessを打つ(-S 長い文字列を改行せず表示)。

less -S myalns.maf

f:id:kazumaxneo:20170623142453j:plain

このような画面が表示される。先頭部分のみ載せた。

 

比較するゲノムがリアレンジメントを起こしており、 途中で配列を切ってアライメントさせる必要がある場合、last-splitに出力をパイプすることが推奨されている(RNAとゲノムの比較など)。

lastal db qurry.fa -P0 | last-split > myalns.maf

 

mafフォーマットからsamへの変換はmaf-convertを使う(last)。このツールでmafからaxt, blast, blasttab, html, psl, sam, tabへ変換できるmamba install -c bioconda last -y )。

maf-convert sam myalns.maf -r ''ID:X SN:Lambda_NEB LN:48502' > last-alignments.sam

 

 

 

LASTは次世代シークエンシングリードのアライメントにも使うことができる。 まずindexを作る。

lastdb -uNEAR -R01 db re.fa

 -uNEARで短いが強い相同性を示すものをアライメントするスキームになる。

 

fastqのアライメントを行う。

lastal -Q1 -D100 db queries.fq | last-split > myalns.maf

ランタイムはあまり短くない。おおよそ10万リードで1分程度の時間がかかる。 

 

 LASTを使ってエラーの多いnanoporeリードをアライメントする設定は、以下の記事を参考にしてください。


 

 

proteinの相同性探索を行うこともできる。流れは塩基配列の場合と同じである。

比較するリファンレスアミノ酸配列(ref.aa)のインデックスを作る。

lastdb -cR01 db ref.aa

 

比較したいゲノム配列(compare.fa)を指定して以下のコマンドを打つ。

lastal db query.aa > myalns.maf

LASTは、長い配列ならば相同性がある程度低くても見つけるスコア設定になっている。40.a.aより短い相同性が非常に高いアミノ酸配列を見つけるなら-pPAM30をつけることが推奨されている (example 4)。

lastal -pPAM30 db query.aa > myalns.maf

 

less -Sで結果を見る。

f:id:kazumaxneo:20170623143459j:plain

 

LASTのレシピは公式サイトに色々紹介されている。

LAST Tutorial

 

 

引用

Parameters for accurate genome alignment
Martin C Frith, Michiaki Hamada, Paul Horton

BMC Bioinformatics. 2010; 11: 80