, Stefano Motta [aut,
cre, cph] , Erik
Lindahl [ctb] (Author of xdrfile C library), David van der
Spoel [ctb] (Author of xdrfile C library)| Title: | Self Organising Maps for the Analysis of Molecular Dynamics Data |
|---|---|
| Description: | Processes data from Molecular Dynamics simulations using Self Organising Maps. Features include the ability to read different input formats. Trajectories can be analysed to identify groups of important frames. Output visualisation can be generated for maps and pathways. Methodological details can be found in Motta S et al (2022) <doi:10.1021/acs.jctc.1c01163>. I/O functions for xtc format files were implemented using the 'xdrfile' library available under open source license. The relevant information can be found in inst/COPYRIGHT. |
| Authors: | Alessandro Pandini [aut, cph]
, Stefano Motta [aut,
cre, cph] , Erik
Lindahl [ctb] (Author of xdrfile C library), David van der
Spoel [ctb] (Author of xdrfile C library) |
| Maintainer: | Stefano Motta <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.2 |
| Built: | 2025-03-03 05:19:04 UTC |
| Source: | https://github.com/cran/SOMMD |
Function to compute the average value of a property for each neuron.
average.neur.property(SOM, P)average.neur.property(SOM, P)
SOM |
the SOM object to cluster |
P |
the property for each frame of the simulation |
The a vector with the per-neuron average of the property.
Stefano Motta [email protected]
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute distance between two atoms in every frame of the simulation Distance <- apply(trj$coord[c(162,1794),,], 3, dist) #Compute average property value for each neuron avg.p <- average.neur.property(som_model, Distance)#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute distance between two atoms in every frame of the simulation Distance <- apply(trj$coord[c(162,1794),,], 3, dist) #Compute average property value for each neuron avg.p <- average.neur.property(som_model, Distance)
Function to compute distances to be used to train the SOM
calc.distances(trj, mol.2 = FALSE, sele = FALSE, atoms = NULL, cap = NULL)calc.distances(trj, mol.2 = FALSE, sele = FALSE, atoms = NULL, cap = NULL)
trj |
contains the trajectory coordinates (array with three dimensions obtained by rioxdr) |
mol.2 |
contains the atom indexes of the second molecule in case only intermolecular distances should be computed |
sele |
contains the selection of distances coming from the native_contacts function |
atoms |
contains a list of atoms indexes on which the distances will be computed |
cap |
If a number is given, distances greater than this value are set at the cap value |
A matrix containing the set of distances computed for all the frames.
Stefano Motta [email protected]
# Read reference structure file with native conformation struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Select only Cbeta atoms to perform the analysis sele_atoms <- which(trj$top$elety=="CB") # Choose only native contacts sele_dists <- native.cont(struct=struct, distance=1.0, atoms=sele_atoms) # Compute distances for SOM training. DIST <- calc.distances(trj, mol.2=FALSE, sele=sele_dists, atoms=sele_atoms)# Read reference structure file with native conformation struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Select only Cbeta atoms to perform the analysis sele_atoms <- which(trj$top$elety=="CB") # Choose only native contacts sele_dists <- native.cont(struct=struct, distance=1.0, atoms=sele_atoms) # Compute distances for SOM training. DIST <- calc.distances(trj, mol.2=FALSE, sele=sele_dists, atoms=sele_atoms)
Compute the pairwise distance matrix of a given set of coordinates, and only retain to some selected distances
calc.dists(coord, mol.1_id = FALSE, mol.2_id = FALSE, sele = FALSE)calc.dists(coord, mol.1_id = FALSE, mol.2_id = FALSE, sele = FALSE)
coord |
matrix of N atomic coordinates (N rows, 3 columns) |
mol.1_id |
vector containing the index of the first molecule for intermolecular distances only |
mol.2_id |
vector containing the index of the second molecule for intermolecular distances only |
sele |
contains the selection of distances coming from the native_contacts function |
A matrix contaning the selected distances for a frame
Stefano Motta[email protected]
Function to concatenate two simulations.
cat.trj(trj1, ...)cat.trj(trj1, ...)
trj1 |
the first trj file |
... |
additional trj files |
A trj object with the simulations concatenated
Stefano Motta [email protected]
# Read the simulations trj1 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) trj2 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD-2.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Concatenate the simulations trj <- cat.trj(trj1, trj2)# Read the simulations trj1 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) trj2 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD-2.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Concatenate the simulations trj <- cat.trj(trj1, trj2)
Cluster pathways according to a time dependent or independent scheme
cluster.pathways( SOM, start, end, time.dep = "independent", method = "complete" )cluster.pathways( SOM, start, end, time.dep = "independent", method = "complete" )
SOM |
a kohonen SOM object. |
start |
the vector specifying the starting frame of each replicas |
end |
the vector specifying the ending frame of each replicas |
time.dep |
choose whether to use time "dependent" or "independent" clustering of pathways |
method |
the method to be passed to hclust for the clustering |
representatives a vector of frames representatives of each neuron
Stefano Motta[email protected]
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Assign length of the replicas trj$start <- seq(1, 25, by=5) trj$end <- seq(5, 25, by=5) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Cluster Pathways using the time dependent algorithm clus.paths.tdep <- cluster.pathways(som_model, start=trj$start, end=trj$end, time.dep="dependent") #Cluster Pathways using the time independent algorithm clus.paths.tindep <- cluster.pathways(som_model, start=trj$start, end=trj$end, time.dep="independent")#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Assign length of the replicas trj$start <- seq(1, 25, by=5) trj$end <- seq(5, 25, by=5) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Cluster Pathways using the time dependent algorithm clus.paths.tdep <- cluster.pathways(som_model, start=trj$start, end=trj$end, time.dep="dependent") #Cluster Pathways using the time independent algorithm clus.paths.tindep <- cluster.pathways(som_model, start=trj$start, end=trj$end, time.dep="independent")
Compute the cluster representatives
cluster.representatives(SOM, clusters)cluster.representatives(SOM, clusters)
SOM |
a kohonen SOM object. |
clusters |
a vector of clusters assignment for each neuron, as returned for example by hclust. |
A vector of frames representatives of each neuron
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) # Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Get representative frames for each cluster cl_repres <- cluster.representatives(som_model, som_cl)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) # Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Get representative frames for each cluster cl_repres <- cluster.representatives(som_model, som_cl)
Compute the transition matrix starting from a vector of subsequent classifications
comp.trans.mat(SOM, start = 1)comp.trans.mat(SOM, start = 1)
SOM |
a kohonen object on which transitions between neurons will be computed |
start |
a vector containing the start frames of each replica (usually contained in trj$start if replicas were merged with cat_trj) |
A matrix of pairwise transitions between neurons
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute transition Matrix tr_mat <- comp.trans.mat(som_model, start = 1)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute transition Matrix tr_mat <- comp.trans.mat(som_model, start = 1)
Coordinate superposition with the Kabsch algorithm. This function make use of the bio3d fit.xyz function to align a SOMMD trj object. If ref is not specified, the trj object is aligned to the first frame of the simulation, otherwise it is aligned to the reference input object.
fit.trj(trj, ref = NULL, trj.inds = NULL, ref.inds = NULL)fit.trj(trj, ref = NULL, trj.inds = NULL, ref.inds = NULL)
trj |
an object with class trj |
ref |
a struct object read with read.struct() to be used as reference |
trj.inds |
a vector of indices that selects the trj atoms upon which fitting should be based. If not specified all atoms will be used. |
ref.inds |
a vector of indices that selects the ref atoms upon which fitting should be based. If not specified all atoms will be used. |
A trj object aligned
Stefano Motta [email protected]
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Fit a trajectory to the first frame based on alpha carbons: ca.inds <- which(trj$top$elety=="CA") trj.fit <- fit.trj(trj, trj.inds=ca.inds)#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Fit a trajectory to the first frame based on alpha carbons: ca.inds <- which(trj$top$elety=="CA") trj.fit <- fit.trj(trj, trj.inds=ca.inds)
Function map a numeric vector of a property to a vector of colors for that property according to that property value.
map.color(x, pal, limits = NULL, na.col = "grey")map.color(x, pal, limits = NULL, na.col = "grey")
x |
a numeric vector |
pal |
a color palette |
limits |
the values of the extremes for the colorscale |
na.col |
the color that will be assigned to the na.values of the vector |
A vector with colors proportional to the values of x
Function to convert a transition matrix to an igraph object
matrix2graph(trans, SOM, SOM.hc, col.set, diag = FALSE)matrix2graph(trans, SOM, SOM.hc, col.set, diag = FALSE)
trans |
a transition matrix (usually obtained from comp.trans.mat) |
SOM |
a kohonen object that form the network |
SOM.hc |
a vector of cluster assignment for SOM neurons |
col.set |
a vector of colors used for the SOM clusters |
diag |
boolean condition to include diagonal elements |
An igraph object, with SOM properties annotated
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Compute transition matrix tr_mat <- comp.trans.mat(som_model, start = 1) #Define a set of colors colors <- c("#1f78b4", "#33a02c", "#e31a1c", "#ffff88", "#6a3d9a") #Create graph object net <- matrix2graph(tr_mat, som_model, som_cl, colors, diag=FALSE)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Compute transition matrix tr_mat <- comp.trans.mat(som_model, start = 1) #Define a set of colors colors <- c("#1f78b4", "#33a02c", "#e31a1c", "#ffff88", "#6a3d9a") #Create graph object net <- matrix2graph(tr_mat, som_model, som_cl, colors, diag=FALSE)
Function to select only distances between residues making contacts in reference file or a frame of the simulation
native.cont( struct = NULL, trj = NULL, trj.frame = 1, distance, mol.2 = FALSE, atoms = NULL )native.cont( struct = NULL, trj = NULL, trj.frame = 1, distance, mol.2 = FALSE, atoms = NULL )
struct |
a struct object read with read.struct() to compute the native.cont |
trj |
a trj object to compute the native.cont |
trj.frame |
The frame of the trj on which the native.cont are computed |
distance |
the distance cut-off |
mol.2 |
can be FALSE (default), use the whole distance matrix, or a vector containing the atomic number of the second molecule (and compute only intermolecular distances) |
atoms |
can be NULL (default), consider all the atoms present in coords, or a vector containing a set of atomic numbers to consider in the calculation (e.g. only CB). atoms can be obtained with the bio3d atom.select function |
A vector containing the index of a subset of selected distances
Stefano Motta [email protected]
# Read reference structure file with native conformation struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Select only Cbeta atoms to perform the analysis sele_atoms <- which(struct$atom$elety=="CB") #Choose only native contacts sele_dists <- native.cont(struct=struct, distance=1.0, atoms=sele_atoms)# Read reference structure file with native conformation struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Select only Cbeta atoms to perform the analysis sele_atoms <- which(struct$atom$elety=="CB") #Choose only native contacts sele_dists <- native.cont(struct=struct, distance=1.0, atoms=sele_atoms)
Function to compute the per-neuron population
neur.population(SOM, start = 1, end = length(SOM$unit.classif), N = 1)neur.population(SOM, start = 1, end = length(SOM$unit.classif), N = 1)
SOM |
the SOM object |
start |
a vector containing the start frames of each replica (usually contained in trj$start if replicas were merged with cat_trj) |
end |
a vector containing the end frames of each replica (usually contained in trj$end if replicas were merged with cat_trj) |
N |
An integer for the portion (replica) of the simulations to be plotted |
A vector containing the per-neuron population
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) pop <- neur.population(som_model)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) pop <- neur.population(som_model)
Compute the representative frame of each neuron (the closest to the neuron codebook)
neur.representatives(SOM)neur.representatives(SOM)
SOM |
a kohonen SOM object. |
A vector containing the index of the representative frames for each neuron
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute representative frame for each neuron neuron_representatives <- neur.representatives(som_model)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute representative frame for each neuron neuron_representatives <- neur.representatives(som_model)
A short description...
## S3 method for class 'struct' print(x, ...)## S3 method for class 'struct' print(x, ...)
x |
trj object |
... |
additional arguments to be passed to further methods |
Called for its effect.
Stefano Motta [email protected]
# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Print basic information print(struct)# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Print basic information print(struct)
A short description...
## S3 method for class 'trj' print(x, ...)## S3 method for class 'trj' print(x, ...)
x |
trj object |
... |
additional arguments to be passed to further methods |
Called for its effect.
Stefano Motta [email protected]
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Print basic informations print(trj)#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Print basic informations print(trj)
Function to read gro files
read.gro(file)read.gro(file)
file |
contains the name and the path to the gro file to be read |
Returns a list of class "gro" with the following components:
atom |
a data frame containing all atomic coordinate with a row per atom and a column per record type. |
xyz |
a numeric matrix of class "xyz" containing the atomic coordinate data. |
box |
a vector of box size. |
call |
the matched call. |
Stefano Motta [email protected]
Function to read pdb and gro files
read.struct(file)read.struct(file)
file |
contains the name and the path to the pdb or gro file to be read |
Returns a list of class "struct" with the following components:
atom |
a data frame containing all atomic coordinate with a row per atom and a column per record type. |
xyz |
a numeric matrix of class "xyz" containing the atomic coordinate data. |
box |
a vector of box size. |
format |
The format of the original file |
call |
the matched call. |
Stefano Motta [email protected]
# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD"))# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD"))
Function to read a trajectory file
read.trj(trjfile, topfile)read.trj(trjfile, topfile)
trjfile |
contains the name and the path to the reference file (pdb or gro files are accepted) |
topfile |
contains the name and the path to the trajectory file (xtc or dcd files are accepted) |
Returns a list of class "trj" with the following components:
topfile |
the input topology file. |
topformat |
the format of the input topology. |
trjfile |
the input trajectory file. |
trjformat |
the format of the input trajectory. |
coord |
a three dimensional array containing atomic coordinates for all the frames. Dimensions are: Natoms:3:Nframes. |
top |
a data.frame containing topological informations with a row per atom and a column per record type (resno, resid, elety, eleno, chain). |
start |
a vector with the first frame of the simulation. When multiple simulations are concatenated with |
end |
a vector with the last frame of the simulation. When multiple simulations are concatenated with |
call |
the matched call. |
Alessandro Pandini
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD"))#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD"))
Assign new data to a pre-trained SOM
remap.data(SOM, X, add = FALSE)remap.data(SOM, X, add = FALSE)
SOM |
a trained SOM |
X |
a data set with the same number of features of the dataset used to train the SOM |
add |
whether to append the new data to the ones used to train the SOM |
An object of class "kohonen" with the new data mapped
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Read a trajectory that was not used to train the som trj_2 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD-2.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read reference structure file gro <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Selection of the same intermolecular distances used to train the SOM protein.sele <- which(gro$atom$resid!="020") ligand.sele <- which(gro$atom$resid=="020") heavy.atoms <- which(startsWith(gro$atom$elety, "H")==FALSE) sele.dists <- native.cont(struct=gro, distance=0.6, mol.2=ligand.sele, atoms=heavy.atoms) # Compute distances on new simulations (the same used for SOM training) dist_2 <- calc.distances(trj_2, mol.2=ligand.sele, sele=sele.dists, atoms=heavy.atoms) # Map new data on the existing SOM som_model_2 <- remap.data(SOM=som_model, X=dist_2)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Read a trajectory that was not used to train the som trj_2 <- read.trj(trjfile = system.file("extdata", "HIF2a-MD-2.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read reference structure file gro <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Selection of the same intermolecular distances used to train the SOM protein.sele <- which(gro$atom$resid!="020") ligand.sele <- which(gro$atom$resid=="020") heavy.atoms <- which(startsWith(gro$atom$elety, "H")==FALSE) sele.dists <- native.cont(struct=gro, distance=0.6, mol.2=ligand.sele, atoms=heavy.atoms) # Compute distances on new simulations (the same used for SOM training) dist_2 <- calc.distances(trj_2, mol.2=ligand.sele, sele=sele.dists, atoms=heavy.atoms) # Map new data on the existing SOM som_model_2 <- remap.data(SOM=som_model, X=dist_2)
Function to read an xtc trajectory file
rio_read_xtc(xtc_filename)rio_read_xtc(xtc_filename)
xtc_filename |
contains the name and the path to the xtc file |
Returns 3D array of cartesian coordinates
Alessandro Pandini
Function to read a xtc trajectory file
rio_read_xtc_natoms(xtc_filename)rio_read_xtc_natoms(xtc_filename)
xtc_filename |
contains the name and the path to the xtc file |
Returns number of atoms in the structure
Alessandro Pandini
Function to read an xtc trajectory file
rio_read_xtc_nframes(xtc_filename)rio_read_xtc_nframes(xtc_filename)
xtc_filename |
contains the name and the path to the xtc file |
Returns number of frames in the trajectory
Alessandro Pandini
Function to read an xtc trajectory file
rio_read_xtc2xyz(xtc_filename)rio_read_xtc2xyz(xtc_filename)
xtc_filename |
contains the name and the path to the xtc file |
Returns bio3d xyz array of Cartesian coordinates
Alessandro Pandini
Function to write an xtc trajectory file
rio_write_xtc(xtc_filename, trj)rio_write_xtc(xtc_filename, trj)
xtc_filename |
contains the name and the path to the xtc file to write |
trj |
trajectory object to save |
Returns status of write execution
Alessandro Pandini
Function to compute the silhouette profile for the Nclus cluster of the SOM neurons
silhouette.profile( SOM, Nclus, dist_clust = "euclidean", clust_method = "complete" )silhouette.profile( SOM, Nclus, dist_clust = "euclidean", clust_method = "complete" )
SOM |
the SOM object to cluster |
Nclus |
the cluster number on which the silhouette profile will be computed |
dist_clust |
the metric for the distance calculation |
clust_method |
the method for the clustering (passed to the hclust function |
A vector of silhouette profile computed with the cluster package
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute the silhouette profile sil_pro <- silhouette.profile(som_model, Nclus=5, clust_method="complete")#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute the silhouette profile sil_pro <- silhouette.profile(som_model, Nclus=5, clust_method="complete")
Function to compute the silhouette score for the clustering of SOM neurons
silhouette.score( SOM, dist_clust = "euclidean", clust_method = "complete", interval = seq(2, 30) )silhouette.score( SOM, dist_clust = "euclidean", clust_method = "complete", interval = seq(2, 30) )
SOM |
the SOM object to cluster |
dist_clust |
the metric for the distance calculation |
clust_method |
the method for the clustering (passed to the hclust function |
interval |
the cluster number on which the silhouette score will be computed |
A vector with the silhouette scores for all the frames
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute the silhouette profile sil_score <- silhouette.score(som_model, clust_method="complete", interval=seq(2,8))#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Compute the silhouette profile sil_score <- silhouette.score(som_model, clust_method="complete", interval=seq(2,8))
Function to add circles to a SOM plot, with dimension proportional to a selected property
som.add.circles(SOM, P, scale = 1, col.circles = "white")som.add.circles(SOM, P, scale = 1, col.circles = "white")
SOM |
the SOM object |
P |
a vector containing the per-neuron property to plot |
scale |
a number to scale up or down the size of the drawn circles |
col.circles |
the background color of the drawn circles |
Called for its effect.
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) # Compute per neuron population pop <- neur.population(som_model) #Plot the som plot(som_model, type = "count", bgcol=c("red", "blue", "yellow", "green"), shape='straight') # Add circles to the SOM plot som.add.circles(som_model, pop, scale=0.9)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) # Compute per neuron population pop <- neur.population(som_model) #Plot the som plot(som_model, type = "count", bgcol=c("red", "blue", "yellow", "green"), shape='straight') # Add circles to the SOM plot som.add.circles(som_model, pop, scale=0.9)
Function to apply a legend of clusters to a SOM map image
som.add.clusters.legend(Nclus, color.scale)som.add.clusters.legend(Nclus, color.scale)
Nclus |
the number of clusters to which put the legent |
color.scale |
the color scale used for the image |
Called for its effect.
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Define a set of colors colors <- c("#1f78b4", "#33a02c", "#e31a1c", "#ffff88", "#6a3d9a") #Plot the som with neurons colored according to clusters plot(som_model, type = "mapping", bgcol=colors[som_cl], col=rgb(0,0,0,0), shape='straight', main="") kohonen::add.cluster.boundaries(som_model, som_cl, lwd=5) #Add legend to the plot som.add.clusters.legend(Nclus=4, color.scale=colors)#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Divide the SOM in the selected number of clusters som_cl <- cutree(hclust(dist(som_model$codes[[1]], method="euclidean"), method="complete"), 4) #Define a set of colors colors <- c("#1f78b4", "#33a02c", "#e31a1c", "#ffff88", "#6a3d9a") #Plot the som with neurons colored according to clusters plot(som_model, type = "mapping", bgcol=colors[som_cl], col=rgb(0,0,0,0), shape='straight', main="") kohonen::add.cluster.boundaries(som_model, som_cl, lwd=5) #Add legend to the plot som.add.clusters.legend(Nclus=4, color.scale=colors)
Add the neuron numbering scheme to the SOM plot
som.add.numbers(SOM, scale = 1, col = "black")som.add.numbers(SOM, scale = 1, col = "black")
SOM |
the SOM object |
scale |
a number to scale up or down the size of the text |
col |
the color of the text |
Called for its effect.
Stefano Motta [email protected]
#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Plot the som plot(som_model, type = "count", bgcol=c("red", "blue", "yellow", "green"), shape='straight') #Add neuron numbers on the som som.add.numbers(som_model, scale=0.5, col="black")#Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #Plot the som plot(som_model, type = "count", bgcol=c("red", "blue", "yellow", "green"), shape='straight') #Add neuron numbers on the som som.add.numbers(som_model, scale=0.5, col="black")
Apply a stride to the frame of a trj object to reduce the number of frames
stride.trj(trj, stride)stride.trj(trj, stride)
trj |
a trj object. |
stride |
the stride to apply to the trajectory |
An object of class trj with a frame every stride
Stefano Motta [email protected]
# Read the simulation trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # keep a frame every 2 frame trj_strd <- stride.trj(trj, 2)# Read the simulation trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # keep a frame every 2 frame trj_strd <- stride.trj(trj, 2)
Convert a struct object into a pdb obtect
struct2pdb(struct)struct2pdb(struct)
struct |
contains the struct object to convert |
Returns an object with class "pdb"
An object of class "pdb"
Stefano Motta [email protected]
# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Convert structure to pdb object pdb <- struct2pdb(struct)# Read structure file struct <- read.struct(system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Convert structure to pdb object pdb <- struct2pdb(struct)
summary method for class struct
## S3 method for class 'struct' summary(object, ...)## S3 method for class 'struct' summary(object, ...)
object |
struct object |
... |
additional arguments to be passed to further methods |
Called for its effect.
Stefano Motta [email protected]
summary method for class trj
## S3 method for class 'trj' summary(object, ...)## S3 method for class 'trj' summary(object, ...)
object |
trajectory object |
... |
additional arguments to be passed to further methods |
Called for its effect.
Stefano Motta [email protected]
Function trace pathway sampled on the SOM
trace.path( SOM, start = 1, end = length(SOM$unit.classif), N = 1, draw.stride = 1, pts.scale = 1, lwd.scale = 1 )trace.path( SOM, start = 1, end = length(SOM$unit.classif), N = 1, draw.stride = 1, pts.scale = 1, lwd.scale = 1 )
SOM |
the SOM object |
start |
a vector containing the start frames of each replica (usually contained in trj$start if replicas were merged with cat_trj) |
end |
a vector containing the end frames of each replica (usually contained in trj$end if replicas were merged with cat_trj) |
N |
The portion of simulation that one want to plot |
draw.stride |
used to plot the pathways with a stride (useful for very complex pathways) |
pts.scale |
a number to scale up or down the size of the circles |
lwd.scale |
a number to scale up or down the size of the lines |
Called for its effect.
Stefano Motta [email protected]
# Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #trace pathway sampled on the SOM trace.path(som_model, start=trj$start, end=trj$end, N=1, pts.scale=0.5)# Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) #Read example SOM data som_model <- readRDS(system.file("extdata", "SOM_HIFa.rds", package = "SOMMD")) #trace pathway sampled on the SOM trace.path(som_model, start=trj$start, end=trj$end, N=1, pts.scale=0.5)
Extract a trj frame to a pdb object
trj2pdb(trj, frame, filename)trj2pdb(trj, frame, filename)
trj |
a trj object. |
frame |
the frame to extract. |
filename |
for the output pdb file |
a pdb object of the selected frame
Called for its effect.
Stefano Motta [email protected]
# Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Write the pdb file for a specific frame trj2pdb(trj = trj, frame=5, filename = tempfile(fileext = '.pdb' ))# Read the trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) # Write the pdb file for a specific frame trj2pdb(trj = trj, frame=5, filename = tempfile(fileext = '.pdb' ))
Convert the trj coordinates 3D-array in a 2D matrix.
trj2xyz(trj, inds = NULL)trj2xyz(trj, inds = NULL)
trj |
an object with class trj |
inds |
indices for the output coordinates |
a xyz matrix with frames on rows and coordinates as columns
Stefano Motta [email protected]
#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) trj2xyz(trj)#Read trajectory trj <- read.trj(trjfile = system.file("extdata", "HIF2a-MD.xtc", package = "SOMMD"), topfile = system.file("extdata", "HIF2a.gro", package = "SOMMD")) trj2xyz(trj)