アライメント。
要するに配列を比較し、並べる操作を指します。
遺伝子の発現量を測定したいときなど、遺伝子解析においてはアライメントを必要とする場面が多々登場します。
配列の類似度を表す基準は複数あり、その計算方法も複数あるため、いろんなアライメントツールが登場してきています。
代表的なアライメントツールとしては、BWA、Bowtie、TopHatなんかがあげられますが、
こんなのを使っていると、全然内部のアルゴリズムを知らないまま、アライメントができてしまいます。
非常にまずい。
と、常日頃から思っているので、EMBOSSのneedleを使ったアライメントを
内部のアルゴリズムに触れつつ、ご紹介します。
アラインメントの種類
アラインメントには、大域アラインメントと、局所アラインメントトの2種類あります。
調べたい配列や、研究目的によって、どちらを行うか決定します。
大域アラインメント（global alignment）
二つの配列の全長にわたってアラインメントを行なったときの最も良いスコアが計算されます。
Needleman-Wunsch algorithm等が挙げられます。
局所アラインメント（local alignment）
各配列の部分配列同士のアラインメントを行なったときの最も良いスコアが計算されます。
Smith-Waterman algorithm、Goad-Kanehisa algorithm等が挙げられます。
Needleman–Wunsch アルゴリズム
Needleman–Wunsch アルゴリズムは、2つの配列の相同性を算出するアラインメント用の（恐らく）最も古いアルゴリズムです。
グローバルアライメント | グローバルアライメントを求める Needleman–Wunsch アルゴリズム
ここを読みましょう。とってもわかりやすいです。
Interactive demo for Needleman–Wunsch algorithm
インタラクティブなデモンストレーションができます。すごい。
https://www.ibm.com/developerworks/jp/java/library/j-seqalign/index.html
Javaによる実装が載ってます。すごい。
要は、配列二つを並べて、
- マッチ
- ミスマッチ
- ギャップ
のどれに該当するかでそれぞれのスコアを決定し、
スコアが最も高くなるような配列を動的計画法によって求めるアルゴリズムです。
EMBOSSによる大域アラインメント（needle）
EMBOSSとは！
- "The European Molecular Biology Open Software Suite" の略。解析用ソフトウェア・パッケージ。
- フリーのオープンソースソフトウェアで、分子生物学のユーザ・コミュニティー (例：EMBnet) のために開発された。
- ウェブからの透過的な配列データ検索、様々なフォーマットを持つデータの自動処理なんかが可能。
- Linux、Digital Unix、Irix、Solaris を含んだ、全ての Unix プラットフォームをサポート。Win、Macでも動作可能。
EMBOSS に含まれる needle は、Needleman-Wunsch algorithmを実装したものです。
Needleman-Wunsch algorithmを使ったアラインメントがお手軽に使えちゃいます。
実践！
GUI版もあるが、コマンドで実行できたほうが何かと便利なので、以下ではCUI版を試します。
EMBOSSのインストール
Ubuntuの場合、
$ sudo apt install emboss
それ以外は以下を参考にインストールしてください。
EMBOSS Downloads
EMBOSS needleの実行
needle コマンドで実行します。
実行の際には、引数を与えてギャップペナルティー等を指定します。
$ needle -asequence JX469991.1.fa \ -bsequence JX205496.1.fa \ -gapopen 10 \ -gapextend 1 \ -endweight True \ -endopen=10 \ -endextend=1 \ -outfile JX469991.1_JX205496.1_needle.txt Needleman-Wunsch global alignment of two sequences
指定可能な他の引数についてはここを参考にどうぞ。
留意すべきこととして、伸長ペナルティーは開始ペナルティーより小さく設定することが多いです。
今回指定した項目は以下の8つ。
| 項目 | 説明 |
|---|---|
| -asequence | 入力ファイル名① |
| -asequence | 入力ファイル名② |
| -gapopen [10.0] | Gapの開始地点でのペナルティースコア。デフォルトは10 |
| -gapextend [0.5] | Gapの伸長地点でのペナルティースコア。デフォルトは0.5 |
| -endweight [false] | 5'末端、3'末端において、Gapペナルティーの重みづけを行うか否か。 |
| -endopen [10.0] | 配列末端におけるGapの開始地点でのペナルティースコア。デフォルトは10 |
| -endextend [0.5] | 配列末端におけるGapの伸長地点でのペナルティースコア。デフォルトは0.5 |
| -outfile | 出力ファイル名。デフォルトは {入力ファイル名①}.needle |
引数を指定しないと、最低限必要な入力ファイルやパラメータについてインタラクティブに聞かれます。
$ needle Input sequence: JX469991.1.fa Second sequence(s): JX205496.1.fa Gap opening penalty [10.0]: 10 Gap extension penalty [0.5]: 0.5 Output alignment [jx469991.needle]: jx469991.needle Needleman-Wunsch global alignment of two sequences
結果ファイル
出力ファイルは以下のようになります。
視覚的にわかりやすい図が作成され、アライメントスコアが算出されます。
今回の場合だと、アライメントスコアは935.0でした。
########################################
# Program: needle
# Rundate: Fri 14 Sep 2018 11:00:48
# Commandline: needle
# -asequence JX469991.1.fa
# -bsequence JX205496.1.fa
# -gapopen 10
# -gapextend 1
# -endweight True
# -endopen 10
# -endextend 1
# -outfile JX469991.1_JX205496.1_needle.txt
# Align_format: srspair
# Report_file: JX469991.1_JX205496.1_needle.txt
########################################
#=======================================
#
# Aligned_sequences: 2
# 1: JX469991.1
# 2: JX205496.1
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 1.0
#
# Length: 2191
# Identity: 780/2191 (35.6%)
# Similarity: 780/2191 (35.6%)
# Gaps: 1141/2191 (52.1%)
# Score: 257.0
#
#
#=======================================
JX469991.1 1 ATGGAAGCCTCCGCCGGCTCGTCGCCACCGCACTCCCAAGAGAACCCGCC 50
|| ||
JX205496.1 1 AT-------TC--------------------------------------- 4
JX469991.1 51 GGAGCACGGTGGCGACATGGGAGGGGCCCCCGCGGAGGAGATCGGAGGGG 100
JX205496.1 4 -------------------------------------------------- 4
JX469991.1 101 AGGCGGCGGATGATTTCATGTTCGCTGAAGACACGTTCCCCTCCCTCCCG 150
||||.|||| ||||||.|.
JX205496.1 5 -------------TTTCCTGTT-------------------TCCCTCTCT 22
JX469991.1 151 GACTTCCCTTGCCTTTCGTCGCCGTCCAGCTCCACCTTCTCGTCCAACTC 200
|.|.||||| ||.|| ||| |||||.||| ||
JX205496.1 23 GCCCTCCCT---CTCTC-------TCC----CCACCATCT--------TC 50
・
・
JX469991.1 1893 GGGAGGCA---AGCACAGGCCCCTCTGTCCAGCAGGTCCAGGGAGAGCCG 1939
|||..|| ||.| ||.|||
JX205496.1 1094 -GGACACATTGAGAA-------CTTTGT---------------------- 1113
JX469991.1 1940 CAGCCGCCGGTGCTCCTGAAGACGCCGTCGTCGACGGGGTCAGCGGCGCC 1989
||....|.|| |||
JX205496.1 1114 ----------TGTATATTAA-------------------TCA-------- 1126
JX469991.1 1990 TGCAAGGGGAGGTCTCCGGAAGGCGTGCGGCGGGTTCGGCAGCAGGGAGC 2039
|| ||| |..|....|||.|
JX205496.1 1127 ---AA--------CTC---------TATGTTTAGTTTG------------ 1144
JX469991.1 2040 CGGCGCCATG----------AGCCAGATGGCGGTGAGCATC 2070
|||| |...|.| |..|..
JX205496.1 1145 ------CATGTAAAAAAAAAAAAAAAA--------AAAAAA 1171
#---------------------------------------
#---------------------------------------
是非お試しくだされ。
