Mar
18

操作pdb文件的类Perl版

Categories: BIO, Perl | Tags: , | Views: 831

 

前面写过一个C++的版本, 功能写得非常的简单, 并且没有考虑到有MODEL的情况.

终于下定决心用Perl来重写这个类. 主要是由于以前没有用过Perl的面向对象写东西,很生疏,怕写不好,但是真正写起来却发现是那么地顺手

perl不用费尽心思一个个地写函数的重载, 正则表达式比substr更好用. 甚至get和set可以写在一个函数里.

包括5个pm,放在目录pdbClass目录下, 所以每个pm里都是use pdbClass::

(pdbClass父目录下的perl文件)调用的时候use pdbClass::Pdb就行了.

Atom.pm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

package Atom;
use strict;
 
sub new {
    my $class    = shift;
    my $self    = {
        SERIAL        => undef,   # id of the atom
        NAME         => undef,    # full name of the atom, like CA
        RESNAME        => undef,    # name of the residue of the atom
        CHAINID        => undef,    # chain label of the atom
        RESSEQ        => undef,    # residue id of the atom
        X            => undef,    # x
        Y            => undef,    # y
        Z            => undef,    # z
        OCCUPANCY    => undef,    # occupancy
        TEMPFACTOR    => undef,    # temperature factor
    };
    bless $self, $class;
    return $self;
}
 
# Parameters: - / $line 
# Get information from an atom text line 
# or return an atom text line
sub line{
    my $self = shift;
    if (@_) { 
        my $line = shift;
        ($self->{SERIAL},$self->{NAME},$self->{RESNAME},$self->{CHAINID},$self->{RESSEQ},$self->{X},$self->{Y},$self->{Z},$self->{OCCUPANCY},$self->{TEMPFACTOR}) = ($line =~ /^ATOM\s+(\d+)\s+([A-Z0-9]+)\s*([A-Z]{3})\s+([A-Z0-9]{1})\s*([a-z0-9]+)\s+([0-9.-]+)\s+([0-9.-]+)\s+([0-9.-]+)\s+([0-9.-]+)\s+([0-9.-]+).+$/); 
        # ATOM   2843  C   ASN C 124     143.867  14.117  34.756  1.00 39.65
    }
    return     sprintf(
                "ATOM  %5d  %-3s %-3s %s%4s    %8.3f%8.3f%8.3f%6.2f%6.2f",
                $self->{SERIAL},
                $self->{NAME},
                $self->{RESNAME},
                $self->{CHAINID},
                $self->{RESSEQ},
                $self->{X},
                $self->{Y},
                $self->{Z},
                $self->{OCCUPANCY},
                $self->{TEMPFACTOR}
            );
}
 
# set or get serial of this atom
sub serial{
    my $self = shift;
    if (@_) { $self->{SERIAL} = shift; }
    return $self->{SERIAL};
}
 
# set or get name of this atom
sub name{
    my $self = shift;
    if (@_) { $self->{NAME} = shift; }
    return $self->{NAME};
}
 
# set or get resName of this atom
sub resName{
    my $self = shift;
    if (@_) { $self->{RESNAME} = shift; }
    return $self->{RESNAME};
}
 
# set or get chainId of this atom
sub chainId{
    my $self = shift;
    if (@_) { $self->{CHAINID} = shift; }
    return $self->{CHAINID};
}
 
# set or get resSeq of this atom
sub resSeq{
    my $self = shift;
    if (@_) { $self->{RESSEQ} = shift; }
    return $self->{RESSEQ};
}
 
#set or get x
sub x{
    my $self = shift;
    if (@_) { $self->{X} = shift; }
    return $self->{X};
}
 
#set or get y
sub y{
    my $self = shift;
    if (@_) { $self->{Y} = shift; }
    return $self->{Y};
}
 
#set or get z
sub z{
    my $self = shift;
    if (@_) { $self->{Z} = shift; }
    return $self->{Z};
}
 
#set or get occupancy
sub occupancy{
    my $self = shift;
    if (@_) { $self->{OCCUPANCY} = shift; }
    return $self->{OCCUPANCY};
}
 
#set or get temperature factor
sub tempFactor{
    my $self = shift;
    if (@_) { $self->{TEMPFACTOR} = shift; }
    return $self->{TEMPFACTOR};
}
 
#calculate the distance from this atom to another 
sub distanceFrom{
    my $self = shift;
    my $atom = shift;
    return sqrt(
        ($self->{X}-$atom->x)*($self->{X}-$atom->x)+
        ($self->{Y}-$atom->y)*($self->{Y}-$atom->y)+
        ($self->{Z}-$atom->z)*($self->{Z}-$atom->z)
    );
}
 
1;

Residue.pm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

package Residue;
use strict;
use pdbClass::Atom;
 
sub new{
    my $class     = shift;
    my $self    = {
        ATOMS    => [],
    };
    bless $self,$class;
    return $self;
}
 
# usage: atoms 
#        atoms($subscript)
#        atoms($serial,1)
# return all atoms or an atom by subscript
# or an atom by its serial/id.
sub atoms{
    my $self = shift;
    if (@_) {
        my $i = shift;
        if (@_) {
            my $j;
            for ($j=0; $j<@{$self->{ATOMS}}; $j++){
                if ($self->{ATOMS}[$j]->id eq $i) {
                    $i = $j; last;
                }
            }
            if ($j == @{$self->{ATOMS}}){
                die("ATOMS: label not exists!\n");
            }
        } elsif ($i<0 || $i>@{$self->{ATOMS}}){
            die("ATOMS: subscript out of range!\n");
        } 
        return $self->{ATOMS}[$i]; 
    }
    return @{$self->{ATOMS}};
}
 
sub name{  #get residue name
    my $self = shift;
    return $self->{ATOMS}[0]->resName;
}
 
sub chainId{ # get chain id
    my $self = shift;
    return $self->{ATOMS}[0]->chainId;
}
 
sub seq{ # get residue id
    my $self = shift;
    return $self->{ATOMS}[0]->resSeq;    
}
 
sub ca{  # get atom CA
    my $self = shift;
    foreach my $atom (@{$self->{ATOMS}}) {
        if($atom->name eq "CA"){ return $atom; }
    }
}
 
sub addAtom{ #add atom to this residue
    my $self = shift;
    push @{$self->{ATOMS}},@_;
}
 
sub size{ # amount of atoms
    my $self = shift;
    return scalar(@{$self->{ATOMS}});
}
 
sub abbr{ # get abbreviation name
    my $self = shift;
    my %aa = (
        'ALA' => 'A',
        'GLY' => 'G',
        'LEU' => 'L',
        'ILE' => 'I',
        'MET' => 'M',
        'SER' => 'S',
        'THR' => 'T',
        'TYR' => 'Y',
        'PHE' => 'F',
        'ASP' => 'D',
        'ASN' => 'N',
        'GLU' => 'E',
        'GLN' => 'Q',
        'ARG' => 'R',
        'HIS' => 'H',
        'CYS' => 'C',
        'TRP' => 'W',
        'VAL' => 'V',
        'LYS' => 'K',
        'PRO' => 'P',
    );
    return $aa{$self->{ATOMS}[0]->resName};
}
1;

Chain.pm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

package Chain;
use strict;
use pdbClass::Residue;
 
sub new{
    my $class     = shift;
    my $self    = {
        RESIDUES=> [],
    };
    bless $self,$class;
    return $self;
}
 
# usage: residues 
#        residues($subscript)
#        residues($id,1)
# return all residues or a residue by subscript
# or a residue by its id.
sub residues{
    my $self = shift;
    if (@_) {
        my $i = shift;
        if (@_) {
            my $j;
            for ($j=0; $j<@{$self->{RESIDUES}}; $j++){
                if ($self->{RESIDUES}[$j]->id eq $i) {
                    $i = $j; last;
                }
            }
            if ($j == @{$self->{RESIDUES}}){
                die("RESIDUES: label not exists!\n");
            }
        } elsif ($i<0 || $i>@{$self->{RESIDUES}}){
            die("RESIDUES: subscript out of range!\n");
        } 
        return $self->{RESIDUES}[$i]; 
    }
    return @{$self->{RESIDUES}};
}
 
sub id{  #get chain id
    my $self = shift;
    return $self->{RESIDUES}[0]->chainId;
}
 
sub addResidue{ #add residue to this chain
    my $self = shift;
    push @{$self->{RESIDUES}},@_;
}
 
sub size{ # how many residues this chain contains
    my $self = shift;
    return scalar(@{$self->{RESIDUES}});
}
 
# get a sequence by abbreviation of residues' name
sub sequence{ 
    my $self = shift;
    my $seq = "";
    foreach my $residue (@{$self->{RESIDUES}}){
        $seq .= $residue->abbr;
    }
    return $seq;
}
 
1;

Model.pm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

package Model;
use strict;
use pdbClass::Chain;
 
sub new{
    my $class     = shift;
    my $self    = {
        CHAINS    => [],
        ID        => undef,
    };
    bless $self,$class;
    return $self;
}
 
# usage: chains 
#        chains($subscript)
#        chains($label,1)
# return all chains or a chain by subscript
# or a chain by its label.
sub chains{
    my $self = shift;
    if (@_) {
        my $i = shift;
        if (@_) {
            my $j;
            for ($j=0; $j<@{$self->{CHAINS}}; $j++){
                if ($self->{CHAINS}[$j]->id eq $i) {
                    $i = $j; last;
                }
            }
            if ($j == @{$self->{CHAINS}}){
                die("CHAINS: label not exists!\n");
            }
        } elsif ($i<0 || $i>@{$self->{CHAINS}}){
            die("CHAINS: subscript out of range!\n");
        } 
        return $self->{CHAINS}[$i]; 
    }
    return @{$self->{CHAINS}};
}
 
sub id{  #get model id
    my $self = shift;
    if (@_) { $self->{ID} = shift; }
    return $self->{ID};
}
 
sub addChain{ #add chain to this model
    my $self = shift;
    push @{$self->{CHAINS}},@_;
}
 
sub size{ # amount of chains
    my $self = shift;
    return scalar(@{$self->{CHAINS}});
}
 
1;

Pdb.pm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

package Pdb;
use strict;
use pdbClass::Model;
 
sub new{
    my $class     = shift;
    my $self    = {
        MODELS    => [],
        ID        => undef,   # like 1A0O
        CLASH    => undef,   # like 1A0O_AB  (comes from 1A0O_AB.pdb)
    };
    bless $self,$class;
    return $self;
}
 
# load pdb from a pdb file
sub load{
    my $self = shift;
    if(!@_){die("Pdb file name needed!\n");}
    my $file = shift;
    $self->{ID} = substr($file,rindex($file,"/")+1,4);
    $self->{CLASH} = substr($file,rindex($file,"/")+1,-4);
    my ($modelFlag, $chainFlag, $residueFlag);
    my ($atom, $residue, $chain, $model);
 
    open(FILE,$file);
    while(my $line = <FILE>){
        if($line =~ m/^MODEL\s+(.+)$/){  # if contains models
            if($modelFlag ne $1){ 
                $modelFlag = $1;
                $model = Model->new;
                $model->id($modelFlag);
                push @{$self->{MODELS}}, $model;
            }
        }
        if($line =~ m/^ATOM.+$/){
            $atom = Atom->new;
            $atom->line($line);
            if(!defined($model)){   # if not contains models,
                $model = Model->new;# define one, and set id 0
                $model->id(0);
                push @{$self->{MODELS}}, $model;
            }
            if($chainFlag ne $atom->chainId){
                $chainFlag = $atom->chainId; 
                $chain = Chain->new;
                $model->addChain($chain); 
            }
            if($residueFlag ne $atom->resSeq){
                $residueFlag = $atom->resSeq;    
                $residue = Residue->new; 
                $chain->addResidue($residue);
            }    
            $residue->addAtom($atom);
        }
    }
    close(FILE);
}
 
# get or set clash
sub clash{
    my $self = shift;
    if (@_) { $self->{CLASH} = @_; }
    return $self->{CLASH};
}
 
# usage: models 
#        models($subscript)
#        models($label,1)
# return all models or a model by subscript
# or a model by its label.
sub models{
    my $self = shift;
    if (@_) {
        my $i = shift;
        if (@_) {
            my $j;
            for ($j=0; $j<@{$self->{MODELS}}; $j++){
                if ($self->{MODELS}[$j]->id eq $i) {
                    $i = $j; last;
                }
            }
            if ($j == @{$self->{MODELS}}){
                die("MODELS: label not exists!\n");
            }
        } elsif ($i<0 || $i>@{$self->{MODELS}}){
            die("MODELS: subscript out of range!\n");
        } 
        return $self->{MODELS}[$i]; 
    }
    return @{$self->{MODELS}};
}
 
sub id{  #get pdb id
    my $self = shift;
    if (@_) { $self->{ID} = @_; }
    return $self->{ID};
}
 
sub addModel{ #add models
    my $self = shift;
    push @{$self->{MODELS}},@_;
}
 
sub size{ # amount of models
    my $self = shift;
    return scalar(@{$self->{MODELS}});
}
 
1;
这篇文章来自 迷途知返(PWWANG.COM), 转载请注明出处。 版权说明

Leave a comment | Trackback
No comments yet.
;) :| :x :twisted: :roll: :oops: :o :mrgreen: :lol: :idea: :evil: :cry: :arrow: :P :D :?: :? :) :( :!: 8O 8)

你可以使用@somebody:开头, 来邮件通知somebody你回复了他的留言(用户名区分大小写).