2018.7.28 星期六 晴 biolearn
通過(guò)blast、hhblits或者其他方式構(gòu)建的多序列比對(duì)（MSA）都存在大量的冗余，這些冗余的序列導(dǎo)致多序列比對(duì)的多樣性較少并相較于其他序列來(lái)說(shuō)提供較少的信息，因?yàn)槿哂嘈蛄械拇嬖谌绻械男蛄袚碛邢嗤臋?quán)重是不合理的，通常采用序列加權(quán)方法用于減少多序列比對(duì)中的冗余序列并維持序列的多樣性，這里介紹常用的Henikoff權(quán)重體系對(duì)序列進(jìn)行加權(quán)。
Henikoff權(quán)重體系原理
Henikoff認(rèn)為可以通過(guò)序列上每一個(gè)位點(diǎn)在MSA中殘基的類(lèi)型和數(shù)目給每個(gè)位點(diǎn)分配權(quán)重，通過(guò)給每個(gè)位置不同的殘基類(lèi)型分配相同的權(quán)重，然后在擁有相同殘基的序列中平均的分配這個(gè)權(quán)重用于代表這個(gè)位置的多樣性，分別給MSA中的每個(gè)位點(diǎn)進(jìn)行權(quán)重最終實(shí)現(xiàn)給每條序列加權(quán)。
如下圖所示，對(duì)于位點(diǎn)1來(lái)說(shuō)，殘基僅有G，出現(xiàn)的次數(shù)卻有4次，所以權(quán)重是1/（1*4），而對(duì)于位點(diǎn)2來(lái)說(shuō)，殘基有Y, F兩種，數(shù)目分別是3次和1次，首先對(duì)殘基給與相同的權(quán)重1/2，然后擁有相同殘基的序列平分這個(gè)權(quán)重，所以對(duì)于殘基Y的權(quán)重是1/（2*3），而殘基F的權(quán)重是1/（2*1），以此計(jì)算每個(gè)位點(diǎn)的權(quán)重，每條序列的權(quán)重即為這條序列所有位點(diǎn)的權(quán)重加和除以所有序列所有位點(diǎn)權(quán)重的加和。
對(duì)于下圖的四條序列來(lái)說(shuō)，可以看到對(duì)序列進(jìn)行權(quán)重之后，序列3的權(quán)重最低，原因是對(duì)于每個(gè)位點(diǎn)來(lái)說(shuō)，序列3中每個(gè)位點(diǎn)的殘基在MSA中每列中出現(xiàn)次數(shù)都是較多的導(dǎo)致其擁有較低的權(quán)重從而導(dǎo)致整個(gè)序列的權(quán)重較低，說(shuō)明通過(guò)這種權(quán)重方式可以很好的消除MSA中序列冗余的情況。
Henikoff權(quán)重體系計(jì)算方法
計(jì)算公式
其中 N 表示MSA中序列的總數(shù)；L 表示序列的長(zhǎng)度；i 表示序列中第 i 個(gè)序列位置；j 表示MSA中第 j 條序列；r_i 表示在 i 位置上不同殘基類(lèi)型的數(shù)目，s_i 表示對(duì)應(yīng)殘基出現(xiàn)的數(shù)目，s_ij 表示第 j 條序列第 i 位置上殘基在第 i 列出現(xiàn)的數(shù)目。
對(duì)應(yīng)上圖的例子，此公式中的分母的值對(duì)應(yīng)的即為每條序列每個(gè)位置權(quán)重的加和，可以看到最終的值即為序列的長(zhǎng)度，因?yàn)樵贛SA中每列的權(quán)重加和都為1，例如上圖的Total為5，所以計(jì)算公式可以進(jìn)一步簡(jiǎn)化為
計(jì)算代碼（perl腳本）
例子MSA
RPYACPVESCDRRFSRSDELTRHIRIHTGQKPFQCRICMRNFSRSDHLTTHIRTHTGEKPFACDICGRKFARSDERKRHTKIHLR
RPYACPVESCDRRFSRSDELTRHIRIHTGQKPFQCRICMRNFSRSDHLTTHIRTHTGEKPFACDICGRKFARSDERKRHTKIHLR
RPYACPVESCDRRFSRSDELTRHIRIHTGQKPFQCRICMRNFSRSDHLTTHIRTHTGEKPFACDICGRKFARSDERKRHTKIHLR
RPYACPVETCVRRFSRSDELTRHIRIHTGQKPFQCRICMRNFSRSDHLTTHIRTHTGEKPFPCDICGRKFARSDERKRHTKIHLR
RPYACPVENCDRRFSRSDELTRHIRIHTGQKPFQCRICMRSFSRSDHLTTHIRTHTGEKPFQCDTCGRKFARSDEKKRHAKVHLK
RPRKYPNRASKTPFSRSDELSRHLRIHTGHKPFQCRICMRSFSRSDHLTTHIRTHTGEKPFSCEQCGRKFARSDERRRHMRIHLR
RPYPCPAEGCDRRFSRSDELTRHIRIHTGHKPFQCRICMRNFSRSDHLTTHIRTHTGEKPFACDFCGRK----------------
KPYHCDI--CGKSFSRSDHLTTHKHIHTRERPYHCDICGKSFSVDSHLTTHKRIHAGEKPYQCDICGKSFSQTNILTTHKRIH--
KPYHCDI--CGKSFIESGQLTRHRRSHTGEKPYPCDICGRPFSYRSTFNKHRRIHTGEMP-------------------------
腳本
#!/usr/bin/perl -w
use strict;
use warnings;
open(IN,"<","F:/msa_test.txt") or die "$!";
chomp(my @msa = <IN>);
close IN;
my %aanum;
foreach my $j(@msa) {
my @tem = split(//,$j);
foreach my $i(0..$#tem) {
$aanum{$i}{$tem[$i]} ++;
}
}
foreach my $t(0..$#msa) {
my @tem = split(//,$msa[$t]);
my $numerator = 0;
foreach my $i(0..$#tem) {
my $aatype = scalar(keys%{$aanum{$i}});
my $numres = $aanum{$i}{$tem[$i]};
$numerator += 1/($aatype * $numres);
}
my $weight = $numerator/scalar@tem;
print "$t\t$weight\n";
}
運(yùn)行結(jié)果
0 0.070844070961718 #第 1 條序列的權(quán)重
1 0.070844070961718 #第 2 條序列的權(quán)重
2 0.070844070961718 #第 3 條序列的權(quán)重
3 0.0762689075630252 #第 4 條序列的權(quán)重
4 0.0856676003734827 #第 5 條序列的權(quán)重
5 0.126055088702147 #第 6 條序列的權(quán)重
6 0.100401493930906 #第 7 條序列的權(quán)重
7 0.183715219421102 #第 8 條序列的權(quán)重
8 0.215359477124183 #第 9 條序列的權(quán)重
可以看到前5條非常相似的序列的權(quán)重都很低，可以有效的排除MSA中序列的冗余并提高多樣性~
參考文獻(xiàn)
Henikoff,S. and Henikoff,J.G. (1994) Position-based sequence weights. J.Mol. Biol., 243, 574–578.
