用bioperl计算比对的相似度

十二月 7, 2014

BioPerl 的比对模块(Bio::Align)中有三个方法用来计算比对相似度,它们分别是 percentage_identity,overall_percentage_identity 与 average_percentage_identity。

下面先以两条序列的比对为例,假设在文本文件 align.fasta 中放如下比对(FASTA 格式):


>seqA
CACAC--CCG
>seqB
CTCACATCCG

这个比对长10bp,第 2 位碱基不同,第 6、7 位还有 gap。下面来看这三个方法计算结果的差异:


use Bio::AlignIO;
my $align_obj = Bio::AlignIO->new(-file=>'align.fasta', -format=>'fasta')->next_aln;
my $i = $align_obj->percentage_identity;
my $j = $align_obj->overall_percentage_identity;
my $k = $align_obj->average_percentage_identity;

计算结果是:


$i: 87.5
$j: 70
$k: 87.5

因此,可以推测:overall_percentage_identity 是计算包括 gap 在内的相似度,gap 位点算作有差异的位点,因此结果是 7/10=70%(注意返回值是百分比的值),而 percentage_identity 与 average_percentage_identity 是计算不包括 gap 在内的相似度,不把 gap 位点计算在内,因此结果是7/8=87.5% 。

如果比对中的序列数目多于 2 条,比如这样时:


>seqA
CACAC--CCG
>seqB
CTCACATCCG
>seqC
CTCACATCCG

这时的计算结果是:


$i: 92.307692
$j: 70
$k: 92.307692

可见,overall_percentage_identity 是计算“总体”的相似度,与两条序列的情况一致。而  percentage_identity 与 average_percentage_identity是在计算“两两之间”的相似度,因此结果是:(7+7+10)/(8+8+10)=92.31

补充一句:如果把 seqB 的序列全部换成小写字母,结果还是一样。因此该方法对大小写是不敏感的。

posted in Perl by billzt

Follow comments via the RSS Feed | Leave a comment | Trackback URL

说点什么

  Subscribe  
提醒
 

Copyright © 2010-2017 | Powered by Wordpress and MySQL. Theme by Shlomi Noach, openark.org