Matlab Tutorials | Examples
Practice 10:
Files: Amino acid bonds
Proteins are made up of amino acids. Every amino acid has one amino group (N), one carboxyl group (C) and one chiral carbon (Cα or CA). The angle between CA-C is called “psi” (ψ) and the angle between N-CA is called “phi” (Φ). Phi and Psi angles range between -180o and 180o.
Given the 3-D space coordinates of N, C, and CA atoms in the from P= [x y z], the following quantities of the protein can be computed:
The operations of the form VT. VU in the above angle formulas represent “dotproducts”, which can be computed by using dot(VT, VU) function in MatLab. You are given a file named “protein.txt”. This file contains the x-y-z coordinates of protein N-CA-C atoms. Our protein is made up of 10 amino acids; so there are 10×3=30 coordinate triplets. The format of the file is as follows:
Write a program that will do the following tasks:
- Get the names of the input and an output files (and their folders) from the user.
- Read atom coordinates into three cell arrays, PN, PCA, and PC. Each cell should contain the x, y, and z coordinates of an atom.
- Calculate the bond vectors for each N-CA, CA-C, and C-N pairs. Determine the average length of each bond type, and write this on the screen, and into the output file.
- Calculate the phi and psi angles (in degrees) for amino acids of the protein. Write this information in tabular form on the screen, and into the output file.
Files: https://drive.google.com/open?id=1roJkdZMWGspbhUgFXpJxjwq74p5a-ii9
Solutions:
|
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 |
inFile= input('Enter input file name: ','s'); Pca = []; Pn = []; Pc = []; Vnca = []; Vcac = []; Vcn = []; [fin,err]= fopen(inFile,'r'); if fin==-1 disp('Could not open input file.'); disp(err); else fout= input('Enter odd output file name: ','s'); [fout,err]= fopen(fout,'w'); if fout==-1 disp('Could not open odd output file.'); disp(err); else while ~feof(fin) tline = fgetl(fin); a = str2num(tline(3:length(tline))); if tline (1:2) == 'CA' disp('CA') Pca = [Pca;a]; elseif tline (1:2) == 'C ' disp('C') Pc = [Pc;a]; elseif tline (1:2) == 'N ' disp('N') Pn = [Pn;a]; end disp(tline) end Vnca = [Vnca (Pca(1:10,:) - Pn(1:10,:))]; Vcac = [Vcac (Pc(1:10,:) - Pca(1:10,:))]; Vcn = [Vcn (Pn(2:10,:) - Pc(1:9,:))]; fprintf('The average bond length for N-Ca bond is %.2f Angstrom\nThe average bond length for Ca-C bond is %.2f Angstrom\nThe average bond length for C-N bond is %.2f Angstrom\n\nPsi angles Phi angles\n---------- ----------\n', mean(mean(Vnca)), mean(mean(Vcac)), mean(mean(Vcn))); fprintf(fout, 'The average bond length for N-Ca bond is %.2f Angstrom\nThe average bond length for Ca-C bond is %.2f Angstrom\nThe average bond length for C-N bond is %.2f Angstrom\n\nPsi angles Phi angles\n---------- ----------\n', mean(mean(Vnca)), mean(mean(Vcac)), mean(mean(Vcn))); for i=1:length(Vcn) dotproduct=dot(Vnca(i),Vcn(i)); dotproduct2=dot(Vcn(i),Vcac(i+1)); X = (Vnca(i,1)^2); Y = (Vnca(i,2)^2); Z = (Vnca(i,3)^2); X1 = (Vcn(i,1)^2); Y1 = (Vcn(i,2)^2); Z1 = (Vcn(i,3)^2); X2 = (Vcac(i,1)^2); Y2 = (Vcac(i,2)^2); Z2 = (Vcac(i,3)^2); f=sqrt(X+Y+Z); f2=sqrt(X1+Y1+Z1); f3=sqrt((X1+Y1+Z1)/(X2+Y2+Z2)); output=radtodeg(acos(dotproduct/(f*f2))); output2=radtodeg(acos(dotproduct2/(f3))); fprintf('%10.2f %10.2f\n', output ,output2); fprintf(fout, '%10.2f %10.2f\n',output ,output2); end fprintf('\n'); fclose(fout); end fclose(fin); end |