EnsEMBL Perl API用例 - mRNA剪接位点上的SNP
取自 PerlChina.org - wiki
EnsEMBL Perl API用例——mRNA剪接位点上的SNP
- 原 名：EnsEMBL Perl API用例——mRNA剪接位点上的SNP
- 作 者：Phil-
- 发 表：2008-1-27
- 审 校：
- 出 处：PerlChina.org
如对本文有任何意见，欢迎到讨论页留言，Emailing to me is also welcome
目录 |
[编辑] 前言
本文试图通过一个实际研究中的生物信息学问题，介绍EnsEMBL Perl API的具体应用方法。(关于EnsEMBL Perl API的基础，请看EnsEMBL Perl API简介)从编程技术上，主要希望介绍EnsEMBL Perl API和BioPERL的一些常用对象及其方法。同时，由于本文涉及实际的生物学问题，在您继续往下阅读之前，建议您先确认自己对下述几个分子生物学概念有所认识。（鉴于文章容量和主题所限，下面关于这几个概念的表述在细节上可能不严谨）
- mRNA ： 信使RNA。在细胞核内，以DNA为模板的转录之后，得到mRNA的前体。mRNA的前体经过剪接和其它形式的加工后成为成熟的mRNA，并被输出到核外，从而成为指导蛋白质合成（“翻译”）的模板。
- mRNA剪接 ： mRNA前体被加工成成熟mRNA的过程之一，主要是将前体中的内含子（intron）去掉，并将剩下的外显子（exon）按顺序拼接起来。
- mRNA剪接位点 ： 在intron的5'端和3'端各有两个保守的碱基，分别是"GT"和"AG"，他们定义了被剪切掉的内含子的边界。（但并不是所有GT/AG都是内含子边界）用SeqLogo表示如：Image:Http://www.lecb.ncifcrf.gov/~toms/gallery/SequenceLogoSculpture.gif（GT/AG原则）。本文中提到的“mRNA剪接位点”均指这四个碱基，而不指intron和exon的交界。
- intron ： 内含子。mRNA剪接中被“丢弃”的部分，但可能还有其他功能。
- exon ： 外显子。mRNA剪接中被保留的部分，内含子被去掉后，外显子会被拼接起来。
- SNP ： 单核甘酸多态性（Single Nucleotide Polymophism）。基因组中最常见的多态性，一般的形式是在基因组中的某个位置，在不同的个体中有不同的碱基，但在群体中一般只有两种可能的碱基。另外，一般要求两种可能的碱基中较稀有的那个在群体中出现的频率>1%，假如<1%则称该位点在某些个体中发生了突变（Mutation），而不会称之为多态。
[编辑] 生物学问题
mRNA剪接的分子过程中，intron边界上的mRNA剪接位点是严格保守的，假如在这些位点上存在SNP，则SNP的其中一种状态（不符合GT/AG原则的那一种）很可能会严重干扰正常的mRNA剪接，进而对基因的功能甚至生物个体的生理功能造成影响。利用EnsEMBL数据库，我们可以找到基因组中的这些SNP。本文以人类基因组为例，简介使用EnsEMBL定位这些SNP的方法。
[编辑] 基本思路
[编辑] 确定mRNA剪接位点
(本文中提到的“mRNA剪接位点”均指"GT/AG"四个碱基，而不指intron和exon的交界。) mRNA剪接位点是mRNA的exon到intron交界处往下游的第1,2个碱基，和intron到exon交界处往上游的第1,2个碱基。因此，我们要确定mRNA剪接位点就必须依次确定exon-intron/intron-exon的交界和mRNA的方向。（实际上，假如忽略接下来要定位SNP这一要求，mRNA的方向是不必确定的。但由于我们需要找SNP，而SNP是定位于基因组上的，因此我们必须确定mRNA剪接位点在基因组上的对应位置）
[编辑] intron-exon/exon-intron交界
在EnsEMBL Core中，所有的mRNA都对应一个Bio::EnsEMBL::Transcript对象，每个Bio::EnsEMBL::Transcript对象又包含若干个Bio::EnsEMBL::Exon对象，这些对象都包含各自在基因组上的坐标信息。据此，我们可以确定任何exon和与之相邻的intron的交界。
[编辑] mRNA的方向
在EnsEMBL Core中，mRNA的方向可以由对应的Bio::EnsEMBL::Transcript对象的"'strand'"方法获得。
[编辑] 寻找mRNA剪接位点上的SNP
SNP（还有其他多态/变异）在基因组上的坐标都记录在EnsEMBL Variation数据库中。一般而言，一个SNP及其在基因组上的位置信息对应一个Bio::EnsEMBL::VariationFeature对象（不考虑其位置信息则是Bio::EnsEMBL::Variation）。获取这一对象的一个方法是指定一个基因组片断（Bio::EnsEMBL::Slice对象），然后交由对应的数据适配器，从EnsEMBL数据库中获取。
[编辑] 代码实现
[编辑] 准备工作
- 首先连接到EnsEMBL的数据库，并获取相关的 DBAdaptor 对象
use Bio::EnsEMBL::Registry; my $registry = 'Bio::EnsEMBL::Registry'; $registry->load_registry_from_db(-host => 'ensembldb.ensembl.org',-user => 'anonymous'); my $transcript_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Transcript' ); #"Transcript":意即转录本，在这个问题里面就是一个mRNA my $slice_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Slice' ); #"Slice":顾名思义，就是基因组上的一个片段 my $exon_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Exon' ); #"Exon":外显子 my $variationfeature_adaptor = $registry->get_adaptor( 'Human', 'Variation', 'VariationFeature'); #"VariationFeature":一个多态及其在基因组上的位置（不一定是SNP）
[编辑] 获取转录本和Exons相关的信息
- 获取转录本，及转录本的strand,chromosome,exons等信息
my $transcript = $transcript_adaptor->fetch_by_stable_id("ENST00000007510");
my $strand = $transcript->strand;
my $chromosome = $transcript->seq_region_name;
my $all_exons = $exon_adaptor->fetch_all_by_Transcript($transcript);
- 获取所有exon的边界（就是相当于所有intron-exon/exon-intron交界，再加上Transcript的开头，末尾位置）
my @borders;
foreach my $this_exon (@$all_exons) {
push(@borders,$this_exon->seq_region_start,$this_exon->seq_region_end);
}
if($strand > 0) { # 因为在正链，所以borders应该从小到大排，其实如果是正链，这一步可以省略，这里写下来只是为了工整一点
@borders = sort {$a <=> $b} @borders;
} else { # 相反地，应该从大到小排
@borders = sort {$b <=> $a} @borders;
}
shift @borders; pop @borders; #去掉Transcript的开头，末尾
[编辑] 获取指定位置的SNP
- 根据exon border和strand，计算出严格保守的剪接位点在基因组上的位置，并获得该基因组位置对应的Slice，然后根据Slice获取VariationFeature
my $exon_cnt = 0;
while(my $this_exon_3_prime_end = shift @borders) { # 一个exon的3'末端
my $next_exon_5_prime_end = shift @borders; # 下一个exon的5'末端
$exon_cnt++;
my ($slice_intron_5_prime,$slice_intron_3_prime);
if($strand > 0 ) { #若是正链，剪接保守位点应该是exon 3'端的碱基+1位和+2位，或是exon 5'端的碱基-2位和-1位
$slice_intron_5_prime = $slice_adaptor->fetch_by_region('chromosome',$chromosome,$this_exon_3_prime_end+1,$this_exon_3_prime_end+2,$strand);
$slice_intron_3_prime = $slice_adaptor->fetch_by_region('chromosome',$chromosome,$next_exon_5_prime_end-2,$next_exon_5_prime_end-1,$strand);
} else { #若是负链，剪接保守位点应该是exon 3'端的碱基-2位和-1位，或是exon 5'端的碱基+1位和+2位
$slice_intron_5_prime = $slice_adaptor->fetch_by_region('chromosome',$chromosome,($this_exon_3_prime_end-2),($this_exon_3_prime_end-1),$strand);
$slice_intron_3_prime = $slice_adaptor->fetch_by_region('chromosome',$chromosome,($next_exon_5_prime_end+1),($next_exon_5_prime_end+2),$strand);
}
foreach my $this_variationfeature (@{$variationfeature_adaptor->fetch_all_by_Slice($slice_intron_5_prime)}) {
next unless ($this_variationfeature->var_class eq 'snp'); # 我们只看SNP
print $this_variationfeature->display_id." is on the 5' end of intron ".$exon_cnt."-".($exon_cnt+1)." of ".$transcript->stable_id."\n"; #输出
}
foreach my $this_variationfeature (@{$variationfeature_adaptor->fetch_all_by_Slice($slice_intron_3_prime)}) {
next unless ($this_variationfeature->var_class eq 'snp'); #同上
print $this_variationfeature->display_id." is on the 3' end of intron ".$exon_cnt."-".($exon_cnt+1)." of ".$transcript->stable_id."\n";
}
}
[编辑] 运行结果
运行后输出为：
rs12975039 is on the 5' end of intron 1-2 of ENST00000007510
[编辑] 讨论
- 显然，上述代码只对一个转录本（Transcript）进行了分析，Bio::EnsEMBL::DBSQL::TranscriptAdaptor（即上面代码中的$transcript_adaptor）可以列出数据库中的所有转录本，把这个方法用上就可以在全转录组范围进行这个分析。不过如果您想做全转录组范围的分析，建议您还是把EnsEMBL数据库本地化，这样你的分析不会由于网络质量而中断，也不会给EnsEMBL服务器带来太大负担以至于影响其他同行的使用。
- 上述代码中计算剪接保守位点的位置并获取相应的 Slice 的过程略嫌繁琐，本人能想到的一个相对快捷的方法是直接把 Exon 转成 Slice ，然后使用 Slice 的 expand 方法，有兴趣的朋友不妨一试。
- 更多关于EnsEMBL PERL API的文档在http://www.ensembl.org/info/using/api/Pdoc/index.html
