using System; using System.IO; public class dihedrals { //Extract Phi, Psi angles from the first chain of the pdb file //V.V. 02/25/2011. vvostri@gmail.com public static void Main() { string filename; do { Console.WriteLine("Enter pdb file name"); filename = Console.ReadLine(); int AAnumber = 0; if (File.Exists(filename)) { StreamReader streamvariable = File.OpenText(filename); string lineinfile; int Nmodels = 0; while ((lineinfile = streamvariable.ReadLine()) != null) { try { string temp = lineinfile.Substring(0, 6); } catch { lineinfile = "XXXXXXXXXXXXXXXXXXXXXXXX"; } if (lineinfile.Substring(0, 6) == "MODEL ") { Nmodels = Nmodels + 1; } } streamvariable.Close(); if (Nmodels == 0) { Nmodels = 1; } //COLUMNS DATA TYPE CONTENTS //-------------------------------------------------------------------------------- // 1 - 6 Record name "ATOM " // 7 - 11 Integer Atom serial number. //13 - 16 Atom Atom name. //17 Character Alternate location indicator. //18 - 20 Residue name Residue name. //22 Character Chain identifier. //23 - 26 Integer Residue sequence number. //27 AChar Code for insertion of residues. //31 - 38 Real(8.3) Orthogonal coordinates for X in Angstroms. //39 - 46 Real(8.3) Orthogonal coordinates for Y in Angstroms. //47 - 54 Real(8.3) Orthogonal coordinates for Z in Angstroms. //55 - 60 Real(6.2) Occupancy. //61 - 66 Real(6.2) Temperature factor (Default = 0.0). //73 - 76 LString(4) Segment identifier, left-justified. //77 - 78 LString(2) Element symbol, right-justified. //79 - 80 LString(2) Charge on the atom. // 1 2 3 4 5 6 7 8 //12345678901234567890123456789012345678901234567890123456789012345678901234567890 //ATOM 145 N VAL A 25 32.433 16.336 57.540 1.00 11.92 A1 N string chainID = null; int current_residue = -1; streamvariable = File.OpenText(filename); lineinfile = streamvariable.ReadLine(); while ((lineinfile.Substring(0, 3) != "END") & (lineinfile.Substring(0, 3) != "TER")) { if ((lineinfile.Substring(0, 4) == "ATOM" || lineinfile.Substring(0, 6) == "HETATM") && Convert.ToInt16(lineinfile.Substring(23, 3)) != current_residue) { current_residue = Convert.ToInt16(lineinfile.Substring(23, 3)); AAnumber = AAnumber + 1; chainID = lineinfile.Substring(21, 1); } lineinfile = streamvariable.ReadLine(); } streamvariable.Close(); string[,] sequence = new string[2, AAnumber]; streamvariable = File.OpenText(filename); lineinfile = streamvariable.ReadLine(); int counter = 0; current_residue = -1; while ((lineinfile.Substring(0, 3) != "END") && (lineinfile.Substring(0, 3) != "TER")) { if ((lineinfile.Substring(0, 4) == "ATOM" || lineinfile.Substring(0, 6) == "HETATM") && Convert.ToInt16(lineinfile.Substring(23, 3)) != current_residue) { //Use "counter" rather than "current_residue" in case the numbering is not sequential sequence[0, counter] = lineinfile.Substring(23, 3); //AA sequence position sequence[1, counter] = lineinfile.Substring(17, 3); //AA name Console.WriteLine("{0}\t{1}", sequence[0, counter], sequence[1, counter]); counter++; current_residue = Convert.ToInt16(lineinfile.Substring(23, 3)); } lineinfile = streamvariable.ReadLine(); } streamvariable.Close(); double[, , ,] coordinates = new double[Nmodels, AAnumber, 3, 3]; streamvariable = File.OpenText(filename); lineinfile = streamvariable.ReadLine(); counter = 0; int m = 0; while (lineinfile != null && m < Nmodels) { if (lineinfile.Substring(0, 3) == "END") { m++; counter = 0; } else if (lineinfile.Substring(0, 3) == "TER") { counter = 0; } else { if (((lineinfile.Substring(0, 4) == "ATOM" || lineinfile.Substring(0, 6) == "HETATM")) && lineinfile.Substring(21, 1) == chainID) { current_residue = Convert.ToInt16(lineinfile.Substring(23, 3)); if (counter < 0.5 * sequence.Length && lineinfile.Substring(23, 3) == sequence[0, counter] && lineinfile.Substring(17, 3) == sequence[1, counter]) { if (lineinfile.Substring(13, 2) == "CA") { coordinates[m, counter, 0, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 0, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 0, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } if (lineinfile.Substring(13, 2) == "N ") { coordinates[m, counter, 1, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 1, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 1, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } if (lineinfile.Substring(13, 2) == "C ") { coordinates[m, counter, 2, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 2, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 2, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } } else { counter++; if (lineinfile.Substring(13, 2) == "CA") { coordinates[m, counter, 0, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 0, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 0, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } if (lineinfile.Substring(13, 2) == "N ") { coordinates[m, counter, 1, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 1, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 1, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } if (lineinfile.Substring(13, 2) == "C ") { coordinates[m, counter, 2, 0] = Convert.ToDouble(lineinfile.Substring(30, 8)); coordinates[m, counter, 2, 1] = Convert.ToDouble(lineinfile.Substring(38, 8)); coordinates[m, counter, 2, 2] = Convert.ToDouble(lineinfile.Substring(46, 8)); } //Atom coordinates: //HA //HN //CA //CB //N //C } } } lineinfile = streamvariable.ReadLine(); } streamvariable.Close(); double[, ,] torsions = new double[Nmodels, AAnumber, 2]; int counter1 = 0; for (counter1 = 0; counter1 < Nmodels; counter1++) { for (counter = 0; counter < AAnumber; counter++) { //Atom coordinates: //0 CA //1 N //2 C //Phi: C-N-CA-C //Psi: N-CA-C-N //AAi - AAj - AAk //General method to get a torsion angle between A-B-C-D //1. Calculate BA, CD, CB bond vectors. //2. Calculate cross products n1 = BA x CB, n2 = CD x CB //3. Calculate dot product cos = n1 . n2 / (|n1| * |n2|) //4. Calculate matrix determinant D = det(n2, n1, CB) //5. Calculate sin = D / (|n2| * |n1| * |CB|) //6. Calculate the torsion tor = atan2(cos, sin) * 180 / Pi //Bond vectors for Phi double[] Nj_Ci = new double[3]; double[] CAj_Cj = new double[3]; double[] CAj_Nj = new double[3]; //Bond vectors for Psi //CAj_Nj - already declared double[] Cj_Nk = new double[3]; double[] Cj_CAj = new double[3]; //Initialize bond vectors and normalize where necessary try { Nj_Ci[0] = coordinates[counter1, counter - 1, 2, 0] - coordinates[counter1, counter, 1, 0]; Nj_Ci[1] = coordinates[counter1, counter - 1, 2, 1] - coordinates[counter1, counter, 1, 1]; Nj_Ci[2] = coordinates[counter1, counter - 1, 2, 2] - coordinates[counter1, counter, 1, 2]; } catch { Nj_Ci[0] = 0; Nj_Ci[1] = 0; Nj_Ci[2] = 0; } CAj_Cj[0] = coordinates[counter1, counter, 2, 0] - coordinates[counter1, counter, 0, 0]; CAj_Cj[1] = coordinates[counter1, counter, 2, 1] - coordinates[counter1, counter, 0, 1]; CAj_Cj[2] = coordinates[counter1, counter, 2, 2] - coordinates[counter1, counter, 0, 2]; double norm_CAj_Cj = norm(CAj_Cj); CAj_Nj[0] = coordinates[counter1, counter, 1, 0] - coordinates[counter1, counter, 0, 0]; CAj_Nj[1] = coordinates[counter1, counter, 1, 1] - coordinates[counter1, counter, 0, 1]; CAj_Nj[2] = coordinates[counter1, counter, 1, 2] - coordinates[counter1, counter, 0, 2]; double norm_CAj_Nj = norm(CAj_Nj); try { Cj_Nk[0] = coordinates[counter1, counter + 1, 1, 0] - coordinates[counter1, counter, 2, 0]; Cj_Nk[1] = coordinates[counter1, counter + 1, 1, 1] - coordinates[counter1, counter, 2, 1]; Cj_Nk[2] = coordinates[counter1, counter + 1, 1, 2] - coordinates[counter1, counter, 2, 2]; double norm_Cj_Nk = norm(Cj_Nk); } catch { Cj_Nk[0] = 0; Cj_Nk[1] = 0; Cj_Nk[2] = 0; double norm_Cj_Nk = 0; } Cj_CAj[0] = coordinates[counter1, counter, 0, 0] - coordinates[counter1, counter, 2, 0]; Cj_CAj[1] = coordinates[counter1, counter, 0, 1] - coordinates[counter1, counter, 2, 1]; Cj_CAj[2] = coordinates[counter1, counter, 0, 2] - coordinates[counter1, counter, 2, 2]; double norm_Cj_CAj = norm(Cj_CAj); //Phi: n1, n2 double[] NjCi_x_CAjNj = crossproduct(Nj_Ci, CAj_Nj); double norm_NjCi_x_CAjNj = norm(NjCi_x_CAjNj); double[] CAjCj_x_CAjNj = crossproduct(CAj_Cj, CAj_Nj); double norm_CAjCj_x_CAjNj = norm(CAjCj_x_CAjNj); //Psi: n1, n2 double[] CAjNj_x_CjCAj = crossproduct(CAj_Nj, Cj_CAj); double norm_CAjNj_x_CjCAj = norm(CAjNj_x_CjCAj); double[] CjNk_x_CjCAj = crossproduct(Cj_Nk, Cj_CAj); double norm_CjNk_x_CjCAj = norm(CjNk_x_CjCAj); double cos1 = dotproduct(NjCi_x_CAjNj, CAjCj_x_CAjNj) / (norm_NjCi_x_CAjNj * norm_CAjCj_x_CAjNj); double cos2 = dotproduct(CAjNj_x_CjCAj, CjNk_x_CjCAj) / (norm_CAjNj_x_CjCAj * norm_CjNk_x_CjCAj); double sin1 = mdet(CAjCj_x_CAjNj, NjCi_x_CAjNj, CAj_Nj) / (norm_CAjCj_x_CAjNj * norm_NjCi_x_CAjNj * norm_CAj_Nj); double sin2 = mdet(CjNk_x_CjCAj, CAjNj_x_CjCAj, Cj_CAj) / (norm_CjNk_x_CjCAj * norm_CAjNj_x_CjCAj * norm_Cj_CAj); double phi = Math.Atan2(sin1, cos1) * 180 / Math.PI; double psi = Math.Atan2(sin2, cos2) * 180 / Math.PI; if (phi.Equals(double.NaN) == true || psi.Equals(double.NaN) == true) { phi = 0; psi = 0; } torsions[counter1, counter, 0] = phi; torsions[counter1, counter, 1] = psi; } } Console.WriteLine("{0} models, {1} AAs", Nmodels, AAnumber); Console.ReadKey(); FileStream file = new FileStream("torsions.txt", FileMode.Create, FileAccess.Write); StreamWriter writer = new StreamWriter(file); for (counter1 = 0; counter1 < Nmodels; counter1++) { for (counter = 0; counter < AAnumber; counter++) { writer.Write(torsions[counter1, counter, 0]); writer.Write(","); writer.Write(torsions[counter1, counter, 1]); writer.Write("\r\n"); } } writer.Close(); file.Close(); } else { Console.WriteLine("File {0} not found\n", filename); } } while (filename != "quit"); } public static void angle(ref double x, ref double y, ref double z, ref double cosine) { cosine = z / (Math.Pow(Math.Pow(x,2) + Math.Pow(y,2) + Math.Pow(z,2),0.5)); } public static double norm(double[] inputarray) { double temp = (Math.Sqrt(Math.Pow(inputarray[0], 2) + Math.Pow(inputarray[1], 2) + Math.Pow(inputarray[2], 2))); return temp; } public static double[] crossproduct(double[] array1, double[] array2) { double[] cp = new double[3]; cp[0] = array1[1] * array2[2] - array1[2] * array2[1]; cp[1] = -(array1[0] * array2[2] - array1[2] * array2[0]); cp[2] = array1[0] * array2[1] - array1[1] * array2[0]; return cp; } public static double dotproduct(double[] array1, double[] array2) { double dp = array1[0] * array2[0] + array1[1] * array2[1] + array1[2] * array2[2]; return dp; } public static double mdet(double[] array1, double[] array2, double[] array3) { double[] cp = crossproduct(array2,array3); double det = array1[0] * cp[0] + array1[1] * cp[1] + array1[2] * cp[2]; return det; } }