2018.7.28 星期六 晴 biolearn
通過blast、hhblits或者其他方式構(gòu)建的多序列比對（MSA）都存在大量的冗余，這些冗余的序列導致多序列比對的多樣性較少并相較于其他序列來說提供較少的信息，因為冗余序列的存在如果所有的序列擁有相同的權重是不合理的，通常采用序列加權方法用于減少多序列比對中的冗余序列并維持序列的多樣性，這里介紹常用的Henikoff權重體系對序列進行加權。
Henikoff權重體系原理
Henikoff認為可以通過序列上每一個位點在MSA中殘基的類型和數(shù)目給每個位點分配權重，通過給每個位置不同的殘基類型分配相同的權重，然后在擁有相同殘基的序列中平均的分配這個權重用于代表這個位置的多樣性，分別給MSA中的每個位點進行權重最終實現(xiàn)給每條序列加權。
如下圖所示，對于位點1來說，殘基僅有G，出現(xiàn)的次數(shù)卻有4次，所以權重是1/（1*4），而對于位點2來說，殘基有Y, F兩種，數(shù)目分別是3次和1次，首先對殘基給與相同的權重1/2，然后擁有相同殘基的序列平分這個權重，所以對于殘基Y的權重是1/（2*3），而殘基F的權重是1/（2*1），以此計算每個位點的權重，每條序列的權重即為這條序列所有位點的權重加和除以所有序列所有位點權重的加和。
對于下圖的四條序列來說，可以看到對序列進行權重之后，序列3的權重最低，原因是對于每個位點來說，序列3中每個位點的殘基在MSA中每列中出現(xiàn)次數(shù)都是較多的導致其擁有較低的權重從而導致整個序列的權重較低，說明通過這種權重方式可以很好的消除MSA中序列冗余的情況。
Henikoff權重體系計算方法
計算公式
其中 N 表示MSA中序列的總數(shù)；L 表示序列的長度；i 表示序列中第 i 個序列位置；j 表示MSA中第 j 條序列；r_i 表示在 i 位置上不同殘基類型的數(shù)目，s_i 表示對應殘基出現(xiàn)的數(shù)目，s_ij 表示第 j 條序列第 i 位置上殘基在第 i 列出現(xiàn)的數(shù)目。
對應上圖的例子，此公式中的分母的值對應的即為每條序列每個位置權重的加和，可以看到最終的值即為序列的長度，因為在MSA中每列的權重加和都為1，例如上圖的Total為5，所以計算公式可以進一步簡化為
計算代碼（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";
}
運行結(jié)果
0 0.070844070961718 #第 1 條序列的權重
1 0.070844070961718 #第 2 條序列的權重
2 0.070844070961718 #第 3 條序列的權重
3 0.0762689075630252 #第 4 條序列的權重
4 0.0856676003734827 #第 5 條序列的權重
5 0.126055088702147 #第 6 條序列的權重
6 0.100401493930906 #第 7 條序列的權重
7 0.183715219421102 #第 8 條序列的權重
8 0.215359477124183 #第 9 條序列的權重
可以看到前5條非常相似的序列的權重都很低，可以有效的排除MSA中序列的冗余并提高多樣性~
參考文獻
Henikoff,S. and Henikoff,J.G. (1994) Position-based sequence weights. J.Mol. Biol., 243, 574–578.
