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
129
|
# Package: Bios2cor
# This file is part of Bios2cor R package.
# 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_circular <- function(dynamic_structure, res_selection= c("C","I","L","M","V","R","H","K","D","E","N","Q","F","Y","W","T","S","P")){
if(missing(dynamic_structure)){
stop("A 'dynamic_structure' object is required")
}
# Importing torsional information
tor <- dynamic_structure$tor
nb_torsions <- length(tor[1,])
# Torsional angles
sidechain.tor <- tor[, 1:nb_torsions]
dihed_names <- colnames(sidechain.tor)
# Transforms the torsion matrix in circular object (circular package)
dihed <- circular(sidechain.tor, units="degrees", type="angles", modulo="2pi")
# Computes circular correlation
dihed_corr <- cor.circular(dihed)
colnames(dihed_corr) <- dihed_names
rownames(dihed_corr) <- dihed_names
diag(dihed_corr) <- 0
# Squares because only absolute value of correlation matters
dihed_corr <- dihed_corr^2
dihed_corr_names <- colnames(dihed_corr)
if(!is.null(res_selection)){
# Position selection
tor.seq <- dynamic_structure$tor.seq
residue_selection.inds <- which(tor.seq %in% res_selection)
dihed_corr <- dihed_corr[, residue_selection.inds]
dihed_corr <- dihed_corr[residue_selection.inds, ]
nb_angles <- length(dihed_corr[1,])
selected_dihed_names <- colnames(dihed_corr)
}
COV2 <- dihed_corr
# Removing "Inf" and "NA" values
COV2[is.infinite(COV2)] <- 0
COV2[is.na(COV2)] <- 0
res <- list()
# Save matrix of scores
res$score <- COV2
# Compute and save matrix of Z_scores
# Mean and stdev must be calculated on off diagonal elements
mean_up <- mean(COV2[upper.tri(COV2)])
stdev <- sd(COV2[upper.tri(COV2)])
COV3 <- (COV2-mean_up)/stdev
diag(COV3) <- 0
res$Zscore <- COV3
# Save matrices of score and Zscores without auto correlation
# Create and save matrix with 0 values for autocorrelation (dihedral angles within the same residue)
COV5 <- COV2
res_num <- c()
res_num <- unlist(lapply(selected_dihed_names, function(x){strsplit(x, "[.]")[[1]][1]}))
# Create matrix with NA that will be used for Zscores
for(i in 1:nb_angles){
for(j in i:nb_angles){
res_i <- res_num[i]
res_j <- res_num[j]
if(res_i == res_j){
COV5[i,j] <- NA
COV5[j,i] <- NA
}
}
}
# Save matrix of score with 0 for autocorrelation
COV4 <-COV5
COV4[is.na(COV4)] <- 0
res$score_noauto <- COV4
# Compute and save matrix of Z_scores without auto correlation
# Mean and stdev are calculated on non NA elements of upper triangle
mean_noauto <- mean(COV5[upper.tri(COV5)],na.rm=TRUE)
stdev_noauto <- sd(COV5[upper.tri(COV5)],na.rm=TRUE)
COV5[is.na(COV5)] <- 0
COV5 <- (COV5-mean_noauto)/stdev_noauto
### Correct diag and autocorrelated elements
diag(COV5) <- 0
for(i in 1:nb_angles){
for(j in i:nb_angles){
res_i <- res_num[i]
res_j <- res_num[j]
if(res_i == res_j){
COV5[i,j] <- 0
COV5[j,i] <- 0
}
}
}
res$Zscore_noauto <- COV5
return(res)
}
|
