LINUX & PERL, 学习和分析生物学信息的电脑工具

# Only we need to change the if where the pattern was found
# Now check the sequence, collected in its entirety,

# for the given pattern

# and check its position in the sequence

  if ($secuencia_total =~ /$patron/) {

   $posicion=length($`)+1;

   print "The sequence query_seq.txt contains the pattern $patron in the following position $posicion\n";
  } else {

  print "The sequence query_seq.txt doesn't contains the pattern $patron\n";

  }

计算氨基酸的频度（Calculus of aminoacid frequences）:

不同蛋白质里，特定的氨基酸出现的频度是不同的，这是因为他们处在不同的环境里面、并且功能不同。下面，我们给出一个例子来展示如何计算给定氨基酸序列里某种氨基酸频度。

#!/usr/bin/perl
# Calculates 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 <STDIN>
if (!$ARGV[0]) {print "The execution line shall be: program.pl file_swissprot\n";}
$fichero = $ARGV[0];
# Initialize the variable $errores
my $errores=0;
# Open the file for reading
open (FICHA, "$fichero") || die "problem opening the file $fichero\n";
# First we check the sequence as did in the example 2
while (<FICHA>) {
chomp $_;
if ($_ =~ /^SQ/) {
$signal_good = 1;
} elsif ($signal_good == 1) {
   last if ($_ =~ /^\/\//);
   $_ =~ s/\s//g;
   $secuencia.=$_;
}
}
close (FICHA);
# Now use a curl that 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' can't be used in foreach, 
# because the array contains the frequencies (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";
}
# Ten it gives the possible errors
# and ends the program
print "\nerrores = $errores\n";
exit;
# Functions
# This one calculates each aminoacid frequency
# from a proteinic sequence
sub comprueba_aa {
# Gets the sequence
my ($secuencia)=@_;
# and runs aminoacid by aminoacid, using a for running
# from 0 until the sequence length
for ($posicion=0 ; $posicion<length $secuencia ; $posicion++ ) {
# Gets the aminoacid
$aa = substr($secuencia, $posicion, 1);
# and checks which one is using if
# when it is checked it aggregates 1 to the correspondant frequency
# in an array using a pointer for each one
# ordered in alphabetic way
if ( $aa eq 'A' ) {
$aa[0]++;
} elsif ( $aa eq 'C' ) {
$aa[1]++;
} elsif ( $aa eq 'D' ) {
$aa[2]++;
} elsif ( $aa eq 'E' ) {
$aa[3]++;
} elsif ( $aa eq 'F' ) {
$aa[4]++;
} elsif ( $aa eq 'G' ) {
$aa[5]++;
} elsif ( $aa eq 'H' ) {
$aa[6]++;
} elsif ( $aa eq 'I' ) {
$aa[7]++;
} elsif ( $aa eq 'K' ) {
$aa[8]++;
} elsif ( $aa eq 'L' ) {
$aa[9]++;
} elsif ( $aa eq 'M' ) {
$aa[10]++;
} elsif ( $aa eq 'N' ) {
$aa[11]++;
} elsif ( $aa eq 'P' ) {
$aa[12]++;
} elsif ( $aa eq 'Q' ) {
$aa[13]++;
} elsif ( $aa eq 'R' ) {
$aa[14]++;
} elsif ( $aa eq 'S' ) {
$aa[15]++;
} elsif ( $aa eq 'T' ) {
$aa[16]++;
} elsif ( $aa eq 'V' ) {
$aa[17]++;
} elsif ( $aa eq 'W' ) {
$aa[18]++;
} elsif ( $aa eq 'Y' ) {
$aa[19]++;
# If the aminoacid is not found
# it aggregates 1 to the errors
} else {
print "ERROR: Aminoacid not found: $aa\n";
$errores++;
}
}
# Finally returns to the frequency array
return @aa;
}

下面就让我们跟着大自然的步伐，看看细胞中的信息流向了何方。其中之一就是转录，RNA 从DNA（基因）中复制出遗传信息，然后又将这些信息传递给蛋白质或者氨基酸序列。为此，我们必须使用与氨基酸对应的基因密码--所谓的RNA／DNA三联密码子。我们要提取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 the file name from the command line
# (SWISS-PROT formatted)
# Also can be asked with print from the <STDIN>
if (!$ARGV[0]) {print "The program line shall be: program.pl ficha_embl\n";}
$fichero = $ARGV[0];
# Open the file for reading
open (FICHA, "$fichero") || die "problem opening the file $fichero\n";
# First we check the sequence as did in the example 2
while (<FICHA>) {
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-9/ /;
$_ =~ s/\s//g;
$secuencia.=$_;
}
}
close (FICHA);
# Now we define 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 program
my(%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',# Aspartic Acid
'GAT' => 'D',# Aspartic Acid
'GAA' => 'E',# Glutamic Acid
'GAG' => 'E',# Glutamic Acid
'GGA' => 'G',# Glicine
'GGC' => 'G',# Glicine
'GGG' => 'G',# Glicine
'GGT' => 'G',# Glicine
);
# Translate every codon in its correspondant aminoacid
# and aggregates to the proteinic sequence
print $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$secuencia\nis the following:\n$protein\n\n";
exit;

Bibliographic References

• http://bioperl.org/
  
  
• http://changjiang.whlib.ac.cn/pylorus/download/book/Beginning%20Perl%20for%20Bioinformatics/contents.html
  
  

• http://www.unix.org.ua/orelly/perl/prog3/
  
  

• Example files :
  - human_kinases_swissprot.txt
  - query_seq.txt
  - ecoli_embl.txt