bio3d-package {bio3d}R Documentation

Biological Structure Analysis

Description

Utilities for the analysis of protein structure and sequence data.

Details

Package: bio3d
Type: Package
Version: 1.0-6
Date: 2009-07-07
License: GPL version 2 or newer
URL: http://mccammon.ucsd.edu/~bgrant/bio3d/

Features include the ability to read and write structure (read.pdb, write.pdb, read.fasta.pdb), sequence (read.fasta, write.fasta) and dynamics trajectory data (read.dcd).

Perform atom summaries (pdb.summary), atom selection (atom.select), re-orientation (orient.pdb), superposition (rot.lsq, fit.xyz), rigid core identification (core.find, plot.core, fit.xyz), torsion/dihedral analysis (torsion.pdb, torsion.xyz), clustering (via hclust) and principal component analysis (pca.xyz, pca.tor, plot.pca, plot.pca.loadings, mktrj.pca) of structure data.

Perform conservation analysis of sequence (conserv, identity, entropy, consensus) and structural (rmsd, core.find) data.

In addition, various utility functions are provided to facilitate manipulation and analysis of biological sequence and structural data (e.g. aa123, aa321, aln2html, atom.select, rot.lsq, fit.xyz, rmsd, is.gap, gap.inspect, orient.pdb and pairwise).

Note

The latest version and further documentation can be obtained from the bio3d website:
http://mccammon.ucsd.edu/~bgrant/bio3d/.

Author(s)

Barry Grant <bgrant@mccammon.ucsd.edu>

References

Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696.

Examples

## Not run: 
lbio3d()  # list the functions within the package

## See the individual functions for further documentation and examples.
#help(read.pdb)

##-- Read a PDB file
pdb <- read.pdb(system.file("examples/1bg2.pdb", package="bio3d"))

##-- Distance matrix
k <- dm(pdb, selection="calpha")
plot(k)

## Extract SEQRES PDB sequence
s <- aa321(pdb$seqres)

## write a FASTA format sequence file
write.fasta(seqs=s, id="seqres", file="eg.fa")

## Reorient and write a new PDB file
xyz <- orient.pdb(pdb)
write.pdb(pdb, xyz = xyz, file = "mov1.pdb")

##-- Torsion angle analysis and basic Ramachandran plot
tor <- torsion.pdb(pdb)
plot(tor$phi, tor$psi)

##-- Secondary structure assignment with DSSP
sse <- dssp(pdb, resno=FALSE)

## Plot of B-factor values along with secondary structure
plot.bio3d(pdb$atom[pdb$calpha,"b"], sse=sse, ylab="B-factor")



##-- Read a FASTA alignment file
aln <- read.fasta( system.file("examples/kinesin_xray.fa",package="bio3d") )

## Entropy score for alignment positions
h <- entropy(aln)    # see also the "conserv()" function

# Alignment consensus
con <- consensus(aln)

## Plot of residue conservation (similarity) with secondary structure
plot.bio3d( conserv(aln)[!is.gap(aln$ali[1,])],  sse=sse )

## Render the alignment as coloured HTML
aln2html(aln, append=FALSE, file="eg.html")


##-- Read an alignment of sequences and their corresponding structures
aln  <- read.fasta( system.file("examples/kif1a.fa", package="bio3d") )
pdbs <- read.fasta.pdb( aln, pdbext = ".ent",
               pdb.path = system.file("examples",package="bio3d") )

##-- DDM: Difference Distance Matrix
a <- dm(pdbs$xyz[2,])
b <- dm(pdbs$xyz[3,])
ddm <- a - b
plot(ddm,key=FALSE, grid=FALSE)

##-- Superpose structures on non gap positions
gaps <- unique(which(is.na(pdbs$xyz) , arr.ind=TRUE)[,2])
gaps <- gap.inspect(pdbs$xyz)

xyz  <- fit.xyz( fixed  = pdbs$xyz[1,],
                mobile = pdbs$xyz,
                fixed.inds  = gaps$f.inds,
                mobile.inds = gaps$f.inds )

##-- RMSD
rmsd(pdbs$xyz[1,gaps$f.inds], pdbs$xyz[,gaps$f.inds])
rmsd(pdbs$xyz[1,gaps$f.inds], pdbs$xyz[,gaps$f.inds], fit=TRUE)


##-- Rigid 'core' identification
aln <- read.fasta(system.file("examples/kinesin_xray.fa",package="bio3d"))
pdb.path = system.file("examples/",package="bio3d")
pdbs <- read.fasta.pdb(aln, pdb.path = pdb.path, pdbext = ".ent")

core <- core.find(pdbs,
                  #write.pdbs=TRUE,
                  rm.island=TRUE,
                  verbose=TRUE)

## Plot core volume vs length
plot(core)

## Core fit the structures (superpose on rigid zones)
xyz <- fit.xyz( fixed = pdbs$xyz[1,],
                mobile = pdbs,
              #  full.pdbs = TRUE,
              #  pdb.path = pdb.path,
              #  pdbext = ".ent",
              #  outpath = "fitlsq/",
                fixed.inds  = core$c0.5A.xyz,
                mobile.inds = core$c0.5A.xyz)

##-- PCA of structures
cut.seqs <- which(pdbs$id %in% c("d1n6mb_","d1ry6a_"))
# Ignore gap containing positions
gaps.res <- gap.inspect(pdbs$ali[-cut.seqs,])
gaps.pos <- gap.inspect(pdbs$xyz[-cut.seqs,])

##-- Do PCA
pc.xray <- pca.xyz(xyz[-cut.seqs, gaps.pos$f.inds])

## Plot results
plot(pc.xray)

#x11()
plot.pca.loadings(pc.xray$U)

## Write PC trajectory
a <- mktrj.pca(pc.xray, pc=1, file="pc1.pdb",
               resno = pdbs$resno[1, gaps.res$f.inds],
               resid = aa123(pdbs$ali[1, gaps.res$f.inds]) )


##-- Read a CHARMM/X-PLOR/NAMD trajectory file
trtfile <- system.file("examples/hivp.dcd", package="bio3d")
trj <- read.dcd(trtfile)

## Read the starting PDB file to determine atom correspondence
pdbfile <- system.file("examples/hivp.pdb", package="bio3d")
pdb <- read.pdb(pdbfile)

## Fit trj on PDB based on residues 23 to 31 and 84 to 87 in both chains
inds <- atom.select(pdb,"///23:31,84:87///CA/")
fit.xyz <- fit.xyz(pdb$xyz, trj, fixed.inds=inds$xyz, mobile.inds=inds$xyz)

##-- RMSD of trj frames from PDB
r <- rmsd(a=pdb, b=fit.xyz)

##-- PCA of trj
pc.trj <- pca.xyz(fit.xyz)

plot.pca.loadings(pc.trj)

## Write PC trajectory for viewing in VMD
a <- mktrj.pca(pc.trj, pc=1, file="pc1_trj.pdb")

## End(Not run)
## other examples include: sdm, alignment, clustering etc...


[Package bio3d version 1.0-6 Index]