Perl 程序使用引用、哈希表和子集来模拟限制性内切酶

Question

我是 Perl 入门课程的学生。我正在寻找有关如何完成任务的建议。我的教授鼓励论坛。任务是：

编写一个 Perl 程序，该程序将从命令行获取两个文件：一个酶文件和一个 DNA 文件。使用限制性内切酶读取文件并将限制性内切酶应用于 DNA 文件。

输出将是按照 DNA 文件中出现的顺序排列的 DNA 片段。输出文件的名称应通过将限制性酶的名称附加到 DNA 文件的名称来构建，并在它们之间使用下划线。

例如，如果酶是 EcoRI 并且 DNA 文件名为 BC161026，则输出文件应命名为 BC161026_EcoRI。

我的方法是创建一个主程序和两个子程序，如下所示：

主要： 不确定如何将我的潜艇绑在一起？

子程序$DNA： 获取 DNA 文件并删除所有新行以生成单个字符串

子程序酶： 读取并存储来自命令行的酶文件中的行 以将酶首字母缩略词与剪切位置分开的方式解析文件。 将剪切位置作为正则表达式存储在哈希表中 将首字母缩略词的名称存储在哈希表中

酶文件格式注意事项： 酶文件遵循称为 Staden 的格式。示例：

AatI/AGG'CCT//
AatII/GACGT'C//
AbsI/CC'TCGAGG//

酶首字母缩略词由第一个斜杠之前的字符组成（AatI，在第一个示例中。识别序列是第一个和第二个斜杠之间的所有内容（AGG'CCT，同样，在第一个示例中）。切点是在识别序列中用撇号表示 酶内 DNA 的标准缩写如下：

R = G 或 A B = 非 A（C 或 G 或 T） 等等...

除了对主要块的建议之外，您是否看到我遗漏的任何缺失部分？您能否推荐您认为对修补该程序有用的特定工具？

输入酶示例：TryII/RRR'TTT//

要读取的示例字符串：CCCCCCGGGTTTCCCCCCCCCCCCAAATTTCCCCCCCCCCCCAGATTTCCCCCCCCCCGAGTTTCCCCC

输出应为：

CCCCCCGGGG

TTTCCCCCCCCCCCCAAA

TTTCCCCCCCCCCCCAGA

TTTCCCCCCCCCCGAG

TTCCCCCC

Answer 1

好吧，我知道我不应该只做你的作业，但是这个有一些有趣的技巧，所以我玩了一下。从中学习，而不是仅仅复制。我评论得不太好，有不明白的地方请追问。这其中有一些微妙的魔力，如果你在课堂上没有讲过，你的教授会知道，所以一定要理解。

#!/usr/bin/env perl

use strict;
use warnings;

use Getopt::Long;

my ($enzyme_file, $dna_file);
my $write_output = 0;
my $verbose = 0;
my $help = 0;
GetOptions(
  'enzyme=s' => \$enzyme_file,
  'dna=s' => \$dna_file,
  'output' => \$write_output,
  'verbose' => \$verbose,
  'help' => \$help
);

$help = 1 unless ($dna_file && $enzyme_file);
help() if $help; # exits

# 'Main'
my $dna = getDNA($dna_file);
my %enzymes = %{ getEnzymes($enzyme_file) }; # A function cannot return a hash, so return a hashref and then store the referenced hash
foreach my $enzyme (keys %enzymes) {
  print "Applying enzyme " . $enzyme . " gives:\n";
  my $dna_holder = $dna;
  my ($precut, $postcut) = ($enzymes{$enzyme}{'precut'}, $enzymes{$enzyme}{'postcut'});

  my $R = qr/[GA]/;
  my $B = qr/[CGT]/;

  $precut =~ s/R/${R}/g;
  $precut =~ s/B/${B}/g;
  $postcut =~ s/R/${R}/g;
  $postcut =~ s/B/${B}/g;
  print "\tPre-Cut pattern: " . $precut . "\n" if $verbose;
  print "\tPost-Cut pattern: " . $postcut . "\n" if $verbose;

  #while(1){
  #  if ($dna_holder =~ s/(.*${precut})(${postcut}.*)/$1/ ) {
  #    print "\tFound section:" . $2 . "\n" if $verbose;
  #    print "\tRemaining DNA: " . $1 . "\n" if $verbose;
  #    unshift @{ $enzymes{$enzyme}{'cut_dna'} }, $2;
  #  } else {
  #    unshift @{ $enzymes{$enzyme}{'cut_dna'} }, $dna_holder;
  #    print "\tNo more cuts.\n" if $verbose;
  #    print "\t" . $_ . "\n" for @{ $enzymes{$enzyme}{'cut_dna'} };
  #    last;
  #  }
  #}
  unless ($dna_holder =~ s/(${precut})(${postcut})/$1'$2/g) {
    print "\tHas no effect on given strand\n" if $verbose;
  }
  @{ $enzymes{$enzyme}{'cut_dna'} } = split(/'/, $dna_holder);
  print "\t$_\n" for @{ $enzymes{$enzyme}{'cut_dna'} };

  writeOutput($dna_file, $enzyme, $enzymes{$enzyme}{'cut_dna'}) if $write_output; #Note that $enzymes{$enzyme}{'cut_dna'} is an arrayref already
  print "\n";
}

sub getDNA {
  my ($dna_file) = @_;

  open(my $dna_handle, '<', $dna_file) or die "Cannot open file $dna_file";
  my @dna_array = <$dna_handle>;
  chomp(@dna_array);

  my $dna = join('', @dna_array);

  print "Using DNA:\n" . $dna . "\n\n" if $verbose;
  return $dna;
}

sub getEnzymes {
  my ($enzyme_file) = @_;
  my %enzymes;

  open(my $enzyme_handle, '<', $enzyme_file) or die "Cannot open file $enzyme_file";
  while(<$enzyme_handle>) {
    chomp;
    if(m{([^/]*)/([^']*)'([^/]*)//}) {
      print "Found Enzyme " . $1 . ":\n\tPre-cut: " . $2 . "\n\tPost-cut: " . $3 . "\n" if $verbose;
      $enzymes{$1} = {
        precut => $2,
        postcut => $3,
        cut_dna => [] #Added to show the empty array that will hold the cut DNA sections
      };
    }
  }

  print "\n" if $verbose;
  return \%enzymes;
}

sub writeOutput {

  my ($dna_file, $enzyme, $cut_dna_ref) = @_;

  my $outfile = $dna_file . '_' . $enzyme;
  print "\tSaving data to $outfile\n" if $verbose; 
  open(my $outfile_handle, '>', $outfile) or die "Cannot open $outfile for writing";

  print $outfile_handle $_ . "\n" for @{ $cut_dna_ref };
}

sub help {

  my $filename = (split('/', $0))[-1];

  my $enzyme_text = <<'END';
AatI/AGG'CCT//
AatII/GACGT'C//
AbsI/CC'TCGAGG//
TryII/RRR'TTT//
Test/AAA'TTT//
END

  my $dna_text = <<'END';
CCCCCCGGGTTTCCCCCCC
CCCCCAAATTTCCCCCCCCCCCCAGATTTC
CCCCCCCCCGAGTTTCCCCC
END

  print <<END;
Usage: 
    $filename --enzyme (-e) <enzyme-filename> --dna (-d) <dna-filename> [options] (files may come in either order)
    $filename -h    (shows this help)

Options: 
    --verbose (-v)  print additional (debugging) information
    --output (-o)   output final data to files


Files:
The DNA file contains one DNA string which may be broken over many lines. E.G.:

$dna_text

The enzymes file constains enzyme definitions, one per line. E.G.:

$enzyme_text
END

exit;
}

编辑：明确添加了 cut_dna 初始化，因为这是每种酶的最终结果持有者，所以我认为更清楚地看到它会很好。

编辑2：添加了输出例程、调用、标志和相应的帮助。

编辑 3：更改了主例程，以结合 Canavanin 的最佳方法，同时删除循环。现在它是一个全局替换，添加临时剪切标记（'），然后将剪切标记拆分到数组中。留下旧方法作为注释，新方法是下面的 5 行。

编辑 4：用于写入多个文件的附加测试用例。 （以下）

my @names = ('cat','dog','sheep'); 
foreach my $name (@names) { #$name is lexical, ie dies after each loop
  open(my $handle, '>', $name); #open a lexical handle for the file, also dies each loop
  print $handle $name; #write to the handle
  #$handles closes automatically when it "goes out of scope"
}

Answer 2

请注意，在酶中，当您将酶存储在哈希中时，酶的名称应该是键，位点应该是值（因为原则上两种酶可以具有相同的位点）。

在 Main 例程中，您可以迭代哈希；对于每种酶产生一个输出文件。最直接的方法是将位点翻译为正则表达式（通过其他正则表达式）并将其应用于 DNA 序列，但还有其他方法。 （可能值得将其分成至少一个其他子部分。）

Answer 3

这是我尝试解决问题的方法（代码如下）。

1) 从参数中获取文件名并创建相应的文件句柄。
2) 为指定格式的输出文件创建一个新的文件句柄
3）从第一个文件中提取“切点”
4) 第二个文件中的 DNA 序列在步骤3中获得的切割点上循环。

#!/usr/bin/perl
use strict;
use warnings;
my $file_enzyme=$ARGV[0];
my $file_dna=$ARGV[1];

open DNASEQ, ">$file_dna"."_"."$file_enzyme";
open ENZYME, "$file_enzyme";
open DNA, "$file_dna";
while (<ENZYME>) {
 chomp;
  if( /'(.*)\/\//) { # Extracts the cut point between ' & // in the enzyme file
    my $pattern=$1;
    while (<DNA>) {
     chomp;
     #print $pattern;
     my @output=split/$pattern/,;
     print DNASEQ shift @output,"\n"; #first recognized sequence being output
     foreach (@output) {
        print DNASEQ "$pattern$_\n"; #prefixing the remaining sequences with the cut point pattern
     }
   }
 }
}
close DNA;
close ENZYME;
close DNASEQ;

Answer 4

我知道已经有几个答案了，但是嘿...我只是想试试运气，所以这是我的建议：

#!/usr/bin/perl

use warnings;
use strict;
use Getopt::Long;

my ($enz_file, $dna_file);

GetOptions( "e=s" => \$enz_file,
            "d=s" => \$dna_file,
          );

if (! $enz_file || ! $dna_file) {
   # some help text 
   print STDERR<<EOF; 

   Usage: restriction.pl -e enzyme_file -d DNA_file

   The enzyme_file should contain one enzyme entry per line.
   The DNA_file may contain the sequence on one single or on
   several lines; all lines will be concatenated to yield a
   single string.
EOF      
   exit();
}

my %enz_and_patterns; # stores enzyme name and corresponding pattern

open ENZ, "<$enz_file" or die "Could not open file $enz_file: $!";
while (<ENZ>) {
   if (m#^(\w+)/([\w']+)//$#) {
      my $enzyme  = $1; # could also remove those two lines and use 
      my $pattern = $2; # the match variables directly, but this is clearer

      $enz_and_patterns{$enzyme} = $pattern;
   }
}
close ENZ;

my $dna_sequence;

open DNA, "<$dna_file" or die "Could not open file $dna_file: $!";
while (my $line = <DNA>) {
   chomp $line;
   $dna_sequence .= $line; # append the current bit to the sequence
                           # that has been read so far
}
close DNA;

foreach my $enzyme (keys %enz_and_patterns) {
   my $dna_seq_processed = $dna_sequence; # local copy so that we retain the original

   # now translate the restriction pattern to a regular expression pattern:
   my $pattern = $enz_and_patterns{$enzyme};
   $pattern    =~ s/R/[GA]/g; # use character classes
   $pattern    =~ s/B/[^A]/g;
   $pattern    =~ s/(.+)'(.+)/($1)($2)/; # remove the ', but due to the grouping
                                         # we "remember" its position

   $dna_seq_processed =~ s/$pattern/$1\n$2/g; # in effect we are simply replacing
                                              # each ' with a newline character
   my $outfile = "${dna_file}_$enzyme";
   open OUT, ">$outfile" or die "Could not open file $outfile: $!";
   print OUT $dna_seq_processed , "\n";
   close OUT;
}

我已经用你的 TryII 示例测试了我的代码，效果很好。

由于这是一项只需编写几行非重复代码即可处理的任务，因此我认为创建单独的子例程是不合理的。我希望我能被原谅...:)