2013-05-03 3 views
1

Это программа для извлечения соседних остатков в пределах 4.5 ангстрема. Я разрешил программу до чисел атомов. Из них я хочу, чтобы извлечьПрограмма Perl для вычисления остатков neibouring из файла pdb

  • остатков, номера,
  • имя остатка,
  • атом № и
  • имя атома.

Я хотел бы выводить эти данные в табличной форме, чтобы я мог напрямую копировать результаты. Но теперь я застрял и нуждаюсь в помощи, чтобы извлечь эти поля для чисел атомов, которые я получил в $close_atomno, а также как использовать программу для нескольких файлов pdb и разных каталитических остатков за один раз.

Любая помощь приветствуется.

#!/usr/bin/perl 
use List::Util qw(sum); 
use Data::Dumper qw(Dumper); 
use 5.010; 

say "enter residue no."; 
chomp(my $r_no = <STDIN>); 

my (@all_pos, @all_data, @matching_pos, @matching_data); 

my $residue_file = "neighbouring_residues.txt"; 
open my $out_fh1, '>', $residue_file or die "Can't open $residue_file: $!"; 

say "enter file name"; 
chomp(my $infile_name = <STDIN>); 
open IN, '<', $infile_name or die "Can't open $infile_name: $!"; 

LINE: while (<IN>) { 
    chomp; 
    /^ATOM/ or next LINE; 
    my @fields = split; 
    push @all_pos, [ @fields[6 .. 8] ]; 
    push @all_data, [ @fields[1 .. 3, 5] ]; 

    if($fields[6] eq $r_no) { 
     say $_; 
     push @matching_pos, [ @fields[6 .. 8] ]; 
     push @matching_data, [ @fields[1 .. 3, 5] ]; 
    } 
} 

say $out_fh1 "Neighbouring residues at position $r_no in the 4.5A region are:"; 
my %close_atoms; 

MATCHING_ATOM: 
for my $i1 (0 .. $#matching_pos) { 
    my $matching_atom = $matching_data[$i1][1]; 
    $matching_atom eq $_ and next MATCHING_ATOM for qw/N CA O C/; 
    for my $i2 (0 .. $#all_pos) { 
     my ($close_atomno, $close_residueno) = @{$all_data[$i2]}[0, 3]; 
     my $dist = distance($matching_pos[$i1], $all_pos[$i2]); 
     if($dist < 4.5 and $close_residueno != $r_no) { 
     $close_atoms{$close_atomno} = 1; 
     } 
    } 
} 

sub distance { sqrt sum map {$_**2} map {$_[0][$_] - $_[1][$_]} 0 .. $#{$_[0]} }; 

my @close_atoms = keys %close_atoms; 

say $out_fh1 "@close_atoms"; 
for my $m (0 .. $#close_atoms) { 
    say $out_fh1 $all_pos[$m];# here is problem i want residue details according to $close_atomno 
} 

say "result in $residue_file"; 

Это будет типичный входной файл:

ATOM 9 N GLU A 1 35.540 1.925 27.662 1.00 19.70 N 
ATOM 10 CA GLU A 1 35.626 1.018 28.802 1.00 20.96 C 
ATOM 11 C GLU A 1 34.264 0.794 29.444 1.00 20.22 C 

с колоннами в таком порядке:

  • атом
  • число атомов
  • имя атома
  • имя остаток
  • цепь
  • номер остаток
  • х коорд
  • у коорд
  • г коорд
  • значения
  • значения
  • значения
+1

Я знаю много о Perl, но ничего о файлах «PDB». Не могли бы вы предоставить короткий тестовый пример (пример ввода + ожидаемый результат для этого ввода), который упрощает понимание происходящего и что происходит неправильно? – amon

+0

Вы можете [изменить] (http://stackoverflow.com/posts/16355323/edit) свой вопрос включить файлы – amon

+0

ATOM 9 N GLU A 1 35.540 1.925 27.662 1.00 19.70 N ATOM 10 CA GLU A 1 35.626 1.018 28.802 1.00 20.96 C ATOM 11 C GLU A 1 34.264 0.794 29.444 1.00 20.22 C 1-й столбец является атомом, второй является номером атома, 3-е имя атома, 4-е имя остатка, 5-е звено 6-е число остатков, 7,8 , 9 - координаты, чтобы взять dist, через атом нет, я хочу извлечь соответствующее имя остатка, нет. и имя атома и нет. столбец 10, 11 не используются. – user2272550

ответ

0
use warnings; 
use strict; 

my $base_r_no = $ARGV[0]; 
open IN, "<$ARGV[1]" or die; 

my @atoms_data = map {[split]} grep {/^ATOM/} <IN>; 
foreach my $base_atom (grep {($$_[2] !~ /^N|CA|O|C$/) && ($$_[5] eq $base_r_no)} @atoms_data) { 
    foreach my $matched_atom (grep {dist($base_atom,$_) < 4.5} @atoms_data) { 
     print join("\t",@{$matched_atom})."\n"; 
    } 
} 

sub dist { 
    return sqrt(($$_[1][5]-$$_[0][5])**2 + ($$_[1][6]-$$_[0][6])**2 + ($$_[1][7]-$$_[0][7])**2); 
} 

Что касается нескольких файлов pdb - я бы просто поместил все файлы в один файл и передал его скрипту.

Смежные вопросы