bioconductor系列教程之一分析基因芯片中(质量控制)

上一节,我们了解了分析基因芯片的预处理的基本知识。其实那只是一个热身。这一节,我们来学习拿到基因芯片数据时更基本的操作:质量控制。只有通过质量检测合格的芯片数据才会真正地进入数据分析的步骤。本节将学习以下内容:

背景

MAS5 标准化

Affymetrix公司制定的内参

教程数据下载

质量控制总览图及报告

使用FitPLM生成权重,残差及NUSE图像

RNA降解曲线及MVA线图

PCA分析

总结

UTF8[……]

Read more

bioconductor系列教程之一分析基因芯片上

我们先用一个图来描述一下基因芯片(通常指Affymetrix microarray chips)实验的流程。

Figure 1基因芯片实验的流程
Figure 1基因芯片实验的流程

在这个流程中,预处理当中图像分析所指的图像就是基因芯片实验完成之后由仪器对杂交结果照像生成的图片。生物芯片实际上就是一堆基因特异探针的集成。当mRNA或cDNA文库标记后与这些探针相杂交,经洗脱,与探针有特异性结合的cDNA就保留了下来,因为它本身带有荧光标计,所以可以被成像系统捕捉下来。为了减少实验产生的噪音,Affy基因芯片在设计时,对每个基因都集成了20组不同的探针,这些特异性探针每个都由25个碱基组成,被称为Perfect Match(PM)探针。与其直接相邻的,也是20组不同探针,这些探针与前面所提的探针序列基本相同,只是正中间的一个碱基被换成了非特异的碱基,理论上来讲,它们应该很难与特异性探针形成相同的阳性结果。而这20组探针被称为MisMatch(MM)探针。有了MM探针,基因芯片就有了阴对内参,可以用于对杂交效果的检测。[……]

Read more

在blast中E值(E value)是什么?

一般的,当我们使用BLAST(是一种用于在数据库当寻找任何蛋白质或者基因序列与你的目标序列一致的程序)时,我们会注意到这里有一个E值。那么这个E value是什么呢?怎么来理解这个值呢?

下面是一个平常的blast结果,

Sequences producing significant alignments:
Score (S)
E

gi|83574104|Moth[……]

Read more

用perl实现细菌蛋白质定位预测批处理

首先声明,这里并不是我用perl来写个算法。如果需要算法相关的东西,得去看文献。
这里,其实只是借对格兰氏阳性细菌蛋白质定位的预测为例,来说明如何利用他人提供的在线查询功能来实现批处理。

假设现在有一个完整的基因组需要您去预测它当中的每一个开放阅读框翻译出来的蛋白质可能的细胞内定位,怎么办呢?手工一个一个提交到网站上去?一共会有四五千个蛋白,等你提交完,你的手和大脑都会不工作了吧?要不自己下载个软件来本地预测吧?我试过去安装那些要求的软件环境,也许是我的系统过新吧,一个c语言库的版本,一个g77让我就头大得不知道该怎么继续下去。于是我还是下定决心,用perl的lwp来伪装成浏览器提交申请,自动批处理吧。也许一觉醒来,全部都做完了。

我们要用到的网站是http://www.psort.org/psortb/。据网站上宣传,Based on a study last performed in 2010, PSORTb v3.0.2 is the most precise bacterial localization prediction tool available. 第二个原因就是can currently submit one or more Gram-positive or Gram-negative bacterial sequences or archaeal sequences in FASTA format。这对于研究格兰氏阳性菌的我来说的确很不错的网站。

最基本的,用lwp来虚拟提交一份表单上去:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#!/usr/bin/perl
use LWP::UserAgent;
$ua = LWP::UserAgent->new;
$ua->agent("Mozilla 8.0 beta");
 
use HTTP::Request::Common qw(POST);
my $req = (POST 'http://www.psort.org/psortb/results.pl',
		["seq" => ">$genes{name}\n$genes{translation}",
		"organism" => "bacteria",
                "gram"=>"positive",
                "advancedgram"=>"none",
                "format"=>"long",
                 "sendresults"=>"display");
 
$request = $ua->request($req);
print $request->as_string;

结果得到的,怎么都是500 Internal Server Error。于是在网上狂google,也没有找到直接的答案。半夜两点,突然想起,为什么不自己想法来解决问题。于是开始用tcpdump抓包,对比从firefox发送出去的包和perl发送出去的包,具体看看有什么不同。[……]

Read more

GO数据库学习笔记(4)

上一次学习了什么DAGs,这次我想实践一下学习的东西。
题目:对于任意给出的属于某属的多个基因,用树状图表示这些基因的功能所属关系。
实例程序如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
#!/usr/bin/perl
#使用DBD-mysql联接gene ontology数据库。
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);
#随机给出一些基因,基因的列表越长,所需时间自然也就越长,显示结果也需要太多的位置,
#所以这里只示例只有五个基因的情况。
#这些基因都来自酿酒酵母。
@gene_list =('AGA2','AIM1','BUD31','BUD32','BUD4');
#sql语句。从gene_product表(这个表内是基因和基因产物相关的记录)出发,
#以gene_product_id为关键字连接association表(gene_product与term的关系表,表中is_not是很有用的信息,如果它的值为真,则表示这个关系没有有效的证据支持。),
#以species_id为关键字连接species表,用species表中的genus和species项来确认数据是属于酿酒酵母的,
#用gene_product表中的symbol项来匹配基因列表中的基因名,用binary来表示基因名区分大小写,
#返回的结果只提取不重复的symbol(基因名),term_id(GO term ID)和product_count(该属中属于该term ID的所有基因的个数)。
$SQL = "SELECT DISTINCT g.symbol,t.term_id,t.product_count FROM gene_product AS g INNER JOIN association AS a ON (g.id=a.gene_product_id) INNER JOIN gene_product_count AS t ON (a.term_id=t.term_id) INNER JOIN species AS s ON (t.species_id=s.id) WHERE s.genus='Saccharomyces' AND s.species='cerevisiae'  AND a.is_not=0 AND BINARY g.symbol IN ('";
$SQL .=join ("','",@gene_list);
$SQL .= "')";
$sth = $dbh->prepare($SQL);
$sth->execute();
#把返回的结果全部都存到db二维数组当中去。
@db = ();
while (@ary = $sth->fetchrow_array()){
my %gene_entry = ("symbol"=>$ary[0],"term_id"=>$ary[1],"count"=>$ary[2]);
push @db, [$gene_entry{symbol},$gene_entry{term_id},$gene_entry{count}];
}
#对db二维数组当中的数据以term_id为索引做统计。
#用哈希表count_list保存某term_id下返回的基因总数(product_count)
%count_list=();
#用哈希表count保存这次所调查的基因中属于某term_id的基因数。
%count=();
#用哈希表gene_list保存属于某term_id的基因名。
%gene_list=();
$sql = "";
for $aref(@db){
$sql.="'@$aref[1]'," if ($sql !~ /@$aref[1]/);
$count_list{@$aref[1]}=@$aref[2];
$count{@$aref[1]}+=1;
$gene_list{@$aref[1]}.=@$aref[0].",";
}
#建立SQL语句,返回所有term_id对应的term_name。
$sql_t = "SELECT term.id, term.name FROM term WHERE term.id IN (".$sql."'1')";
#建立SQL语句,返回所有term_id的所属关系和步长。
$sql = "SELECT t.term1_id,t.term2_id,t.distance FROM graph_path AS t WHERE t.term1_id IN (".$sql."'1') AND t.term2_id IN (".$sql."'1') AND t.term1_id<>t.term2_id";
$sth_p = $dbh->prepare($sql);
$sth_p->execute();
#将返回的所属关系和步长保存到graph_path二维数组当中。
@graph_path=();
while(@ary=$sth_p->fetchrow_array()){
my %term2term=("term1_id"=>$ary[0],"term2_id"=>$ary[1],"distance"=>$ary[2]);
push @graph_path, [$term2term{term1_id},$term2term{term2_id},$term2term{distance}];
#  print join ("\t",@ary),"\n";
}
#对graph_path二维数组以步长升序排序。
@path_sort = sort {@$a[2] <=> @$b[2]} @graph_path;
$sth_t = $dbh->prepare($sql_t);
$sth_t->execute();
#将返回的term_id和term_name对应列表保存到哈希表term_list当中去。
%term_list=();
while(@ary=$sth_t->fetchrow_array()){
$term_list{$ary[0]}=$ary[1];
}
#对排序后的path_sort列表以term_id为关键字增加term_name,某term_id的基因数,某term_id下返回的基因总数(product_count),某term_id下的基因名,和一个脏标记。
for($i=0;$i<@path_sort;$i++){
$path_sort[$i][3]=$term_list{$path_sort[$i][0]};
$path_sort[$i][4]=$term_list{$path_sort[$i][1]};
$path_sort[$i][5]=$count{$path_sort[$i][1]};
$path_sort[$i][6]=$count_list{$path_sort[$i][1]};
$path_sort[$i][7]=$gene_list{$path_sort[$i][1]};
$path_sort[$i][8]=0;
}
#for $aref(@path_sort){
#   print "@$aref,\n";
#}
#对path_sort二维数组中的数据依据所属关系重新保存到tree_view表当中去
print "\nroot.\n";
&draw_path(1,0,\@path_sort,\@tree_view);
#合并统计相关数据
for $hit(@tree_view){
@daughter_gene_ary = split(/,/,$gene_list{$hit->{this_id}});
@parent_gene_ary = split(/,/,$gene_list{$hit->{parent_id}}) if($hit->{parent_id}!=1);
@parent_gene_ary = (@parent_gene_ary,@daughter_gene_ary);
%screen = ();
@parent_ary = grep (!$screen{$_}++, @parent_gene_ary);
$gene_list{$hit->{parent_id}}= join (",",@parent_ary);
for $par(@tree_view){
if ($par->{this_id}==$hit->{parent_id}){
$par->{genes_name}=$gene_list{$par->{this_id}};
$count{$par->{this_id}}=@parent_ary;
$par->{this_test}=$count{$par->{this_id}};
}
}
}
#画出树状图。
for $hit(@tree_view){
print " " x $hit->{pre_distance};
$temp_l=length($hit->{parent_name})-1;
print " " x $temp_l;
print "|";
print "-" x $hit->{this_distance};
print "$hit->{term_name}(step:$hit->{this_distance})($hit->{this_test}/$hit->{total})($hit->{genes_name})\n";
}
 
$sth_p->finish();
$sth_t->finish();
$sth->finish();
$dbh->disconnect();
print "\n\nthis is the END!";
 
#层进递归生成tree_view数组,用以生成树状图
sub draw_path{
my ($parent,$distance,$path_t,$tree_v)=@_;
my ($sep) = "-";
my ($pre) = $sep x $distance;
for(my $i=0;$i<@$path_t;$i++){
if($$path_t[$i][0]==$parent){
my $flag = 1;
$flag = 0 if(&get_step($$path_t[$i][0],$$path_t[$i][1],0,$$path_t[$i][2],$path_t));
if($flag){
$$path_t[$i][8]=1;
my %tree = ("parent_id"=>$$path_t[$i][0],"this_id"=>$$path_t[$i][1],"pre_distance"=>$distance,"this_distance"=>$$path_t[$i][2],"parent_name"=>$$path_t[$i][3],"term_name"=>$$path_t[$i][4],"this_test"=>$$path_t[$i][5],"total"=>$$path_t[$i][6],"genes_name"=>$$path_t[$i][7]);
push @$tree_v,\%tree;
my $out ="$pre$$path_t[$i][3]";
$out .=  $sep x $$path_t[$i][2];
draw_path($$path_t[$i][1],length($out),$path_t,$tree_v);
}
}
}
}
 
#防止tree_view当中产生重复数据的递归函数。
sub get_step{
my($parent,$daughter,$step,$step_t,$path_t)=@_;
my $flag;
for(my $i=0; $i<@$path_t;$i++){
if($$path_t[$i][8]&&$$path_t[$i][1]==$daughter){
$step += $$path_t[$i][2];
if($$path_t[$i][0]==$parent){
if ($step==$step_t){
return $flag=1;
}#if
else{
$step=0;
next;
}#else
}#if
elsif($$path_t[$i][0]==1){
$step=0;
next;
}#elsif
else{
if(get_step($parent,$$path_t[$i][0],$step,$step_t,$path_t)){
return $flag = 1;
}
else{
$step=0;
next;
}
}#else
}#if
}#for
return $flag;
}

两个递归没把我的头想破了……递归怎么就这么绕呢?从一开始C的时候学递归就学不利索。
现在看一下结果吧:

root.
|-molecular_function(s[……]

Read more

GO数据库学习笔记(3)

为了搞清楚gene ontology database的数据结构,必须得先了解清楚一个概念,什么是Directed Acyclic Graphs (DAG)。这个DAG被翻译成有向无环图,具体是个啥意思呢?
对于一个DAG来讲,有以下两个特征:

  1. 有向。数据间有指向性关系,A属于B,或者说A就是B的一种形式之类的。
  2. 无环。当你从任何一点开始,顺着指向走的话永远不可能走回原点。

下[……]

Read more

GO数据库学习笔记(2)

这次要对数据进行处理了。原代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
#!/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);
}

运行以后,得到结果:

please input the file name
c:\\testfile.txt
there are 5716 entries
there are 7 columns
please set the columns you want to com[……]

Read more

GO数据库学习笔记(1)

前几天学习了KEGG API之后,感觉这个数据虽然好,但是还是有些局限性,比如它比较重视代谢途径而会遗漏其它信息。为了找到更完整的基因注解信息,同学为我介绍了GO(the Gene Ontology)数据库。如果想快速了解这个数据库能提供些什么有用的信息,可以去http://www.geneontology.org/GO.teaching.resources.shtml?all下载一些它的PPT来[……]

Read more

KEGG API学习笔记(3)

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

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

$file[……]

Read more