3

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

编写一个 Perl 程序,它将从命令行获取两个文件,一个酶文件和一个 DNA 文件。使用限制酶读取文件并将限制酶应用于 DNA 文件。

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

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

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

Main:不确定如何将我的潜艇绑定在一起?

子程序 $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

输出应该是:

中交CCGGGG

TTTCCCCCCCCCCCCAAA

TTTCCCCCCCCCCCCAGA

TTTCCCCCCCCCCGAG

TTTCCCCC

4

4 回答 4

3

好的,我知道我不应该只做你的作业,但是这个有一些有趣的技巧,所以我玩了一下。从中学习,而不仅仅是复制。我没有很好地评论,所以如果你有什么不明白的地方,请提问。这有一些小魔法,如果你没有在课堂上讲过,你的教授会知道的,所以一定要明白。

#!/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"
}
于 2010-11-28T17:36:50.783 回答
3

这是我尝试解决问题的方法(下面的代码)。
1) 从参数中提取文件名并filehandles创建相应的文件名。
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;
于 2010-11-28T06:50:01.577 回答
3

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

在 Main 例程中,您可以遍历哈希;为每种酶产生一个输出文件。最直接的方法是将站点翻译成正则表达式(通过其他正则表达式)并将其应用于 DNA 序列,但还有其他方法。(可能值得将其拆分为至少一个其他子。)

于 2010-11-28T06:24:00.593 回答
2

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

#!/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 示例测试了我的代码,效果很好。

由于这是一项可以通过编写几行非重复代码来处理的任务,我认为创建单独的子例程是不合理的。我希望我会被原谅... :)

于 2010-11-28T20:29:27.030 回答