Side chain selection only for b-factor analysis

[Former user]

Hi there

I would like to make a b-factor plot of the side-chain atoms only, what is the best way to do this? Also can I label specific residues in the plot involved in the binding?

Thanks a lot

[Lars Skjærven]

Absolutely. Take a look at functions plotb3 and atom.select. There are some examples that might be relevant for you. Here is one way:

# read PDB
pdb <- read.pdb("1rx2")

# select side chain atoms
sele <- atom.select(pdb, "side")
b <- pdb$atom$b[ sele$atom ]

# First plot - side chain b-factors
plotb3(b)

# atom labels
labs <- paste(pdb$atom$resid, pdb$atom$resno, pdb$atom$elety, sep="-")
labs <- labs[sele$atom]

# identify enables labeling selected data points (right click on plot to exit identification mode)
identify(1:length(labs), b, labs)

# highlight binding site residues
bs <- binding.site(pdb, a.inds=atom.select(pdb, "protein"), b.inds=atom.select(pdb, resid="FOL"))
b <- pdb$atom$b

# make a vector containing data points to highlight in plot
hl <- rep(NA, length(b))
hl[bs$inds$atom] <-  pdb$atom$b[bs$inds$atom]

# all atoms
plot(b, type="h")
lines(hl, col=2)

# side-chain atoms only
plot(b[sele$atom], type="h")
lines(hl[sele$atom], col=2)

Hope it helps!

[chdo]

Hi Lars Thanks a lot for the quick and helpful answer.. The example works very nicely, the labeling option is great!

I'm new to R, so I have to get familiar For 4wnu, its quite large structure, but I only need chain A (I trim the pdb after readin) When I look at the side-chain plot only I would like to improve it still bit

Could it start at the number of the residue in the pdb instead of 1, and with b-factor average value of the side chains? Also the secondary elements would be nice to include.

Last question how is are the binding residues calculated in bio3d? I got an error when I tried to calculate bs... Error in colSums((a[i, ] - t(b))^2) : 'x' must be an array of at least two dimensions

[Lars Skjærven]

Here is one option- but using all atoms

# read and trim protein
pdb <- read.pdb("1rx2")
pdb <- trim(pdb, "protein", resid="FOL", operator="OR")
sele <- atom.select(pdb, "protein")
b <- pdb$atom$b[sele$atom]

# binding site 
bs <- binding.site(pdb, a.inds=sele, b.inds=atom.select(pdb, resid="FOL"))
bs.inds <- combine.select(bs$inds, atom.select(pdb, "calpha"))

# calcualte mean b-values
unq <-  paste(pdb$atom$resid[sele$atom], pdb$atom$resno[sele$atom], sep="-")
b.mean <- sapply(unique(unq), 
                            function(x) { 
                                mean(b[unq==x])
                            })

# make a vector containing data points to highlight in plot
hl <- rep(NA, length(b.mean))
names(hl) <- unique(unq)
hl[ unq[bs.inds$atom] ] <-  b.mean[unq[bs.inds$atom]]

# plot mean b-values
sse <- dssp(trim(pdb, sele))
plotb3(b.mean, sse=sse)
points(hl, pch=16, col=2)

The binding site is defined from the distance between the protein and ligand. Make sure you make the correct selection of ligand:

pdb <- read.pdb("4wnu")
pdb <- trim(pdb, chain="A")
pdb <- trim(pdb, "protein", resid="HEM", operator="OR")
bs <- binding.site(pdb)

[chdo]

All clear, great thanks a lot for the option and explanation!!