How to read and extract X Y Z coordinates of pdb file without using pdbstruct and pdbread but using file handling commands?

pdb='C:\Users\Priya\Documents\MATLAB\1ubq.pdb'
fileID=fopen(pdb)
x = fileread('C:\Users\Priya\Documents\MATLAB\1ubq.pdb');
x.Model.Atom(:).AtomName;
x.Model.Atom(:).X
x.Model.Atom(:).Y
x.Model.Atom(:).Z
pos = find(strcmp({x.Model.Atom(:).AtomName},'CA'));
X=[x.Model.Atom(pos).X];
Y=[x.Model.Atom(pos).Y];
Z=[x.Model.Atom(pos).Z];

[Walter Roberson]

Why? What is the restriction against using pdbstruct() or pdbread() ? The reason why you have been forbidden to use those will inform the answers we give. For example is the reason just that you do not have the Biosystems Toolbox ?

[Walter Roberson]

https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html

describes the PDB file format.

Answer 1

[priyanka ilango]

I am doing internship.My professor asked me to use file handling commands.After reading a file I have to extract coordinate and calculate angle between successive calpha atom.

Answer 2

[Walter Roberson]

S = fileread(pdb);
info = regexp(S, '(?<=^ATOM.{27})(?<X>.{8})(?<Y>.{8})(?<Z>.{8})', 'names', 'lineanchors');

If all went well, then info is now a struct array with fields X, Y, and Z, giving the coordinates of each atom as character vectors. struct2cell() and str2double() should then produce numeric arrays.

[Walter Roberson]

S = fileread(pdb);
info = regexp(S, '(?<=^ATOM.{5})(?<AN>.{4})(?<AL>.).{15}(?<X>.{8})(?<Y>.{8})(?<Z>.{8})', 'names', 'lineanchors');

The AN field of the struct will now contain the 4-character right-justified atom "name" information, which is not the same as the element information.

The AL field will now contain the "alternate location" information. Atom "name" C together with Alternate Location 'A' indicate C alpha atoms.

You can find out which entries are CA by using

mask = strcmp({info.AN}, '   C') && strcmp({info.AL}, 'A');
info(mask)

However I still do not know what you want to take the angle with respect to.

[priyanka ilango]

Can you please tell me (? <=^ATOM. {27})that 27 represents what?

[Walter Roberson]

Do not use a space between the period and the {27} : that would mean something different.

https://www.mathworks.com/help/matlab/ref/regexp.html

https://www.mathworks.com/matlabcentral/fileexchange/48930-interactive-regular-expression-tool

In short, the period means "any one character" (which, by default, can include carriage returns and new lines, and any binary character that might happen to be present. any character, once.

The {27} right after that means that the pattern immediately before that is to be matched exactly 27 times. So .{27} means any 27 characters -- even if that would require going to the next line.

(There is an option to regexp, 'dotexceptnewline', which tells regexp that . is not permitted to match the newline at the end of a line. There are cases where that is very useful.)

The ^ATOM.{27 part means the beginning of a line is to be matched, followed by the four literal characters A T O M, followed by any 27 characters. Wrapping (?<= ) around that means to the pattern that follows the ) is to be matched only when it is immediately preceeded by ATOM followed by any 27 characters. This particular construct "positions" you to somewhere in the file that meets the pattern requirements, without holding on to any of the content of what was inside the (?<=) expression.

(?<X>.{8}) at that point means match any one character, exactly 8 times -- so exactly 8 characters -- and that the matched characters are to be put into a field named X in the returned data. We do the same thing (?<Y>.{8}) to grab exactly 8 characters and put them into the Y field, and the (?<Z>.{8}) to grab exactly 8 characters and put them into the Z field. Note these are character vectors, unintepreted, not yet numbers. You can do things like

str2double({info.X})

to ask that those character vectors be converted to numbers.

[Guillaume]

There is an option to regexp, 'dotexceptnewline'

One of the strangest peculiarity of matlab regular expression engine, which tripped me up a few times initially. All other regular expression engines that I know (.Net, Java, Python, Perl, PHP, C++, Ruby) have the opposite default: dot does not match newlines.

[priyanka ilango]

Thank you.That angle calculation is using x y z coordinate values of calpha atom to calculate angle using basic geometry.I can see some formulas for 2 coordinate only.Use three successive calpha atom to calculate like this I have to calculate for all calpha atom.

[priyanka ilango]

There are 602 atoms among them 76 atoms are CAatoms how to take CA atom coordinate values to calculate angle .I saw some formula based on giving x1 y1 z1 X2 y2 z2 x3 y3 z3 then calculate using atan2

[Walter Roberson]

Two-coordinate form:

acos(dot([x2-x1, y2-y2, z2-z1]) ./ sqrt(sum([x2-x1, y2-y1, z2-z1].^2)))

However when you do this, you are implicitly creating a vector between the coordinates and the origin (0,0,0) . The angle you would get out would change if you were to rotate the structure in 3-space.

Unless you choose a coordinate such as (0,0,0) as your implicit starting point, two points in 3-space form a line, not an angle.

[priyanka ilango]

Thank you.

Answer 3

[priyanka ilango]

Thank you so much.I got the coordinate values.How to calculate angle between successive calpha atoms from these values?

[Walter Roberson]

Are the calpha atoms the ones with atom name 'C' and "alternate location" 'A", and which therefore appear to have atom name 'CA' ?

You need at least three points to create an angle, and I do know that what you would be taking the angle relative to.

Answer 4

[priyanka ilango]

CA is calpha atom of a protein.Here I extract coordinates values for all atoms.But I need to calculate angle for all calpha atom using these coordinate values.I am using 1UBQ protein.

[Walter Roberson]

Please click on "Comment on this Answer" to respond to things, instead of creating new Answer each time. Answer are intended for substantial steps towards solutions to the original question, not for back and forth dialog.

Answer 5

[priyanka ilango]

If there is a values in linear line like X= 27.567 24.789 45.879 27.568 26.897... upto 100 values How to extract or take particular values from this?