perl - 如何解析文件、创建记录并对记录执行操作,包括术语频率和距离计算

标签 perl hash bioinformatics

我是 Perl 入门类(class)的一名学生,正在寻找有关我编写一个分析原子数据的小型(但棘手)程序的方法的建议和反馈。我的教授鼓励论坛。我对 Perl 子程序或模块(包括 Bioperl)并不熟悉,因此请将回复限制在适当的“初学者水平”,以便我可以理解并从您的建议和/或代码中学习(也请限制“Magic”)。

该计划的要求如下:

  1. Read a file (containing data about Atoms) from the command line & create an array of atom records (one record/atom per newline). For each record the program will need to store:

    • The atom's serial number (cols 7 - 11)
    • The three-letter name of the amino acid to which it belongs (cols 18 - 20)
    • The atom's three coordinates (x,y,z) (cols 31 - 54 )
    • The atom's one- or two-letter element name (e.g. C, O, N, Na) (cols 77-78 )

  2. 提示输入三个命令之一:频率、长度、密度 d(d 是某个数字):

    • freq - how many of each type of atom is in the file (example Nitrogen, Sodium, etc would be displayed like this: N: 918 S: 23
    • length - The distances among coordinates
    • density d (where d is a number) - program will prompt for the name of a file to save computations to and will containing the distance between that atom and every other atom. If that distance is less than or equal to the number d, it increments the count of the number of atoms that are within that distance, unless that count is zero into the file. The output will look something like:
    1: 5
    2: 3
    3: 6
    ... (very big file) and will close when it finishes.

我正在寻找有关我在下面的代码中编写(和需要编写)的内容的反馈。我特别感谢任何有关如何编写我的字幕的反馈。我在底部提供了示例输入数据。

我看到的程序结构和功能描述:

$^W = 1; # turn on warnings
use strict; # behave!

my @fields;
my @recs;

while ( <DATA> ) {
 chomp;
 @fields = split(/\s+/);
 push @recs, makeRecord(@fields);
}

for (my $i = 0; $i < @recs; $i++) {
 printRec( $recs[$i] );
}
    my %command_table = (
 freq => \&freq,
 length => \&length,
 density => \&density,
 help => \&help, 
 quit => \&quit
 );

print "Enter a command: ";
while ( <STDIN> ) {
 chomp; 
 my @line = split( /\s+/);
 my $command = shift @line;
 if ($command !~ /^freq$|^density$|length|^help$|^quit$/ ) {
    print "Command must be: freq, length, density or quit\n";
    }
  else {
    $command_table{$command}->();
    }
 print "Enter a command: ";
 }

sub makeRecord 
    # Read the entire line and make records from the lines that contain the 
    # word ATOM or HETATM in the first column. Not sure how to do this:
{
 my %record = 
 (
 serialnumber => shift,
 aminoacid => shift,
 coordinates => shift,
 element  => [ @_ ]
 );
 return\%record;
}

sub freq
    # take an array of atom records, return a hash whose keys are 
    # distinct atom names and whose values are the frequences of
    # these atoms in the array.  

sub length
    # take an array of atom records and return the max distance 
    # between all pairs of atoms in that array. My instructor
    # advised this would be constructed as a for loop inside a for loop. 

sub density
    # take an array of atom records and a number d and will return a
    # hash whose keys are atom serial numbers and whose values are 
    # the number of atoms within that distance from the atom with that
    # serial number. 

sub help
{
    print "To use this program, type either\n",
          "freq\n",
          "length\n",
          "density followed by a number, d,\n",
          "help\n",
          "quit\n";
}

sub quit
{
 exit 0;
}

# truncating for testing purposes. Actual data is aprox. 100 columns 
# and starts with ATOM or HETATM.
__DATA__
ATOM   4743  CG  GLN A 704      19.896  32.017  54.717  1.00 66.44           C  
ATOM   4744  CD  GLN A 704      19.589  30.757  55.525  1.00 73.28           C  
ATOM   4745  OE1 GLN A 704      18.801  29.892  55.098  1.00 75.91           O 

最佳答案

看起来您的 Perl 技能正在进步——使用引用和复杂的数据结构。以下是一些提示和一般建议。

  • 使用use warnings启用警告,而不是$^W = 1。前者是自记录的,其优点是位于封闭 block 的本地而不是全局设置。

  • 使用命名良好的变量,这将有助于记录程序的行为,而不是依赖 Perl 的特殊 $_。例如:

    while (my $input_record = <DATA>){
    }
    
  • 在用户输入场景中,无限循环提供了一种避免重复指令(例如“输入命令”)的方法。见下文。

  • 您的正则表达式可以简化以避免重复 anchor 的需要。见下文。

  • 一般来说,肯定测试比否定测试更容易理解。请参阅下面修改后的 if-else 结构。

  • 将程序的每个部分包含在其自己的子例程中。出于多种原因,这是一个很好的一般做法,所以我就开始养成这个习惯。

  • 一个相关的良好实践是尽量减少全局变量的使用。作为练习,您可以尝试编写程序,使其完全不使用全局变量。相反,任何需要的信息都将在子例程之间传递。对于小程序,不一定需要严格避免全局变量,但牢记理想并不是一个坏主意。

  • 为您的length 子例程指定一个不同的名称。该名称已被内置 length 函数使用。

  • 关于您关于makeRecord的问题,一种方法是忽略makeRecord内部的过滤问题。相反,makeRecord 可以包含一个额外的哈希字段,并且过滤逻辑将驻留在其他地方。例如:

    my $record = makeRecord(@fields);
    push @recs, $record if $record->{type} =~ /^(ATOM|HETATM)$/;
    

上述一些要点的说明:

use strict;
use warnings;

run();

sub run {
    my $atom_data = load_atom_data();
    print_records($atom_data);
    interact_with_user($atom_data);
}

...

sub interact_with_user {
    my $atom_data = shift;
    my %command_table = (...);

    while (1){
        print "Enter a command: ";
        chomp(my $reply = <STDIN>);

        my ($command, @line) = split /\s+/, $reply;

        if ( $command =~ /^(freq|density|length|help|quit)$/ ) {
            # Run the command.
        }
        else {
            # Print usage message for user.
        }
    }
}

...

关于perl - 如何解析文件、创建记录并对记录执行操作,包括术语频率和距离计算,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/4351080/

相关文章:

perl - Perl Tk 模块的缺点是什么?

javascript - 是否有一种类似 JS 的语言允许对本地文件进行写 Access ?

php - password_verify 正在返回 true

ruby - 如何使用ruby更改小写的哈希键

Python 和 Matplotlib : characters as the x axis

gzip - 有没有办法将 vcf.gz 文件直接转换为 vcf.bgz 文件?

regex - 对 Perl 中正则表达式基本规则的困惑

perl - 为什么我用 Perl 编写的图像下载 CGI 脚本不起作用?

将哈希值生成为数字的 Python 库

r - 按相似行折叠数据框