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
|
# Bios2cor is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# See the GNU General Public License at:
# http://www.gnu.org/licenses/
dynamic_structure <- function(pdb, trj, frames=NULL){
if (missing(pdb)){
stop("Missing pdb file")
}
if (missing(trj)){
stop("Missing trajectory file")
}
## Transform the xyz trajectory into dihedral trajectory
# Reading pdb file
pdb <- read.pdb(pdb)
trj <- read.dcd(trj)
# Selecting only "CA" atoms for superposition
ca.inds <- atom.select(pdb, elety = "CA")
# Getting xyz coordinates using fit.xyz from bio3D
xyz <- fit.xyz(fixed = pdb$xyz, mobile = trj, fixed.inds = ca.inds$xyz, mobile.inds = ca.inds$xyz)
nb_frames <- length(xyz[,1])
# Creating torsion object using the xyz2torsion function from bio3D
if(is.null(frames)) {
tor <- xyz2torsion(pdb, xyz, tbl = "sidechain", ncore= 1)
} else {
tor <- xyz2torsion(pdb, xyz[frames,], tbl = "sidechain", ncore= 1)
}
nb_torsions <- length(tor[1,])
nb_frames <- length(tor[,1])
# Associating one letter AA code to each residue, using the pdbseq function from bio3D
prot.seq <- pdbseq(pdb)
tor.names <- colnames(tor)
tor.resno <- sub("\\..*$", "", tor.names)
tor.angle <- sub("^[0-9]*\\.","", tor.names)
tor.seq <- prot.seq[tor.resno]
# Filling torsional structure
res <- list()
res$pdb <- pdb
res$trj <- trj
res$ca.inds <- ca.inds
res$xyz <- xyz
res$nb_frames <- nb_frames
if(!is.null(frames)) {
res$frames <- frames
} else {
res$frames <- c(1:nb_frames)
}
res$tor <- tor
res$nb_torsions <- nb_torsions
res$prot.seq <- prot.seq
res$tor.names <- tor.names
res$tor.resno <- tor.resno
res$tor.angle <- tor.angle
res$tor.seq <- tor.seq
class (res) <- c("structure")
return (res)
}
|
