Mercurial > repos > bitlab > plidflow
comparison PLIDflow/scripts/pdbcenter_npts_finder.R @ 0:6fcfa4756040 draft
Uploaded
| author | bitlab |
|---|---|
| date | Tue, 14 Jan 2020 06:09:42 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:6fcfa4756040 |
|---|---|
| 1 #pdbcenter_npts_finder.R calculates the geometric center of a PDB file | |
| 2 | |
| 3 #!/usr/bin/env Rscript | |
| 4 args = commandArgs(trailingOnly=TRUE) | |
| 5 | |
| 6 if(length(args) < 1){ | |
| 7 stop("USE: Rscript pdbcenter_npts_finder.R <receptor_name> <session_dir>") | |
| 8 } | |
| 9 | |
| 10 #Arguments definition | |
| 11 receptor_name <- args[1] | |
| 12 session_dir <- args[2] | |
| 13 | |
| 14 #Read PDB file | |
| 15 archivo <- scan(args[1], what = character(), quiet = TRUE) | |
| 16 | |
| 17 #Make a table where positions of carbon (C) atoms and their x, y, z coordenates will be placed | |
| 18 salto <- 0 | |
| 19 for(i in 1:(length(archivo)-2)){ | |
| 20 if(archivo[i] == "ATOM" && archivo[i+2] == "C"){ | |
| 21 salto <- salto + 1 | |
| 22 } | |
| 23 } | |
| 24 tablaC <- matrix(1, nrow = salto, ncol = 4) | |
| 25 | |
| 26 #Header. Note: first column: C atom position (will be filled in when the number of C will be known) | |
| 27 #Row labels | |
| 28 colnames(tablaC) <- c("C Position","X", "Y", "z") | |
| 29 rownames(tablaC) <- 1:salto | |
| 30 | |
| 31 #Fill in tablaC with x, y,z coordenates of the C atoms as numeric | |
| 32 salto <- 0 | |
| 33 posicionC <- c() | |
| 34 for(i in 1:(length(archivo)-2)){ | |
| 35 if(archivo[i] == "ATOM" && archivo[i+2] == "C"){ | |
| 36 salto <- salto + 1 | |
| 37 posicionC <- c(posicionC, archivo[i+1]) | |
| 38 | |
| 39 tablaC[salto,2] <- as.numeric(archivo[i+6]) | |
| 40 tablaC[salto,3] <- as.numeric(archivo[i+7]) | |
| 41 tablaC[salto,4] <- as.numeric(archivo[i+8]) | |
| 42 } | |
| 43 } | |
| 44 | |
| 45 #Add C atom positions in the column 1 | |
| 46 for(i in 1:salto){ | |
| 47 tablaC[i,1] <- as.numeric(posicionC[i]) | |
| 48 } | |
| 49 | |
| 50 #Geometric center calculation as grid center x, y, z to be introduced in the GPF file | |
| 51 pto_x_medio <- (max(tablaC[,2]) + min(tablaC[,2]))/2 | |
| 52 pto_y_medio <- (max(tablaC[,3]) + min(tablaC[,3]))/2 | |
| 53 pto_z_medio <- (max(tablaC[,4]) + min(tablaC[,4]))/2 | |
| 54 | |
| 55 pto_x_medio_r <- round(pto_x_medio,digits=3) | |
| 56 pto_y_medio_r <- round(pto_y_medio,digits=3) | |
| 57 pto_z_medio_r <- round(pto_z_medio,digits=3) | |
| 58 | |
| 59 #print(paste(pto_x_medio_r, pto_y_medio_r, pto_z_medio_r, sep = ";")) | |
| 60 | |
| 61 #Calculate distances in each coordenate to have de number of points (npts) x, y, z to be introduced in the GPF file | |
| 62 eu_x <- max(tablaC[,2])-min(tablaC[,2]) | |
| 63 eu_y <- max(tablaC[,3])-min(tablaC[,3]) | |
| 64 eu_z <- max(tablaC[,4])-min(tablaC[,4]) | |
| 65 | |
| 66 eu_x_r <- round(eu_x,digits=0) | |
| 67 eu_y_r <- round(eu_y,digits=0) | |
| 68 eu_z_r <- round(eu_z,digits=0) | |
| 69 | |
| 70 #print(paste(eu_x, eu_y, eu_z, sep = ";")) | |
| 71 #print(paste(eu_x_r, eu_y_r, eu_z_r, sep = ";")) | |
| 72 | |
| 73 if(eu_x > 126 || eu_y > 126 || eu_x > 126){ | |
| 74 print("At least one dimension in number of points is out of range (>126) according AutoDock software") | |
| 75 break | |
| 76 } | |
| 77 | |
| 78 #Write PDB geometric center and eu_x_r, eu_y_r, eu_z_r in a file | |
| 79 #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") | |
| 80 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") | |
| 81 | |
| 82 | |
| 83 | |
| 84 |