Merge branch 'kabsch'

src/KABSCH/runtests.jl

@@ -1,7 +1,8 @@
 include(Pkg.dir("BiomolecularStructures", "src/KABSCH", "kabsch.jl")) 
+include(Pkg.dir("BiomolecularStructures", "src/PDB", "pdb.jl")) 
 using Base.Test
 using Kabsch
-
+using PDB
 # julia uses MATLAB-style syntax for matrices
 #P = [0.1 0.2 0.3; 0.5 -0.5 -0.3; 0.2 -0.3 8.3]
 
@@ -20,4 +21,26 @@ Q = [9 8 7; 6 5 4; 3 2 1]
 
 @test_approx_eq_eps rmsd(P,Q) 8.94427191 1e-4
 
-@test_approx_eq_eps kabsch_rmsd(P,Q) 2.4323767778e-15 1e-6
+@test_approx_eq_eps kabsch_rmsd(P,Q) 2.4323767778e-15 1e-6
+
+structure = get_structure("examples/data/2HHB.pdb")
+chains = get_chains(structure)
+
+alpha_1 = chains[1]
+
+alpha_2 = chains[3]
+
+beta_1 = chains[2]
+beta_2 = chains[4]
+
+println("RMSD for Alpha chains of HUMAN DEOXYHAEMOGLOBIN (PDB ID: 2HHB)")
+println("Naive RMSD:")
+println(rmsd(alpha_1,alpha_2))
+println("RMSD after Kabsch:")
+println(kabsch_rmsd(alpha_1,alpha_2))
+
+println("RMSD for Beta chains of HUMAN DEOXYHAEMOGLOBIN (PDB ID: 2HHB)")
+println("Naive RMSD:")
+println(rmsd(beta_1,beta_2))
+println("RMSD after Kabsch:")
+println(kabsch_rmsd(beta_1,beta_1))

src/PDB/pdb.jl

@@ -3,13 +3,57 @@
 # sudo pip install biopython
 # sudo apt-get install numpy | unless you want to compile it..
 # julia: Pkg.add("PyCall")
+module PDB
+export get_structure, get_chains, structure_to_matrix
+
 using PyCall
 @pyimport Bio.PDB as pdb
 
-pdbparser = pdb.PDBParser()
+# get a structure from a PDB File
+function get_structure(filename::String)
+	pdbparser = pdb.PDBParser(PERMISSIVE = 1)
+
+	# parse the file
+	structure = pdbparser[:get_structure](filename, filename)
+	return structure
+end
+
+# get chains from structure
+function get_chains(structure)
+	chains = structure[:get_chains]()
+
+	chainMatrices = Any[]
+
+	for c in chains
+		push!(chainMatrices,structure_to_matrix(c))
+	end
+
+	return chainMatrices
+end
+
+# Read in a PDB file and return a matrix of C_Alpha atom coordinates
+function structure_to_matrix(structure)
+
+
+	atoms = structure[:get_atoms]()
+	# Filter out C_Alpha atoms
+	atoms = filter(x -> x[:get_name]() == "CA", atoms) 
+	# Get the coordinates
+	atoms = map(x -> x[:get_coord](), atoms)
+
+	n = length(atoms)
+
+	# initialize an empty Nx3 matrix
+	matrix = zeros(n,3)
+
+	# fill it
+	for i in 1:n
+		matrix[i,:] = [atoms[i][1] atoms[i][2] atoms[i][3]]
+	end
+
+	matrix = reshape(matrix,n,3)
 
-structure = pdbparser[:get_structure]("test", "test.pdb")
+	return matrix
+end
 
-for atom in structure[:get_atoms]()
-	println(atom[:get_coord]())
 end