diff PLIDflow/scripts/pdbcenter_npts_finder.R @ 0:6fcfa4756040 draft

Uploaded
author bitlab
date Tue, 14 Jan 2020 06:09:42 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/PLIDflow/scripts/pdbcenter_npts_finder.R	Tue Jan 14 06:09:42 2020 -0500
@@ -0,0 +1,84 @@
+#pdbcenter_npts_finder.R calculates the geometric center of a PDB file
+
+#!/usr/bin/env Rscript      
+args = commandArgs(trailingOnly=TRUE) 
+
+if(length(args) < 1){
+  stop("USE: Rscript pdbcenter_npts_finder.R <receptor_name> <session_dir>")
+}
+
+#Arguments definition 
+receptor_name <- args[1]
+session_dir <- args[2]
+
+#Read PDB file 
+archivo <- scan(args[1], what = character(), quiet = TRUE)
+
+#Make a table where positions of carbon (C) atoms and their x, y, z coordenates will be placed
+salto <- 0
+for(i in 1:(length(archivo)-2)){
+  if(archivo[i] == "ATOM" && archivo[i+2] == "C"){
+    salto <- salto + 1
+  }
+}
+tablaC <- matrix(1, nrow = salto, ncol = 4)
+
+#Header. Note: first column: C atom position (will be filled in when the number of C will be known)
+#Row labels
+colnames(tablaC) <- c("C Position","X", "Y", "z")
+rownames(tablaC) <- 1:salto
+
+#Fill in tablaC with x, y,z coordenates of the C atoms as numeric
+salto <- 0
+posicionC <-  c()
+for(i in 1:(length(archivo)-2)){
+  if(archivo[i] == "ATOM" && archivo[i+2] == "C"){
+    salto <- salto + 1
+    posicionC <- c(posicionC, archivo[i+1])
+    
+    tablaC[salto,2] <- as.numeric(archivo[i+6])
+    tablaC[salto,3] <- as.numeric(archivo[i+7])
+    tablaC[salto,4] <- as.numeric(archivo[i+8])
+ }
+}
+
+#Add C atom positions in the column 1
+for(i in 1:salto){
+  tablaC[i,1] <- as.numeric(posicionC[i]) 
+}
+
+#Geometric center calculation as grid center x, y, z to be introduced in the GPF file 
+pto_x_medio <- (max(tablaC[,2]) + min(tablaC[,2]))/2
+pto_y_medio <- (max(tablaC[,3]) + min(tablaC[,3]))/2
+pto_z_medio <- (max(tablaC[,4]) + min(tablaC[,4]))/2
+
+pto_x_medio_r <- round(pto_x_medio,digits=3)
+pto_y_medio_r <- round(pto_y_medio,digits=3)
+pto_z_medio_r <- round(pto_z_medio,digits=3)
+
+#print(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, sep = ";"))
+
+#Calculate distances in each coordenate to have de number of points (npts) x, y, z to be introduced in the GPF file 
+eu_x <- max(tablaC[,2])-min(tablaC[,2])
+eu_y <- max(tablaC[,3])-min(tablaC[,3])
+eu_z <- max(tablaC[,4])-min(tablaC[,4])
+
+eu_x_r <- round(eu_x,digits=0)
+eu_y_r <- round(eu_y,digits=0)
+eu_z_r <- round(eu_z,digits=0)
+
+#print(paste(eu_x, eu_y, eu_z, sep = ";"))
+#print(paste(eu_x_r, eu_y_r, eu_z_r, sep = ";"))
+
+if(eu_x > 126 || eu_y > 126 || eu_x > 126){
+  print("At least one dimension in number of points is out of range (>126) according AutoDock software")
+  break
+}
+
+#Write PDB geometric center and eu_x_r, eu_y_r, eu_z_r in a file
+#write(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, eu_x_r, eu_y_r, eu_z_r,sep= "\n"), file = "/home/galaxy/galaxy/tools/proteindocking/scripts/pdbcenter_npts.txt")
+write(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, eu_x_r, eu_y_r, eu_z_r,sep= "\n"), file = "pdbcenter_npts.txt")
+
+
+
+