遺伝子やタンパクの機能・構造を予測するために配列類似性を検討することは有効な手段である。
配列どうしの類似性を評価するためにギャップを適宜配して対応文字を揃えて並べる操作を配列アライメントという。
並列対象が2つの配列である場合はペアワイズアライメント、3つ以上の配列である場合は多重配列アライメント(multiple sequence alignment, MSA)という。
並列範囲が配列全体に渡る場合は大域的アライメントといい、類似部分のみに限局する場合は局所的アラインメントという。
次節では生物配列解析の基本となるペアワイズアライメントについて述べる。
ペアワイズアライメント
2つの配列の類似度は文字の類似度スコアを定義し、並列された配列の各文字の類似度の和をとることで評価することができる*1。
最適な配列アライメントは配列類似度の最大値を与える場合として得られる。
その統計的評価は、ランダム配列対における類似度スコア分布に基づき、スコアが配列類似度の最大値以上となる\( p \)値を計算し、これが有意に小さいかどうかを判断する。
配列類似度の計算において、最も単純なモデルでは文字類似度を一致、不一致、ギャップ状態の観点でのみスコア化するが、一般的には配列アライメントの背後に進化モデルを考慮する。
その場合、アライメントの不一致は置換変異、ギャップは挿入・欠失変異に相当する。
アミノ酸配列の場合では配列アライメントの一致・不一致に対して、共通祖先に由来する相同配列間で観察される置換頻度に基づき様々なアミノ酸置換スコアが定義されている。
ギャップはアライメントの評価ではペナルティーとなるよう設定されるが、挿入・欠失が断続的か連続的かを考慮したペナルティー値の設定がなされる*2。
最適配列アライメントは動的計画法により得られる。
動的計画法は問題を分割して部分的な解を求め、それを利用することで漸次範囲拡大した問題を解決する手続きの総称である。
配列アライメントは、２つの配列を縦横に配置した格子構造において、左上端から右下端に向けて縦・横（ギャップ）、斜め（一致または不一致）のいずれかの辺を辿る経路として表現することができる。
動的計画法に従い、左上端の原点から始めて右下端に至るまで、各頂点にそこまでの類似度スコアの最大値を格納し順次計算していくことで、配列類似度を最大化するための経路履歴が得られる。
その経路復元(traceback)をしたものが最適配列アライメントとなる。
上述の動的計画法に基づき、最適大域的アライメントを求める手続きをNeedleman-Wunsch法という。
一方、この手続きの一部を変更することで最適局所的アライメントを求めることができるようにしたものにSmith-Waterman法がある。
*1
２つの配列\( x, y \)のアライメントを\( x', y'\)としアライメント長を\( L \)とすると、配列類似度\( S(x', y')\)は\( x', y' \)の\( i \)番目の文字\( x_{i}', y_{i}' \)の類似度\( s(x_{i}', y_{i}') \)の和
\( S(x', y') = \sum_{i=1}^{L} s(x_{i}', y_{i}') \)
として得られる。
*2
一文字あたりのギャップペナルティを(\( -d,\ d \gt 0\))とした場合、\( k \)個の連続したギャップペナルティには主に以下のいずれかの値が用いられる。
- \( -dk \) (linear gap penaltyという)
- \( -d \)
- \( -d-ke \) where \( e \gt 0, \ d \gt e \) (affine gap penaltyという)
2つの配列類似性の直感的な表現方法にdotplotがある。
2つの配列を縦横に配置した格子構造を考え、同じ文字（塩基/アミノ酸）に対して打点する。
全てを打点した後、２つの配列の類似部分には斜線が現れる。
塩基の場合は４種類しかなくノイズが多くみられるため、複数文字をまとめたwindow単位で比較し、全て一致する場合のみ打点することがある。
dotplotの実装としてはEMBOSSのdottupがある。
相同性検索
進化的類縁性を検討するために、配列データベースに対して局所的に類似した配列の検索がなされる。
これは相同性検索と称するが、その名称が示す目的以外にも遺伝子やタンパクの構造・機能の予測など幅広い用途に利用されている。
Smith-Waterman法をデータベース検索に用いると計算コストが膨大なために実用的でない。
データベース検索を高速化するために、元の質問配列のk-merに対して完全に一致するデータベース配列部位をスクリーニング後*1, *2、その周辺でギャップを許容したアライメントを求めることが行われる。
代表的な相同性検索プログラムにはFASTAや
BLASTがあり、FASTAよりもBLASTの方が現在頻用されている。
*1 k-merが短いと偽陽性数が増えて、長いと検出感度が低下する。
*2
FASTAではハッシュ表、BLASTでは有限オートマトン(Aho Corasickアルゴリズム)を利用して質問配列とデータベース配列のk-merどうしを高速照合している。
ハッシュ表は「キー」と「値」の組を表として格納したデータ構造であり、「キー」に対応する「値」を定数時間で参照できる。
k-merを4進数（塩基配列の場合）または20進数（タンパク配列の場合）の数値とみてハッシュ表の「キー」とし、その「値」としてk-merの出現位置をリスト表現する。
BLASTを利用するにはデータベース配列の索引を作成する。
$ makeblastdb -in sbj.fa -dbtype nucl -hash_index -parse_seqids # 塩基配列の場合
$ makeblastdb -in sbj.fa -dbtype prot -hash_index -parse_seqids # アミノ酸配列の場合
質問配列及びデータベース配列が塩基配列であるかタンパク配列であるかの相違により、使用するプログラムが異なる（下記参照）。
blastnプログラムでは塩基配列間の配列アライメントをするが、他のプログラムではタンパク配列に翻訳後に配列アライメントする。
翻訳は相補鎖も含めた6通りの読み枠で行われる
核酸配列の比較は主に近縁種間の比較で用いられ、遠縁種間の比較ではタンパク配列を用いることが多い。
PSI-BLASTではPSSMを用いた検索を反復実行することで非常に遠い類縁関係にある相同配列でも検出することができる。
質問配列中における低複雑性領域はマスク(存在しないアミノ酸Xまたは塩基Nで置換)することで検索対象から除外する。
質問配列中のアミノ酸組成が平均的な組成から大きく逸脱している場合は偽陽性や偽陰性が生じやすくなるため、BLASTには質問配列のアミノ酸組成に合わせてアミノ酸置換スコア行列の補正を行うオプションがある。
コマンドの基本書式は以下の通りである。
$ command -query qry.fa -db sbj.fa
NGSリード配列のアライメント
NGSショートリードをゲノム配列に対して検索するためにはsuffix arrayに基づく手法が頻用されている。
suffix arrayは文字列の接頭辞の開始位置を要素とし、それらを辞書順に整序した配列である。
任意長の部分文字列の出現位置情報を保持しており、二分探索を用いて効率的に探索できる。
suffix arrayの発展技術にBurrows-Wheeler Transformation (BWT)後の文字列を利用したFM indexがある。
BWTでは接尾辞の代わりに1文字ずつ巡回シフトした文字列を用い、それらを辞書順に整序することで得られた文字列の最後の文字だけをつなげる。
BWT後の文字列中では同じ文字列が連続する傾向があるため圧縮性に優れるとともに、BWT後文字列のみから元の文字列を復元することも可能である。
FM indexではこの性質を利用して質問文字列と一致する接尾辞配列の範囲を高速に求めることができる。
多重配列アライメント
多重配列アライメント(multiple sequence alignment, MSA)に動的計画法を適用することは多大な計算コストを要するために実用的でない。
累進法（Progressive alignment）では、まず全配列間でペアワイズアライメントを実施し、配列間の距離行列を求めて、案内木を構築する。
次に案内木で示された類似度が高いものからペアワイズアライメントを順次行うことでMSAを構築する*1。
実装プログラムにはClustal Omegaがある。
反復改善法（Iterative approaches）では構築済みのMSAを任意の２つのグループに分割し、再度アライメントを実行する操作を反復することでアライメントエラーを修正する仕組みが採られている。
実装プログラムにはMAFFT (Consistency-based method)がある。
MSAのスコアはアライメントの各位置における類似度スコアの和で表される*2。
N本のMSAの場合、N個の文字がその位置に揃う確率を算出するのが本来適切ではあるが、計算量が膨大となり現実的でない。
そのため、全てのペアの類似度スコアを求め、それらの和をMSAの類似度スコアとしたsum of pairs (SP)と称する方式が一般的に採用されている*3。
MSAの可視化・編集プログラムにはJalviewがある。
Clustal Omegaのclustaloコマンドを用いてMSAを構築する。
$ clustalo -i in.fa > out.fa
*1
累進法では各段階で生じたアライメントエラーが固定され維持され続ける難点がある。
類似度の低い配列どうしのアライメントの場合にアライメントエラーが生じやすく、アライメント初期段階で生じたエラーは以後の全てのアライメントに悪影響を及ぼすことを考慮すると、初めに案内木を構築することで、適切なアライメント順序を指定する必要がある。
*2
\( N \)本の配列アライメント\( M \)のスコア\( S(M) \)は、アライメント長をLとすると、各位置\( i \)の類似度スコア\( S(Mi) \)の和で算出可能である。
\[
S(M) = \sum^{L}_{i=1} S(M_{i})
\]
*3
k', l'を配列k, lがMSAで並置された状態の配列であるとすると
\[
S(M_{i}) = \sum_{k \lt l} s(k'_{i}, l'_{i}), \ s(-,-) = 0
\]
プロファイルを用いた検索
MSAにより得られた配列特徴はコンセンサス配列、正規表現、出現確率行列、位置特異的スコア行列(position specific scoring matrix, PSSM)、隠れマルコフモデル(HMM)のいずれかで表現される。
後三者を利用する場合、確率論的な検索が可能である。
出現確率行列やPSSMは重み行列と総称されるもので、出現確率行列では各位置の文字の出現確率を行列要素とするが、PSSMではこれをバックグラウンド(配列全体での各文字の出現頻度)を考慮した相対値として定義される*1。
PSSMは挿入・欠失の処理に不適であるが、HMMではこれらを考慮した配列尤度の計算が可能である。
*1
PSSMの要素\( PSSM(i, j) \)は位置\( i \)における文字jの出現頻度\( f(i, j) \)を、文字\( j \)のbackgroundにおける出現頻度\( g(j) \)で除した値の対数値(対数尤度比)
\( PSSM(i, j) = log \frac{f(i, j)}{g(j)} \)
で定義される。
\( f(i, j)=0 \)のときの発散を防ぐために、
\( PSSM(i, j) = log \frac{f(i, j) + 1}{g(j) + n} \)という計算式
（ただし塩基配列では\( n=4 \)、タンパク配列では\( n=20 \)）で定義する場合もある。
タンパク配列解析の場合では位置\( i \)における文字\( j \)の出現しやすさを文字\( k \)から\( j \)への置換スコア\( s \)を考慮して
\( PSSM(i, j) = \sum^{n}_{i=1} \frac{f(i, k)}{g(j)} s(k,j) \)
のように定義する場合もある。
タンパクドメイン解析
InterProはモチーフやドメインなどのタンパクの配列特徴に関連するデータベースを統合し、それらを横断的に検索することができるようにした統合データベースである。
タンパク配列特徴の重み行列情報はPROSITEのMATRIXエントリー、HMM情報はPfamからそれぞれ入手可能ある。
機能的に保存されたタンパクドメインは単なるアミノ酸配列の類似性からは検索することができない。
HMMによるタンパクドメイン検索プログラムにHMMERがある。
対象配列群sbj.fa.gzに対して指定ドメインdomain_A.hmmを検索するにはhmmsearchを用いる。
$ hmmsearch domain_A.hmm sbj.fa.gz
指定タンパクのドメイン構成を明らかにするためにはhmmscanを用いる。
前段階としてタンパクドメインのHidden Markov Model (HMM)プロファイルをPfamから入手し、hmmpressを用いて索引を作成後、hmmscanを用いてドメイン構成を出力する。
$ curl -O ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz
$ gunzip Pfam-A.hmm.gz
$ hmmpress Pfam-A.hmm
$ hmmscan -o out Pfam-A.hmm in.fa
構成情報の可視化プログラムは
こちらの一覧から確認できる。