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

云之南

风声,雨声,读书声,声声入耳;家事,国事,天下事,事事关心

 
 
 

日志

 
 
关于我

专业背景:计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

网易考拉推荐

KEGG API学习笔记(3)(转)  

2010-11-05 16:18:42|  分类: 生物信息学 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
KEGG API学习笔记(3)

http://www.qiuworld.com/blog/archives/83

今天接着昨天的任务继续。昨天已经实现了文件读取了表达差异基因的筛选工作,今天的任务就是把这些基因在代谢途径中显示出来。原代码如下:


#!/usr/bin/perl
use SOAP::Lite;
$wsdl = ‘http://soap.genome.jp/KEGG.wsdl’;
$soap_ser = SOAP::Lite->service($wsdl);

$filename = “C:\\TPS1NS.xls”;
open(MICROARRAY, $filename)||die (“Could not open file”);
$up_reg = 0;
$down_reg = 0;
@map_hash = ();
@path_hash = ();
$j = 0;
$times = 1;
print “gene_name\texp_ratio\tpath\n”;
while($line = <MICROARRAY>){
$pattern = “\x09″;
@words = split(/$pattern/, $line);
if ($words[2] eq “P” && $words[4] eq “P”){
$ratio=log($words[1])-log($words[3]);
if(abs($ratio)>=$times){
if($ratio>=$times){
$up_reg += 1;
}
if($ratio<=-$times){
$down_reg +=1;
}
if($words[5]=~/DB_XREF=(\D+)(\d+)\s/){
$db_pre=$1;
$db_id=$2;
if($db_pre=~/GI/){
@kegg_id = split(/\s+/,$soap_ser->bconv(“ncbi-gi:$db_id”));
for($i=0;$i<3;$i++){
$map_hash[$j][$i] = $kegg_id[$i];
}
$map_hash[$j][3] = $ratio;
#通过KEGG API获取基因所在的代谢途径,并存贮起来。
$map_hash[$j][4]=$soap_ser->get_pathways_by_genes([$map_hash[$j][1]]);
#如果返回的代谢途径不为空,则将其打印出来。
if(@{$map_hash[$j][4]}){
foreach $hit(@{$map_hash[$j][4]}){
my %path = (“gene_name”=>$map_hash[$j][1],
“exp_ratio”=>$map_hash[$j][3],
“path”=>$hit);
push @path_hash, \%path;
print “$path{gene_name}\t$path{exp_ratio}\t$path{path}\n”;
}
}
$j++;
}
}
}
}
}
close(MICROARRAY);



#将相关的代谢途径打印一遍
print “\nrelated pathways\n”;
for($i=0;$i<scalar(@path_hash);$i++){
$tmp = $path_hash[$i];
%path_t = %$tmp;
print $path_t{path}.”\n”;
}

#以代谢途径名排序后再打印一遍
print “\nsorted related pathways\n”;
@path_sort = sort {$a->{path} cmp $b->{path}} @path_hash;
for($i=0;$i<scalar(@path_sort);$i++){
$tmp = $path_sort[$i];
%path_t = %$tmp;
print “$path_t{gene_name}\t$path_t{exp_ratio}\t$path_t{path}\n”;
}

#通过KEGG API获取代谢图,并在图上以颜色标注差异表达的基因。
@url_list = ();
$path_tmp = “none”;
@obj_list = ();
@fg_list = ();
@bg_list = ();
for($i=0;$i<scalar(@path_sort);$i++){
$tmp = $path_sort[$i];
%path_t = %$tmp;
if($path_tmp eq “none”){
$path_tmp = $path_t{path};
push @obj_list, $path_t{gene_name};
push @fg_list, “#ff0000″;
push @bg_list, “#ffff00″;
next;
}
elsif($path_t{path} eq $path_tmp){
push @obj_list, $path_t{gene_name};
push @fg_list, “#ff0000″;
push @bg_list, “#ffff00″;
next;
}
else{
$object_url=$soap_ser->color_pathway_by_objects($path_tmp,
\@obj_list,\@fg_list,\@bg_list);
#    print “$path_tmp,@obj_list,@fg_list,@bg_list\n” ;
push @url_list, $object_url;
$path_tmp = $path_t{path};
@obj_list=($path_t{gene_name});
@fg_list=(“#ff0000″);
@bg_list=(“#ffff00″);
if($i==scalar(@path_sort)-1){
$object_url=$soap_ser->color_pathway_by_objects($path_tmp,
\@obj_list,\@fg_list,\@bg_list);
push @url_list, $object_url;
}
}
}

#用于显示图片
foreach $url(@url_list){
print “<img src=$url>\n”;
}

print $up_reg.” genes up regulated. And “.$down_reg.” genes down regulated.”;
print “\nThis is the end!”;

#将Perl数组转换成SOAP数组
sub SOAP::Serializer::as_ArrayOfstring{
my ($self, $value, $name, $type, $attr) = @_;
return [$name, {'xsi:type' => 'array', %$attr}, $value];
}

sub SOAP::Serializer::as_ArrayOfint{
my ($self, $value, $name, $type, $attr) = @_;
return [$name, {'xsi:type' => 'array', %$attr}, $value];
}


运行结果如下:


gene_name    exp_ratio    path
sce:YDR468C    1.349167127091    path:sce04130
sce:YGL201C    1.06268290373857    path:sce03030
sce:YGL201C    1.06268290373857    path:sce04111
sce:YLR340W    2.44280064588198    path:sce03010

related pathways
path:sce04130
path:sce03030
path:sce04111
path:sce03010

sorted related pathways
sce:YLR340W    2.44280064588198    path:sce03010
sce:YGL201C    1.06268290373857    path:sce03030
sce:YGL201C    1.06268290373857    path:sce04111
sce:YDR468C    1.349167127091    path:sce04130
<img src=http://soap.genome.jp/tmp/mark_pathway_www_api.125119056924788/sce03010.png>

<img src=http://soap.genome.jp/tmp/mark_pathway_www_api.125119057124812/sce03030.png>

<img src=http://soap.genome.jp/tmp/mark_pathway_www_api.125119057324832/sce04111.png>

<img src=http://soap.genome.jp/tmp/mark_pathway_www_api.125119057524854/sce04130.png>

31 genes up regulated. And 0 genes down regulated.
This is the end!


作业完成。下一步,应该探索如何对基因进行功能分类。

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

历史上的今天

在LOFTER的更多文章

评论

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

页脚

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