Delete residues with vmd atomselect command

Question

I am trying to use the atomselect command to delete lipids within an area of proteins in my .gro file. Since it is a martini coarse grained file I can just use the keywords resname for residuenames, and name or type for bead-types (my pseudoatoms). So the default singlewords are not defined.

I tried it with the command:

atomselect 0 "all not resname DPPE DOPE POPE POPG within 1 of resname ALA ARG ASN ASP CYS GLN GLU GLY HIS ILE LEU LYS MET PHE PRO SER THR TRP TYR VAL"

and got the following error:

ERROR) Selection terminated too early

ERROR) syntax error atomselect: cannot parse selection text: all not resname DPPE DOPE POPE POPG within 1 of resname ALA ARG ASN ASP CYS GLN GLU GLY HIS ILE LEU LYS MET PHE PRO SER THR TRP TYR VAL

so apparently I didn't get the syntax right. I tried a few different versions and to store the resname selections in variables, but nothing worked. How could I solve this problem?

Is there a line break after ASP in the code? Tcl has no problems with line breaks within quotes, but atomselect might have. — Peter Lewerin
I found another solution instead: I use the bio3d package in R, helps a lot! If I will try it again with vmd for integrity reasons I will comment back :) — gwendolinese

gwendolinese gwendolinese · Accepted Answer · 2015-09-15T14:54:58

So, I wrote a little code in R using bio3d to do the task:

library(bio3d)
library (gdata)

#read pdb input
insane <- read.pdb("nice.pdb")


#vector with residuenumbers to delete 
select_resi <- c(1000, 1136, 1026, 1252, 1449, 970, 1067, 1298, 
                 1287, 1357, 1051, 993, 1241, 1282, 1341, 1344, 
                 1048, 1154, 1205, 1274, 1465, 1322, 1418, 992)

select_resi <- sort(select_resi)

#select residues to delete
selected_resi <-  atom.select(insane, resno = select_resi)


#delete selected residues
insane$atom <- insane$atom[-selected_resi$atom,]
insane$calpha <- insane$calpha[-selected_resi$atom]
insane$xyz <- insane$xyz[-selected_resi$xyz]

#renumber residuenumbers and convert to gromacs type pdb
printstuff <- convert.pdb(insane, type = "gromacs", 
                          renumber = TRUE, first.resno = 1, first.eleno = 1)

#write pdb file
write.pdb(pdb = printstuff, file = "memb2R.pdb", xyz = printstuff$xyz, 
          type = printstuff$atom$type, resno = printstuff$atom$resno, 
          resid = printstuff$atom$resid, eleno = printstuff$atom$eleno, 
          elety = printstuff$atom$elety, end = TRUE, verbose = TRUE)

if you need the .gro format afterwards (that's why gdata is loaded), you can continue with this:

#preparing xyz to match with .gro format
xyz_vector <- insane$xyz/10 

pos_x <- c()
pos_y <- c()
pos_z <- c()

#sorting loop for x, y, z coordinates
for (i in seq(along=xyz_vector)) {
  pos_x <- c(pos_x, xyz_vector[i])
  xyz_vector <- xyz_vector[-i]
    pos_y <- c(pos_y, xyz_vector[i])
    xyz_vector <- xyz_vector[-i]
      pos_z <- c(pos_z, xyz_vector[i])
      print(i)
}

#delete redundant entries
pos_x <- pos_x[1:length(xyz_vector)]
pos_y <- pos_y[1:length(xyz_vector)]
pos_z <- pos_z[1:length(xyz_vector)]

#prepare other passing vectors for .gro vector
resi_numb <- insane$atom$resno
resi_name <- insane$atom$resid
atom_name <- insane$atom$elety
atom_numb <- insane$atom$eleno


#prepare .gro vector
gro_vec <- sprintf ("%5d%-5s%5s%5d%8.3f%8.3f%8.3f",

                      #for velocity fields add %8.4f%8.4f%8.4f to string and create 3 velocity vectors similar to pos_x etc.

                      resi_numb, resi_name, atom_name, atom_numb,
                      pos_x, pos_y, pos_z)

#transform gro_vec in matrix
gro_matrix <- as.matrix(gro_vec, byrow = TRUE)


#write output
write.fwf (gro_matrix, file = "yourfile.txt", sep = "")

remember to rename the .txt file in .gro. Write in the first line of your file the systemname, the second line the atom entity and in the last line the box vectors, if available. Not the most beautiful piece of code, but I'm still new and it gets the job done.
Improvement suggestions are always welcome :)

Delete residues with vmd atomselect command

1 Answers