BioStructures documentation

The BioStructures.jl package provides functionality to manipulate macromolecular structures, and in particular to read and write Protein Data Bank (PDB), mmCIF and MMTF files. It is designed to be used for standard structural analysis tasks, as well as acting as a platform on which others can build to create more specific tools.

It compares favourably in terms of performance to other PDB parsers - see some benchmarks online and the benchmark suite. The PDB, mmCIF and MMTF parsers currently read in the whole PDB without explicit errors (with a handful of known exceptions). Help on individual functions can be found in the BioStructures API section or by using ?function_name from within Julia.

Basics

To download a PDB file:

using BioStructures

# Stored in the current working directory by default
downloadpdb("1EN2")

To parse a PDB file into a Structure-Model-Chain-Residue-Atom framework:

julia> struc = read("/path/to/pdb/file.pdb", PDB)
ProteinStructure 1EN2.pdb with 1 models, 1 chains (A), 85 residues, 754 atoms

mmCIF files can be read into the same data structure with read("/path/to/cif/file.cif", MMCIF). The keyword argument gzip, default false, determines if the file is gzipped. If you want to read an mmCIF file into a dictionary to query yourself (e.g. to access metadata fields), use MMCIFDict:

julia> mmcif_dict = MMCIFDict("/path/to/cif/file.cif")
mmCIF dictionary with 716 fields

julia> mmcif_dict["_entity_src_nat.common_name"]
1-element Array{String,1}:
 "great nettle"

A MMCIFDict can be accessed in similar ways to a standard dictionary, and if necessary the underlying dictionary of MMCIFDict d can be accessed with d.dict. Note that the values of the dictionary are always an Array{String,1}, even if only one value was read in or the data is numerical.

MMTF files can be read into the same data structure with read("/path/to/mmtf/file.mmtf", MMTF). The keyword argument gzip, default false, determines if the file is gzipped. In a similar manner to mmCIF dictionaries, a MMTF file can be read into a dictionary with MMTFDict. The values of the dictionary are a variety of types depending on the MMTF specification. To convert a MMCIFDict or MMTFDict to the Structure-Model-Chain-Residue-Atom framework, use the ProteinStructure constructor, e.g. ProteinStructure(mmcif_dict).

The elements of struc can be accessed as follows:

You can use begin and end to access the first and last elements, for example struc[1][begin] gets the first chain in model 1. Disordered atoms are stored in a DisorderedAtom container but calls fall back to the default atom, so disorder can be ignored if you are not interested in it. Disordered residues (i.e. point mutations with different residue names) are stored in a DisorderedResidue container. The idea is that disorder will only bother you if you want it to. See the Biopython discussion for more.

Properties can be retrieved as follows:

The strip keyword argument determines whether surrounding whitespace is stripped for atomname, element, charge, resname and atomnames (default true).

The coordinates of an atom can be changed using x!, y!, z! and coords!. The chain ID of a chain or residue can be changed using chainid!. Currently these are the only setter functions available.

Manipulating structures

Elements can be looped over to reveal the sub-elements in the correct order:

for mod in struc
    for ch in mod
        for res in ch
            for at in res
                # Do something
            end
        end
    end
end

Models are ordered numerically; chains are ordered by chain ID character ordering, except the empty chain is last; residues are ordered by residue number and insertion code with hetero residues after standard residues; atoms are ordered by atom serial. Since the ordering of elements is defined you can use the sort function. For example sort(res) sorts a list of residues as described above, or sort(res, by=resname) will sort them alphabetically by residue name.

collect can be used to get arrays of sub-elements. collectatoms, collectresidues, collectchains and collectmodels return arrays of a particular type from a structural element or element array. Since most operations should use a single version of an atom or residue, disordered entities are not expanded by default and only one entity is present in the array. This can be changed by setting expand_disordered to true in collectatoms or collectresidues.

Selectors are functions passed as additional arguments to these functions. Only elements that return true when passed to all the selector are retained. For example:

The selectors available are:

To create a new atomnameselector, resnameselector or sscodeselector:

cdeltaselector(at::AbstractAtom) = atomnameselector(at, ["CD"])

It is easy to define your own atom, residue, chain or model selectors. The below will collect all atoms with x coordinate less than 0:

xselector(at) = x(at) < 0
collectatoms(struc, xselector)

Selectors can be inverted with !, e.g. collectatoms(struc, !xselector). Alternatively, you can use an anonymous function:

collectatoms(struc, at -> x(at) < 0)

countatoms, countresidues, countchains and countmodels can be used to count elements with the same selector API. For example:

julia> countatoms(struc)
754

julia> countatoms(struc, calphaselector)
85

julia> countresidues(struc, standardselector)
85

julia> countatoms(struc, expand_disordered=true)
819

The amino acid sequence of a protein can be retrieved by passing an element to LongAA with optional residue selectors:

julia> LongAA(struc['A'], standardselector)
85aa Amino Acid Sequence:
RCGSQGGGSTCPGLRCCSIWGWCGDSEPYCGRTCENKCW…RCGAAVGNPPCGQDRCCSVHGWCGGGNDYCSGGNCQYRC

The gaps keyword argument determines whether to add gaps to the sequence based on missing residue numbers (default true). threeletter_to_aa provides a lookup table of amino acid codes should that be required. See BioSequences.jl and BioAlignments.jl for more on how to deal with sequences. LongAA is an alias for LongSequence{AminoAcidAlphabet}. For example, to see the alignment of CDK1 and CDK2 (this example also makes use of Julia's broadcasting):

julia> struc1, struc2 = retrievepdb.(["4YC6", "1HCL"])
2-element Array{ProteinStructure,1}:
 ProteinStructure 4YC6.pdb with 1 models, 8 chains (A,B,C,D,E,F,G,H), 1420 residues, 12271 atoms
 ProteinStructure 1HCL.pdb with 1 models, 1 chains (A), 294 residues, 2546 atoms

julia> seq1, seq2 = LongAA.([struc1["A"], struc2["A"]], standardselector, gaps=false)
2-element Vector{LongAA}:
 MEDYTKIEKIGEGTYGVVYKGRHKTTGQVVAMKKIRLES…SHVKNLDENGLDLLSKMLIYDPAKRISGKMALNHPYFND
 MENFQKVEKIGEGTYGVVYKARNKLTGEVVALKKIRTEG…RSLLSQMLHYDPNKRISAKAALAHPFFQDVTKPVPHLRL

julia> using BioAlignments

julia> scoremodel = AffineGapScoreModel(BLOSUM62, gap_open=-10, gap_extend=-1);

julia> alres = pairalign(GlobalAlignment(), seq1, seq2, scoremodel)
PairwiseAlignmentResult{Int64, LongAA, LongAA}:
  score: 945
  seq:   1 MEDYTKIEKIGEGTYGVVYKGRHKTTGQVVAMKKIRLESEEEGVPSTAIREISLLKELRH  60
           ||   | ||||||||||||| | | || ||| |||| |    |||||||||||||||| |
  ref:   1 MENFQKVEKIGEGTYGVVYKARNKLTGEVVALKKIRTE----GVPSTAIREISLLKELNH  56

  seq:  61 PNIVSLQDVLMQDSRLYLIFEFLSMDLKKYLD-SIPPGQYMDSSLVKSYLYQILQGIVFC 119
           |||| | ||      ||| ||||  ||||  | |   |      | |||| | |||  ||
  ref:  57 PNIVKLLDVIHTENKLYLVFEFLHQDLKKFMDASALTG--IPLPLIKSYLFQLLQGLAFC 114

  seq: 120 HSRRVLHRDLKPQNLLIDDKGTIKLADFGLARAFGV----YTHEVVTLWYRSPEVLLGSA 175
           || ||||||||||||||   | ||||||||||||||    ||||||||||| || |||  
  ref: 115 HSHRVLHRDLKPQNLLINTEGAIKLADFGLARAFGVPVRTYTHEVVTLWYRAPEILLGCK 174

  seq: 176 RYSTPVDIWSIGTIFAELATKKPLFHGDSEIDQLFRIFRALGTPNNEVWPEVESLQDYKN 235
            ||| ||||| | ||||  |   || ||||||||||||| ||||   ||| | |  ||| 
  ref: 175 YYSTAVDIWSLGCIFAEMVTRRALFPGDSEIDQLFRIFRTLGTPDEVVWPGVTSMPDYKP 234

  seq: 236 TFPKWKPGSLASHVKNLDENGLDLLSKMLIYDPAKRISGKMALNHPYFND---------- 285
            ||||        |  ||| |  ||| || ||| |||| | || || | |          
  ref: 235 SFPKWARQDFSKVVPPLDEDGRSLLSQMLHYDPNKRISAKAALAHPFFQDVTKPVPHLRL 294

In fact, pairalign is extended to carry out the above steps and return the alignment by calling pairalign(struc1["A"], struc2["A"], standardselector) in this case. scoremodel and aligntype are keyword arguments with the defaults shown above.

Spatial calculations

Various functions are provided to calculate spatial quantities for proteins:

The omegaangle, phiangle and psiangle functions can take either a pair of residues or a chain and a position. The omegaangle and phiangle functions measure the angle between the residue at the given index and the one before. The psiangle function measures between the given index and the one after. For example:

julia> distance(struc['A'][10], struc['A'][20])
10.782158874733762

julia> rad2deg(bondangle(struc['A'][50]["N"], struc['A'][50]["CA"], struc['A'][50]["C"]))
110.77765846083398

julia> rad2deg(dihedralangle(struc['A'][50]["N"], struc['A'][50]["CA"], struc['A'][50]["C"], struc['A'][51]["N"]))
-177.38288114072924

julia> rad2deg(psiangle(struc['A'][50], struc['A'][51]))
-177.38288114072924

julia> rad2deg(psiangle(struc['A'], 50))
-177.38288114072924

ContactMap takes in a structural element or a list, such as a Chain or Vector{Atom}, and returns a ContactMap object showing the contacts between the elements for a specified distance. ContactMap can also be given two structural elements as arguments, in which case a non-symmetrical 2D array is returned showing contacts between the elements. The underlying BitArray{2} for ContactMap contacts can be accessed with contacts.data if required.

julia> contacts = ContactMap(collectatoms(struc['A'], cbetaselector), 8.0)
Contact map of size (85, 85)

A plot recipe is defined for this so it can shown with Plots.jl:

using Plots
plot(contacts)

For a quick, text-based representation of a ContactMap use showcontactmap.

DistanceMap works in an analogous way to ContactMap and gives a map of the distances. It can also be plotted:

dists = DistanceMap(collectatoms(struc['A'], cbetaselector))
using Plots
plot(dists)

Structural elements can be superimposed, and superposition-dependent properties such as the RMSD can be calculated. To carry out superimposition, BioStructures.jl carries out a sequence alignment and superimposes aligned residues using the Kabsch algorithm. For example:

# Change the coordinates of element 1 to superimpose it onto element 2
# Do sequence alignment with standard residues and calculate the transformation with Cα atoms (the default)
superimpose!(el1, el2, standardselector)

# The transformation object for the above superimposition
Transformation(el1, el2, standardselector)

# Calculate the transformation with backbone atoms
superimpose!(el1, el2, standardselector, alignatoms=backboneselector)

# Calculate RMSD on Cα atoms (the default) after superimposition
rmsd(el1, el2, standardselector)

# Superimpose based on backbone atoms and calculate RMSD based on Cβ atoms
rmsd(el1, el2, standardselector, alignatoms=backboneselector, rmsdatoms=cbetaselector)

# Do not do a superimposition - assumes the elements are already superimposed
rmsd(el1, el2, standardselector, superimpose=false)

displacements is used in a similar way to rmsd but returns the vector of distances for each superimposed atom. dispatoms selects the atoms to calculate the displacements on.

These transformation functions may be useful beyond the context of protein structures. For example, Transformation(c1, c2) calculates the transformation to map one set of coordinates to another. The coordinate sets must be the same size and have the number of dimensions in the first axis and the number of points in the second axis.

The contacting elements in a molecular structure form a graph, and this can be retrieved using MetaGraph. This extends MetaGraph from MetaGraphs.jl, allowing you to use all the graph analysis tools in Graphs.jl. For example:

julia> using Graphs, MetaGraphs

julia> mg = MetaGraph(collectatoms(struc["A"], cbetaselector), 8.0)
{85, 423} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)

julia> nv(mg)
85

julia> ne(mg)
423

julia> get_prop(mg, :contactdist)
8.0

julia> mg[10, :element]
Atom CB with serial 71, coordinates [-3.766, 4.031, 23.526]

See the Graphs docs for details on how to calculate properties such as shortest paths, centrality measures, community detection and more. Similar to ContactMap, contacts are found between any element type passed in. So if you wanted the graph of chain contacts in a protein complex you could give a Model as the first argument.

Assigning secondary structure

The secondary structure code of a residue or atom can be accessed after assigning the secondary structure using DSSP or STRIDE. To assign secondary structure when reading the structure:

# Assign secondary structure using DSSP
read("/path/to/pdb/file.pdb", PDB, run_dssp=true)

# Assign secondary structure using STRIDE
read("/path/to/pdb/file.pdb", PDB, run_strude=true)

rundssp!, runstride!, rundssp and runstride can also be used to assign secondary structure to a ProteinStructure or Model:

rundssp!(struc)
runstride!(struc)

The assignment process may fail if the structure is too large, since we use an intermediate PDB file where the atom serial cannot exceed 99999 and the chain ID must be a single character.

To get access to the secondary structure code of a residue or atom as a Char:

sscode(res)
sscode(at)

The secondary structure code of a residue can be changed using sscode!.

rundssp and runstride can also be run directly on structure files:

rundssp("/path/to/pdb/file.pdb", "out.dssp") # Also works with mmCIF files
runstride("/path/to/pdb/file.pdb", "out.stride")

See the documentation for DSSP_jll and STRIDE_jll, and also ProteinSecondaryStructures.jl, for other ways to run these programs.

Downloading PDB files

To download a PDB file to a specified directory:

downloadpdb("1EN2", dir="path/to/pdb/directory")

To download multiple PDB files to a specified directory:

downloadpdb(["1EN2", "1ALW", "1AKE"], dir="path/to/pdb/directory")

To download a PDB file in PDB, XML, mmCIF or MMTF format use the format argument:

# To get mmCIF
downloadpdb("1ALW", dir="path/to/pdb/directory", format=MMCIF)

# To get XML
downloadpdb("1ALW", dir="path/to/pdb/directory", format=PDBXML)

To apply a function to a downloaded file and delete the file afterwards:

downloadpdb(f, "1ALW")

Or, using Julia's do syntax:

downloadpdb("1ALW") do fp
    s = read(fp, PDB)
    # Do something
end

Note that some PDB entries, e.g. large viral assemblies, are not available as PDB format files. In this case download the mmCIF file or MMTF file instead.

Reading PDB files

To parse an existing PDB file into a Structure-Model-Chain-Residue-Atom framework:

julia> struc = read("/path/to/pdb/file.pdb", PDB)
ProteinStructure 1EN2.pdb with 1 models, 1 chains (A), 85 residues, 754 atoms

Read a mmCIF/MMTF file instead by replacing PDB with MMCIF/MMTF. Various options can be set through optional keyword arguments when parsing PDB/mmCIF/MMTF files:

Use retrievepdb to download and parse a PDB file into a Structure-Model-Chain-Residue-Atom framework in a single line:

julia> struc = retrievepdb("1ALW", dir="path/to/pdb/directory")
INFO: Downloading PDB: 1ALW
ProteinStructure 1ALW.pdb with 1 models, 2 chains (A,B), 346 residues, 2928 atoms

If you prefer to work with data frames rather than the data structures in BioStructures, the DataFrame constructor from DataFrames.jl has been extended to construct relevant data frames from lists of atoms or residues:

julia> using DataFrames

julia> df = DataFrame(collectatoms(struc));

julia> first(df, 3)
3×17 DataFrame. Omitted printing of 5 columns
│ Row │ ishetero │ serial │ atomname │ altlocid │ resname │ chainid │ resnumber │ inscode │ x       │ y       │ z       │ occupancy │
│     │ Bool     │ Int64  │ String   │ Char     │ String  │ String  │ Int64     │ Char    │ Float64 │ Float64 │ Float64 │ Float64   │
├─────┼──────────┼────────┼──────────┼──────────┼─────────┼─────────┼───────────┼─────────┼─────────┼─────────┼─────────┼───────────┤
│ 1   │ false    │ 1      │ N        │ ' '      │ GLU     │ A       │ 94        │ ' '     │ 15.637  │ -47.066 │ 18.179  │ 1.0       │
│ 2   │ false    │ 2      │ CA       │ ' '      │ GLU     │ A       │ 94        │ ' '     │ 14.439  │ -47.978 │ 18.304  │ 1.0       │
│ 3   │ false    │ 3      │ C        │ ' '      │ GLU     │ A       │ 94        │ ' '     │ 14.141  │ -48.183 │ 19.736  │ 1.0       │

julia> df = DataFrame(collectresidues(struc));

julia> first(df, 3)
3×8 DataFrame
│ Row │ ishetero │ resname │ chainid │ resnumber │ inscode │ countatoms │ modelnumber │ isdisorderedres │
│     │ Bool     │ String  │ String  │ Int64     │ Char    │ Int64      │ Int64       │ Bool            │
├─────┼──────────┼─────────┼─────────┼───────────┼─────────┼────────────┼─────────────┼─────────────────┤
│ 1   │ false    │ GLU     │ A       │ 94        │ ' '     │ 9          │ 1           │ false           │
│ 2   │ false    │ GLU     │ A       │ 95        │ ' '     │ 9          │ 1           │ false           │
│ 3   │ false    │ VAL     │ A       │ 96        │ ' '     │ 7          │ 1           │ false           │

As with file writing disordered entities are expanded by default but this can be changed by setting expand_disordered to false.

Reading multiple mmCIF data blocks

You can read and write files containing multiple mmCIF data blocks (equivalent to a MMCIFDict in this package) with the readmultimmcif and writemultimmcif functions. An example of such a file is the PDB's chemical component dictionary.

julia> ccd = readmultimmcif("components.cif.gz"; gzip=true);

julia> ccd["2W4"]
mmCIF dictionary with 64 fields

Writing PDB files

PDB format files can be written:

writepdb("1EN2_out.pdb", struc)

Any element type can be given as input to writepdb. The first argument can also be a stream. Atom selectors can also be given as additional arguments:

# Only backbone atoms are written out
writepdb("1EN2_out.pdb", struc, backboneselector)

To write mmCIF format files, use the writemmcif function with similar arguments. The gzip keyword argument, default false, determines whether to gzip the written file. A MMCIFDict can also be written using writemmcif:

writemmcif("1EN2_out.dic", mmcif_dict)

To write out a MMTF file, use the writemmtf function with any element type or a MMTFDict as an argument. The gzip keyword argument, default false, determines whether to gzip the written file.

Unlike for the collection functions, expand_disordered is set to true when writing files as it is usually desirable to retain all entities. Set expand_disordered to false to not write out more than one atom or residue at each location.

Multi-character chain IDs can be written to mmCIF and MMTF files but will throw an error when written to a PDB file as the PDB file format only has one character allocated to the chain ID.

If you want the PDB record line for an AbstractAtom, use pdbline. For example:

julia> pdbline(at)
"HETATM  101  C  A  X B  20      10.500  20.123  -5.123  0.50 50.13           C1+"

If you want to generate a PDB record line from values directly, do so using an AtomRecord:

julia> pdbline(AtomRecord(false, 669, "CA", ' ', "ILE", "A", 90, ' ', [31.743, 33.11, 31.221], 1.00, 25.76, "C", ""))
"ATOM    669  CA  ILE A  90      31.743  33.110  31.221  1.00 25.76           C  "

This can be useful when writing PDB files from your own data structures.

RCSB PDB utility functions

To get the list of all PDB entries:

l = pdbentrylist()

To download the entire RCSB PDB database in your preferred file format:

downloadentirepdb(dir="path/to/pdb/directory", format=MMTF)

This operation takes a lot of disk space and time to complete, depending on internet connection.

To update your local PDB directory based on the weekly status list of new, modified and obsolete PDB files from the RCSB server:

updatelocalpdb(dir="path/to/pdb/directory", format=MMTF)

Obsolete PDB files are stored in the auto-generated obsolete directory inside the specified local PDB directory.

To maintain a local copy of the entire RCSB PDB database, run the downloadentirepdb function once to download all PDB files and set up a CRON job or similar to run updatelocalpdb function once a week to keep the local PDB directory up to date with the RCSB server.

There are a few more functions that may be useful:

Visualising structures

The Bio3DView.jl package can be used to visualise molecular structures. For example:

using Bio3DView
using Blink
viewpdb("1CRN")

Use the mouse to move and zoom the images. You can view BioStructures data types:

struc = retrievepdb("1AKE")
viewstruc(struc['A'], surface=Surface(Dict("colorscheme"=> "greenCarbon")))

See the Bio3DView.jl tutorial for more information. BioMakie.jl can also be used to visualise BioStructures objects.

Related software

Other packages in the Julia ecosystem that deal with structural bioinformatics or related fields include:

• MIToS.jl - protein sequence and structure analysis.
• PDBTools.jl - read and write PDB files.
• PdbTool.jl - read and work with PDB files.
• Chemfiles.jl - read and write various chemistry trajectory files.
• Bio3DView.jl - view molecular structures (see Visualising structures).
• BioMakie.jl - view molecular structures and sequence alignments.
• ProteinSecondaryStructures.jl - parses DSSP and STRIDE outputs.
• MMTF.jl - read and write MMTF files. BioStructures.jl builds on top of MMTF.jl.
• ProteinEnsembles.jl - model ensembles of protein structures.
• Molly.jl - molecular dynamics of biomolecules.