摘要:
生物信息學(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。http://www.cpan.org網站上有安裝和使用Perl的很多實用信息。
在linux下,運行Perl程序,是將這個程序的文件名作為perl這個命令的一個參數,然后perl會依次解釋執行這個程序里的命令。
另一種常用的方法,不需要運行perl這個命令,為此,我們需要做以下兩件事:(a)在程序的文件里加入一行特殊的注釋:
#!/usr/bin/envperl
(b)保存此文件并給它加上可執行的屬性:
chmod xgreetings.pl這樣,我們就可以直接通過文件名來運行這個程序:
./greetings.pl
用Perl來文件管理:
當我們有了文本格式的分子序列,我們可以用Perl寫一個序列搜索工具。下面的例子我們可以看到如何在SWISS-PROT(db_human_swissprot)格式的數據庫中用id碼來查找蛋白質序列。
#!/usr/bin/perl
#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)
#!/usr/bin/perl#Searcherforaminoacidpatterns#Asktheuserthepatternsforsearchprint"Please,introducethepatterntosearchinquery.seq:";$patron=<STDIN>;chomp$patron;#Openthedatabasefile#butifitcan'titendstheprogramopen(query,"query_seq.txt")||die"problemopeningthefilequery_seq.txt/n";#LooklinebylinetheSWISS-PROTsequencewhile(<query>){chomp$_;#WhenarrivestotheSQfield,putthemarkin1
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也是非常有用的,它會給出一串數據的長度。
#Onlyweneedtochangetheifwherethepatternwasfound#Nowcheckthesequence,collectedinitsentirety,
#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://changjiang.whlib.ac.cn/pylorus/download/book/Beginning Perl for Bioinformatics/contents.html
ecoli_embl.txt