1、摘 要 : 这篇文章介绍了对 DNA、RNA 和蛋白质序列数据库的生物信息提取时,在 Unix 上的 Perl 程序的一些优点。这些 Perl 程序可用来作数据比对处理和分析。人类基因组计划和 DNA 克隆技术的发展加速了这个领域的进步。这些领域每天产生的大量信息使得我们处理这些信息的方式不得不有所改进。 不同基因组(生物有机体上的一组基因)上激增的生物信息推动着生物信息学成为处理分析这些数据的基础手段。 生 物 信 息 学 ( Bioinformatics)生物信息学开始于科学家们将生物学数据以数字格式存放并且用程序来处理这些数据。很长一段时间以来,生物信息学都限制在序列的分析上。然而,随着
2、构建分子的结构模型的重要性开始显现,电子计算机也开始成为理论生物化学的重要工具。每天都不断有关于分子 3D 信息和数据被采集,人们对基因的认识和研究也从单个的基因研究转变为从整体上或者扩展式的研究。由于生物信息学的发展,现在更容易理解蛋白质之间为何相互之间那样作用、又是如何通过新陈代谢来组织相互的。而且我们现在也越来越清醒的认识到组织好这些数据的重要性。生物信息学里面至少有两点特征使她变得非常有趣。其一,生物信息学的研究目标是找出各种生命分子的关系;而这个目标恰恰是一个有趣的程序设计问题,因为这就需要我们联合并整合我们得到的那些信息,然后从中得到对生命活动的一些整体的和有效的一些认识。我们还发
3、现,将计算机科学中的不同领域的知识结合起来是非常必要的,比如数据的管理和整合、高效可靠的算法和强劲的硬件格点技术、多处理器的使用等。 PerlLarry Wall 于 1Array86 年开始开发 Perl。 Perl 是一种解释型的语言,是处理文本、文件和进程的强大的工具。Perl 使得我们能够很快的开发出小程序。可以说,Perl 是高级编程语言(例如 C)和脚本语言(如 bash)的一种有效组合。 Perl 程序可以运行在多种操作系统平台上,尽管 Perl 是在 Unix 上诞生并且快速发展的。由于 Perl 广泛的用于 web 程序设计,其发展很快便超出了其预想。在 Perl 之前,人们
4、使用 awk,thirst 和 grep 来分析文件并提取信息。Perl 将这些 UNIX 上广泛使用的工具统一在一个程序里面,并将这些功能扩展和现代化以适应各种需求。Perl 是一种免费自由的程序语言,可以运行在现代生物实验室里使用的各种操作系统上。在 UNIX 和MacOSX 上,它是预安装好的,在其他系统上,得先安装好 Perl。http:/www.cpan.org 网站上有安装和使用 Perl 的很多实用信息。在 Linux 下,运行 Perl 程序,是将这个程序的文件名作为 perl 这个命令的一个参数,然后 perl 会依次解释执行这个程序里的命令。另一种常用的方法,不需要运行 p
5、erl 这个命令,为此,我们需要做以下两件事: (a)在程序的文件里加入一行特殊的注释: #!/usr/bin/env perlprint “Hin“;(b) 保存此文件并给它加上可执行的属性:% chmod +x greetings.pl这样,我们就可以直接通过文件名来运行这个程序: % ./greetings.pl用 Perl 来 文 件 管 理 :当我们有了文本格式的分子序列,我们可以用 Perl 写一个序列搜索工具。下面的例子我们可以看到如何在 SWISS-PROT(db_human_swissprot)格式的数据库中用 id 码来查找蛋白质序列。#!/usr/bin/perl# Lo
6、ok for aminoacid sequence in a database# SWISS-PROT formated, with a given id code# Ask for the code in the ID field# and it assigns it from the input(STDIN)to a variableprint “Enter the ID to search: “;$id_query=;chomp $id_query;# We open the database file# but if it isnt possible the program endso
7、pen (db, “human_kinases_swissprot.txt“) |die “problem opening the file human_kinases_swissprot.txtn“;# Look line by line in the databasewhile () chomp $_;# Check if we are in the ID fieldif ($_ = /ID/) # If it is possitive we gather the information# breaking the line by spaces($a1,$id_db) = split (/
8、s+/,$_);# but if there is no coincidence of ID we continue to the followingnext if ($id_db ne $id_query);# When they coincide, we put a mark$signal_good=1;# Then we check the sequence field# and if the mark is 1 (chosen sequence)# If possitive, we change the mark to 2,to collect the sequence elsif (
9、$_ = /SQ/) # Finally, if the mark is 2, we present each line# of the sequence, until the line begins with /# is such case we broke the while elsif ($signal_good = 2) last if ($_ = /);print “$_n“;# When we left the while instruction we check the mark# if negative that means that we dont find the chos
10、en sequence# that will give us an errorif (!$signal_good) print “ERROR: “.“Sequence not foundn“;# Finally, we close the file# that still si openclose (db);exit;查 找 氨 基 酸 的 模 式 ( Search for aminoacid patterns)#!/usr/bin/perl# Searcher for aminoacid patterns# Ask the user the patterns for searchprint
11、“Please, introduce the pattern to search in query.seq: “;$patron = ;chomp $patron;# Open the database file# but if it cant it ends the programopen (query, “query_seq.txt“) | die “problem opening the file query_seq.txtn“;# Look line by line the SWISS-PROT sequencewhile () chomp $_;# When arrives to t
12、he SQ field,put the mark in 1if ($_ = /SQ/) $signal_seq = 1;# When arrive to the end of sequence, leave the curl# Check that this expression is put before to check# the mark=1,because this line doesnt belong to the aminoacid sequence elsif ($_ = /) last;# Check the mark if it is equal to 1, if possi
13、tive# eliminate the blank spaces in the sequence line# and join every line in a new variable# To concatenate, we also can do:# $secuencia_total.=$_; elsif ($signal_seq = 1) $_ = s/ /g;$secuencia_total=$secuencia_total.$_;# Now check the sequence, collected in its entirety,# for the given patternif (
14、$secuencia_total = /$patron/) print “The sequence query.seq contains the pattern $patronn“; else print “The sequence query.seq doesnt contains the pattern $patronn“;# Finally we close the file# and leave the programclose (query);exit;如果想知道数据库里模式的具体位置,我们必须使用特殊变量$print “The sequence query_seq.txt cont
15、ains the pattern $patron in the following position $posicionn“; else print “The sequence query_seq.txt doesnt contains the pattern $patronn“;计 算 氨 基 酸 的 频 度 ( Calculus of aminoacid frequences) :不同蛋白质里,特定的氨基酸出现的频度是不同的,这是因为他们处在不同的环境里面、并且功能不同。下面,我们给出一个例子来展示如何计算给定氨基酸序列里某种氨基酸频度。#!/usr/bin/perl# Calculate
16、s the frequency of aminoacid in a proteinic sequence# Gets the file name from the command line# (SWISS-PROT formatted)# Also can be asked with print from the if (!$ARGV0) print “The execution line shall be: program.pl file_swissprotn“;$fichero = $ARGV0;# Initialize the variable $erroresmy $errores=0
17、;# Open the file for readingopen (FICHA, “$fichero“) | die “problem opening the file $ficheron“;# First we check the sequence as did in the example 2while () chomp $_;if ($_ = /SQ/) $signal_good = 1; elsif ($signal_good = 1) last if ($_ = /);$_ = s/s/g;$secuencia.=$_;close (FICHA);# Now use a curl t
18、hat checks every position of the aminoacid# in the sequence (from a funcion of its own,that can be used after in other# programs)comprueba_aa ($secuencia);# Print the results to the screen# First the 20 aminoacids and then the array with their frequencies# In this case sort cant be used in foreach,#
19、 because the array contains the frequencies (numbers)print“AtCtDtEtFtGtHtItKtLtMtNtPtQtRtStTtVtWtYn“;foreach $each_aa (aa) print “$each_aat“;# Ten it gives the possible errors# and ends the programprint “nerrores = $erroresn“;exit;# Functions# This one calculates each aminoacid frequency# from a pro
20、teinic sequencesub comprueba_aa # Gets the sequencemy ($secuencia)=_;# and runs aminoacid by aminoacid, using a for running# from 0 until the sequence lengthfor ($posicion=0 ; $posicion下面就让我们跟着大自然的步伐,看看细胞中的信息流向了何方。其中之一就是转录,RNA 从 DNA(基因)中复制出遗传信息,然后又将这些信息传递给蛋白质或者氨基酸序列。为此,我们必须使用与氨基酸对应的基因密码-所谓的 RNADNA 三
21、联密码子。我们要提取 Escherichia coli(一种埃 舍利 希氏杆菌属的大肠杆菌) 的基因所对应的氨基酸序列,而这些信息都是以 EMBL(European Molecular Biology Laboratory)要求的格式。做完这些转换之后,我们将与已有的转录信息校验。对这个例子,非常有必要引进数组的关联变量(associative variables of arrays)和哈希表。#!/usr/bin/perl# Translates an ADN sequence from an EMBL fiche# to the aminoacid correspondant# Gets
22、the file name from the command line# (SWISS-PROT formatted)# Also can be asked with print from the if (!$ARGV0) print “The program line shall be: program.pl ficha_embln“;$fichero = $ARGV0;# Open the file for readingopen (FICHA, “$fichero“) | die “problem opening the file $ficheron“;# First we check
23、the sequence as did in the example 2while () chomp $_;if ($_ = /FT CDS/) $_ = tr/ /;($a1,$a2,$a3,$a4) = split (“ “,$_);elsif ($_ = /SQ/) $signal_good = 1; elsif ($signal_good = 1) last if ($_ = /);# Eliminate numbers and spaces$_ = tr/0-Array/ /;$_ = s/s/g;$secuencia.=$_;close (FICHA);# Now we defin
24、e an associate array with the correpondence# of every aminoacids with their nucleotide# correspondants (also in an own function,# for if the same genetic code is used in other programmy(%codigo_genetico) = (TCA = S,# SerineTCC = S,# SerineTCG = S,# SerineTCT = S,# SerineTTC = F,# FenilalanineTTT = F
25、,# FenilalanineTTA = L,# LeucineTTG = L,# LeucineTAC = Y,# TirosineTAT = Y,# TirosineTAA = *,# StopTAG = *,# StopTGC = C,# CysteineTGT = C,# CysteineTGA = *,# StopTGG = W,# TryptofaneCTA = L,# LeucineCTC = L,# LeucineCTG = L,# LeucineCTT = L,# LeucineCCA = P,# ProlineCCC = P,# ProlineCCG = P,# Proli
26、neCCT = P,# ProlineCAC = H,# HystidineCAT = H,# HystidineCAA = Q,# GlutamineCAG = Q,# GlutamineCGA = R,# ArginineCGC = R,# ArginineCGG = R,# ArginineCGT = R,# ArginineATA = I,# IsoLeucineATC = I,# IsoLeucineATT = I,# IsoLeucineATG = M,# MethioninaACA = T,# TreoninaACC = T,# TreoninaACG = T,# Treonin
27、aACT = T,# TreoninaAAC = N,# AsparaginaAAT = N,# AsparaginaAAA = K,# LisinaAAG = K,# LisinaAGC = S,# SerineAGT = S,# SerineAGA = R,# ArginineAGG = R,# ArginineGTA = V,# ValineGTC = V,# ValineGTG = V,# ValineGTT = V,# ValineGCA = A,# AlanineGCC = A,# AlanineGCG = A,# AlanineGCT = A,# AlanineGAC = D,#
28、 Aspartic AcidGAT = D,# Aspartic AcidGAA = E,# Glutamic AcidGAG = E,# Glutamic AcidGGA = G,# GlicineGGC = G,# GlicineGGG = G,# GlicineGGT = G,# Glicine);# Translate every codon in its correspondant aminoacid# and aggregates to the proteinic sequenceprint $a3;for($i=$a3 - 1; $i $a4 - 3 ; $i += 3) $codon = substr($secuencia,$i,3);# Pass the codon from subcase (EMBL format) to uppercase$codon = tr/a-z/A-Z/;$protein.= codon2aa($codon);print “This proteinic sequence of the gen:n$secuencianis the following:n$proteinnn“;exit;