model.cpp 7.71 KB
Newer Older
Pilar Cossio's avatar
License    
Pilar Cossio committed
1
2
3
4
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        < BioEM software for Bayesian inference of Electron Microscopy images>
            Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
            Max Planck Institute of Biophysics, Frankfurt, Germany.
5

Pilar Cossio's avatar
License    
Pilar Cossio committed
6
7
8
9
                See license statement for terms of distribution.

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

10
11
12
13
14
15
16
17
18
19
20
21
22
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstring>

#include "model.h"

using namespace std;


int bioem_model::readModel()
{
David Rohr's avatar
David Rohr committed
23
24
25
	// **************************************************************************************
	// ***************Reading reference Models either PDB or x,y,z,r,d format****************
	// **************************************************************************************
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

	ofstream exampleReadCoor;
	exampleReadCoor.open ("COORDREAD");

	std::ifstream input(filemodel);
	if(readPDB)
	{
		ifstream input(filemodel);
		if (!input.good())
		{
			cout << "PDB Failed to open file" << endl ; // pdbfilename << " ("<<filename<<")\n";
			exit(0);
		}

		char line[512] = {' '};
		char tmpLine[512] = {' '};
42
43
		int numres = 0;
		NormDen = 0.0;
44
45
46
47
48

		//  cout << " HERE	" << filemodel ;
		// for eachline in the file
		while (!input.eof())
		{
49
			input.getline(line, 512);
50

51
52
			strncpy(tmpLine, line, strlen(line));
			char *token = strtok(tmpLine, " ");
53

54
			if (strcmp(token, "ATOM") == 0) // Optional,Mandatory if standard residues exist
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
			{
				/*
				1-6 			"ATOM  "
				7 - 11		 Integer		 serial		Atom serial number.
				13 - 16		Atom			name		  Atom name.
				17			 Character	   altLoc		Alternate location indicator.
				18 - 20		Residue name	resName	   Residue name.
				22			 Character	   chainID	   Chain identifier.
				23 - 26		Integer		 resSeq		Residue sequence number.
				27			 AChar		   iCode		 Code for insertion of residues.
				31 - 38		Real(8.3)	   x			 Orthogonal coordinates for X in
				39 - 46		Real(8.3)	   y			 Orthogonal coordinates for Y in
				47 - 54		Real(8.3)	   z			 Orthogonal coordinates for Z in
					 */

				char name[5] = {"PPP"};
				char resName[3] = {' '};
72
73
74
				float x = 0.0;
				float y = 0.0;
				float z = 0.0;
75
76
77
				char tmp[6] = {' '};

				// parse name
78
79
				strncpy(tmp, line + 12, 4);
				sscanf (tmp, "%s", name);
80
81

				// parse resName
82
83
				strncpy(tmp, line + 17, 3);
				sscanf (tmp, "%s", resName);
84
85
86
87

				// parse x, y, z
				char tmpVals[36]  = {' '};

88
89
				strncpy (tmpVals, line + 30, 8);
				sscanf (tmpVals, "%f", &x);
90

91
92
				strncpy (tmpVals, line + 38, 8);
				sscanf (tmpVals, "%f", &y);
93

94
95
				strncpy (tmpVals, line + 46, 8);
				sscanf (tmpVals, "%f", &z);
96

97
				if (strcmp(name, "CA") == 0)
98
				{
99
100
101
102
103
					if (numres >= BIOEM_MODEL_SIZE)
					{
						cout << "BIOEM_MODEL_SIZE too small\n";
						exit(1);
					}
104
					//Getting residue Radius and electron density
105
106
107
					radiusPointsModel[numres] = getAminoAcidRad(resName);
					densityPointsModel[numres] = getAminoAcidDensity(resName);
					NormDen += densityPointsModel[numres];
108
109

					//Getting the coordinates
110
111
112
					PointsModel[numres].pos[0] = (myfloat_t) x;
					PointsModel[numres].pos[1] = (myfloat_t) y;
					PointsModel[numres].pos[2] = (myfloat_t) z;
113
					exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres]  << "\n";
114
115
116
117
118
119
					numres++;
				}
			}


		}
120
		nPointsModel = numres;
David Rohr's avatar
David Rohr committed
121
		cout << "Protein structure read from PDB\n";
122
123
124
125
	}
	else //Reading model from FILE FORMAT x,y,z,rad,density
	{
		char line[128];
126
127
		int numres = 0;
		NormDen = 0.0;
128
		FILE *file = fopen ( filemodel , "r" );
David Rohr's avatar
David Rohr committed
129
		if (file == NULL)
130
		{
131
			cout << "Error opening file " << filemodel << "\n";
David Rohr's avatar
David Rohr committed
132
133
134
			exit(1);
		}
		while ( fgets ( line, sizeof line, file ) != NULL )
135
		{
David Rohr's avatar
David Rohr committed
136
			if (numres >= BIOEM_MODEL_SIZE)
137
			{
David Rohr's avatar
David Rohr committed
138
139
				cout << "BIOEM_MODEL_SIZE too small\n";
				exit(1);
140
			}
David Rohr's avatar
David Rohr committed
141
142
143
144
145
146
147
148
149
150
151
			float tmpval[5];
			sscanf(line, "%f %f %f %f %f", &tmpval[0], &tmpval[1], &tmpval[2], &tmpval[3], &tmpval[4]);
			PointsModel[numres].pos[0] = (myfloat_t) tmpval[0];
			PointsModel[numres].pos[1] = (myfloat_t) tmpval[1];
			PointsModel[numres].pos[2] = (myfloat_t) tmpval[2];
			radiusPointsModel[numres] = (myfloat_t) tmpval[3];
			densityPointsModel[numres] = (myfloat_t) tmpval[4];

			exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
			NormDen += densityPointsModel[numres];
			numres++;
152
		}
153
		fclose(file);
David Rohr's avatar
David Rohr committed
154
155
		nPointsModel = numres;
		cout << "Protein structure read from Standard File\n";
156
	}
David Rohr's avatar
David Rohr committed
157
158
	cout << "Total Number of Voxels " << nPointsModel ;
	cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
159
160
161
162
163
	exampleReadCoor.close();

	//Moving to Model to its center of mass:
	myfloat3_t r_cm;

164
	for(int n = 0; n < 3; n++)r_cm.pos[n] = 0.0;
165

166
	for(int n = 0; n < nPointsModel; n++)
167
168
169
170
171
	{
		r_cm.pos[0] += PointsModel[n].pos[0];
		r_cm.pos[1] += PointsModel[n].pos[1];
		r_cm.pos[2] += PointsModel[n].pos[2];
	}
172
173
174
175
	r_cm.pos[0] = r_cm.pos[0] / (myfloat_t) nPointsModel;
	r_cm.pos[1] = r_cm.pos[1] / (myfloat_t) nPointsModel;
	r_cm.pos[2] = r_cm.pos[2] / (myfloat_t) nPointsModel;
	for(int n = 0; n < nPointsModel; n++)
176
	{
177
178
179
		PointsModel[n].pos[0] -= r_cm.pos[0];
		PointsModel[n].pos[1] -= r_cm.pos[1];
		PointsModel[n].pos[2] -= r_cm.pos[2];
180
181
182
183

	}

	return(0);
184
185
186
187
}

myfloat_t bioem_model::getAminoAcidRad(char *name)
{
David Rohr's avatar
David Rohr committed
188
	// *************** Function that gets the radius for each amino acid ****************
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
	myfloat_t iaa = 0;

	if(std::strcmp(name, "CYS") == 0)iaa = 2.75;
	if(std::strcmp(name, "PHE") == 0)iaa = 3.2;
	if(std::strcmp(name, "LEU") == 0)iaa = 3.1;
	if(std::strcmp(name, "TRP") == 0)iaa = 3.4;
	if(std::strcmp(name, "VAL") == 0)iaa = 2.95;
	if(std::strcmp(name, "ILE") == 0)iaa = 3.1;
	if(std::strcmp(name, "MET") == 0)iaa = 3.1;
	if(std::strcmp(name, "HIS") == 0)iaa = 3.05;
	if(std::strcmp(name, "TYR") == 0)iaa = 3.25;
	if(std::strcmp(name, "ALA") == 0)iaa = 2.5;
	if(std::strcmp(name, "GLY") == 0)iaa = 2.25;
	if(std::strcmp(name, "PRO") == 0)iaa = 2.8;
	if(std::strcmp(name, "ASN") == 0)iaa = 2.85;
	if(std::strcmp(name, "THR") == 0)iaa = 2.8;
	if(std::strcmp(name, "SER") == 0)iaa = 2.6;
	if(std::strcmp(name, "ARG") == 0)iaa = 3.3;
	if(std::strcmp(name, "GLN") == 0)iaa = 3.0;
	if(std::strcmp(name, "ASP") == 0)iaa = 2.8;
	if(std::strcmp(name, "LYS") == 0)iaa = 3.2;
	if(std::strcmp(name, "GLU") == 0)iaa = 2.95;
211
212
213
214
215
216
217

	if(iaa == 0)
	{
		cout << "PROBLEM WITH AMINO ACID " << name << endl;
		exit(0);
	}
	return iaa;
218
219
220
221
222

}

myfloat_t bioem_model::getAminoAcidDensity(char *name)
{
David Rohr's avatar
David Rohr committed
223
	// *************** Function that gets the number of electrons for each amino acid ****************
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
	myfloat_t iaa = 0.0;

	if(std::strcmp(name, "CYS") == 0)iaa = 64.0;
	if(std::strcmp(name, "PHE") == 0)iaa = 88.0;
	if(std::strcmp(name, "LEU") == 0)iaa = 72.0;
	if(std::strcmp(name, "TRP") == 0)iaa = 108.0;
	if(std::strcmp(name, "VAL") == 0)iaa = 64.0;
	if(std::strcmp(name, "ILE") == 0)iaa = 72.0;
	if(std::strcmp(name, "MET") == 0)iaa = 80.0;
	if(std::strcmp(name, "HIS") == 0)iaa = 82.0;
	if(std::strcmp(name, "TYR") == 0)iaa = 96.0;
	if(std::strcmp(name, "ALA") == 0)iaa = 48.0;
	if(std::strcmp(name, "GLY") == 0)iaa = 40.0;
	if(std::strcmp(name, "PRO") == 0)iaa = 62.0;
	if(std::strcmp(name, "ASN") == 0)iaa = 66.0;
	if(std::strcmp(name, "THR") == 0)iaa = 64.0;
	if(std::strcmp(name, "SER") == 0)iaa = 56.0;
	if(std::strcmp(name, "ARG") == 0)iaa = 93.0;
	if(std::strcmp(name, "GLN") == 0)iaa = 78.0;
	if(std::strcmp(name, "ASP") == 0)iaa = 59.0;
	if(std::strcmp(name, "LYS") == 0)iaa = 79.0;
	if(std::strcmp(name, "GLU") == 0)iaa = 53.0;
246
247
248
249
250
251
252

	if(iaa == 0.0)
	{
		cout << "PROBLEM WITH AMINO ACID " << name << endl;
		exit(0);
	}
	return iaa;
253
254
}