注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Bioinformatics home

 
 
 

日志

 
 

Blat的"percent identity"的Perl代码 (已经测试通过)  

2014-04-20 03:42:41|  分类: 生物信息编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

网上有个代码,这兄弟问了怎么我的代码跑出来是100,

好了,拿出UCSC上原始代码,重新用perl写了下。测试通过!放心使用。
代码中只能使用于mRNA-DNA的比对

use strict;

my ($file,$line,$pid,$query,$subject);

open $file,"alignment_up_gar.psl";

<$file>; <$file>;  <$file>; <$file>; <$file>;

print join(' ',qw(query subject p_identify(%))),"\n";   

while($line=<$file>)

{

$line=trim($line);

($query,$subject,$pid) = pslCalcMilliBad($line);

print join(' ',($query,$subject,$pid)),"\n";   


 

sub trim

{

    my $value=shift;

chomp($value);

$value=~s/^\s*//;

$value=~s/\s*$//;

return $value;

}

 


sub pslCalcMilliBad 

#/* Calculate badness in parts per thousand. */

{

    my $line=shift;   

my $sizeMul =1; #not protein

my $isMrna=1;  # mRNA

my @cols=split(/\t/,$line);

my ($qAliSize, $tAliSize, $aliSize);

my $milliBad = 0;

my $sizeDif;

my $insertFactor;

my $total;


# cols[0]  matches 

# cols[1]  misMatches 

# cols[2]  repMaches 

# cols[4]  qNumInsert 

# cols[6]  tNumInsert 

# cols[11] qStart 

# cols[12] qEnd 

# cols[15] tStart 

# cols[16] tEnd  

# cols[13] subject

# cols[9] query

my ($matches, $misMatches,$repMaches,$qNumInsert,$tNumInsert,$qStart,$qEnd,$tStart,$tEnd,$subject,$query)=($cols[0],$cols[1],$cols[2],$cols[4],$cols[6],$cols[11],$cols[12],$cols[15],$cols[16],$cols[13],$cols[9]);

$qAliSize = $sizeMul * ($qEnd - $qStart);

$tAliSize = $tEnd - $tStart;

my $aliSize; 

if($qAliSize < $tAliSize)

{

  $aliSize = $qAliSize;

}

else

{

$aliSize = $tAliSize; 

}


if ($aliSize <= 0)

{

return ($query,$subject,100);

}

$sizeDif = $qAliSize - $tAliSize;

if ($sizeDif < 0)

{

if ($isMrna)

{

$sizeDif = 0;

}

else

{

$sizeDif = -($sizeDif);

}

}


$insertFactor = $qNumInsert;

if (!$isMrna)

{

$insertFactor += $tNumInsert;

}


$total = ($sizeMul * ($matches + $repMaches + $misMatches));

if ($total != 0)

{

$milliBad = (1000 * ($misMatches*$sizeMul + $insertFactor + 

round(3*log(1+$sizeDif)))) / $total;

}  

return ($query,$subject,(100.0 - ($milliBad * 0.1 ) ));

}

 

sub round {

    my $number = shift;

    return int( $number + .5 );

}

  评论这张
 
阅读(546)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017