Package 'bio3d'

Description

Utilities for the analysis of protein structure and sequence data.

Details

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, read.ncdf, write.ncdf).

Perform sequence and structure database searches (blast.pdb, hmmer), atom summaries (summary.pdb), atom selection (atom.select), alignment (pdbaln, seqaln, mustang) superposition (rot.lsq, fit.xyz), pdbfit), rigid core identification (core.find, plot.core, fit.xyz), dynamic domain analysis (geostas), torsion/dihedral analysis (torsion.pdb, torsion.xyz), clustering (via hclust), principal component analysis (pca.xyz, pca.pdbs, pca.tor, plot.pca, plot.pca.loadings, mktrj.pca), dynamical cross-correlation analysis (dccm, plot.dccm) and correlation network analysis (cna, plot.cna, cnapath) of structure data.

Perform conservation analysis of sequence (seqaln, conserv, seqidentity, entropy, consensus) and structural (pdbaln, rmsd, rmsf, core.find) data.

Perform normal mode analysis (nma, build.hessian), ensemble normal mode analysis (nma.pdbs), mode comparison (rmsip) and (overlap), atomic fluctuation prediction (fluct.nma), cross-correlation analysis (dccm.nma), cross-correlation visualization (pymol.dccm), deformation analysis (deformation.nma), and mode visualization (pymol.modes, mktrj.nma).

In addition, various utility functions are provided to facilitate manipulation and analysis of biological sequence and structural data (e.g. get.pdb, get.seq, aa123, aa321, pdbseq, aln2html, atom.select, rot.lsq, fit.xyz, is.gap, gap.inspect, orient.pdb, pairwise, plot.bio3d, plot.nma, plot.blast, biounit, etc.).

Note

The latest version, package vignettes and documentation with worked example outputs can be obtained from the bio3d website:
http://thegrantlab.org/bio3d/.
http://thegrantlab.org/bio3d/reference/.
https://bitbucket.org/Grantlab/bio3d/.

Author(s)

Barry Grant <[email protected]> Xin-Qiu Yao <[email protected]> Lars Skjaerven <[email protected]> Julien Ide <[email protected]>

References

Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696. Skjaerven, L. et al. (2014) BMC Bioinformatics 15, 399. Grant, B.J. et al. (2021) Protein Science 30, 20–30.

Examples

help(package="bio3d")     # list the functions within the package
#lbio3d()                  # list bio3d function names only

## Or visit:
##   http://thegrantlab.org/bio3d/reference/

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

## Or online:
##    http://thegrantlab.org/bio3d/reference/read.pdb.html

## Not run: 
##-- See the list of Bio3D demos
demo(package="bio3d")

## Try some out, e.g:
demo(pdb) # PDB Reading, Manipulation, Searching and Alignment
demo(pca) # Principal Component Analysis
demo(md)  # Molecular Dynamics Trajectory Analysis
demo(nma) # Normal Mode Analysis

## See package vignettes and tutorals online:
##   http://thegrantlab.org/bio3d/articles/

## End(Not run)

Description

A collection of published indices, or scales, of numerous physicochemical and biological properties of the 20 standard aminoacids (Release 9.1, August 2006).

Usage

data(aa.index)

Format

A list of 544 named indeces each with the following components:

1. H character vector: Accession number.

2. D character vector: Data description.

3. R character vector: LITDB entry number.

4. A character vector: Author(s).

5. T character vector: Title of the article.

6. J character vector: Journal reference.

7. C named numeric vector: Correlation coefficients of similar indeces (with coefficients of 0.8/-0.8 or more/less). The correlation coefficient is calculated with zeros filled for missing values.

8. I named numeric vector: Amino acid index data.

Source

‘AAIndex’ was obtained from:
https://www.genome.jp/aaindex/
For a description of the ‘AAindex’ database see:
https://www.genome.jp/aaindex/aaindex_help.html.

References

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

‘AAIndex’ is the work of Kanehisa and co-workers:
Kawashima and Kanehisa (2000) Nucleic Acids Res. 28, 374;
Tomii and Kanehisa (1996) Protein Eng. 9, 27–36;
Nakai, Kidera and Kanehisa (1988) Protein Eng. 2, 93–100.

Examples

## Load AAindex data
data(aa.index)

## Find all indeces described as "volume"
ind <- which(sapply(aa.index, function(x)
                    length(grep("volume", x$D, ignore.case=TRUE)) != 0))

## find all indeces with author "Kyte"
ind <- which(sapply(aa.index, function(x) length(grep("Kyte", x$A)) != 0))

## examine the index
aa.index[[ind]]$I

## find indeces which correlate with it
all.ind <- names(which(Mod(aa.index[[ind]]$C) >= 0.88))

## examine them all
sapply(all.ind, function (x) aa.index[[x]]$I)

Description

This data set provides the atomic masses of a selection of amino acids regularly occuring in proteins.

Usage

aa.table

Format

A data frame with the following components.

aa3

a character vector containing three-letter amino acid code.

aa1

a character vector containing one-letter amino acid code.

mass

a numeric vector containing the mass of the respective amino acids.

formula

a character vector containing the formula of the amino acid in which the mass calculat was based.

name

a character vector containing the full names of the respective amino acids.

References

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

See Also

aa2mass, aa.index, atom.index, elements,

Examples

data(aa.table)
aa.table

## table look up
aa.table["HIS", ]

## read PDB, and fetch residue masses
pdb <- read.pdb(system.file("examples/1hel.pdb", package="bio3d"))
aa2mass(pdb)

Description

Convert between one-letter IUPAC aminoacid codes and three-letter PDB style aminoacid codes.

Usage

aa123(aa)
aa321(aa)

Arguments

Details

Standard conversions will map ‘A’ to ‘ALA’, ‘G’ to ‘GLY’, etc. Non-standard codes in aa will generate a warning and return ‘UNK’ or ‘X’.

Value

A character vector of aminoacid codes.

Author(s)

Barry Grant

References

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

For a description of IUPAC one-letter codes see:
https://www.insdc.org/documents/feature_table.html#7.4.3

For more information on PDB residue codes see:
http://ligand-expo.rcsb.org/ld-search.html

See Also

read.pdb, read.fasta, pdbseq

Examples

# Simple conversion
aa123(c("D","L","A","G","S","H"))
aa321(c("ASP", "LEU", "ALA", "GLY", "SER", "HIS"))

## Not run: 
# Extract sequence from a PDB file's ATOM and SEQRES cards
pdb <- read.pdb("1BG2") 
s <- aa321(pdb$seqres)                   # SEQRES
a <- aa321(pdb$atom[pdb$calpha,"resid"]) # ATOM

# Write both sequences to a fasta file
write.fasta(alignment=seqbind(s,a), id=c("seqres","atom"), file="eg2.fa")

# Alternative approach for ATOM sequence extraction
pdbseq(pdb)
pdbseq(pdb, aa1=FALSE )

## End(Not run)

Description

Converts sequences to aminoacid indeces from the ‘AAindex’ database.

Usage

aa2index(aa, index = "KYTJ820101", window = 1)

Arguments

Details

By default, this function simply returns the index values for each amino acid in the sequence. It can also be set to perform a crude sliding window average through the window argument.

Value

Returns a numeric vector.

Author(s)

Ana Rodrigues

References

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

‘AAIndex’ is the work of Kanehisa and co-workers: Kawashima and Kanehisa (2000) Nucleic Acids Res. 28, 374; Tomii and Kanehisa (1996) Protein Eng. 9, 27–36; Nakai, Kidera and Kanehisa (1988) Protein Eng. 2, 93–100.

For a description of the ‘AAindex’ database see:
https://www.genome.jp/aaindex/ or the aa.index documentation.

See Also

aa.index, read.fasta

Examples

## Residue hydropathy values
seq <- c("R","S","D","X","-","X","R","H","Q","V","L")
aa2index(seq)

## Not run: 
## Use a sliding window average
aa2index(aa=seq, index=22, window=3)

## Use an alignment

aln  <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))
prop <- t(apply(aln$ali, 1, aa2index, window=1))

## find and use indices for volume calculations
i <- which(sapply(aa.index,
       function(x) length(grep("volume", x$D, ignore.case=TRUE)) != 0))
sapply(i, function(x) aa2index(aa=seq, index=x, window=5)) 

## End(Not run)

Description

Convert a sequence of amino acid residue names to mass.

Usage

aa2mass(pdb, inds=NULL, mass.custom=NULL, addter=TRUE, mmtk=FALSE)

Arguments

Details

This function converts amino acid residue names to their corresponding masses. In the case of a non-standard amino acid residue name mass.custom can be used to map the residue to the correct mass. User-defined amino acid masses (with argument mass.custom) will override mass entries obtained from the database.

See examples for more details.

Value

Returns a numeric vector of masses.

Note

When object of type pdb is provided, non-calpha atom records are omitted from the selection.

Author(s)

Lars Skjaerven

References

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

See Also

atom.index, atom2mass, aa.index

Examples

resi.names <- c("LYS", "ALA", "CYS", "HIS")
masses <- aa2mass(resi.names, addter=FALSE)

## Not run: 
## Fetch atomic masses in a PDB object
pdb <- read.pdb("3dnd")
masses <- aa2mass(pdb)

## or
masses <- aa2mass(pdb$atom[1:10,"resid"])

## Dealing with unconventional residues
#pdb <- read.pdb("1xj0")

#mass.cust <- list("CSX"=122.166)
#masses <- aa2mass(pdb, mass.custom=mass.cust)

## End(Not run)

Description

Perform all-atom elastic network model normal modes calculation of a protein structure.

Usage

aanma(...)

## S3 method for class 'pdb'
aanma(pdb, pfc.fun = NULL, mass = TRUE, temp = 300,
  keep = NULL, hessian = NULL, outmodes = "calpha", rm.wat = TRUE,
  reduced = FALSE, rtb = FALSE, nmer = 1, ...)

rtb(hessian, pdb, mass = TRUE, nmer = 1, verbose = TRUE)

Arguments

Details

This function builds an elastic network model (ENM) based on all heavy atoms of input pdb, and performs subsequent normal mode analysis (NMA) in various manners. By default, the ‘aaenm2’ force field (defining of the spring constants between atoms) is used, which was obtained by fitting to a local energy minimum of a crambin model derived from the AMBER99SB force field. It employs a pair force constant function which falls as r^-6, and specific force constants for covalent and intra-residue atom pairs. See also load.enmff for other force field options.

The outmodes argument controls the type of output modes. There are two standard types of output modes: ‘noh’ and ‘calpha’. outmodes='noh' invokes regular all-atom based ENM-NMA. When outmodes='calpha', an effective Hessian with respect to all C-alpha atoms will be first calculated using the same formula as in Hinsen et al. NMA is then performed on this effective C-alpha based Hessian. In addition, users can provide their own atom selection (see atom.select) as the value of outmodes for customized output modes generation.

When reduced=TRUE, only a selection of all heavy atoms is used to build the ENM. More specifically, three to five atoms per residue constitute the model. Here the N, CA, C atoms represent the protein backbone, and zero to two selected side chain atoms represent the side chain (selected based on side chain size and the distance to CA). This coarse-grained ENM has significantly improved computational efficiency and similar prediction accuracy with respect to the all-atom ENM.

When rtb=TRUE, rotation-translation block (RTB) based approximate modes will be calculated. In this method, each residue is assumed to be a rigid body (or ‘block’) that has only rotational and translational degrees of freedom. Intra-residue deformation is thus ignored. (See Durand et al 1994 and Tama et al. 2000 for more details). N residues per block is also supported, where N=1, 2, 3, etc. (See argument nmer). The RTB method has significantly improved computational efficiency and similar prediction accuracy with respect to the all-atom ENM.

By default the function will diagonalize the mass-weighted Hessian matrix. The resulting mode vectors are moreover scaled by the thermal fluctuation amplitudes.

Value

Returns an object of class ‘nma’ with the following components:

Author(s)

Lars Skjaerven & Xin-Qiu Yao

References

Hinsen, K. et al. (2000) Chem. Phys. 261, 25. Durand, P. et al. (1994) Biopolymers 34, 759. Tama, F. et al. (2000) Proteins 41, 1.

See Also

nma.pdb for C-alpha based NMA, aanma.pdbs for ensemble all-atom NMA, load.enmff for available ENM force fields, and fluct.nma, mktrj.nma, and dccm.nma for various post-NMA calculations.

Examples

## Not run: 
   # All-atom NMA takes relatively long time - Don't run by default.

   ## Fetch stucture
   pdb <- read.pdb( system.file("examples/1hel.pdb", package="bio3d") )

   ## Calculate all-atom normal modes
   modes.aa <- aanma(pdb, outmodes='noh')

   ## Calculate all-atom normal modes with RTB approximation
   modes.aa.rtb <- aanma(pdb, outmodes='noh', rtb=TRUE)

   ## Compare the two modes
   rmsip(modes.aa, modes.aa.rtb)

   ## Calculate C-alpha normal modes.
   modes <- aanma(pdb)

   ## Calculate C-alpha normal modes with reduced ENM.
   modes.cg <- aanma(pdb, reduced=TRUE)

   ## Calculate C-alpha normal modes with RTB approximation
   modes.rtb <- aanma(pdb, rtb=TRUE)

   ## Compare modes
   rmsip(modes, modes.cg)
   rmsip(modes, modes.rtb)

   ## Print modes
   print(modes)

   ## Plot modes
   plot(modes)

   ## Visualize modes
   #m7 <- mktrj.nma(modes, mode=7, file="mode_7.pdb", pdb=pdb)

## End(Not run)

Description

Perform normal mode analysis (NMA) on an ensemble of aligned protein structures using all-atom elastic network model (aaENM).

Usage

## S3 method for class 'pdbs'
aanma(pdbs, fit = TRUE, full = FALSE, subspace = NULL,
  rm.gaps = TRUE, ligand = FALSE, outpath = NULL, gc.first = TRUE,
  ncore = NULL, ...)

Arguments

Details

This function builds elastic network model (ENM) using all heavy atoms and performs subsequent normal mode analysis (NMA) on a set of aligned protein structures obtained with function read.all. The main purpose is to automate ensemble normal mode analysis using all-atom ENMs.

By default, the effective Hessian for all C-alpha atoms is calculated based on the Hessian built from all heavy atoms (including ligand atoms if ligand=TRUE). Returned values include aligned mode vectors and (when full=TRUE) a list containing the full ‘nma’ objects one per each structure. When ‘rm.gaps=TRUE’ the unaligned atoms are ommited from output. With default arguments ‘rmsip’ provides RMSIP values for all pairwise structures.

When outmodes is provided and is not ‘calpha’ (e.g. ‘noh’. See aanma for more details), the function simply returns a list of ‘nma’ objects, one per each structure, and no aligned mode vector is returned. In this case, the arguments full, subspace, and rm.gaps are ignored. This is equivalent to a wrapper function repeatedly calling aanma.

Value

Returns a list of ‘nma’ objects (outmodes is provided and is not ‘calpha’) or an ‘enma’ object with the following components:

Author(s)

Xin-Qiu Yao & Lars Skjaerven

See Also

For normal mode analysis on single structure PDB: aanma

For conventional C-alpha based normal mode analysis: nma, nma.pdbs.

For the analysis of the resulting ‘eNMA’ object: mktrj.enma, dccm.enma, plot.enma, cov.enma.

Similarity measures: sip, covsoverlap, bhattacharyya, rmsip.

Related functionality: read.all.

Examples

# Needs MUSCLE installed - testing excluded
  if(check.utility("muscle")) {

    ## Fetch PDB files and split to chain A only PDB files
    ids <- c("1a70_A", "1czp_A", "1frd_A", "1fxi_A", "1iue_A", "1pfd_A")
    files <- get.pdb(ids, split = TRUE, path = tempdir())
    
    ## Sequence Alignement
    aln <- pdbaln(files, outfile = tempfile())
    
    ## Read all pdb coordinates
    pdbs <- read.all(aln)

    ## Normal mode analysis on aligned data
    modes <- aanma(pdbs, rm.gaps=TRUE)
    
    ## Plot fluctuation data
    plot(modes, pdbs=pdbs)
    
    ## Cluster on Fluctuation similariy
    sip <- sip(modes)
    hc <- hclust(dist(sip))
    col <- cutree(hc, k=3)
    
    ## Plot fluctuation data
    plot(modes, pdbs=pdbs, col=col)
    
    ## RMSIP is pre-calculated
    heatmap(1-modes$rmsip)
    
    ## Bhattacharyya coefficient
    bc <- bhattacharyya(modes)
    heatmap(1-bc)

  }

Description

Renders a sequence alignment as coloured HTML suitable for viewing with a web browser.

Usage

aln2html(aln, file="alignment.html", Entropy=0.5, append=TRUE,
         caption.css="color: gray; font-size: 9pt",
         caption="Produced by <a href=http://thegrantlab.org/bio3d/>Bio3D</a>",
         fontsize="11pt", bgcolor=TRUE, colorscheme="clustal")

Arguments

Value

Called for its effect.

Note

Your web browser should support style sheets.

Author(s)

Barry Grant

References

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

See Also

read.fasta, write.fasta, seqaln

Examples

## Not run: 
## Read an example alignment
aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))

## Produce a HTML file for this alignment
aln2html(aln, append=FALSE, file=file.path("eg.html"))
aln2html(aln, colorscheme="ent", file="eg.html")
## View/open the file in your web browser
#browseURL("eg.html")

## End(Not run)

Description

A function for basic bond angle determination.

Usage

angle.xyz(xyz, atm.inc = 3)

Arguments

Value

Returns a numeric vector of angles.

Note

With atm.inc=1, angles are calculated for each set of three successive atoms contained in xyz (i.e. moving along one atom, or three elements of xyz, between sucessive evaluations). With atm.inc=3, angles are calculated for each set of three successive non-overlapping atoms contained in xyz (i.e. moving along three atoms, or nine elements of xyz, between sucessive evaluations).

Author(s)

Barry Grant

References

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

See Also

torsion.pdb, torsion.xyz, read.pdb, read.dcd.

Examples

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

## Angle between N-CA-C atoms of residue four
inds <- atom.select(pdb, resno=4, elety=c("N","CA","C"))
angle.xyz(pdb$xyz[inds$xyz])

## Basic stats of all N-CA-C bound angles
inds <- atom.select(pdb, elety=c("N","CA","C"))
summary( angle.xyz(pdb$xyz[inds$xyz]) )
#hist( angle.xyz(pdb$xyz[inds$xyz]), xlab="Angle" )

Description

Convert alignment/sequence in matrix/vector format to FASTA object.

Usage

as.fasta(x, id=NULL, ...)

Arguments

Details

This function provides basic functionality to convert a sequence character matrix/vector to a FASTA object.

Value

Returns a list of class "fasta" with the following components:

Author(s)

Lars Skjaerven

References

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

See Also

get.seq, seqaln, seqbind, pdbaln

Examples

as.fasta(c("A", "C", "D"))

Description

Convert Tripos Mol2 format, or Amber parameter/topology and coordinate data to PDB format.

Usage

as.pdb(...)

## S3 method for class 'mol2'
as.pdb(mol, ...)

## S3 method for class 'prmtop'
as.pdb(prmtop, crd=NULL, inds=NULL, inds.crd=inds, ncore=NULL, ...)

## Default S3 method:
as.pdb(pdb=NULL, xyz=NULL, type=NULL, resno=NULL,
                   resid=NULL, eleno=NULL, elety=NULL, chain=NULL, 
                   insert=NULL, alt=NULL, o=NULL, b=NULL, segid=NULL, 
                   elesy=NULL, charge=NULL, verbose=TRUE, ...)

Arguments

Details

This function converts Tripos Mol2 format, Amber formatted parameter/topology (PRMTOP) and coordinate objects, and vector data to a PDB object.

While as.pdb.mol2 and as.pdb.prmtop converts specific objects to a PDB object, as.pdb.default provides basic functionality to convert raw data such as vectors of e.g. residue numbers, residue identifiers, Cartesian coordinates, etc to a PDB object. When pdb is provided the returned PDB object is built from the input object with fields replaced by any input vector arguments. e.g. as.pdb(pdb, xyz=crd) will return the same PDB object, with only the Cartesian coordinates changed to crd.

Value

Returns a list of class "pdb" with the following components:

Author(s)

Lars Skjaerven

References

Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696. http://ambermd.org/FileFormats.php

See Also

read.crd, read.ncdf, atom.select, read.pdb

Examples

## Vector(s) to PDB object
pdb <- as.pdb(resno=1:6, elety="CA", resid="ALA", chain="A")
pdb

## Not run: 
## Read a PRMTOP file
prmtop <- read.prmtop(system.file("examples/crambin.prmtop", package="bio3d"))

## Read Amber coordinates
crds <- read.crd(system.file("examples/crambin.inpcrd", package="bio3d"))

## Atom selection
ca.inds <- atom.select(prmtop, "calpha")

## Convert to PDB format
pdb <- as.pdb(prmtop, crds, inds=ca.inds)


## Read a single entry MOL2 file
## (returns a single object)
mol <- read.mol2( system.file("examples/aspirin.mol2", package="bio3d") )

## Convert to PDB
pdb <- as.pdb(mol)

## End(Not run)

Description

Convert atomic indices to a select object with ‘atom’ and ‘xyz’ components.

Usage

as.select(x, ...)

Arguments

Details

Convert atomic indices to a select object with ‘atom’ and ‘xyz’ components.

Value

Returns a list of class "select" with the following components:

Author(s)

Lars Skjaerven

References

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

See Also

atom.select, read.pdb

Examples

as.select(c(1,2,3))

Description

This data set gives for various atom names/types the corresponding atomic symbols.

Usage

atom.index

Format

A data frame with the following components.

name

a character vector containing atom names/types.

symb

a character vector containing atomic symbols.

References

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

See Also

elements, atom.index, atom2ele

Examples

data(atom.index)
atom.index

# Get the atomic symbol of some atoms
atom.names <- c("CA", "O", "N", "OXT")
atom.index[match(atom.names, atom.index$name), "symb"]

Description

Return the ‘atom’ and ‘xyz’ coordinate indices of ‘pdb’ or ‘prmtop’ structure objects corresponding to the intersection of a hierarchical selection.

Usage

atom.select(...)

## S3 method for class 'pdb'
atom.select(pdb, string = NULL,
                          type  = NULL, eleno = NULL, elety = NULL,
                          resid = NULL, chain = NULL, resno = NULL,
                          insert = NULL, segid = NULL, 
                          operator = "AND", inverse = FALSE,
                          value = FALSE, verbose=FALSE, ...)
## S3 method for class 'pdbs'
atom.select(pdbs, string = NULL, 
                           resno = NULL, chain = NULL, resid = NULL,
                           operator="AND", inverse = FALSE,
                           value = FALSE, verbose=FALSE, ...)
## S3 method for class 'mol2'
atom.select(mol, string=NULL,
                           eleno = NULL, elena = NULL, elety = NULL,
                           resid = NULL, chain = NULL, resno = NULL,
                           statbit = NULL,
			   operator = "AND", inverse = FALSE,
                           value = FALSE, verbose=FALSE,  ...)

## S3 method for class 'prmtop'
atom.select(prmtop, ...)

## S3 method for class 'select'
print(x, ...)

Arguments

Details

This function allows for the selection of atom and coordinate data corresponding to the intersection of various input criteria.

Input selection criteria include selection string keywords (such as "calpha", "backbone", "sidechain", "protein", "nucleic", "ligand", etc.) and individual named selection components (including ‘chain’, ‘resno’, ‘resid’, ‘elety’ etc.).

For example, atom.select(pdb, "calpha") will return indices for all C-alpha (CA) atoms found in protein residues in the pdb object, atom.select(pdb, "backbone") will return indices for all protein N,CA,C,O atoms, and atom.select(pdb, "cbeta") for all protein N,CA,C,O,CB atoms.

Note that keyword string shortcuts can be combined with individual selection components, e.g. atom.select(pdb, "protein", chain="A") will select all protein atoms found in chain A.

Selection criteria are combined according to the provided operator argument. The default operator AND (or &) will combine by intersection while OR (or |) will take the union.

For example, atom.select(pdb, "protein", elety=c("N", "CA", "C"), resno=65:103) will select the N, CA, C atoms in the protein residues 65 through 103, while atom.select(pdb, "protein", resid="ATP", operator="OR") will select all protein atoms as well as any ATP residue(s).

Other string shortcuts include: "calpha", "back", "backbone", "cbeta", "protein", "notprotein", "ligand", "water", "notwater", "h", "noh", "nucleic", and "notnucleic".

In addition, the combine.select function can further combine atom selections using ‘AND’, ‘OR’, or ‘NOT’ logical operations.

Value

Returns a list of class "select" with the following components:

Note

Protein atoms are defined as any atom in a residue matching the residue name in the attached aa.table data frame. See aa.table$aa3 for a complete list of residue names.

Nucleic atoms are defined as all atoms found in residues with names A, U, G, C, T, I, DA, DU, DG, DC, DT, or DI.

Water atoms/residues are defined as those with residue names H2O, OH2, HOH, HHO, OHH, SOL, WAT, TIP, TIP, TIP3, or TIP4.

Author(s)

Barry Grant, Lars Skjaerven

References

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

See Also

read.pdb, as.select, combine.select, trim.pdb, write.pdb, read.prmtop, read.crd, read.dcd, read.ncdf.

Examples

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

# Select protein atoms of chain A
atom.select(pdb, "protein", chain="A")

# Select all atoms except from the protein
atom.select(pdb, "protein", inverse=TRUE, verbose=TRUE)

# Select all C-alpha atoms with residues numbers between 43 and 54
sele <- atom.select(pdb, "calpha", resno=43:54, verbose=TRUE)

# Access the PDB data with the selection indices
print( pdb$atom[ sele$atom, "resid" ] )
print( pdb$xyz[ sele$xyz ] )

# Trim PDB to selection
ca.pdb <- trim.pdb(pdb, sele)

## Not run: 

##- PRMTOP example
prmtop <- read.prmtop(system.file("examples/crambin.prmtop", package="bio3d"))

## Atom selection
ca.inds <- atom.select(prmtop, "calpha")


## End(Not run)

Description

Convert atom names/types into atomic symbols

Usage

atom2ele(...)

## Default S3 method:
atom2ele(x, elety.custom=NULL, rescue=TRUE, ...)

## S3 method for class 'pdb'
atom2ele(pdb, inds=NULL, ...)

Arguments

Details

The default method searchs for the atom names/types in the atom.index data set and returns their corresponding atomic symbols. If elety.custom is specified it is combined with atom.index (using rbind) before searching. Therefore, elety.custom must contains columns named name and symb.

The S3 method for object of class ‘pdb’, pass pdb$atom[,"elety"] to the default method.

Value

Return a character vector of atomic symbols

Author(s)

Julien Ide, Lars Skjaerven

See Also

atom.index, elements, read.pdb, atom2mass, formula2mass

Examples

atom.names <- c("CA", "O", "N", "OXT")
atom2ele(atom.names)


# PDB server connection required - testing excluded

## Get atomic symbols from a PDB object with a customized data set
pdb <- read.pdb("3RE0",verbose=FALSE)
lig <- trim(pdb, "ligand")

## maps PT1 to Pt, CL2 to Cl, C4A to C
atom2ele(lig)

## map atom name to element manually
myelety <- data.frame(name = "CL2", symb = "Cl")
atom2ele(lig, elety.custom = myelety)

Description

Convert atom names/types into atomic masses.

Usage

atom2mass(...)
## Default S3 method:
atom2mass(x, mass.custom=NULL, elety.custom=NULL,
                            grpby=NULL, rescue=TRUE, ...) 
## S3 method for class 'pdb'
atom2mass(pdb, inds=NULL, mass.custom=NULL,
                        elety.custom=NULL, grpby=NULL, rescue=TRUE, ...)

Arguments

Details

The default method first convert atom names/types into atomic symbols using the atom2ele function. Then, atomic symbols are searched in the elements data set and their corresponding masses are returned. If mass.custom is specified it is combined with elements (using rbind) before searching. Therefore, mass.custom must have columns named symb and mass (see examples). If grpby is specified masses are splitted (using split) to compute the mass of groups of atoms defined by grpby.

The S3 method for object of class ‘pdb’, pass pdb$atom$elety to the default method.

Value

Return a numeric vector of masses.

Author(s)

Julien Ide, Lars Skjaerven

See Also

elements, atom.index, atom2ele, read.pdb

Examples

atom.names <- c("CA", "O", "N", "OXT")
atom2mass(atom.names)


# PDB server connection required - testing excluded

## Get atomic symbols from a PDB object with a customized data set
pdb <- read.pdb("3RE0", verbose=FALSE)
inds <- atom.select(pdb, resno=201, verbose=FALSE)

## selected atoms
print(pdb$atom$elety[inds$atom])

## default will map CL2 to C
atom2mass(pdb, inds)

## map element CL2 correctly to Cl
myelety  <- data.frame(name = c("CL2","PT1","N1","N2"), symb = c("Cl","Pt","N","N"))
atom2mass(pdb, inds, elety.custom = myelety)

## custom masses
mymasses <- data.frame(symb = c("Cl","Pt"), mass = c(35.45, 195.08))
atom2mass(pdb, inds, elety.custom = myelety, mass.custom = mymasses)

Description

Basic functions to convert between xyz and their corresponding atom indices.

Usage

atom2xyz(num)
xyz2atom(xyz.ind)

Arguments

Value

A numeric vector of either xyz or atom indices.

Author(s)

Barry Grant

References

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

See Also

atom.select, read.pdb

Examples

xyz.ind <- atom2xyz(c(1,10,15))
xyz2atom( xyz.ind )

Description

Removes all of the path up to and including the last path separator (if any) and the final ‘.pdb’ extension.

Usage

basename.pdb(x, mk4 = FALSE, ext=".pdb")

Arguments

Details

This is a simple utility function for the common task of PDB file name manipulation. It is used internally in several bio3d functions and van be thought of as basename for PDB files.

Value

A character vector of the same length as the input ‘x’.

Paths not containing any separators are taken to be in the current directory.

If an element of input is ‘x’ is ‘NA’, so is the result.

Author(s)

Barry Grant

References

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

See Also

basename, dirname

Examples

basename.pdb("/somedir/somewhere/1bg2_myfile.pdb")
	basename.pdb("/somedir/somewhere/1bg2_myfile.pdb", TRUE)

Description

Calculate the Bhattacharyya Coefficient as a similarity between two modes objects.

Usage

bhattacharyya(...)

## S3 method for class 'enma'
bhattacharyya(enma, covs=NULL, ncore=NULL, ...)

## S3 method for class 'array'
bhattacharyya(covs, ncore=NULL, ...)

## S3 method for class 'matrix'
bhattacharyya(a, b, q=90, n=NULL, ...)

## S3 method for class 'nma'
bhattacharyya(...)

## S3 method for class 'pca'
bhattacharyya(...)

Arguments

Details

Bhattacharyya coefficient provides a means to compare two covariance matrices derived from NMA or an ensemble of conformers (e.g. simulation or X-ray conformers).

Value

Returns the similarity coefficient(s).

Author(s)

Lars Skjaerven

References

Skjaerven, L. et al. (2014) BMC Bioinformatics 15, 399. Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696. Fuglebakk, E. et al. (2013) JCTC 9, 5618–5628.

See Also

Other similarity measures: sip, covsoverlap, rmsip.

Description

Determines the interacting residues between two PDB entities.

Usage

binding.site(a, b=NULL, a.inds=NULL, b.inds=NULL, cutoff=5,
            hydrogens=TRUE, byres=TRUE, verbose=FALSE)

Arguments

Details

This function reports the residues of a closer than a cutoff to b. This is a wrapper function calling the underlying function dist.xyz.

If b=NULL then b.inds should be elements of a upon which the calculation is based (typically chain A and B of the same PDB file).

If b=a.inds=b.inds=NULL the function will use atom.select with arguments "protein" and "ligand" to determine receptor and ligand, respectively.

Value

Returns a list with the following components:

Author(s)

Lars Skjaerven

References

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

See Also

read.pdb, atom.select, dm

Examples

# PDB server connection required - testing excluded

     pdb <- read.pdb('3dnd')

     ## automatically identify 'protein' and 'ligand'
     bs <- binding.site(pdb)

     bs$resnames
     #pdb$atom[bs$inds$atom, ]

     # provide indices
     rec.inds <- atom.select(pdb, chain="A", resno=1:350)
     lig.inds <- atom.select(pdb, chain="A", resno=351)
     bs <- binding.site(pdb, a.inds=rec.inds, b.inds=lig.inds)
   


   ## Not run:   
     # Interaction between peptide and protein
     rec.inds <- atom.select(pdb, chain='A', resno=c(1:350))
     lig.inds <- atom.select(pdb, chain='I', resno=c(5:24))
     bs <- binding.site(pdb, a.inds=rec.inds, b.inds=lig.inds)
   
## End(Not run)

   
    # Redundant testing excluded

     # Interaction between two PDB entities   
     #rec <- read.pdb("receptor.pdb")
     #lig <- read.pdb("ligand.pdb")
     rec <- trim.pdb(pdb, inds=rec.inds)
     lig <- trim.pdb(pdb, inds=lig.inds)
     bs <- binding.site(rec, lig, hydrogens=FALSE)

Description

Construct biological assemblies/units based on a 'pdb' object.

Usage

biounit(pdb, biomat = NULL, multi = FALSE, ncore = NULL)

Arguments

Details

A valid structural/simulation study should be performed on the biological unit of a protein system. For example, the alpha2-beta2 tetramer form of hemoglobin. However, canonical PDB files usually contain the asymmetric unit of the crystal cell, which can be:

1. One biological unit

2. A portion of a biological unit

3. Multiple biological units

The function performs symmetry operations to the coordinates based on the transformation matrices stored in a 'pdb' object returned by read.pdb, and returns biological units stored as a list of pdb objects.

Value

a list of pdb objects with each representing an individual biological unit.

Author(s)

Xin-Qiu Yao

See Also

read.pdb

Examples

# PDB server connection required - testing excluded

   pdb <- read.pdb("2dn1")
   biounit <- biounit(pdb)
   pdb
   biounit

## Not run: 
   biounit <- biounit(read.pdb("2bfu"), multi=TRUE)
   write.pdb(biounit[[1]], file="biounit.pdb")
   # open the pdb file in VMD to have a look on the biological unit

## End(Not run)

Description

Run NCBI blastp, on a given sequence, against the PDB, NR and swissprot sequence databases. Produce plots that facilitate hit selection from the match statistics of a BLAST result.

Usage

blast.pdb(seq, database = "pdb", time.out = NULL, chain.single=TRUE)

get.blast(urlget, time.out = NULL, chain.single=TRUE)

## S3 method for class 'blast'
plot(x, cutoff = NULL, cut.seed=NULL, cluster=TRUE, mar=c(2, 5, 1, 1), cex=1.5, ...)

Arguments

Details

The blast.pdb function employs direct HTTP-encoded requests to the NCBI web server to run BLASTP, the protein search algorithm of the BLAST software package.

BLAST, currently the most popular pairwise sequence comparison algorithm for database searching, performs gapped local alignments via a heuristic strategy: it identifies short nearly exact matches or hits, bidirectionally extends non-overlapping hits resulting in ungapped extended hits or high-scoring segment pairs(HSPs), and finally extends the highest scoring HSP in both directions via a gapped alignment (Altschul et al., 1997)

For each pairwise alignment BLAST reports the raw score, bitscore and an E-value that assess the statistical significance of the raw score. Note that unlike the raw score E-values are normalized with respect to both the substitution matrix and the query and database lengths.

Here we also return a corrected normalized score (mlog.evalue) that in our experience is easier to handle and store than conventional E-values. In practice, this score is equivalent to minus the natural log of the E-value. Note that, unlike the raw score, this score is independent of the substitution matrix and and the query and database lengths, and thus is comparable between BLASTP searches.

Examining plots of BLAST alignment lengths, scores, E-values and normalized scores (-log(E-Value) from the blast.pdb function can aid in the identification sensible hit similarity thresholds. This is facilitated by the plot.blast function.

If a ‘cutoff’ value is not supplied then a basic hierarchical clustering of normalized scores is performed with initial group partitioning implemented at a hopefully sensible point in the vicinity of ‘h=cut.seed’. Inspection of the resultant plot can then be use to refine the value of ‘cut.seed’ or indeed ‘cutoff’. As the ‘cutoff’ value can vary depending on the desired application and indeed the properties of the system under study it is envisaged that ‘plot.blast’ will be called multiple times to aid selection of a suitable ‘cutoff’ value. See the examples below for further details.

Value

The function blast.pdb returns a list with three components, hit.tbl, raw, and url. The function plot.blast produces a plot on the active graphics device and returns a list object with four components, hits, pdb.id, acc, and inds. See below:

Note

Online access is required to query NCBI blast services.

Author(s)

Barry Grant

References

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

‘BLAST’ is the work of Altschul et al.: Altschul, S.F. et al. (1990) J. Mol. Biol. 215, 403–410.

Full details of the ‘BLAST’ algorithm, along with download and installation instructions can be obtained from:
https://www.ncbi.nlm.nih.gov/BLAST/.

See Also

plot.blast, hmmer, seqaln, get.pdb

Examples

## Not run: 
pdb <- read.pdb("4q21")
blast <- blast.pdb( pdbseq(pdb) )

head(blast$hit.tbl)
top.hits <- plot(blast)
head(top.hits$hits)

## Use 'get.blast()' to retrieve results at a later time.
#x <- get.blast(blast$url)
#head(x$hit.tbl)

# Examine and download 'best' hits
top.hits <- plot.blast(blast, cutoff=188)
head(top.hits$hits)
#get.pdb(top.hits)

## End(Not run)

Description

Find the ‘bounds’ (i.e. start, end and length) of consecutive numbers within a larger set of numbers in a given vector.

Usage

bounds(nums, dup.inds=FALSE, pre.sort=TRUE)

Arguments

Details

This is a simple utility function useful for summarizing the contents of a numeric vector. For example: find the start position, end position and lengths of secondary structure elements given a vector of residue numbers obtained from a DSSP secondary structure prediction.

By setting ‘dup.inds’ to TRUE then the indices of the first (start) and last (end) duplicated elements of the vector are returned. For example: find the indices of atoms belonging to a particular residue given a vector of residue numbers (see below).

Value

Returns a three column matrix listing starts, ends and lengths.

Author(s)

Barry Grant

References

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

Examples

test <- c(seq(1,5,1),8,seq(10,15,1))
 bounds(test)

 test <- rep(c(1,2,4), times=c(2,3,4))
 bounds(test, dup.ind=TRUE)

Description

Inverse process of the funciton pdb2sse.

Usage

bounds.sse(x, pdb = NULL)

Arguments

Details

call for its effects.

Value

a 'sse' object.

Note

In both $helix and $sheet, an additional $id component is added to indicate the original numbering of the sse. This is particularly useful in e.g. trim.pdb() function.

Author(s)

Xin-Qiu Yao & Barry Grant

See Also

pdb2sse

Examples

# PDB server connection required - testing excluded

   pdb <- read.pdb("1a7l")
   sse <- pdb2sse(pdb)
   sse.ind <- bounds.sse(sse)
   sse.ind

Description

Create a vector of ‘n’ “contiguous” colors forming either a Blue-White-Red or a White-Gray-Black color palette.

Usage

bwr.colors(n)
mono.colors(n)

Arguments

Details

The function bwr.colors returns a vector of n color names that range from blue through white to red.

The function mono.colors returns color names ranging from white to black. Note: the first element of the returned vector will be NA.

Value

Returns a character vector, cv, of color names. This can be used either to create a user-defined color palette for subsequent graphics with palette(cv), or as a col= specification in graphics functions and par.

Author(s)

Barry Grant

References

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

The bwr.colors function is derived from the gplots package function colorpanel by Gregory R. Warnes.

See Also

vmd_colors, cm.colors, colors, palette, hsv, rgb, gray, col2rgb

Examples

# Redundant testing excluded

# Color a distance matrix
pdb <- read.pdb( system.file("examples/1hel.pdb", package="bio3d") )
d <- dm(pdb,"calpha")

plot(d, color.palette=bwr.colors)

plot(d,
     resnum.1 = pdb$atom[pdb$calpha,"resno"],
     color.palette = mono.colors,
     xlab="Residue Number", ylab="Residue Number")

Description

Produce a new concatenated PDB object from two or more smaller PDB objects.

Usage

cat.pdb(..., renumber=FALSE, rechain=TRUE)

Arguments

Details

This is a basic utility function for creating a concatenated PDB object based on multipe smaller PDB objects.

Value

Returns an object of class "pdb". See read.pdb for further details.

Author(s)

Lars Skjaerven

References

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

See Also

read.pdb, atom.select, write.pdb, trim.pdb

Examples

## Not run: 
## Read a PDB file from the RCSB online database
pdb1 <- read.pdb("1etl")
pdb2 <- read.pdb("1hel")

## Concat
new.pdb <- cat.pdb(pdb1, pdb2, pdb1, rechain=TRUE, renumber=TRUE)

## Write to file
write.pdb(new.pdb, file="concat.pdb")

## End(Not run)

Description

Find possible chain breaks based on connective Calpha or peptide bond (C-N) atom separation.

Usage

chain.pdb(pdb, ca.dist = 4, bond=TRUE, bond.dist=1.5, blank = "X", rtn.vec = TRUE)

Arguments

Details

This is a basic function for finding possible chain breaks in PDB structure files, i.e. connective Calpha atoms that are further than ca.dist apart or peptide bond (C-N) atoms separated by at least bond.dist.

Value

Prints basic chain information and if rtn.vec is TRUE returns a character vector of chain ids consisting of the 26 upper-case letters of the Roman alphabet plus possible blank entries for non-protein atoms.

Author(s)

Barry Grant

References

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

See Also

read.pdb, atom.select, trim.pdb, write.pdb

Examples

# PDB server connection required - testing excluded

full.pdb <- read.pdb( get.pdb("5p21", URLonly=TRUE) )
inds <- atom.select(full.pdb, resno=c(10:20,30:33))
cut.pdb <- trim.pdb(full.pdb, inds)
chain.pdb(cut.pdb)

Description

Internally used in examples, tests, or vignettes.

Usage

check.utility(x = c("muscle", "clustalo", "dssp", "stride", 
                "mustang", "makeup"), quiet = TRUE)

Arguments

Details

Check if requested utility programs are availabe or not.

Value

logical, TRUE if programs are available and FALSE if any one of them is missing.

Examples

check.utility(c("muscle", "dssp"), quiet=FALSE)
   if(!check.utility("mustang")) 
      cat(" The utility program, MUSTANG, is missing on your system\n")

Description

Inspect alternative coordinates, chain breaks, bad residue numbering, non-standard/unknow amino acids, etc. Return a 'clean' pdb object with fixed residue numbering and optionally relabeled chain IDs, corrected amino acid names, removed water, ligand, or hydrogen atoms. All changes are recorded in a log in the returned object.

Usage

clean.pdb(pdb, consecutive = TRUE, force.renumber = FALSE,
  fix.chain = FALSE, fix.aa = FALSE, rm.wat = FALSE, rm.lig = FALSE,
  rm.h = FALSE, verbose = FALSE)

Arguments

Details

call for its effects.

Value

a 'pdb' object with an additional $log component storing all the processing messages.

Author(s)

Xin-Qiu Yao & Barry Grant

See Also

read.pdb

Examples

# PDB server connection required - testing excluded

   pdb <- read.pdb("1a7l")
   clean.pdb(pdb)

Description

Construct a Contact Map for Given Protein Structure(s).

Usage

cmap(...)

## Default S3 method:
cmap(...)

## S3 method for class 'xyz'
cmap(xyz, grpby = NULL, dcut = 4, scut = 3, pcut=1, binary=TRUE,
            mask.lower = TRUE, collapse=TRUE, gc.first=FALSE, ncore=1, nseg.scale=1, ...)

## S3 method for class 'pdb'
cmap(pdb, inds = NULL, verbose = FALSE, ...)

## S3 method for class 'pdbs'
cmap(pdbs, rm.gaps=FALSE, all.atom=FALSE, ...)

Arguments

Details

A contact map is a simplified distance matrix. See the distance matrix function dm for further details.

Function "cmap.pdb" is a wrapper for "cmap.xyz" which selects all ‘notwater’ atoms and calculates the contact matrix grouped by residue number.

Value

Returns a N by N numeric matrix composed of zeros and ones, where one indicates a contact between selected atoms.

Author(s)

Barry Grant

References

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

See Also

dm, dccm, dist, dist.xyz

Examples

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

## Atom Selection indices
inds <- atom.select(pdb, "calpha")

## Reference contact map
ref.cont <- cmap( pdb$xyz[inds$xyz], dcut=6, scut=3 )
plot.cmap(ref.cont)

## Not run: 
##- Read Traj file
trj <- read.dcd( system.file("examples/hivp.dcd", package="bio3d") )
## For each frame of trajectory
sum.cont <- NULL
for(i in 1:nrow(trj)) {

  ## Contact map for frame 'i'
  cont <- cmap(trj[i,inds$xyz], dcut=6, scut=3)

  ## Product with reference
  prod.cont <- ref.cont * cont
  sum.cont <- c(sum.cont, sum(prod.cont,na.rm=TRUE))
}

plot(sum.cont, typ="l")

## End(Not run)

Description

This function builds both residue-based and community-based undirected weighted network graphs from an input correlation matrix, as obtained from the functions ‘dccm’, ‘dccm.nma’, and ‘dccm.enma’. Community detection/clustering is performed on the initial residue based network to determine the community organization and network structure of the community based network.

Usage

cna(cij, ...)
  ## S3 method for class 'dccm'
cna(cij, cutoff.cij=0.4, cm=NULL,  vnames=colnames(cij),
      cluster.method="btwn", collapse.method="max", 
      cols=vmd_colors(), minus.log=TRUE, ...)
  ## S3 method for class 'ensmb'
cna(cij, ..., ncore = NULL)

Arguments

Details

The input to this function should be a correlation matrix as obtained from the ‘dccm’, ‘dccm.mean’ or ‘dccm.nma’ and related functions. Optionally, a contact map ‘cm’ may also given as input to filter the correlation matrix resulting in the exclusion of network edges between non-contacting atom pairs (as defined in the contact map).

Internally this function calls the igraph package functions ‘graph.adjacency’, ‘edge.betweenness.community’, ‘walktrap.community’, ‘fastgreedy.community’, and ‘infomap.community’. The first constructs an undirected weighted network graph. The second performs Girvan-Newman style clustering by calculating the edge betweenness of the graph, removing the edge with the highest edge betweenness score, calculates modularity (i.e. the difference between the current graph partition and the partition of a random graph, see Newman and Girvan, Physical Review E (2004), Vol 69, 026113), then recalculating edge betweenness of the edges and again removing the one with the highest score, etc. The returned community partition is the one with the highest overall modularity value. ‘walktrap.community’ implements the Pons and Latapy algorithm based on the idea that random walks on a graph tend to get "trapped" into densely connected parts of it, i.e. a community. The random walk process is used to determine a distance between nodes. Nodes with low distance values are joined in the same community. ‘fastgreedy.community’ instead determines the community structure based on the optimization of the modularity. In the starting state each node is isolated and belongs to a separated community. Communities are then joined together (according to the network edges) in pairs and the modularity is calculated. At each step the join resulting in the highest increase of modularity is chosen. This process is repeated until a single community is obtained, then the partitioning with the highest modularity score is selected. ‘infomap.community’ finds community structure that minimizes the expected description length of a random walker trajectory.

Value

Returns a list object that includes igraph network and community objects with the following components:

If an ensemble of correlation matrices is provided, a list of ‘cna’ object, of the ‘ecna’ class, will be returned.

Author(s)

Guido Scarabelli and Barry Grant

See Also

plot.cna, summary.cna, vmd.cna, graph.adjacency, edge.betweenness.community, walktrap.community, fastgreedy.community, infomap.community

Examples

# PDB server connection required - testing excluded

if (!requireNamespace("igraph", quietly = TRUE)) {
   message('Need igraph installed to run this example')
} else {

##-- Build a correlation network from NMA results
## Read example PDB
pdb <- read.pdb("4Q21")

## Perform NMA
modes <- nma(pdb)
#plot(modes, sse=pdb)

## Calculate correlations 
cij <- dccm(modes)
#plot(cij, sse=pdb)

## Build, and betweenness cluster, a network graph
net <- cna(cij, cutoff.cij=0.35)
#plot(net, pdb)

## within VMD set 'coloring method' to 'Chain' and 'Drawing method' to Tube
#vmd.cna(net, trim.pdb(pdb, atom.select(pdb,"calpha")), launch=TRUE )

##-- Build a correlation network from MD results
## Read example 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)

## select residues 24 to 27 and 85 to 90 in both chains
inds <- atom.select(pdb, resno=c(24:27,85:90), elety='CA')

## lsq fit of trj on pdb
xyz <- fit.xyz(pdb$xyz, trj, fixed.inds=inds$xyz, mobile.inds=inds$xyz)

## calculate dynamical cross-correlation matrix
cij <- dccm(xyz)

## Build, and betweenness cluster, a network graph
net <- cna(cij)

# Plot coarse grained network based on dynamically coupled communities
xy <- plot.cna(net)
plot.dccm(cij, margin.segments=net$communities$membership)


##-- Begin to examine network structure - see CNA vignette for more details
net
summary(net)
attributes(net)
table( net$communities$members )

}

Description

Find k shortest paths between a pair of nodes, source and sink, in a correlation network.

Usage

cnapath(cna, from, to=NULL, k=10, collapse=TRUE, ncore=NULL, ...)
## S3 method for class 'cnapath'
summary(object, ..., pdb = NULL, label = NULL, col = NULL,
   plot = FALSE, concise = FALSE, cutoff = 0.1, normalize = TRUE, weight = FALSE)
## S3 method for class 'cnapath'
print(x, ...)
## S3 method for class 'cnapath'
plot(x, ...)
## S3 method for class 'ecnapath'
plot(x, ...)

Arguments

Value

The function cnapath returns a (or a list of) ‘cnapath’ class of list containing following three components:

The function summary.cnapath returns a matrix of (normalized) node degeneracy for ‘on path’ residues.

Author(s)

Xin-Qiu Yao

References

Yen, J.Y. (1971) Management Science 17, 712–716.

See Also

cna, cna.dccm, vmd.cna, vmd.cnapath, get.shortest.paths.

Examples

# Redundant testing excluded

if (!requireNamespace("igraph", quietly = TRUE)) {
   message('Need igraph installed to run this example')
} else {

attach(transducin)
inds = match(c("1TND_A", "1TAG_A"), pdbs$id)

npdbs <- trim(pdbs, row.inds=inds)
gaps.res <- gap.inspect(npdbs$ali)

modes <- nma(npdbs)
cij <- dccm(modes)
net <- cna(cij, cutoff.cij=0.3)

# get paths
pa1 <- cnapath(net[[1]], from = 314, to=172, k=50)
pa2 <- cnapath(net[[2]], from = 314, to=172, k=50)

# print the information of a path
pa1

# print two paths simultaneously
pas <- list(pa1, pa2)
names(pas) <- c("GTP", "GDP")
print.cnapath(pas)

# Or, for the same effect,
# summary(pa1, pa2, label=c("GTP", "GDP"))

# replace node numbers with residue name and residue number in the PDB file
pdb <- read.pdb("1tnd")
pdb <- trim.pdb(pdb, atom.select(pdb, chain="A", resno=npdbs$resno[1, gaps.res$f.inds]))
print.cnapath(pas, pdb=pdb)

# plot path length distribution and node degeneracy
print.cnapath(pas, pdb = pdb, col=c("red", "darkgreen"), plot=TRUE)

# View paths in 3D molecular graphic with VMD
#vmd.cnapath(pa1, pdb, launch = TRUE)
#vmd.cnapath(pa1, pdb, colors = 7, launch = TRUE)
#vmd.cnapath(pa1, pdb, spline=TRUE, colors=c("pink", "red"), launch = TRUE)
#pdb2 <- read.pdb("1tag")
#pdb2 <- trim.pdb(pdb2, atom.select(pdb2, chain="A", resno=npdbs$resno[2, gaps.res$f.inds]))
#vmd.cnapath(pa2, pdb2, launch = TRUE)

detach(transducin)

}

Description

Calculate the center of mass of a PDB object.

Usage

com(...)

## S3 method for class 'pdb'
com(pdb, inds=NULL, use.mass=TRUE, ...)

## S3 method for class 'xyz'
com(xyz, mass=NULL, ...)

Arguments

Details

This function calculates the center of mass of the provided PDB structure / Cartesian coordiantes. Atom names found in standard amino acids in the PDB are mapped to atom elements and their corresponding relative atomic masses.

In the case of an unknown atom name elety.custom and mass.custom can be used to map an atom to the correct atomic mass. See examples for more details.

Alternatively, the atom name will be mapped automatically to the element corresponding to the first character of the atom name. Atom names starting with character H will be mapped to hydrogen atoms.

Value

Returns the Cartesian coordinates at the center of mass.

Author(s)

Lars Skjaerven

References

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

See Also

read.pdb, atom2mass

Examples

# PDB server connection required - testing excluded

## Stucture of PKA:
pdb <- read.pdb("3dnd")

## Center of mass:
com(pdb)

## Center of mass of a selection
inds <- atom.select(pdb, chain="I")
com(pdb, inds)

## using XYZ Cartesian coordinates
xyz <- pdb$xyz[, inds$xyz]
com.xyz(xyz)

## with mass weighting
com.xyz(xyz, mass=atom2mass(pdb$atom[inds$atom, "elety"]) )



## Not run: 
## Unknown atom names
pdb <- read.pdb("3dnd")
inds <- atom.select(pdb, resid="LL2")
mycom <- com(pdb, inds, rescue=TRUE)
#warnings()


## Map atom names manually
pdb <- read.pdb("3RE0")
inds <- atom.select(pdb, resno=201)

myelety  <- data.frame(name = c("CL2","PT1","N1","N2"), symb = c("Cl","Pt","N","N"))
mymasses <- data.frame(symb = c("Cl","Pt"), mass = c(35.45, 195.08))
mycom    <- com(pdb, inds, elety.custom=myelety, mass.custom=mymasses)

## End(Not run)

Description

Do "and", "or", or "not" set operations between two or more atom selections made by atom.select

Usage

combine.select(sel1=NULL, sel2=NULL, ..., operator="AND", verbose=TRUE)

Arguments

Details

The value of operator should be one of following: (1) "AND", "and", or "&" for set intersect, (2) "OR", "or", "|", or "+" for set union, (3) "NOT", "not", "!", or "-" for set difference sel1 - sel2 - sel3 ....

Value

Returns a list of class "select" with components:

Author(s)

Xin-Qiu Yao

References

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

See Also

atom.select, as.select read.pdb, trim.pdb

Examples

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

## - Build atom selections to be operated
# Select C-alpha atoms of entire system
ca.global.inds <- atom.select(pdb, "calpha")

# Select C-beta atoms of entire protein
cb.global.inds <- atom.select(pdb, "protein", elety="CB")

# Select backbone atoms of entire system
bb.global.inds <- atom.select(pdb, "backbone")

# Select all atoms with residue number from 46 to 50
aa.local.inds <- atom.select(pdb, resno=46:50)


# Do set intersect:
# - Return C-alpha atoms with residue number from 46 to 50
ca.local.inds <- combine.select(ca.global.inds, aa.local.inds)
print( pdb$atom[ ca.local.inds$atom, ] )

# Do set subtract:
# - Return side-chain atoms with residue number from 46 to 50
sc.local.inds <- combine.select(aa.local.inds, bb.global.inds, operator="-")
print( pdb$atom[ sc.local.inds$atom, ] )

# Do set union:
# - Return C-alpha and side-chain atoms with residue number from 46 to 50
casc.local.inds <- combine.select(ca.local.inds, sc.local.inds, operator="+")
print( pdb$atom[ casc.local.inds$atom, ] )

# More than two selections:
# - Return side-chain atoms (but not C-beta) with residue number from 46 to 50
sc2.local.inds <- combine.select(aa.local.inds, bb.global.inds, cb.global.inds, operator="-")
print( pdb$atom[ sc2.local.inds$atom, ] )

Description

Find equivalent communities from two or more networks and re-assign colors to them in a consistent way across networks. A ‘new.membership’ vector is also generated for each network, which maps nodes to community IDs that are renumbered according to the community equivalency.

Usage

community.aln(x, ..., aln = NULL)

Arguments

Details

This function facilitates the inspection on the variance of the community partition in a group of similar networks. The original community numbering (and so the colors of communities in the output of plot.cna and vmd.cna) can be inconsistent across networks, i.e. equivalent communities may display different colors, impeding network comparison. The function calculates the dissimilarity between all communities and clusters communities with ‘hclust’ funciton. In each cluster, 0 or 1 community per network is included. The color attribute of communities is then re-assigned according to the clusters through all networks. In addition, a ‘new.membership’ vector is generated for each network, which mapps nodes to new community IDs that are numbered consistently across networks.

Value

Returns a list of updated cna objects.

See Also

cna, plot.cna, vmd.cna

Examples

# Needs MUSCLE installed - testing excluded
  if(check.utility("muscle")) {

    if (!requireNamespace("igraph", quietly = TRUE)) {
      message('Need igraph installed to run this example')
    } else {

    ## Fetch PDB files and split to chain A only PDB files
    ids <- c("1tnd_A", "1tag_A")
    files <- get.pdb(ids, split = TRUE, path = tempdir())
    
    ## Sequence Alignement
    pdbs <- pdbaln(files, outfile = tempfile())
    
    ## Normal mode analysis on aligned data
    modes <- nma(pdbs, rm.gaps=TRUE)
    
    ## Dynamic Cross Correlation Matrix
    cijs <- dccm(modes)$all.dccm
 
    ## Correlation Network
    nets <- cna(cijs, cutoff.cij=0.3)

    ## Align network communities
    nets.aln <- community.aln(nets)

    ## plot all-residue and coarse-grained (community) networks
    pdb <- pdbs2pdb(pdbs, inds=1, rm.gaps=TRUE)[[1]]
    op <- par(no.readonly=TRUE)

    # before alignment
    par(mar=c(0.1, 0.1, 0.1, 0.1), mfrow=c(2,2))
    invisible( lapply(nets, function(x) 
       plot(x, layout=layout.cna(x, pdb=pdb, k=3, full=TRUE)[, 1:2], 
               full=TRUE)) )
    invisible( lapply(nets, function(x) 
       plot(x, layout=layout.cna(x, pdb=pdb, k=3)[, 1:2])) )

    # after alignment
    par(mar=c(0.1, 0.1, 0.1, 0.1), mfrow=c(2,2))
    invisible( lapply(nets.aln, function(x) 
       plot(x, layout=layout.cna(x, pdb=pdb, k=3, full=TRUE)[, 1:2], 
               full=TRUE)) )
    invisible( lapply(nets.aln, function(x) 
       plot(x, layout=layout.cna(x, pdb=pdb, k=3)[, 1:2])) )

    par(op)     

    }
  }

Description

This function reconstructs the community tree of the community clustering analysis performed by the ‘cna’ function. It allows the user to explore different network community partitions.

Usage

community.tree(x, rescale=FALSE)

Arguments

Details

The input of this function should be a ‘cna’ class object containing ‘network’ and ‘communities’ attributes.

This function reconstructs the community residue memberships for each modularity value. The purpose is to facilitate inspection of alternate community partitioning points, which in practice often corresponds to a value close to the maximum of the modularity, but not the maximum value itself.

Value

Returns a list object that includes the following components:

Author(s)

Guido Scarabelli

See Also

cna, network.amendment, summary.cna

Examples

# PDB server connection required - testing excluded

if (!requireNamespace("igraph", quietly = TRUE)) {
   message('Need igraph installed to run this example')
} else {

###-- Build a CNA object
pdb <- read.pdb("4Q21")
modes <- nma(pdb)
cij <- dccm(modes)
net <- cna(cij, cutoff.cij=0.2)


##-- Reconstruct the community membership vector for each clustering step.
tree <- community.tree(net, rescale=TRUE)

## Plot modularity vs number of communities
plot( tree$num.of.comms, tree$modularity )

## Inspect the maximum modularity value partitioning
max.mod.ind <- which.max(tree$modularity)

## Number of communities (k) at max modularity
tree$num.of.comms[ max.mod.ind ]

## Membership vector at this partition point 
tree$tree[max.mod.ind,]

# Should be the same as that contained in the original CNA network object
net$communities$membership == tree$tree[max.mod.ind,]

# Inspect a new membership partitioning (at k=7)
memb.k7 <- tree$tree[ tree$num.of.comms == 7, ]

## Produce a new k=7 community network  
net.7 <- network.amendment(net, memb.k7)
plot(net.7, pdb)
#view.cna(net.7, trim.pdb(pdb, atom.select(pdb,"calpha")), launch=TRUE )

}

Description

Determines the consensus sequence for a given alignment at a given identity cutoff value.

Usage

consensus(alignment, cutoff = 0.6)

Arguments

Value

A vector containing the consensus sequence, where ‘-’ represents positions with no consensus (i.e. under the cutoff)

Author(s)

Barry Grant

References

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

See Also

read.fasta

Examples

#-- Read HIV protease alignment
aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))

# Generate consensus
con <- consensus(aln)
print(con$seq)

# Plot residue frequency matrix
##png(filename = "freq.png", width = 1500, height = 780)
col <- mono.colors(32)
aa  <- rev(rownames(con$freq))

image(x=1:ncol(con$freq),
      y=1:nrow(con$freq),
      z=as.matrix(rev(as.data.frame(t(con$freq)))),
      col=col, yaxt="n", xaxt="n",
      xlab="Alignment Position", ylab="Residue Type")

# Add consensus along the axis
axis(side=1, at=seq(0,length(con$seq),by=5))
axis(side=2, at=c(1:22), labels=aa)
axis(side=3, at=c(1:length(con$seq)), labels =con$seq)
axis(side=4, at=c(1:22), labels=aa)
grid(length(con$seq), length(aa))
box()

# Add consensus sequence
for(i in 1:length(con$seq)) {
  text(i, which(aa==con$seq[i]),con$seq[i],col="white")
}

# Add lines for residue type separation
abline(h=c(2.5,3.5, 4.5, 5.5, 3.5, 7.5, 9.5,
         12.5, 14.5, 16.5, 19.5), col="gray")

Description

Quantifies residue conservation in a given protein sequence alignment by calculating the degree of amino acid variability in each column of the alignment.

Usage

conserv(x, method = c("similarity","identity","entropy22","entropy10"),
        sub.matrix = c("bio3d", "blosum62", "pam30", "other"),
        matrix.file = NULL, normalize.matrix = TRUE)

Arguments

Details

To assess the level of sequence conservation at each position in an alignment, the “similarity”, “identity”, and “entropy” per position can be calculated.

The “similarity” is defined as the average of the similarity scores of all pairwise residue comparisons for that position in the alignment, where the similarity score between any two residues is the score value between those residues in the chosen substitution matrix “sub.matrix”.

The “identity” i.e. the preference for a specific amino acid to be found at a certain position, is assessed by averaging the identity scores resulting from all possible pairwise comparisons at that position in the alignment, where all identical residue comparisons are given a score of 1 and all other comparisons are given a value of 0.

“Entropy” is based on Shannons information entropy. See the entropy function for further details.

Note that the returned scores are normalized so that conserved columns score 1 and diverse columns score 0.

Value

Returns a numeric vector of scores

Note

Each of these conservation scores has particular strengths and weaknesses. For example, entropy elegantly captures amino acid diversity but fails to account for stereochemical similarities. By employing a combination of scores and taking the union of their respective conservation signals we expect to achieve a more comprehensive analysis of sequence conservation (Grant, 2007).

Author(s)

Barry Grant

References

Grant, B.J. et al. (2006) Bioinformatics 22, 2695–2696. Grant, B.J. et al. (2007) J. Mol. Biol. 368, 1231–1248.

See Also

read.fasta, read.fasta.pdb

Examples

## Read an example alignment
aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))

## Score conservation
conserv(x=aln$ali, method="similarity", sub.matrix="bio3d")
##conserv(x=aln$ali,method="entropy22", sub.matrix="other")

Description

Renumber and convert between CHARMM, Amber, Gromacs and Brookhaven PDB formats.

Usage

convert.pdb(pdb, type=c("original", "pdb", "charmm", "amber", "gromacs"),
                renumber = FALSE, first.resno = 1, first.eleno = 1,
                consecutive=TRUE, rm.h = TRUE, rm.wat = FALSE,
                verbose=TRUE)

Arguments

Details

Convert atom names and residue names, renumber atom and residue records, strip water and hydrogen atoms from pdb objects.

Format type can be one of “ori”, “pdb”, “charmm”, “amber” or “gromacs”.

Value

Returns a list of class "pdb", with the following components:

Note

For both atom and het list components the column names can be used as a convenient means of data access, namely: Atom serial number “eleno” , Atom type “elety”, Alternate location indicator “alt”, Residue name “resid”, Chain identifier “chain”, Residue sequence number “resno”, Code for insertion of residues “insert”, Orthogonal coordinates “x”, Orthogonal coordinates “y”, Orthogonal coordinates “z”, Occupancy “o”, and Temperature factor “b”. See examples for further details.

Author(s)

Barry Grant

References

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

For a description of PDB format (version3.3) see:
http://www.wwpdb.org/documentation/format33/v3.3.html.

See Also

atom.select, write.pdb, read.dcd, read.fasta.pdb, read.fasta

Examples

## Not run: 

# Read a PDB file
pdb <- read.pdb("4q21")
pdb
head( pdb$atom[pdb$calpha,"resno"] )

# Convert to CHARMM format
new <- convert.pdb(pdb, type="amber", renumber=TRUE, first.resno=22 )
head( new$atom[new$calpha,"resno"] )

# Write a PDB file
#write.pdb(new, file="tmp4amber.pdb")


## End(Not run)

Description

Find core positions that have the largest number of contact with neighboring residues.

Usage

core.cmap(pdbs, write.pdb = FALSE, outfile="core.pdb",
          cutoff = NULL, refine = FALSE, ncore = NULL, ...)

Arguments

Details

This function calculates eigenvector centrality of the weighted contact network built based on input structure data and uses it to determine the core positions.

In this context, core positions correspond to the most invariant C-alpha atom positions across an aligned set of protein structures. Traditionally one would use the core.find function to for their identification and then use these positions as the basis for improved structural superposition. This more recent function utilizes a much faster approach and is thus preferred in time sensitive applications such as shiny apps.

Value

Returns a list of class "select" containing ‘atom’ and ‘xyz’ indices.

Author(s)

Xin-Qiu Yao

References

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

See Also

core.find, read.fasta.pdb, fit.xyz

Examples

## Not run: 
##-- Generate a small kinesin alignment and read corresponding structures
pdbfiles <- get.pdb(c("1bg2","2ncd","1i6i","1i5s"), URLonly=TRUE)
pdbs <- pdbaln(pdbfiles)

##-- Find 'core' positions
core <- core.cmap(pdbs)
xyz <- pdbfit(pdbs, core, outpath="corefit_structures")

## End(Not run)

Description

Perform iterated rounds of structural superposition to identify the most invariant region in an aligned set of protein structures.

Usage

core.find(...)

## S3 method for class 'pdbs'
core.find(pdbs, shortcut = FALSE, rm.island = FALSE,
          verbose = TRUE, stop.at = 15, stop.vol = 0.5,
          write.pdbs = FALSE, outpath="core_pruned",
          ncore = 1, nseg.scale = 1, progress = NULL, ...)

## Default S3 method:
core.find(xyz, ...)

## S3 method for class 'pdb'
core.find(pdb, verbose=TRUE, ...)

Arguments

Details

This function attempts to iteratively refine an initial structural superposition determined from a multiple alignment. This involves iterated rounds of superposition, where at each round the position(s) displaying the largest differences is(are) excluded from the dataset. The spatial variation at each aligned position is determined from the eigenvalues of their Cartesian coordinates (i.e. the variance of the distribution along its three principal directions). Inspired by the work of Gerstein et al. (1991, 1995), an ellipsoid of variance is determined from the eigenvalues, and its volume is taken as a measure of structural variation at a given position.

Optional “core PDB files” containing core positions, upon which superposition is based, can be written to a location specified by outpath by setting write.pdbs=TRUE. These files are useful for examining the core filtering process by visualising them in a graphics program.

Value

Returns a list of class "core" with the following components:

Note

The relevance of the ‘core positions’ identified by this procedure is dependent upon the number of input structures and their diversity.

Author(s)

Barry Grant

References

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

Gerstein and Altman (1995) J. Mol. Biol. 251, 161–175.

Gerstein and Chothia (1991) J. Mol. Biol. 220, 133–149.

See Also

read.fasta.pdb, plot.core, fit.xyz

Examples

## Not run: 
##-- Generate a small kinesin alignment and read corresponding structures
pdbfiles <- get.pdb(c("1bg2","2ncd","1i6i","1i5s"), URLonly=TRUE)
pdbs <- pdbaln(pdbfiles)

##-- Find 'core' positions
core <- core.find(pdbs)
plot(core)

##-- Fit on these relatively invarient subset of positions 
#core.inds <- print(core, vol=1)
core.inds <- print(core, vol=0.5)
xyz <- pdbfit(pdbs, core.inds, outpath="corefit_structures")

##-- Compare to fitting on all equivalent positions
xyz2 <- pdbfit(pdbs)

## Note that overall RMSD will be higher but RMSF will
##  be lower in core regions, which may equate to a
##  'better fit' for certain applications
gaps <- gap.inspect(pdbs$xyz)
rmsd(xyz[,gaps$f.inds])
rmsd(xyz2[,gaps$f.inds])

plot(rmsf(xyz[,gaps$f.inds]), typ="l", col="blue", ylim=c(0,9))
points(rmsf(xyz2[,gaps$f.inds]), typ="l", col="red")

## End(Not run)

## Not run: 
##-- Run core.find() on a multimodel PDB file
pdb <- read.pdb('1d1d', multi=TRUE)
core <- core.find(pdb)

##-- Run core.find() on a trajectory
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)

## select calpha coords from a manageable number of frames
ca.ind <- atom.select(pdb, "calpha")$xyz
frames <- seq(1, nrow(trj), by=10)

core <- core.find( trj[frames, ca.ind], write.pdbs=TRUE )

## have a look at the various cores "vmd -m core_pruned/*.pdb"

## Lets use a 6A^3 core cutoff
inds <- print(core, vol=6)
write.pdb(xyz=pdb$xyz[inds$xyz],resno=pdb$atom[inds$atom,"resno"], file="core.pdb")


##- Fit trj onto starting structure based on core indices
xyz <- fit.xyz( fixed = pdb$xyz,
               mobile = trj,
               fixed.inds  = inds$xyz,
               mobile.inds = inds$xyz)

##write.pdb(pdb=pdb, xyz=xyz, file="new_trj.pdb")
##write.ncdf(xyz, "new_trj.nc")


## End(Not run)