|
|
1#
大 中
小 发表于 07-12-13 10:14 只看该作者
MATLAB 7.X生物信息工具箱的应用
MATLAB 7.X生物信息工具箱的应用
——基因序列分析(一)
刘新星 李红燕 杨英杰
（中南大学资源加工与生物工程学院湖南 长沙 410083）
近年来，生物信息学的发展对生物学的发展起了巨大的推动作用。MathWorks公司推出的MATLAB7.0以上版本为广大科研工作者提供了一个用于生物信息分析的生物信息工具箱,它具有简单易学，操作方便而且功能强大等特点,即使不懂编程也能用它进行生物信息的分析研究[1]，用户可利用网络数据库资源，使研究工作更加事半功倍。本刊将分期系列的介绍MATLAB 7.X生物信息工具箱在基因序列分析、微阵列分析和系统发生分析等几个方面的应用，希望对从事基因组学、蛋白质组学和药物研究的科研人员和相关专业的师生有一些帮助。本文主要介绍了其在基因序列分析中的应用，包括确定核苷酸组成，密码子组成，氨基酸转化和组成等内容。
1.2MATLAB及生物信息工具箱简介
1.2.1MATLAB简介
MATLAB，全称为“Matrix Laboratory”，即矩阵实验室之意，它是由MathWorks公司推出的一款工具软件，功能强大，适合多学科，多工作平台的数据处理。
MATLAB被称为第四代计算机语言，它的一大特色是含有许多功能强大的工具箱。其工具箱分为两大类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能，图示建模仿真功能，文字处理功能以及硬件实时交互功能，功能性工具箱用于多种学科；而学科性工具箱的专业性比较强，这些工具箱都是由该领域内学术水平很高的专家编写的，生物信息工具箱就是学科性工具箱中的一种。生物信息工具箱提供了一个开放的环境，用户可以按照自己的目的来设计和利用分析工具。
1.2.2MATLAB生物信息工具箱
MATLAB生物信息工具箱（Bioinformatics Toolbox）提供了许多用于基因组学和蛋白质组学分析的函数，其中绝大多数是M代码的工具（MATLAB程序语言），并有可供查看的源代码；它扩展了MATLAB在生物信息方面的应用，提供了一个用于基因组和蛋白质组分析的综合软件环境。科研人员可以利用生物信息工具箱解决药物研究和设计，基因工程和生物信息方面的问题，并可利用这个工具箱提供的基本生物信息函数来创建更复杂的算法和应用程序。
1.2.3MATLAB生物信息工具箱特点及功能
生物信息工具箱可以进入许多网络数据库和其他的网络数据资源，它支持许多普通的基因组文件格式，用户可以直接复制序列和基因表达信息到MATLAB中。目前主要的序列数据库有国际基因库GenBank，蛋白质序列库GenPept，欧洲分子生物数据库EMBL，蛋白质数据库PSD 。用户也可以从同源蛋白家族数据库PFAM中得到多重比对序列，隐藏Markov模型图谱，系统发生树数据。此外，MATLAB还可以从基因测序仪、质谱仪和Agilent微阵列扫描仪上读取数据。
利用生物信息工具箱，用户可以选择一系列的分析方法来执行成对的或多重的序列比对,确定序列的一系列统计量，在一个序列中查找特定模型，或是查找开放阅读框。此外，还可以创建随机序列，从一组多序列比对的氨基酸、核苷酸序列中找出共同序列或序列文件，对序列格式化或用频率数据图解显示一个序列比对。MATLAB的附加功能还可以用正则表达式有效地处理字符串，以便在一个序列中查找特定模型，并查找串匹配的信息库。通过查找回文序列在DNA/RNA序列中查找可能的断点。
用户可以使用一系列的蛋白质分析方法从数据中选取信息。工具箱提供了用于计算蛋白质序列性质的多个函数，比如基本组成，分子量和等电点。用户可以用酶裂解一个蛋白质并为PDB数据创建距离和拉氏图。工具箱中包含了用于蛋白质分析和使用PDB数据库绘制蛋白质3D结构图的图形用户界面。用户可以预测氨基酸序列的统计量并获得有关特征的编码信息；计算两生物序列间的距离和计算替换率；利用距离数据构建系统发生树，在交互性的图形用户界面里观察系统发生树和编辑数据；也可以在这个图形用户界面中修剪分支、重排、改名和测量距离。
总之，MATLAB可用来执行成对序列或多序列的比对，进行序列转换，绘制序列图谱，进行蛋白质分析和氨基酸序列分析，创建并分析系统发生树，进行微阵列数据分析等。用户还可以创建自己的算法和应用程序，并与其他用户分享。
2序列分析
序列分析是利用计算机方法来寻找有关核苷酸或氨基酸序列的信息。序列分析的一般工作是基因识别，确定两个基因的相似性，确定一个基因的蛋白质编码以及研究另一有机体中相似基因的功能[2]。
2.1 序列分析
在分析完一段DNA序列之后，首要任务就是研究序列中的核苷酸含量。本节使用序列统计函数来确定核苷酸含量，并找出开放阅读框。
如无特别说明，均以人类线粒体基因组为例，逐个介绍序列分析的相关函数。
2.1.1搜索网络数据库资源
首先需要搜索网络数据库资源，查找有关人类线粒体的信息，找出基因组的核苷酸序列。
1） 连接网络：可用web 函数连接到网络，下面的命令以一个独立的浏览器窗口打开NCBI网站的主页。
web('http://www.ncbi.nlm.nih.gov/')；
2） 查找信息：如在NCBI网站上查找人类线粒体基因组，可在搜索列表【Search】选择基因
【Genome】，在目的栏【for】输入人科线粒体【mitochondrion homo sapiens】后进行搜索，如图2-1所示。
图 2-1查找相关信息 NCBI网站查找并返回一系列相关页的链接，如图2-2。
图 2-2查找结果
3） 显示信息：点击链接标签NC_001807即显示人类线粒体基因组信息页面，图2-3给出了显示页面的主要部分。
图 2-3线粒体基因信息
2.1.2获取序列信息
MATLAB提供了一个读取序列信息的getgenbank函数。由于国际基因库GenBank 的条目数量非常巨大，而用户也许只对特定的序列感兴趣，所以利用getgenbank函数可以通过登陆号码查找特定序列的信息。
比如，人类线粒体基因组共同序列的GenBank的登陆号码为NC_001807，可用getgenbank函数从GenBank网络数据库中搜索人类线粒体基因组的序列信息，并读入MATLAB工作区。
mitochondria= getgenbank('NC_001807','SequenceOnly',true)
mitochondria =gatcacaggtctatcaccctattaaccactcacgggagctctccatgcatttggtattttcgtctggggggtgtgcacgcgatagcattgcgagacgctggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctactaaagt . . .
2.1.3确定核苷酸的组成
一个含A+T丰富的核苷酸DNA片段通常是序列的一部分，而含A+T低，含G+C丰富的核苷酸则是潜在的基因。通常，一个基因的CG二核苷酸含量都是已经确定的。在读取一段序列到MATLAB中以后，用户可以使用序列统计函数确定这个序列是否含有蛋白质编码域的特征。
1） 绘制密度图：可用ntdensity函数绘制单体密度和联合体密度图，图2-4显示这个基因组富含A+T。
ntdensity(mitochondria)；
图 2-4单体和联合体密度图 2）计算核苷酸数目：可用basecount函数计算5’-3’链中的核苷酸数目。
basecount(mitochondria)
ans =
A: 5113
C: 5192
G: 2180
T: 4086
3）计算互补核苷酸数目：可用seqrcomplement函数计算互补5’-3’链中的核苷酸数目。
basecount(seqrcomplement(mitochondria))
ans =
A: 4086
C: 2180
G: 5192
T: 5113
4）显示核苷酸分布：可用basecount函数显示核苷酸分布的饼状图，如图2-5所示。
basecount(mitochondria,'chart','pie')；
图 2-5核苷酸分布饼状图 图 2-6二聚体数目条形图
5） 计算二聚体个数：可用dimercount函数计算一个序列中的二聚体个数,并在一个条形图中显示出来，见图2-6。
dimercount(mitochondria,'chart','bar')；
2.1.4确定密码子组成
三核苷酸（密码子）编码一个氨基酸，在一个核苷酸序列中有64个可能的密码子。知道序列中密码子的百分比有助于用户比较预设的密码子排列。
1）计算密码子数目：可用codoncount函数计算一个核苷酸序列中的密码子数目。
codoncount(mitochondria)
MATLAB即显示出第一个阅读框的密码子数目。
AAA
| 172 AAC
| 157 AAG
| 67 AAT
| 123
| ACA
| 153 ACC
| 163 ACG
| 42 ACT
| 130
| AGA
| 58 AGC
| 90 AGG
| 50 AGT
| 43
| ATA
| 132 ATC
| 103 ATG
| 57 ATT
| 96
| CAA
| 166 CAC
| 167 CAG
| 68 CAT
| 135
| CCA
| 146 CCC
| 215 CCG
| 50 CCT
| 182
| CGA
| 33 CGC
| 60 CGG
| 18 CGT
| 20
| CTA
| 187 CTC
| 126 CTG
| 52 CTT
| 98
| GAA
| 68 GAC
| 62 GAG
| 47 GAT
| 39
| GCA
| 67 GCC
| 87 GCG
| 23 GCT
| 61
| GGA
| 53 GGC
| 61 GGG
| 23 GGT
| 25
| GTA
| 61 GTC
| 49 GTG
| 26 GTT
| 36
| TAA
| 136 TAC
| 127 TAG
| 82 TAT
| 107
| TCA
| 143 TCC
| 126 TCG
| 37 TCT
| 103
| TGA
| 64 TGC
| 35 TGG
| 27 TGT
| 25
| TTA
| 115 TTC
| 113 TTG
| 37 TTT
| 99
|
2）绘制热红外分布图：可用下列程序绘制热红外分布图显示出6个阅读框中的所有的64个密码子，见图2-7(a-c)。
for frame = 1:3
figure('color',[1 1 1])
subplot(2,1,1);
codoncount(mitochondria,'frame',frame,'figure',true);
title(sprintf('Codons for frame %d',frame));
subplot(2,1,2);
codoncount(mitochondria,'reverse',true,'frame',frame,'figure',true);
title(sprintf('Codons for reverse frame %d',frame));
end
(a) (b) (c)
图 2-7密码子分布图
2.1.5开放阅读框
为一个真核基因确定蛋白质编码序列是一项困难的工作，因为内含子和外显子是间杂的。但是，原核基因通常是没有内含子的，mRNA序列的内含子是移动的。通过识别翻译的起始密码子和终止密码子，可以确定序列中蛋白质编码段，即开放阅读框ORF。一旦用户知道基因或者mRNA的ORF，就可以将一个核苷酸序列转化成相应的氨基酸序列。
1）显示核苷酸序列的ORF：可用seqshoworfs函数显示核苷酸序列的ORF。
seqshoworfs(mitochondria)
将命令的执行结果与NABI网页上NC_001807的基因作比较，会发现基因比想像中的要少。这是由于脊椎动物线粒体的遗传密码与标准遗传密码稍有不同。
2） 显示脊椎动物线粒体编码的ORF：
orfs= seqshoworfs(mitochondria,...
'GeneticCode','Vertebrate Mitochondrial',...
'alternativestart',true)
执行结果中，第一个阅读框有两个比较大的ORF。一个的开始位置是4471，另一个的开始位置是5905，它们分别对应于基因ND2（人类NADH脱氢酶亚基2）和COX1（细胞色素C氧化酶亚基1）。
3）查找终止密码子：可用find函数找出相应的终止密码子，ORF的起始和终止位置与起始和终止域的开始位置都有一样的引物。
ND2Start = 4471;
StartIndex = find(orfs(1).Start == ND2Start);
ND2Stop = orfs(1).Stop(StartIndex)
ND2Stop =
5512
4） 摘录子序列：利用基因的起始和终止位置的序列引物，从序列中摘录子序列。下面语句将子序列（蛋白质编码区域）摘录到ND2Seq中，并被显示在屏幕上。
ND2Seq = mitochondria(ND2Start:ND2Stop)
attaatcccctggcccaacccgtcatctactctaccatctttgcaggcacactcatcacagcgctaagctcgcactgattttttacctgagtaggcctagaaataaacatgctagcttttattccagttctaaccaaaaaaataaaccctcgttccacagaagctgccatcaagtatttcctcacgcaagcaaccgcatccataatccttc . . .
5） 确定密码子分布：下面的密码子计算结果显示ACC,ATA,CAT,ATC的含量较丰富。
codoncount (ND2Seq)
AAA
| 10 AAC
| 14 AAG
| 2 AAT
| 6
| ACA
| 11 ACC
| 24 ACG
| 3 ACT
| 5
| AGA
| 0 AGC
| 4 AGG
| 0 AGT
| 1
| ATA
| 22 ATC
| 24 ATG
| 2 ATT
| 8
| CAA
| 8 CAC
| 3 CAG
| 2 CAT
| 1
| CCA
| 4 CCC
| 12 CCG
| 2 CCT
| 5
| CGA
| 0 CGC
| 3 CGG
| 0 CGT
| 1
| CTA
| 26 CTC
| 18 CTG
| 4 CTT
| 7
| GAA
| 5 GAC
| 0 GAG
| 1 GAT
| 0
| GCA
| 8 GCC
| 7 GCG
| 1 GCT
| 4
| GGA
| 5 GGC
| 7 GGG
| 0 GGT
| 1
| GTA
| 3 GTC
| 2 GTG
| 0 GTT
| 3
| TAA
| 0 TAC
| 8 TAG
| 0 TAT
| 2
| TCA
| 7 TCC
| 11 TCG
| 1 TCT
| 4
| TGA
| 10 TGC
| 0 TGG
| 1 TGT
| 0
| TTA
| 8 TTC
| 7 TTG
| 1 TTT
| 8
|
6）查找编码的氨基酸：可用aminolookup函数查找给定密码子编码的氨基酸，下面是查找密码子ATA,CTA,ACC和ATC编码的氨基酸。
aminolookup('code',nt2aa('ATA'))
aminolookup('code',nt2aa('CTA'))
aminolookup('code',nt2aa('ACC'))
aminolookup('code',nt2aa('ATC'))
Ile isoleucine
Leu leucine
Thr threonine
Ile isoleucine
2.1.6氨基酸转化和组成
确定蛋白质相关的氨基酸组成可以提供给用户蛋白质特征图谱。通常，这个图谱含有足够的用来识别蛋白质的信息。利用氨基酸组成、基本成分和分子量，用户可以在公共数据库中查找类似的蛋白质。在用户定位基因的一个ORF之后，就可以将它转换成一个氨基酸序列，并确定它的氨基酸组成。
1）核苷酸序列转氨基酸序列：可用nt2aa函数将核苷酸序列转换为氨基酸序列。在下面的例子中，利用脊椎动物线粒体遗传密码将ND2Seq序列转换，仅起始密码子和终止密码子间的蛋白质编码序列被转换。由于默认情况下，交替的起始密码子性质被设定为“真”，所以第一个密码子ATT被转化为M而非I。
ND2AASeq = nt2aa(ND2Seq,'geneticcode','Vertebrate Mitochondrial')
MNPLAQPVIYSTIFAGTLITALSSHWFFTWVGLEMNMLAFIPVLTKKMNPRSTEAAIKYFLTQATASMILLMAILFNNMLSGQWTMTNTTNQYSSLMIMMAMAMKLGMAPFHFWVPEVTQGTPLTSGLLLLTWQKLAPISIMYQISPSLNVSLLLTLSILSIMAGSWGGLNQTQLRKILAYSSITHMGWMMAVLPYNPNMTILNLTIYIILTTTAFLLLNLNSSTTTLLLSRTWNKLTWLTPLIPSTLLSLGGLPPLTGFLPKWAIIEEFTKNNSLIIPTIMATITLLNLYFYLRLIYSTSITLLPMSNNVKMKWQFEHTKPTPFLPTLIALTTLLLPISPFMLMIL
2）转化结果比较：可用getgenpept函数把转化结果与GenPept中已公布的转化结果作比较，下面的语句从NCBI数据库中获取已公布的转化结果并读取入MATLAB工作区。
ND2protein = getgenpept('NP_536844','sequenceonly',true);
3）计算氨基酸数目：可用aacount函数计算蛋白质序列中的氨基酸数目，并绘制出条形图。图2-8显示亮氨酸，苏氨酸和异亮氨酸的含量较高，而半胱氨酸和天冬氨酸的含量较低。
aacount(ND2AASeq, 'chart','bar');
图 2-8氨基酸数目条形图
4） 确定氨基酸组成和分子量：可用atomiccomp函数和molweight函数确定蛋白质的氨基酸组成和分子量。
atomiccomp(ND2AASeq)
ans =
C: 1818
H: 3574
N: 420
O: 817
S: 25
molweight (ND2AASeq)
ans =
3.8960e+004
如果此序列是未知的，用户还可以通过将其基本组成与数据库中其他蛋白质作比较来识别蛋白质。
3小 结
本文简单介绍了生物信息学的基本情况和生物信息工具箱的功能特点，初步介绍了序列分析的基本操作，关于序列比对等进一步的内容将在后续文章中逐一阐述。为了便于读者查阅，表3-1给出了本文介绍过的序列分析的生物信息学函数。
表3-1 序列分析的生物信息学函数（1）
函
数
| 功
能
| 函
数
| 功
能
| web
| 打开浏览器显示网页或文件
| seqshoworfs
| 显示序列中的开放阅读框
| getgenbank
| 返回NCBI的基因数据库中的序列信息
| codoncount
| 计算序列中的密码子数目
| ntdensity
| 绘制序列中核苷酸分布的密度图
| aminolookup
| 查找特定氨基酸名称或密码子等
| basecount
| 计算序列中的核苷酸数目
| nt2aa
| 将核苷酸序列转换为氨基酸序列
| seqrcomplement
| 计算DNA互补链中核苷酸数目
| getgenpept
| 返回NCBI的蛋白质数据库中的序列信息
| molweight
| 计算蛋白质的分子量
| aacount
| 计算序列中氨基酸数目
| dimercount
| 计算序列中二聚体的数目
| atomiccomp
| 计算蛋白质序列的成分组成
|
参考文献
[1]罗军辉, 冯平, 哈力旦·A等. MATLAB 7.0在图像处理中的应用[M]. 北京：机械工业出版社，2005.
[2] （美）Mount, David W.Bioinformatics: sequence and genome analysis = 生物信息学 ：序列与基因组分析.北京：科学出版社,2006.
附件: 您所在的用户组无法下载或查看附件
搜索更多相关主题的帖子:
MATLAB 生物 工具箱 应用 MATLAB 生物 工具箱 应用
生命的精彩在于大家共同的努力～～
|
