model.cpp 7.01 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstring>

#include "model.h"

using namespace std;


int bioem_model::readModel()
{
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
	/**************************************************************************************/
	/***************Reading reference Models either PDB or x,y,z,r,d format****************/
	/**************************************************************************************/

	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] = {' '};
		int numres=0;
		NormDen=0.0;

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

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

			if (strcmp(token,"ATOM")==0) // Optional,Mandatory if standard residues exist
			{
				/*
				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] = {' '};
63
64
65
				float x = 0.0;
				float y = 0.0;
				float z = 0.0;
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
				char tmp[6] = {' '};

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

				// parse resName
				strncpy(tmp,line+17,3);
				sscanf (tmp,"%s",resName);

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

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

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

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

				if (strcmp(name,"CA") == 0)
				{
90
91
92
93
94
					if (numres >= BIOEM_MODEL_SIZE)
					{
						cout << "BIOEM_MODEL_SIZE too small\n";
						exit(1);
					}
95
96
97
98
99
100
					//Getting residue Radius and electron density
					radiusPointsModel[numres]=getAminoAcidRad(resName);
					densityPointsModel[numres]=getAminoAcidDensity(resName);
					NormDen+=densityPointsModel[numres];

					//Getting the coordinates
101
102
103
					PointsModel[numres].pos[0] = (myfloat_t) x;
					PointsModel[numres].pos[1] = (myfloat_t) y;
					PointsModel[numres].pos[2] = (myfloat_t) z;
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
					exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]  << "\n";
					numres++;
				}
			}


		}
		nPointsModel=numres;
		cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
	}
	else //Reading model from FILE FORMAT x,y,z,rad,density
	{
		char line[128];
		int numres=0;
		NormDen=0.0;
		FILE *file = fopen ( filemodel , "r" );
		if ( file != NULL )
		{
			while ( fgets ( line, sizeof line, file ) != NULL )
			{
125
126
127
128
129
				if (numres >= BIOEM_MODEL_SIZE)
				{
					cout << "BIOEM_MODEL_SIZE too small\n";
					exit(1);
				}
130
131
132
133
134
135
136
137
				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];

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
				exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]<< "\n";
				NormDen+=densityPointsModel[numres];
				numres++;
			}
			nPointsModel=numres;
			cout << "Protein structure read from Standard File \nTotal Number of Residues " << nPointsModel ;
			cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
		}
	}
	exampleReadCoor.close();

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

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

	for(int n=0; n < nPointsModel; n++)
	{
		r_cm.pos[0] += PointsModel[n].pos[0];
		r_cm.pos[1] += PointsModel[n].pos[1];
		r_cm.pos[2] += PointsModel[n].pos[2];
	}
160
161
162
	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;
163
164
165
166
167
168
169
170
171
	for(int n=0; n < nPointsModel; n++)
	{
		PointsModel[n].pos[0]-= r_cm.pos[0];
		PointsModel[n].pos[1]-= r_cm.pos[1];
		PointsModel[n].pos[2]-= r_cm.pos[2];

	}

	return(0);
172
173
174
175
}

myfloat_t bioem_model::getAminoAcidRad(char *name)
{
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
	/*************** Function that gets the radius for each amino acid ****************/
	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;

	if(iaa == 0)
	{
		cout << "PROBLEM WITH AMINO ACID " << name << endl;
		exit(0);
	}
	return iaa;
206
207
208
209
210

}

myfloat_t bioem_model::getAminoAcidDensity(char *name)
{
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
	/*************** Function that gets the number of electrons for each amino acid ****************/
	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;

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