LINUX
2006-03-05 11:27:47 来源:WEB开发网摘要:
生物信息学(Bioinformatics)
生物信息学开始于科学家们将生物学数据以数字格式存放并且用程序来处理这些数据。很长一段时间以来,生物信息学都限制在序列的分析上。然而,随着构建分子的结构模型的重要性开始显现,电子计算机也开始成为理论生物化学的重要工具。每天都不断有关于分子3D信息和数据被采集,人们对基因的认识和研究也从单个的基因研究转变为从整体上或者扩展式的研究。由于生物信息学的发展,现在更容易理解蛋白质之间为何相互之间那样作用、又是如何通过新陈代谢来组织相互的。而且我们现在也越来越清醒的认识到组织好这些数据的重要性。
生物信息学里面至少有两点特征使她变得非常有趣。其一,生物信息学的研究目标是找出各种生命分子的关系;而这个目标恰恰是一个有趣的程序设计问题,因为这就需要我们联合并整合我们得到的那些信息,然后从中得到对生命活动的一些整体的和有效的一些认识。我们还发现,将计算机科学中的不同领域的知识结合起来是非常必要的,比如数据的管理和整合、高效可靠的算法和强劲的硬件--格点技术、多处理器的使用等。
Perl
LarryWall于1986年开始开发Perl。Perl是一种解释型的语言,是处理文本、文件和进程的强大的工具。Perl使得我们能够很快的开发出小程序。可以说,Perl是高级编程语言(例如C)和脚本语言(如bash)的一种有效组合。
Perl程序可以运行在多种操作系统/平台上,尽管Perl是在Unix上诞生并且快速发展的。由于Perl广泛的用于web程序设计,其发展很快便超出了其预想。在Perl之前,人们使用awk,thirst和grep来分析文件并提取信息。
Perl将这些UNIX上广泛使用的工具统一在一个程序里面,并将这些功能扩展和现代化以适应各种需求。
Perl是一种免费/自由的程序语言,可以运行在现代生物实验室里使用的各种操作系统上。在UNIX和MacOSX上,它是预安装好的,在其他系统上,得先安装好Perl。
在linux下,运行Perl程序,是将这个程序的文件名作为perl这个命令的一个参数,然后perl会依次解释执行这个程序里的命令。
另一种常用的方法,不需要运行perl这个命令,为此,我们需要做以下两件事:(a)在程序的文件里加入一行特殊的注释:
PRint"Hi\n";
(b)保存此文件并给它加上可执行的属性:
这样,我们就可以直接通过文件名来运行这个程序:
用Perl来文件管理:
当我们有了文本格式的分子序列,我们可以用Perl写一个序列搜索工具。下面的例子我们可以看到如何在SWISS-PROT(db_human_swissprot)格式的数据库中用id码来查找蛋白质序列。
#Lookforaminoacidsequenceinadatabase
#SWISS-PROTformated,withagivenidcode
#AskforthecodeintheIDfield
#anditassignsitfromtheinput(STDIN)toavariable
print"EntertheIDtosearch:";$id_query=<STDIN>;chomp$id_query;#Weopenthedatabasefile
#butifitisn'tpossibletheprogramends
open(db,"human_kinases_swissprot.txt")||die"problemopeningthefilehuman_kinases_swissprot.txt\n";#Looklinebylineinthedatabase
while(<db>){chomp$_;#CheckifweareintheIDfieldif($_=~/^ID/){#Ifitispossitivewegathertheinformation
#breakingthelinebyspaces
($a1,$id_db)=split(/\s /,$_);#butifthereisnocoincidenceofIDwecontinuetothefollowing
nextif($id_dbne$id_query);#Whentheycoincide,weputamark
$signal_good=1;#Thenwecheckthesequencefield
#andifthemarkis1(chosensequence)#Ifpossitive,wechangethemarkto2,tocollectthesequence
}elsif(($_=~/^SQ/)&&($signal_good==1)){$signal_good=2;#Finally,ifthemarkis2,wepresenteachline
#ofthesequence,untilthelinebeginswith//#issuchcasewebrokethewhile}elsif($signal_good==2){lastif($_=~/^\/\//);print"$_\n";}}#Whenweleftthewhileinstructionwecheckthemark
#ifnegativethatmeansthatwedon'tfindthechosensequence
#thatwillgiveusanerror
if(!$signal_good){print"ERROR:"."Sequencenotfound\n";}#Finally,weclosethefile#thatstillsiopen
close(db);exit;
查找氨基酸的模式(Searchforaminoacidpatterns)
if($_=~/^SQ/){
$signal_seq=1;#Whenarrivetotheendofsequence,leavethecurl
#Checkthatthisexpressionisputbeforetocheck
#themark=1,becausethislinedoesn'tbelongtotheaminoacidsequence
}elsif($_=~/^\/\//){
last;#Checkthemarkifitisequalto1,ifpossitive
#eliminatetheblankspacesinthesequenceline
#andjoineverylineinanewvariable
#Toconcatenate,wealsocando:
#$secuencia_total.=$_;
}elsif($signal_seq==1){
$_=~s///g;
$secuencia_total=$secuencia_total.$_;
}
}#Nowcheckthesequence,collectedinitsentirety,
#forthegivenpattern
if($secuencia_total=~/$patron/){
print"Thesequencequery.seqcontainsthepattern$patron\n";
}else{
print"Thesequencequery.seqdoesn'tcontainsthepattern$patron\n";
}#Finallyweclosethefile
#andleavetheprogram
close(query);
exit;
如果想知道数据库里模式的具体位置,我们必须使用特殊变量`$&',这个变量在对正则表达式求值后仍然保存着找到的模式(应该将它放在`if($$secuencia_total>=~/$$patron>/一句的后面)。另外,可以将变量`$`'和`$´'组合起来使用,它们会将找到的模式的左右位置的信息保存。将这些变量正确的加入前面的程序中,我们就可以给出模式的确切位置。注意:length也是非常有用的,它会给出一串数据的长度。
#forthegivenpattern
#andcheckitspositioninthesequence
if($secuencia_total=~/$patron/){
$posicion=length($`) 1;
print"Thesequencequery_seq.txtcontainsthepattern$patroninthe
followingposition$posicion\n";}else{
print"Thesequencequery_seq.txtdoesn'tcontainsthepattern$patron\n";
}
计算氨基酸的频度(Calculusofaminoacidfrequences):
不同蛋白质里,特定的氨基酸出现的频度是不同的,这是因为他们处在不同的环境里面、并且功能不同。下面,我们给出一个例子来展示如何计算给定氨基酸序列里某种氨基酸频度。
#!/usr/bin/perl#Calculatesthefrequencyofaminoacidinaproteinicsequence#Getsthefilenamefromthecommandline#(SWISS-PROTformatted)#Alsocanbeaskedwithprintfromthe<STDIN>if(!$ARGV[0]){print"Theexecutionlineshallbe:program.plfile_swissprot\n";}$fichero=$ARGV[0];#Initializethevariable$erroresmy$errores=0;#Openthefileforreadingopen(FICHA,"$fichero")||die"problemopeningthefile$fichero\n";#Firstwecheckthesequenceasdidintheexample2while(<FICHA>){chomp$_;if($_=~/^SQ/){$signal_good=1;}elsif($signal_good==1){lastif($_=~/^\/\//);$_=~s/\s//g;$secuencia.=$_;}}close(FICHA);#Nowuseacurlthatcheckseverypositionoftheaminoacid#inthesequence(fromafuncionofitsown,thatcanbeusedafterinother#programs)comprueba_aa($secuencia);#Printtheresultstothescreen#Firstthe20aminoacidsandthenthearraywiththeirfrequencies#Inthiscase'sort'can'tbeusedinforeach,#becausethearraycontainsthefrequencies(numbers)print"A\tC\tD\tE\tF\tG\tH\tI\tK\tL\tM\tN\tP\tQ\tR\tS\tT\tV\tW\tY\n";foreach$each_aa(@aa){print"$each_aa\t";}#Tenitgivesthepossibleerrors#andendstheprogramprint"\nerrores=$errores\n";exit;#Functions#Thisonecalculateseachaminoacidfrequency#fromaproteinicsequencesubcomprueba_aa{#Getsthesequencemy($secuencia)=@_;#andrunsaminoacidbyaminoacid,usingaforrunning#from0untilthesequencelengthfor($posicion=0;$posicion<length$secuencia;$posicion ){#Getstheaminoacid$aa=substr($secuencia,$posicion,1);#andcheckswhichoneisusingif#whenitischeckeditaggregates1tothecorrespondantfrequency#inanarrayusingapointerforeachone#orderedinalphabeticwayif($aaeq'A'){$aa[0] ;}elsif($aaeq'C'){$aa[1] ;}elsif($aaeq'D'){$aa[2] ;}elsif($aaeq'E'){$aa[3] ;}elsif($aaeq'F'){$aa[4] ;}elsif($aaeq'G'){$aa[5] ;}elsif($aaeq'H'){$aa[6] ;}elsif($aaeq'I'){$aa[7] ;}elsif($aaeq'K'){$aa[8] ;}elsif($aaeq'L'){$aa[9] ;}elsif($aaeq'M'){$aa[10] ;}elsif($aaeq'N'){$aa[11] ;}elsif($aaeq'P'){$aa[12] ;}elsif($aaeq'Q'){$aa[13] ;}elsif($aaeq'R'){$aa[14] ;}elsif($aaeq'S'){$aa[15] ;}elsif($aaeq'T'){$aa[16] ;}elsif($aaeq'V'){$aa[17] ;}elsif($aaeq'W'){$aa[18] ;}elsif($aaeq'Y'){$aa[19] ;#Iftheaminoacidisnotfound#itaggregates1totheerrors}else{print"ERROR:Aminoacidnotfound:$aa\n";$errores ;}}#Finallyreturnstothefrequencyarrayreturn@aa;}
下面就让我们跟着大自然的步伐,看看细胞中的信息流向了何方。其中之一就是转录,RNA从DNA(基因)中复制出遗传信息,然后又将这些信息传递给蛋白质或者氨基酸序列。为此,我们必须使用与氨基酸对应的基因密码--所谓的RNA/DNA三联密码子。我们要提取Escherichiacoli(一种埃[舍利]希氏杆菌属的大肠杆菌)的基因所对应的氨基酸序列,而这些信息都是以EMBL(EuropeanMolecularBiologyLaboratory)要求的格式。做完这些转换之后,我们将与已有的转录信息校验。对这个例子,非常有必要引进数组的关联变量(associativevariablesofarrays)和哈希表。
#!/usr/bin/perl#TranslatesanADNsequencefromanEMBLfiche#totheaminoacidcorrespondant#Getsthefilenamefromthecommandline#(SWISS-PROTformatted)#Alsocanbeaskedwithprintfromthe<STDIN>if(!$ARGV[0]){print"Theprogramlineshallbe:program.plficha_embl\n";}$fichero=$ARGV[0];#Openthefileforreadingopen(FICHA,"$fichero")||die"problemopeningthefile$fichero\n";#Firstwecheckthesequenceasdidintheexample2while(<FICHA>){chomp$_;if($_=~/^FTCDS/){$_=~tr/..//;($a1,$a2,$a3,$a4)=split("",$_);}elsif($_=~/^SQ/){$signal_good=1;}elsif($signal_good==1){lastif($_=~/^\/\//);#Eliminatenumbersandspaces$_=~tr/0-9//;$_=~s/\s//g;$secuencia.=$_;}}close(FICHA);#Nowwedefineanassociatearraywiththecorrepondence#ofeveryaminoacidswiththeirnucleotide#correspondants(alsoinanownfunction,#forifthesamegeneticcodeisusedinotherprogrammy(codigo_genetico)=('TCA'=>'S',#Serine'TCC'=>'S',#Serine'TCG'=>'S',#Serine'TCT'=>'S',#Serine'TTC'=>'F',#Fenilalanine'TTT'=>'F',#Fenilalanine'TTA'=>'L',#Leucine'TTG'=>'L',#Leucine'TAC'=>'Y',#Tirosine'TAT'=>'Y',#Tirosine'TAA'=>'*',#Stop'TAG'=>'*',#Stop'TGC'=>'C',#Cysteine'TGT'=>'C',#Cysteine'TGA'=>'*',#Stop'TGG'=>'W',#Tryptofane'CTA'=>'L',#Leucine'CTC'=>'L',#Leucine'CTG'=>'L',#Leucine'CTT'=>'L',#Leucine'CCA'=>'P',#Proline'CCC'=>'P',#Proline'CCG'=>'P',#Proline'CCT'=>'P',#Proline'CAC'=>'H',#Hystidine'CAT'=>'H',#Hystidine'CAA'=>'Q',#Glutamine'CAG'=>'Q',#Glutamine'CGA'=>'R',#Arginine'CGC'=>'R',#Arginine'CGG'=>'R',#Arginine'CGT'=>'R',#Arginine'ATA'=>'I',#IsoLeucine'ATC'=>'I',#IsoLeucine'ATT'=>'I',#IsoLeucine'ATG'=>'M',#Methionina'ACA'=>'T',#Treonina'ACC'=>'T',#Treonina'ACG'=>'T',#Treonina'ACT'=>'T',#Treonina'AAC'=>'N',#asparagina'AAT'=>'N',#Asparagina'AAA'=>'K',#Lisina'AAG'=>'K',#Lisina'AGC'=>'S',#Serine'AGT'=>'S',#Serine'AGA'=>'R',#Arginine'AGG'=>'R',#Arginine'GTA'=>'V',#Valine'GTC'=>'V',#Valine'GTG'=>'V',#Valine'GTT'=>'V',#Valine'GCA'=>'A',#Alanine'GCC'=>'A',#Alanine'GCG'=>'A',#Alanine'GCT'=>'A',#Alanine'GAC'=>'D',#AsparticAcid'GAT'=>'D',#AsparticAcid'GAA'=>'E',#GlutamicAcid'GAG'=>'E',#GlutamicAcid'GGA'=>'G',#Glicine'GGC'=>'G',#Glicine'GGG'=>'G',#Glicine'GGT'=>'G',#Glicine);#Translateeverycodoninitscorrespondantaminoacid#andaggregatestotheproteinicsequenceprint$a3;for($i=$a3-1;$i<$a4-3;$i =3){$codon=substr($secuencia,$i,3);#Passthecodonfromsubcase(EMBLformat)touppercase$codon=~tr/a-z/A-Z/;$protein.=codon2aa($codon);}print"Thisproteinicsequenceofthegen:\n$secuencia\nisthefollowing:\n$protein\n\n";exit;
BibliographicReferences
http://bioperl.org/ http://changjiang.whlib.ac.cn/pylorus/download/book/Beginning Perl for Bioinformatics/contents.html http://www.unix.org.ua/orelly/perl/prog3/ - Examplefiles:
-human_kinases_swissprot.txt
-query_seq.txt
-ecoli_embl.txt
- ››linux下两台服务器文件实时同步方案设计和实现
- ››Linux文件描述符中的close on exec标志位
- ››Linux下管道使用的一些限制
- ››Linux 误删/usr/bin 解决方法
- ››linux 添加新用户并赋予sudo执行权限
- ››linux常用软件安装方法
- ››Linux的分区已经被你从Windows中删除,系统启动后...
- ››linux enable命令大全
- ››Linux实现基于Loopback的NVI(NAT Virtual Interfa...
- ››Linux远程访问windows时,出现"连接被对端重...
- ››linux中使用head命令和tail命令查看文件中的指定行...
- ››linux swap 分区调控(swap分区 lvm管理)
更多精彩
赞助商链接