【问题标题】:Translating a cDNA to amino acids using Perl使用 Perl 将 cDNA 翻译成氨基酸
【发布时间】:2014-02-27 19:25:28
【问题描述】:

所以我试图将一条互补的 DNA 链翻译成它各自的氨基酸。到目前为止,我有这个代码:

#!/usr/bin/perl

open (INFILE, "sumaira2.out");
open (OUTFILE3, ">>sumaira3.out");

%aacode = (
  TTT => "F", TTC => "F", TTA => "L", TTG => "L",
  TCT => "S", TCC => "S", TCA => "S", TCG => "S",
  TAT => "Y", TAC => "Y", TAA => "STOP", TAG => "STOP",
  TGT => "C", TGC => "C", TGA => "STOP", TGG => "W",
  CTT => "L", CTC => "L", CTA => "L", CTG => "L",
  CCT => "P", CCC => "P", CCA => "P", CCG => "P",
  CAT => "H", CAC => "H", CAA => "Q", CAG => "Q",
  CGT => "R", CGC => "R", CGA => "R", CGG => "R",
  ATT => "I", ATC => "I", ATA => "I", ATG => "M",
  ACT => "T", ACC => "T", ACA => "T", ACG => "T",
  AAT => "N", AAC => "N", AAA => "K", AAG => "K",
  AGT => "S", AGC => "S", AGA => "R", AGG => "R",
  GTT => "V", GTC => "V", GTA => "V", GTG => "V",
  GCT => "A", GCC => "A", GCA => "A", GCG => "A",
  GAT => "D", GAC => "D", GAA => "E", GAG => "E",
  GGT => "G", GGC => "G", GGA => "G", GGG => "G",
); # this is the hash table for the amino acids

while ($line=<INFILE>){
  $codon = $codon.$line;
  @array = split "",$codon;
} # splits all the characters in the text

for ($count = 0; $count<scalar@array; $count= $count + 3) {
  $codon = $codon.$array[$count].$array[$count+1].$array[$count+2];
  $aminoacid = $aacode{$codon};
} # tells how to read the codon and execute the hash table

$protein = $protein.$aminoacid; #catenate the string

print OUTFILE3 $protein;

我的 infile 已经有反向互补 DNA,我只是想翻译它。出于某种原因,我的输出中没有任何内容。我不知道出了什么问题,因为终端也没有给我任何错误。任何帮助将不胜感激。

这是我要翻译的文件示例:

TCGTCGCCTCCCCAACCTAGGTAGTCCGTTGCTGCCCGACGACGGCCGGTAGTCGCCT GCGTCCCTCCTGAAAGGCGTTGGCCGGCAAGCTACGCCGTGGCTACCGGAAGCGCGTCCCCATCAC GCGGTCCTAACTGAACGCGACGGGATGGAGAGTGATCACTCCCCGCCGTCGCGTAGTTCGCCACTC

它会继续运行 17 行。

【问题讨论】:

    标签: arrays perl hashtable dna-sequence


    【解决方案1】:

    好吧,伙计们,

    所以我问了我的教授,我的代码有很多问题。首先,我使用 $codon 两次,同时希望它做两件不同的事情(我在 while 循环中使用了一次,在 for 循环中使用了一次)。所以它将整个 infile 视为一个 $codon,然后在它之后执行哈希表。第二个错误(正如其他人之前提到的)是 $protein 不在 for 循环中,因此只会给我最后一个氨基酸。无论如何,这是更正后的有效代码:

    open (INFILE, "sumaira2.out");
    open (OUTFILE3, ">sumaira3.out");
    
    %aacode = (
    TTT => "F", TTC => "F", TTA => "L", TTG => "L",
    TCT => "S", TCC => "S", TCA => "S", TCG => "S",
    TAT => "Y", TAC => "Y", TAA => "STOP", TAG => "STOP",
    TGT => "C", TGC => "C", TGA => "STOP", TGG => "W",
    CTT => "L", CTC => "L", CTA => "L", CTG => "L",
    CCT => "P", CCC => "P", CCA => "P", CCG => "P",
    CAT => "H", CAC => "H", CAA => "Q", CAG => "Q",
    CGT => "R", CGC => "R", CGA => "R", CGG => "R",
    ATT => "I", ATC => "I", ATA => "I", ATG => "M",
    ACT => "T", ACC => "T", ACA => "T", ACG => "T",
    AAT => "N", AAC => "N", AAA => "K", AAG => "K",
    AGT => "S", AGC => "S", AGA => "R", AGG => "R",
    GTT => "V", GTC => "V", GTA => "V", GTG => "V",
    GCT => "A", GCC => "A", GCA => "A", GCG => "A",
    GAT => "D", GAC => "D", GAA => "E", GAG => "E",
    GGT => "G", GGC => "G", GGA => "G", GGG => "G",
    ); # this is the hash table for the amino acids
    
    while ($line=<INFILE>){
    $line =~ s/\s+$//;
    $sequence = $sequence.$line;
    @array = split "",$sequence;
     } # splits all the characters in the text
    
    for ($count = 0; $count<=scalar @array-3; $count= $count + 3) {
    $codon = $array[$count].$array[$count+1].$array[$count+2];
    $aminoacid = $aacode{$codon};
    $protein = $protein.$aminoacid; #catenate the string
    
    } # tells how to read the codon and execute the hash table
    
    
    print OUTFILE3 $protein;
    

    再次感谢大家的帮助,很抱歉我花了这么长时间才回来!

    【讨论】:

    • 这可能有效,但此代码在技术上存在很多问题。对于初学者,请始终将use strict;use warnings; 放在脚本顶部并使用3-argument version of open。这些东西会为你省去很多麻烦。此外,我强烈建议您使用我在回答中描述的 BioPerl(出于多种原因),除非这是家庭作业。
    • 是的,我确实注意到每个代码都以 use strict 开头;使用警告;。但是,我的教授没有向我介绍这些命令,我​​不确定将它们用于作业是否合适。所以我的问题是,这些命令的目的是什么?就 Bioperl 而言,我认为我无法在本课程中使用它。但希望在我找到一个研究实验室后,我会使用 BioPerl。感谢您的建议!
    • Perl strictwarnings 是通过警告您不安全代码或在出现错误时停止来强制执行更好习惯的编译指示。这些在最新版本的 Perl 中默认启用,这表明确实有必要理解(或至少使用)这些 pragma 并使用词法变量。您的代码可能在没有这些做法的情况下工作(使用某些版本的 Perl),但强烈建议不要这样做。
    【解决方案2】:

    我强烈建议使用 BioPerl 或其他一些库/工具包来解决这类任务。原因是除了有3个阅读框外,还有16个密码子表。在我看来,人们已经在这个问题上花费了太多精力(我也没有看到任何正确的解决方案),并且做任何不平凡的事情都需要更多的工作和代码。这是一个使用标准密码子表进行翻译的简单示例。

    #!/usr/bin/env perl
    
    use strict;
    use warnings;
    use Bio::SeqIO;
    
    my $usage = "$0 nt.fasta";
    my $file  = shift or die $usage;
    my $seqio = Bio::SeqIO->new(-file => $file); 
    
    my $seqobj = $seqio->next_seq;   # create a Bio::Seq object
    my $trans  = $seqobj->translate; # call the translate method 
                                     # on the Bio::Seq object
    
    print $trans->seq;               # $trans is a Bio::Seq object, 
                                     # so we call the seq method to get the sequence
    

    您可以针对多个序列稍作修改,或者使用不同的密码子表。您还可以包含自定义密码子表。 BioPerl HOWTO 页面上有一个很好的翻译序列的教程。

    编辑:我尝试的另外两个解决方案仅适用于序列,但不要像我假设的那样解析 Fasta 格式。一个主要的实际考虑是您应该在翻译中插入一个符号(默认是 BioPerl 的星号,但您可以将其更改为任何您喜欢的)而不是单词“STOP”,因为它不会被任何其他工具识别.肉眼也很难辨别。

    【讨论】:

      【解决方案3】:

      也许以下内容会有所帮助:

      use strict;
      use warnings;
      
      my %aacode = (
        TTT => "F", TTC => "F", TTA => "L", TTG => "L",
        TCT => "S", TCC => "S", TCA => "S", TCG => "S",
        TAT => "Y", TAC => "Y", TAA => "STOP", TAG => "STOP",
        TGT => "C", TGC => "C", TGA => "STOP", TGG => "W",
        CTT => "L", CTC => "L", CTA => "L", CTG => "L",
        CCT => "P", CCC => "P", CCA => "P", CCG => "P",
        CAT => "H", CAC => "H", CAA => "Q", CAG => "Q",
        CGT => "R", CGC => "R", CGA => "R", CGG => "R",
        ATT => "I", ATC => "I", ATA => "I", ATG => "M",
        ACT => "T", ACC => "T", ACA => "T", ACG => "T",
        AAT => "N", AAC => "N", AAA => "K", AAG => "K",
        AGT => "S", AGC => "S", AGA => "R", AGG => "R",
        GTT => "V", GTC => "V", GTA => "V", GTG => "V",
        GCT => "A", GCC => "A", GCA => "A", GCG => "A",
        GAT => "D", GAC => "D", GAA => "E", GAG => "E",
        GGT => "G", GGC => "G", GGA => "G", GGG => "G",
      ); # this is the hash table for the amino acids
      
      my $compDNA = uc do { local $/; <> };
      $compDNA =~ s/\s+//g;
      
      my @codons = unpack '(A3)*', $compDNA;
      my @aminoAcids = map { exists $aacode{$_} ? $aacode{$_} : "?$_?" } @codons;
      print join '', @aminoAcids;
      

      用法:perl script.pl compDNA_File [&gt;aminoAcid_File]

      最后一个可选参数将输出定向到文件。

      首先,整个文件被 slurped(并转换为全部大写)到一个变量中。接下来,删除所有空格。 unpack 用于创建三字符元素(密码子)的列表。 map 用于使用您提供的哈希将密码子翻译成氨基酸。 (请注意,如果没有密码子的键,则插入密码子,并用问号括起来。)最后,将这些氨基酸joined 形成一个字符串,结果为printed。

      【讨论】:

      • 很好地使用解包和地图
      • @MatthewLock - 感谢您的评论。谢谢。
      • 大家好,非常感谢您的帮助。我尝试了每个人的更正,但是,我仍然无法使其正常工作。我想我要在 "$aminoacid = $aacode{$codon}" 中创建一个 while 循环。我只是不确定如何执行这个哈希,我会用 map 和 join 命令尝试不同的变体!再次感谢您的帮助!
      • @user3268152 - 它可以帮助您发布您需要翻译的数据样本;将其添加到您的原始问题中。
      • 我添加了几行需要翻译的数据并添加到我的原始问题中。
      【解决方案4】:

      尝试以scriptname &lt; sumaira2.out &gt;&gt; sumaira3.out 执行以下脚本。
      $DEBUG 设置为零以删除调试输出(如果它按预期工作)。

      #!/usr/bin/perl
      use strict; use warnings;
      
      my $DEBUG = 2;
      
      my %aacode = (
      TTT => "F", TTC => "F", TTA => "L", TTG => "L",
      TCT => "S", TCC => "S", TCA => "S", TCG => "S",
      TAT => "Y", TAC => "Y", TAA => "STOP", TAG => "STOP",
      TGT => "C", TGC => "C", TGA => "STOP", TGG => "W",
      CTT => "L", CTC => "L", CTA => "L", CTG => "L",
      CCT => "P", CCC => "P", CCA => "P", CCG => "P",
      CAT => "H", CAC => "H", CAA => "Q", CAG => "Q",
      CGT => "R", CGC => "R", CGA => "R", CGG => "R",
      ATT => "I", ATC => "I", ATA => "I", ATG => "M",
      ACT => "T", ACC => "T", ACA => "T", ACG => "T",
      AAT => "N", AAC => "N", AAA => "K", AAG => "K",
      AGT => "S", AGC => "S", AGA => "R", AGG => "R",
      GTT => "V", GTC => "V", GTA => "V", GTG => "V",
      GCT => "A", GCC => "A", GCA => "A", GCG => "A",
      GAT => "D", GAC => "D", GAA => "E", GAG => "E",
      GGT => "G", GGC => "G", GGA => "G", GGG => "G",
      ); # this is the hash table for the amino acids
      
      my ($codon, $protein)  = ('','');
      while (<STDIN>){
        chomp; # remove end of line characters
        s/\s//g; # remove whitespaces
        $codon .= $_;
      }
      
      print STDERR "DBG Codon: ", $codon, "\n" if $DEBUG >= 1;
      
      my @aminoacids = ( $codon =~ /(...)/sg );
      
      print STDERR "Aminoacids: ", join(" ", @aminoacids), "\n" if $DEBUG >= 2;
      
      for my $aminoacid (@aminoacids) {
        die "Unknown aminoacid: $aminoacid\n" unless exists $aacode{$aminoacid};
        $protein .=  $aacode{$aminoacid};
      }
      
      print STDERR "DBG Protein: ", $protein, "\n" if $DEBUG >= 1;
      
      print $protein, "\n";
      

      【讨论】:

      • 感谢您的帮助!但是,它仍然没有正确执行代码。我想我必须在 "$aminoacid = $aacode{$codon}" 中添加一个 while 循环并使用 map 和 join 命令。再次感谢您的帮助!
      • 您能否显示 SHORT 序列(3-5 个氨基酸)和调试输出产生了什么问题? $codon 字符串被单条指令分割成 @aminoacids 数组 my @aminoacids = ( $codon =~ /(...)/sg );
      【解决方案5】:

      你不想放

      print OUTFILE3 $protein;
      

      在你的 for 循环中,这样你就可以打印出你正在处理的每一个蛋白质,而不是你在 for 循环完成运行后剩下的最后一个蛋白质,就像这样?

      for ($count = 0; $count<scalar@array; $count= $count + 3) {
        $codon = $codon.$array[$count].$array[$count+1].$array[$count+2];
        $aminoacid = $aacode{$codon};
      
        print OUTFILE3 $aminoacid;
      
      } # tells how to read the codon and execute the hash table
      

      【讨论】:

        猜你喜欢
        • 2017-09-29
        • 2016-07-07
        • 1970-01-01
        • 1970-01-01
        • 2014-03-28
        • 2021-12-07
        • 1970-01-01
        • 2010-11-07
        • 2014-04-08
        相关资源
        最近更新 更多