GO数据库学习笔记(2)

#!/usr/bin/perl
#从键盘输入数据文件名并打开文件。其中<STDIN>指的就是标准输入设备。
print "please input the file name\n";
$filename = <STDIN>;
open(MICROARRAY, $filename)||die ("Could not open file");
#将文件一次性读入内存。
@whole_file = <MICROARRAY>;
#统计文件一共有多少行
$rows = scalar(@whole_file);
print "there are $rows entries\n";
#统计文件一共有多少以TAB为间隔的列。
$pattern ="\t";
@cells = split(/$pattern/, $whole_file[0]);
$cols = scalar(@cells);
print "there are $cols columns\n";
#设置计入统计的表达差异的比例。
$times = 1;
#设置脏标记用来退出程序。
$flag = 1;
while($flag==1){
#用first_col和second_col来设置需要比较的列。
$first_col=0;
$second_col=0;
print "please set the columns you want to compare\n";
while ($first_col<=0 || $first_col>=$cols){
print "the first column is:\n";
$first_col = <STDIN>;
}
$first_col--;
 
while ($second_col<=0 || $second_col>=$cols){
print "the second column is:\n";
$second_col = <STDIN>;
}
$second_col--;
#
$up_reg = 0;
$down_reg = 0;
@total_hash = ();
@go_hash = ();
$j = 0;
print "gene_name\texp_ratio\n";
#对文件进行逐行处理。
for($i=0;$i<$rows;$i++){
#对每一行的内容按TAB间隔存到数组words当中去。
@words = split(/$pattern/, $whole_file[$i]);
next if($words[$first_col]<=0 || $words[$second_col]<=0);
$ratio = log($words[$first_col])- log( $words[$second_col]);
#比较需要对比的两数据的差异，当差异大于两倍时进行处理。
if(abs($ratio)>=$times){
#对最后一列的注释进行处理，以空隔为分隔把内容存到description数组里去。
#当数据差异大于两倍时分为上调和下调分别进行统计。
@description = split(/\s+/,$words[$cols-1]);
my %reg_gene = ("gene_name"=>$description[2],"exp_ratio"=>$ratio);
push @total_hash, \%reg_gene;
printf ("%s\t%8.4f\n",$reg_gene{gene_name},$reg_gene{exp_ratio});
$j++;
$up_reg++ if($ratio>=$times);
$down_reg++ if($ratio<=-$times);
#查找GO信息，提取goid，并把每一个goid和对应的基因名存入哈希表go_hash当中去。
while($words[$cols-1]=~/\[goid\s(GO:\w+)\]/g){
$go_id = $1;
my %go = ("gene_name"=>$description[2],"exp_ratio"=>$ratio,"go_id"=>$go_id);
push @go_hash, \%go;
}
}
}
print "totally $j genes affected\n";
print $up_reg." genes up regulated. And ".$down_reg." genes down regulated.\n";
 
print "the GO info list as follows\n";
print "go_id\tgene_name\n";
#对存有go_id的哈希表go_hash进行排序，便于统计每个go_id一共有多少基因关联。
@go_sort = sort {$a->{go_id} cmp $b->{go_id}} @go_hash;
#统计一共有多少个go_id。
$sample_n = scalar(@go_sort);
#从GO数据库调取库中酿酒酵母的基因总数
$sample_N = &get_N();
#准备上调基因和下调基因哈希表。
%go_list_up = ();
%go_list_down = ();
#对每一个go_id进行统计
for($i=0;$i<$sample_n;$i++){
$tmp = $go_sort[$i];
%path_t = %$tmp;
#如果是上调，则把它压入上调的哈希统计表当中去，反之亦然。
$go_list_up{$path_t{go_id}} +=1 if ($path_t{exp_ratio}>=$times);
$go_list_down{$path_t{go_id}} +=1 if ($path_t{exp_ratio}<=-$times);
print "$path_t{go_id}\t$path_t{gene_name}\n";
}
#计算并显示所有上调和下调的GO统计。
&go_sql(1,$sample_N,$j,%go_list_up);
&go_sql(0,$sample_N,$j,%go_list_down);
 
print "\n\n0 for exit\n1 for continue\n";
$flag = <STDIN>;
}
 
close(MICROARRAY);
 
print "\nThis is the end!";
#从GO数据库当中调取属于酿酒酵母的基因总数
sub get_N{
use DBI;
$username = 'go_select';$password = 'amigo';$database = 'go_latest';
$hostname = 'mysql.ebi.ac.uk';$port='4085';
$dbh = DBI->connect("dbi:mysql:database=$database;" .
"host=$hostname;port=$port", $username, $password) or exit(1);
#count(*)统计出有多少条记录，把gene_product数据表和species数据表以species_id相连，查找到所有符合种属为酿酒酵母的记录条数。
$SQL = "SELECT count(*) AS n".
" FROM gene_product AS p".
" INNER JOIN species AS s ON (p.species_id=s.id)".
" WHERE s.genus='Saccharomyces'".
" AND s.species='cerevisiae'";
$sth = $dbh->prepare($SQL);
$sth->execute();
$sth->bind_columns(\$sample_no);
while ($sth->fetch()){
$sample_No=$sample_no;
}
$sth->finish();
$dbh->disconnect();
return $sample_No;
}
#查询，计算，并显示统计信息
sub go_sql{
use DBI;
my ($up_or_down,$N,$n,%list_sql)=@_;
$username = 'go_select';$password = 'amigo';$database = 'go_latest';
$hostname = 'mysql.ebi.ac.uk';$port='4085';
$dbh = DBI->connect("dbi:mysql:database=$database;" .
"host=$hostname;port=$port", $username, $password) or exit(1);
if($up_or_down){
print "\n\nup regulated genes\n";
}
else {
print "\n\ndown regulated genes\n";
}
print "\n\ngo_id\tm\tk\tN\tn\tp_value\tterm_type\n";
#对于记录中每条数据，都可以提取到两个数据，一是go_id，二是该go_id下有多少表达差异的基因。
while (($go_id_holder,$record) = each(%list_sql)){
#调取该go_id下属于酿酒酵母的基因总数，同时调取该go_id的注释。
$SQL= "SELECT gene_product_count.product_count, term.name".
" FROM term ".
" INNER JOIN gene_product_count ON (term.id=gene_product_count.term_id)".
" INNER JOIN species ON (gene_product_count.species_id=species.id)".
" WHERE term.acc='$go_id_holder'".
" AND species.genus='Saccharomyces'".
" AND species.species='cerevisiae'";
 
$sth = $dbh->prepare($SQL);
$sth->execute();
$sth->bind_columns(\$product_count,\$term_type);
while($sth->fetch()){
#计算p_value。这个值的计算方法还不确认是否正确，但我们实验室的人以前都是用的这个公式，不过与文献方法还是有差异。实验室的人用的matlab内的库函数实现的，而且调用的是本地数据库，这两方面我都必须换成自己的方式。
$p_value=&phyper($N,$product_count, $n, $record);
print "$go_id_holder\t$product_count\t$record\t$N\t$n\t$p_value\t$term_type\n";
}
$sth->finish();
}
$dbh->disconnect();
}
#超几何统计。超几何统计的原理可以百度google一下。
sub phyper{
($N,$m,$n,$k) = @_;
@x = (0, $m+1, $k+1, $N-$m+1, $n-$k+1, $N+1, $n+1);
if($x[2]>$x[1]){
$e=0;
}
else{
$e = &gammln($x[1]) - &gammln($x[2]) - &gammln($x[1]-$x[2]+1)
+ &gammln($x[3]) - &gammln($x[4]) - &gammln($x[3]-$x[4]+1)
- &gammln($x[5]) + &gammln($x[6]) + &gammln($x[5]-$x[6]+1);
$e = exp($e);
}
return $e;
}
#用于计算大数阶乘的近似算法，原程序来自numerical recipes in C第六单第一节。
sub gammln
{
my ($xx) = @_;
my ($x,$y,$tmp,$ser,$eye)=(0,0,0,0,0);
my @cof = (76.18009172947146,-86.50532032941677,
24.01409824083091,-1.231739572450155,
0.1208650973866179e-1,-0.5395239384953e-5);
$y=$x=$xx;
$tmp = $x+5.5;
$tmp -= ($x+0.5)*log($tmp);
$ser=1.000000000190015;
for ($eye=0;$eye<=5;$eye++){
$ser += $cof[$eye]/++$y;
}
return -$tmp+log(2.506628746310005*$ser/$x);
}