Skip to content

Commit 1d33cf5

Browse files
jkleinjjkleinj
committed
CA distance matrix between chains is implemented
1 parent 3203ea5 commit 1d33cf5

8 files changed

Lines changed: 92 additions & 50 deletions

File tree

‎POPSC/src/getmmcif.c‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ __inline__ static char aacode(char *code3)
3232
/* three-letter code of amino acid residues, exception HET (X) */
3333
char *aa3[] = {"ALA","---","CYS","ASP","GLU","PHE","GLY","HIS","ILE","---","LYS","LEU","MET","ASN","---","PRO","GLN","ARG","SER","THR","UNL","VAL","TRP","HET","TYR","UNK"};
3434
/* nucleotide residues */
35-
char *nuc[] = {" A"," DA"," C"," DC","---","---"," G"," DG"," I"," DI","---","---","---"," N"," DN","---","---","---"," DT"," T"," U"," DU","---","---","---","---"};
35+
char *nuc[] = {"A","DA","C","DC","---","---","G","DG","I"," DI","---","---","---","N","DN","---","---","---","DT","T","U","DU","---","---","---","---"};
3636

3737
/* match against amino acid residues */
3838
residue = scan_array(code3, aa3, 65);
@@ -100,7 +100,8 @@ int map_structure_mmcif(Arg *arg, Argpdb *argpdb, Str *str, Structure *s) {
100100
regex_t *regexPattern = 0; /* regular atom patterns */
101101
/* allowed HETATM atom types (standard N,CA,C,O) and elements (any N,C,O,P,S) */
102102
const int nHetAtom = 9;
103-
char hetAtomPattern[9][32] = {{" N "},{" CA "},{" C "},{" O "},{".{1}C[[:print:]]{1,3}"},{".{1}N[[:print:]]{1,3}"},{".{1}O[[:print:]]{1,3}"},{".{1}P[[:print:]]{1,3}"},{".{1}S[[:print:]]{1,3}"}};
103+
104+
char hetAtomPattern[9][32] = {{"N"},{"CA"},{"C"},{"O"},{".{1}C[[:print:]]{1,3}"},{".{1}N[[:print:]]{1,3}"},{".{1}O[[:print:]]{1,3}"},{".{1}P[[:print:]]{1,3}"},{".{1}S[[:print:]]{1,3}"}};
104105

105106
/*____________________________________________________________________________*/
106107
/* initialise/allocate memory for set of (64) selected (CA) atom entries */
@@ -209,8 +210,8 @@ int map_structure_mmcif(Arg *arg, Argpdb *argpdb, Str *str, Structure *s) {
209210
}
210211

211212
/* detect CA/N3 atoms */
212-
if ((strncmp(str->atom[i].atomName, " CA ", 4) == 0) ||
213-
(strncmp(str->atom[i].atomName, " N3 ", 4) == 0)) {
213+
if ((strncmp(str->atom[i].atomName, "CA", 4) == 0) ||
214+
(strncmp(str->atom[i].atomName, "N3", 4) == 0)) {
214215

215216
str->resAtom[k] = i;
216217
str->sequence.res[k++] = aacode(str->atom[i].residueName);

‎POPSC/src/getpdb.c‎

Lines changed: 25 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,27 @@
11
/*==============================================================================
22
getpdb.c : routines for reading PDB structures
3-
Copyright (C) 2004 Jens Kleinjung
3+
Copyright (C) 2004-2026 Jens Kleinjung
44
Read the COPYING file for license information.
55
==============================================================================*/
66

77
#include "getpdb.h"
88

9+
/*____________________________________________________________________________*/
10+
/* remove spaces from atom and residue names */
11+
void remove_spaces(char *s)
12+
{
13+
char *src = s;
14+
char *dst = s;
15+
16+
while (*src) {
17+
if (*src != ' ')
18+
*dst++ = *src;
19+
src++;
20+
}
21+
22+
*dst = '\0';
23+
}
24+
925
/*____________________________________________________________________________*/
1026
/* match PDB residue name against constant residue name array */
1127
__inline__ static char scan_array(char *code3, char *residue_array[], int shift)
@@ -54,8 +70,8 @@ __inline__ static char aacode(char *code3)
5470
__inline__ static int standardise_name(char *residueName, char *atomName)
5571
{
5672
/* GRO 'ILE CD' to PDB 'ILE CD1' */
57-
if ((strcmp(residueName, "ILE") == 0) && (strcmp(atomName, " CD ") == 0))
58-
strcpy(atomName, " CD1");
73+
if ((strcmp(residueName, "ILE") == 0) && (strcmp(atomName, "CD") == 0))
74+
strcpy(atomName, "CD1");
5975

6076
return 0;
6177
}
@@ -155,8 +171,7 @@ int read_pdb(FILE *pdbInFile, gzFile *pdbgzInFile, Arg *arg, Argpdb *argpdb, Str
155171
regex_t *regexPattern = 0; /* regular atom patterns */
156172
/* allowed HETATM atom types (standard N,CA,C,O) and elements (any N,C,O,P,S) */
157173
const int nHetAtom = 9;
158-
char hetAtomPattern[9][32] = {{" N "},{" CA "},{" C "},{" O "},{".{1}C[[:print:]]{1,3}"},{".{1}N[[:print:]]{1,3}"},{".{1}O[[:print:]]{1,3}"},{".{1}P[[:print:]]{1,3}"},{".{1}S[[:print:]]{1,3}"}};
159-
/*char hetAtomNewname[9][32] = {{" N "},{" CA "},{" C "},{" O "},{" C_ "},{" N_ "},{" O_ "},{" P_ "},{" S_ "}};*/
174+
char hetAtomPattern[9][32] = {{"N"},{"CA"},{"C"},{"O"},{".{1}C[[:print:]]{1,3}"},{".{1}N[[:print:]]{1,3}"},{".{1}O[[:print:]]{1,3}"},{".{1}P[[:print:]]{1,3}"},{".{1}S[[:print:]]{1,3}"}};
160175

161176
/*____________________________________________________________________________*/
162177
/* initialise/allocate memory for set of (64) selected (CA) atom entries */
@@ -262,12 +277,14 @@ int read_pdb(FILE *pdbInFile, gzFile *pdbgzInFile, Arg *arg, Argpdb *argpdb, Str
262277
str->atom[str->nAtom].recordName[j++] = line[i++];
263278
}
264279
str->atom[str->nAtom].recordName[j] = '\0';
280+
265281

266282
/* atom name */
267283
for (i = 12, j = 0; i < 16; ) {
268284
str->atom[str->nAtom].atomName[j++] = line[i++];
269285
}
270286
str->atom[str->nAtom].atomName[j] = '\0';
287+
remove_spaces(str->atom[str->nAtom].atomName);
271288

272289
/* alternative location */
273290
str->atom[str->nAtom].alternativeLocation[0] = line[16];
@@ -278,10 +295,12 @@ int read_pdb(FILE *pdbInFile, gzFile *pdbgzInFile, Arg *arg, Argpdb *argpdb, Str
278295
str->atom[str->nAtom].residueName[j++] = line[i++];
279296
}
280297
str->atom[str->nAtom].residueName[j] = '\0';
298+
remove_spaces(str->atom[str->nAtom].residueName);
281299

282300
/* chain identifier */
283301
str->atom[str->nAtom].chainIdentifier[0] = line[21];
284302
str->atom[str->nAtom].chainIdentifier[1] = '\0';
303+
remove_spaces(str->atom[str->nAtom].chainIdentifier);
285304

286305
/* residue number */
287306
str->atom[str->nAtom].residueNumber = atoi(&line[22]);
@@ -318,6 +337,7 @@ int read_pdb(FILE *pdbInFile, gzFile *pdbgzInFile, Arg *arg, Argpdb *argpdb, Str
318337
str->atom[str->nAtom].element[j++] = line[i++];
319338
}
320339
str->atom[str->nAtom].element[j] = '\0';
340+
remove_spaces(str->atom[str->nAtom].element);
321341

322342
/* charge */
323343
/*for (i = 78, j = 0; i < 80; )

‎POPSC/src/getpdb.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ Read the COPYING file for license information.
2727

2828
/*____________________________________________________________________________*/
2929
/* prototypes */
30+
void remove_spaces(char *s);
3031
int read_pdb(FILE *pdbfile, gzFile *pdbgzInFile, Arg *arg, Argpdb *argpdb, Str *str);
3132
void read_structure(Arg *arg, Argpdb *argpdb, Str *pdb);
3233

‎POPSC/src/pops‎

1.38 KB
Binary file not shown.

‎POPSC/src/pops.c‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
135135
/*____________________________________________________________________________*/
136136
/** compute molecular topology */
137137
if (! arg.silent) fprintf(stdout, "Topology\n");
138-
init_topology(&pdb, &topol);
138+
init_topology(&arg, &pdb, &topol);
139139
get_topology(&pdb, &type, &topol, constant_sasa, &argpdb, &arg);
140140

141141
/*____________________________________________________________________________*/
@@ -204,7 +204,7 @@ int main(int argc, char *argv[])
204204
assert(traj.frame[i].nAtom == pdb.nAllAtom);
205205
copy_coordinates(&pdb, &traj, i);
206206
/* topology */
207-
init_topology(&pdb, &topol);
207+
init_topology(&arg, &pdb, &topol);
208208
get_topology(&pdb, &type, &topol, constant_sasa, &argpdb, &arg);
209209
/* SASA */
210210
init_sasa(&pdb, &type, &molSasa, constant_sasa, &arg);
@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
228228
free(pdb.resAtom);
229229
free(pdb.atomMap);
230230
free(pdb.sequence.res);
231-
free(pdb.sequence.name);
231+
/*free(pdb.sequence.name);*/
232232

233233
/* trajectory */
234234
if (arg.trajInFileName) {

‎POPSC/src/putDistMatCA.c‎

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -12,17 +12,8 @@ Read the COPYING file for license information.
1212
/** print Calpha distance matrix */
1313
void print_distMatCA(Arg *arg, Topol *topol)
1414
{
15-
unsigned int i, j;
16-
1715
arg->distMatCAOutFile = safe_open(arg->distMatCAOutFileName, "w");
18-
19-
for (i = 0; i < topol->nCA; ++ i) {
20-
for (j = 0; j < topol->nCA; ++ j) {
21-
print_mat2D_float(arg->distMatCAOutFileName,
22-
topol->distMatCA, i, j);
23-
}
24-
}
25-
16+
print_mat2D_float(arg->distMatCAOutFileName, topol->distMatCA, topol->nCA1, topol->nCA2);
2617
fclose(arg->distMatCAOutFile);
2718
}
2819

‎POPSC/src/topol.c‎

Lines changed: 54 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,11 @@ __inline__ static void print_torsion(Str *pdb, int t1, int t2, int t3, int t4)
4242

4343
/*___________________________________________________________________________*/
4444
/** init topology */
45-
void init_topology(Str *pdb, Topol *topol)
45+
void init_topology(Arg *arg, Str *pdb, Topol *topol)
4646
{
4747
unsigned int i;
48+
int *nCA = NULL;
49+
char *chain1 = NULL;
4850
/* assuming an upper limit of 63 bonded interactions per atom */
4951
topol->bondState = alloc_mat2D_int(topol->bondState, pdb->nAtom, 64);
5052
/* assuming an upper limit of 255 non-bonded interactions per atom */
@@ -56,16 +58,23 @@ void init_topology(Str *pdb, Topol *topol)
5658
topol->interfaceNn = safe_malloc(pdb->nAtom * sizeof(int));
5759
topol->interfaceNnDist = safe_malloc(pdb->nAtom * sizeof(float));
5860

59-
/* Calpha distance matrix */
60-
/* determine number of Calpha atoms */
61-
topol->nCA = 0;
62-
for (i = 0; i < pdb->nAtom; ++ i) {
63-
if (strncmp(pdb->atom[i].atomName, " CA ", 4) == 0) {
64-
++ topol->nCA;
61+
/* between-chain or between domain Calpha distance matrix */
62+
/* determine number of Calpha atoms per chain or domain */
63+
nCA = &(topol->nCA1);
64+
for (i = 0, topol->nCA1 = 0, topol->nCA2 = 0; i < pdb->nAtom; ++ i) {
65+
if (strcmp(pdb->atom[i].atomName, "CA") == 0) {
66+
if ((*nCA) == 0) {
67+
chain1 = &(pdb->atom[i].chainIdentifier[0]);
68+
}
69+
if ((*nCA) > 0 && strcmp(pdb->atom[i].chainIdentifier, chain1) != 0) {
70+
nCA = &(topol->nCA2);
71+
}
72+
++ (*nCA);
6573
}
6674
}
67-
topol->distMatCA = alloc_mat2D_float(topol->distMatCA, topol->nCA, topol->nCA);
68-
init_mat2D_float(topol->distMatCA, topol->nCA, topol->nCA, 0.);
75+
printf("CA distance matrix has dimensions %d x %d \n", topol->nCA1, topol->nCA2);
76+
topol->distMatCA = alloc_mat2D_float(topol->distMatCA, topol->nCA1, topol->nCA2);
77+
init_mat2D_float(topol->distMatCA, topol->nCA1, topol->nCA2, 0.);
6978

7079
for (i = 0; i < pdb->nAtom; ++ i) {
7180
topol->bondState[i][0] = 0; /* no bonded pairs recorded */
@@ -100,7 +109,7 @@ void free_topology(Str *pdb, Topol *topol)
100109
free_mat2D_float(topol->neighbourPar, dimMat2D); /* POPS parameters of neighbours */
101110
free(topol->interfaceNn);
102111
free(topol->interfaceNnDist);
103-
free_mat2D_float(topol->distMatCA, topol->nCA);
112+
free_mat2D_float(topol->distMatCA, topol->nCA1);
104113
}
105114

106115
/*___________________________________________________________________________*/
@@ -643,31 +652,50 @@ int get_topology(Str *pdb, Type *type, Topol *topol, ConstantSasa *constant_sasa
643652
}
644653

645654
/*____________________________________________________________________________*/
646-
/** Calpha distances */
655+
/** Calpha distances between different chains */
647656
int calpha_distances(Arg *arg, Str *pdb, Topol *topol) {
648657
unsigned int i, j;
649658
int n_cai, n_caj;
659+
char chain1 = '\0';
660+
661+
n_cai = 0;
650662

651-
for (i = 0, n_cai = 0; i < pdb->nAtom; ++ i) {
652-
if (strncmp(pdb->atom[i].atomName, " CA ", 4) == 0) {
653-
654-
for (j = i+1, n_caj = n_cai + 1; j < pdb->nAtom; ++ j) {
655-
if (strncmp(pdb->atom[j].atomName, " CA ", 4) == 0) {
663+
/* find first CA chain */
664+
for (i = 0; i < pdb->nAtom; ++i) {
665+
if (strcmp(pdb->atom[i].atomName, "CA") == 0) {
666+
chain1 = pdb->atom[i].chainIdentifier[0];
667+
break;
668+
}
669+
}
656670

657-
topol->distMatCA[n_cai][n_caj] =
658-
topol->distMatCA[n_caj][n_cai] =
659-
atom_distance(pdb, i, j);
671+
for (i = 0; i < pdb->nAtom; ++i) {
672+
if (strcmp(pdb->atom[i].atomName, "CA") != 0)
673+
continue;
660674

661-
/*printf("Distance between CA_i %d (%d) and CA_j %d (%d): %f\n",
662-
n_cai, i, n_caj, j, topol->distMatCA[n_cai][n_caj]);*/
675+
/* only rows from chain 1 */
676+
if (pdb->atom[i].chainIdentifier[0] != chain1)
677+
continue;
663678

664-
++ n_caj;
665-
}
666-
}
667-
++ n_cai;
668-
}
679+
n_caj = 0;
680+
681+
for (j = 0; j < pdb->nAtom; ++j) {
682+
if (strcmp(pdb->atom[j].atomName, "CA") != 0)
683+
continue;
684+
685+
/* only columns from other chain(s) */
686+
if (pdb->atom[j].chainIdentifier[0] == chain1)
687+
continue;
688+
689+
topol->distMatCA[n_cai][n_caj] = atom_distance(pdb, i, j);
690+
++n_caj;
691+
}
692+
693+
++n_cai;
669694
}
670695

696+
assert(n_cai == topol->nCA1);
697+
assert(n_caj == topol->nCA2);
698+
671699
return(0);
672700
}
673701

‎POPSC/src/topol.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,8 @@ typedef struct
4848
float **neighbourPar; /* records the POPS parameters 'p_ij * b_ij' of neighbours */
4949
int *interfaceNn; /* nearest neighbour on separate chain */
5050
float *interfaceNnDist; /* distance to nearest neighbour on separate chain */
51-
int nCA; /* number of Calpha atoms */
51+
int nCA1; /* number of Calpha atoms in chain or domain 1 */
52+
int nCA2; /* number of Calpha atoms in chain or domain 2 */
5253
float **distMatCA; /* Calpha distance matrix */
5354
} Topol;
5455

@@ -60,7 +61,7 @@ int get_bonds(Str *pdb, Type *type, Topol *topol, ConstantSasa *constant_sasa, A
6061
int get_angles(Str *pdb, Topol *topol); /* calculate angles (from bonds) */
6162
int get_torsions(Str *pdb, Type *type, Topol *topol, ConstantSasa *constant_sasa); /* calculate torsions (from angles) */
6263
int nonbonded_overlaps(Str *pdb, Type *type, Topol *topol, ConstantSasa *constant_sasa, Arg *arg); /* calculate overlapping atoms */
63-
void init_topology(Str *pdb, Topol *topol);
64+
void init_topology(Arg *arg, Str *pdb, Topol *topol);
6465
void free_topology(Str *pdb, Topol *topol);
6566
int get_topology(Str *pdb, Type *type, Topol *topol, ConstantSasa *constant_sasa, Argpdb *argpdb, Arg *arg); /* call topology routines */
6667
int calpha_distances(Arg *arg, Str *pdb, Topol *topol);

0 commit comments

Comments
 (0)