Using neighbor segmentation in voxelspace for fast and consistent spatial and temporal segmentation.
This software has been developed to allow for a higher selection syntax than atom and or residue index, such as abstract complex particles (e.g. lipid monolayers). MDVoxelSegmentation combines multiple layers of neighbor segmentation with a voxel mask, and makes tracking segments over time possible at high quality.
This code was published in an open access article.
- Open software: Apache 2 license
- v1.1.7 is a stable build and segmentation should be of high quality. The code is usable for high throughput with minimal effort with optimization. By default no more than a GRO and XTC (or equivalent) are required for successful segmentation of Martini systems.
- Voxel based neighbour segmentation under all periodic boundary conditions
- Fast contour segmentation
- Compatible with most MD file formats due to its tight link to MDAnalysis
- Consistent segmentation over time on trajectories
- Compatible with VMD using standard visualization files (python compiled VMD)
- A wide range of examples (bottom of this readme)
- Membrane leaflet assignment of lipids of most topologies
- Bilayers
- Vesicles
- (Inverted) hexagonal
- Cubic
- Membrane tethers
- Stalks
- Pores
- Complex lipids formulations including cholesterol
- Proteins
- Up to millions of beads per frame (possibly billions)
experimental features
- Make non-bonded densities whole over PBC (mdvwhole, currently a separate repository that will be merged upon completion). This feature is stable, tested and useful, but its integration in MDVoxelSegmentation has to be completed.
@article{Bruininks2021,
doi = {10.1021/acs.jctc.1c00446},
url = {https://doi.org/10.1021/acs.jctc.1c00446},
year = {2021},
month = oct,
publisher = {American Chemical Society ({ACS})},
volume = {17},
number = {12},
pages = {7873--7885},
author = {Bart M. H. Bruininks and Albert S. Thie and Paulo C. T. Souza and Tsjerk A. Wassenaar and Shirin Faraji and Siewert J. Marrink},
title = {Sequential Voxel-Based Leaflet Segmentation of Complex Lipid Morphologies},
journal = {Journal of Chemical Theory and Computation}
}
pip install mdvoxelsegmentation
mdvseg -h (run in terminal)
To perform default segmentation on a GRO and XTC file containing a coarse grain Martini system, you have to specify the GRO and XTC file. The final segmentation assignment will be written to clusters.npy. This file can be used using numpy in python to perform the required analysis.
mdvseg -f path_to_your.gro -x path_to_your.xtc
MDAnalysis will probably throw some warnings stating that it cannot estimate the masses for you coarse grain particles. This cannot be suppressed but is harmless. A useful graph of your segmentation is automatically generated in your folder. The plotting can be manually controlled by altering the plot script, which is also place in your active folder:
python plotting.py # alter this and rerun
The graph is created by default and the plotter is written to the folder. Therefore adjusting the plotting script to make exactly what you need should be rather straight forward.
It is mainly the force segmentation flag (-fs) and it associated recursion depth (-rd) which have a big impact on performance. Turning force segmentation off (-fs 0) is fine if perfect final quality is not needed. By default mdvseg generates its own selections.inp which should cover all basic lipids in Martini. However, if some definitions are missing, you can always manually add them to the selections.inp (mdvseg does NOT overwrite an already present selection.inp).
For visualization with VMD you need to make sure that you have a version of VMD compiled against a python version supporting numpy. A compatible VMD compilation will be distributed in the future, for now you have to figure this out yourself, I did ask the developers to support anyone asking for such compilation and they said yes! So just send an e-mail to the VMD mailing list if you would need it. If you are running Ubuntu 20.04, you are in luck and I can supply you with a custom version of VMD 1.9.4 which should be relatively easy to install. Open a ticket and I'll see what I can do for you.
If you have the right flavor of VMD, all you need to type next is:
vmd -e vmd_clusters_visualization.vmd
The first 32 segments will automatically be assigned a color and material/style. They can be used to make selections using 'user your_segment' in the VMD selection syntax. These representations should automatically be set to update every frame. By typing 'hide' in the VMD terminal, we can easily turn off all segment representations. Segment 0 always contains everything which was not assigned a segment and is hidden. 'user' 32 always shows segment 32 to 1000, to show all segments which might have a very high index. The downside is that all segments from 32 onward have the same color.
Using MDVoxelSegmentation on coarse grain Martini lipid/protein systems should work without needing much prior knowledge. However, to make the most out of the created clusters.npy it is useful to know some python (numpy, MDAnalysis, Matplotlib). If you are working with atomistic systems and have to specify your own headgroups/linkers/tails, you need to known what the relevant names are from your PDB/GRO and make your own selection entries in the selections.inp. The selections.inp uses the MDAnalysis selection syntax (very close to the VMD selection syntax). Below are some basic lines of code to help you on your way with using the segmentation data. First we will give an example for some basic plotting, fetching segment information for residues and/or complete selections, followed by an example for an atomistic CHARMM selections.inp for DOPE lipids and how to segments it. Finally we show how to use MDVoxelSegmentation for non-amphipathic density segmentation (normal segmentation).
A basic python example to plot the number of segments over time
## Importing numpy and matplotlib.
import numpy as np
import matplotlib.pyplot as plt
## Loading the segmentation data.
segments_over_time = np.load('clusters.npy')
## Calculating the amount of segments in each frame.
# Make an empty array which has one int32 for every frame.
segments_per_frame = np.zeros(segments_over_time.shape[0], dtype='int32')
# Fill the array with the amount of non-zero segments in each frame.
segments_per_frame[:] = [len(np.unique(frame)) - 1 for frame in segments_over_time]
## Plotting the results.
# Making an empty plot.
fig, ax = plt.subplots()
# Adding required data to plot.
ax.plot(segments_per_frame)
# Set ticks to a sensible regime.
start, end = [round(limit) for limit in ax.get_ylim()]
ax.yaxis.set_ticks(np.arange(start, end + 1, 1))
# Add labels to axes.
ax.set_xlabel('Frame count')
ax.set_ylabel('Number of segments')
# Save the plot.
fig.savefig('amount_of_segments_over_time.png', dpi=300)
# Usually people don't like it if you pop windows, however
# if you would like to automatically show the result uncomment
# the following line.
#fig.show()Basic fetching of segment ID for residues
class Container():
"A simple container class for our universe and segmentation data."
def __init__(self, universe, segmentation):
self.u = universe
self.segmentation = segmentation
def get_segment_from_resid(self, resid):
"""Returns the residue segment id in the current frame."""
residue_index = self.u.residues[resid].atoms[0].ix
current_frame = self.u.trajectory.frame
residue_segment = self.segmentation[current_frame, residue_index]
return residue_segment
def get_segments_from_selection(self, selection):
"""Returns an array of lipid segments for the given selection
in the current frame. The selection should adhere
to the basic mda selection syntax."""
selection = self.u.select_atoms(selection)
resids = selection.residues.ix
segments = [container.get_segment_from_resid(resid)
for resid in resids]
return np.asarray(segments, dtype=int)
# File paths
gro = 'your.gro'
xtc = 'your.xtc'
segmentation_data = 'clusters.npy'
# Creating universe and loading segmentation
u = mda.Universe(gro, xtc)
segmentation = np.load(segmentation_data)
# Create our container
container = Container(u, segmentation)
# Segment from resid
resid1_segmentID = container.get_segment_from_resid(1)
# Segment from selection
all_cholesterol_segmentsIDs = container.get_segments_from_selection('resname CHOL')An example file for flip-flop analysis is added under 'mdvoxelsegmentation/templates/lipid_flip-flop.ipynb'
An atomistic segmentation example for DOPE lipids with the CHARMM force field
$ vi selections.inp ## Create an empty `selections.inp` and add the following lines, the selection ## syntax should always be one line and directly follow its header description. # It is not bad to include the linkers in the headgroups as well, but this is often # not important. [charmm_heads] (name N P C12 C11 O11 O12 O13 O14) [charmm_linkers] (name C1 C2 O21 C21 C3 O31 C31) # Not adding the first carbons of the tail can improve performance, but again, # this is usually not important. [charmm_tails] (name C22 C23 C24 C25 C26 C27 C28 C29 C210 C211 C212 C213 C214 C215 C216 C217 C218 C32 C33 C34 C35 C36 C37 C38 C39 C310 C311 C312 C313 C314 C315 C316 C317 C318) ## Run the mdvseg, hyper resolution can be turned off for there are more than ## enough particles due to the atomistic resolution! $ mdvseg -f your.gro -x your.xtc -hg charmm_heads -lg charmm_linkers -tg charmm_tails -hres 0
Segmenting non-amphipathic densities (normal segmentation)
# Define the density selection $ vi selections.inp [density] (not resname W WF ION) # Basic non-solvent Martini density [none] False # Run mdvseg with the selections.inp without an exclusion group nor forced segmenation. mdvseg -f your.gro -x your.xtc -hg density -tg none -eg none -fs 0
We would love to see how you used MDVoxelSegmentation in your projects and are always happy to see cool systems and screenshots. Just upload them to an issue of this repository. There is also place for comments on usability and bugs.
If you are interested in joining this project after its initial release just post an issue, or better yet, send an email to [email protected]. We are currently still filled with ideas to be implemented and all hands are welcome. Some of the open topics are in the issues. Anybody who contributes for a fair share will off course be included in future publications.
Segmentation of the inverted hexagonal phase with four inner channels connected to a bilayer with a fusion stalk
Inside the channels is a fragment of dsDNA. The leaflet segmentation was performed using a resolution of 0.5 and hyperesolution turned on. This to allow for the correct segmentation of the tight geometry of the channels in coarse grain data (Martini, we used hyper resolution for all CG data!), also force segmentation was turned on to have (almost?) every lipid assigned up to a distance of 2 nm.
Leaflet segmentation of a complex plasmamembrane tether
The two leaflets of the plasmamembrane are clearly assigned correctly and depicted as a transparent surface. The cholesterol inside the two leaflets is drawn in VDW spheres and their headgroups have a slightly altering color. All cholesterol seems to be assigned correctly. Segmentation was performed with a 0.5 nm resolution and iterative forced segmentation to assign the diving cholesterol up to a distance of 2 nm.
Leaflet segmentation of a plasma membrane including multiple proteins
Less than 30 lipids remain unassigned of the roughly 1 million present. The leaflet assignment seemed to have worked correctly. For segmentation a resolution of 0.5 nm and iterative forced segmentation within 2 nm was used. The protein was used as exclusion to prevent them acting as pores in our segmentation. In total 1.3 millions beads were segmented in 30 minutes on a desktop. Mainly the force segmentation to assign all diving leaflets took a while. Keep in mind that this will never change the amount of segments present, so forced segmentation could be skipped in many situations. This is also only making use of a single core (a single frame cannot be hypterthreaded in the current code).
A collection of notoriously hard bilayer problems for segmentation
For the cholesterol flip-flopping we use non iterative forced segmentation with a cutoff of 1.5 nm to act as a deadzone of 1 nm (A, B, C; recursion depth set to 1). We see that intercalating close contact leaflets do not cause faulty segmentation (D, E). Pores are also handled correctly and the minimum pore size at a resolution of 0.5 nm is actually the pore itself (F/G). If the pore is only a water channel, but the lipids do not reorient, its not considered a pore. Since the leaflets are not even continuous. In short we are able to detect all toroidal pores in a membrane. Water pores are a different game which we might solve in the future with a similar set based approach (ohh yhea we got something nice brewing, if only we had time :D).
Luis S. Mayorga and Diego Masone wrote some heartwarming words in their manuscript (The Secret Ballet Inside Multivesicular Bodies) work about MDVoxelSegmentation. We cannot express how motivating it is for us to read these words of appreciation. If you are reading this README, you are probably into interesting looking lipid morphologies and we can highly recommend reading their work. In their manuscript they write the following in the Results section:
The key feature in this work has been the application to the model autophagosome trajectories of a convenient computational tool by Bruininks et al. (18) that breaks down very complex lipid superstructures in clearly delimited monolayers, otherwise undistinguishable by normal means. This powerful tool tracks the individual lipid leaflets (named segments) in a molecular dynamics simulation unveiling topological changes produced by spontaneous fusion and fission events. Even for the trained eye, without segmentation it is difficult, if not impossible, to figure out the shapes that are confined in the structures that result from the dynamics of the model autophagosome (Figure 1). Despite the use of powerful tools that allow for cross-sectional visualizations (e.g., VMD (35)), inner vesicles are highly irregular and cannot be understood in most cases. Therefore, leaflet segmentation of the model autophagosome has revealed to be a game-changing technique for identifying lipid superstructures and their interactions.
Bart M. H. Bruininks, Albert Thie, Paulo C. T. de Souza, Tsjerk A. Wassenaar, Shirin Faraji & Siewert J. Marrink