Differences

This shows you the differences between two versions of the page.

--- pairwise_2d-rdf_and_clustering [2020/07/26 12:32] – [Pairwise 2D-RDF and Clustering Analysis] edit
+++ pairwise_2d-rdf_and_clustering [2022/10/23 21:33] (current) – [Pairwise 2D-RDF and Clustering Analysis] edit
@@ Line 4: / Line 4: @@
 Author: Yalun Yu (alanyu17@terpmail.umd.edu; yalun.research@gmail.com)
+Bug fix log:
+/24/2020: An asymmetric membrane only bug found in clustering.py (line # 757); status: code fixed and updated, example not updated.
+{{ :cluster.tar.gz |}}
+{{ :cluster-example.tar.gz |}}
 This new code allows pair-specific "Dcut" in clustering analysis. To find the proper pair/lipid-specific "Dcut", please use the 2D-RDF code included in the same mini-package.
@@ Line 10: / Line 17: @@
 **Background**
-Different lipids have different sizes, so that we want to allow pair/lipid-specific Dcuts. You can compare this to LJ radius, where Rmin = Rmin(a)/2 + Rmin(b)/2 (we call this combination rule. Apparently, this is **atom-specific**). Or, if you are using a NBFIX (**pair-specific** interaction) for the LJ interactions between atom type a & b, Rmin=Rmin (NBFIX between a & b). The concept of NBFIX can be found at:
+Different lipids have different sizes, so that we want to allow pair/lipid-specific Dcuts. You can compare this to LJ radius, where Rmin = Rmin(a)/2 + Rmin(b)/2 (we call this combination rule. Apparently, this is **atom-specific**). Or, if you are using a NBFIX (**pair-specific** interaction) for the LJ interaction between atom type a & b, Rmin=Rmin (NBFIX between a & b). The concept of NBFIX can be found at:
 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2853947/
@@ Line 20: / Line 27: @@
 The DBSCAN algorithm used in the old clustering code uses the XY coordinates as inputs. However, DBSCAN also allows distance matrix as inputs, which means we can scale the elements in the real distance matrix as we want. For example, we have a matrix:
-  [Particle] A           B          C     ...
+  [Particle] A           B         C     ...
   A          0           1.5       1.5
   B          1.5         0         2
@@ Line 32: / Line 39: @@
 numpy, mdtraj, scikit-learn, openmm, matplotlib
+I recommend using **bash** instead of csh to avoid glitches when using conda! On DT2, you can change to bash by just typing "bash", you can also change your default shell following the instructions on Deepthought2 website.
 A quick one step command could be (after installing conda):
@@ Line 37: / Line 46: @@
 conda create -n cluster -c omnia python=3.6 openmm mdtraj numpy scikit-learn matplotlib
-Using command, we created a conda environment (env) called "cluster" and installed several packages into that environment. To use this env, type "source activate cluster" in the command line, and you should see your command line changing from **[login-1 ~]$ ...** to **(cluster) [login-1 ~]$ ...**
+Using this command, we created a conda environment (env) called "cluster" and installed several packages into that environment. To use this env, type "conda activate cluster" in the command line, and you should see your command line changing from **[login-1 ~]$ ...** to **(cluster) [login-1 ~]$ ...**
 (2) The RDF code
-Two versions of the RDF code are in the "rdf" and "rdf-explicit" subfolders. The two are basically the same, and the only difference between them is the method to separate upper and lower leaflets. In "rdf-explicit", a coordinate file extracted from a re-centered trajectory is needed as input, while the other one ("rdf") needs the user to find the last atom index of the upper leaflet from the psf file. My experience is that the "rdf" template is good for membrane containing only phospholipids built from CHARMM-GUI. If you have any lipid built by patching two residues (for example, the glycolipids), then you should use the "rdf-explicit" folder as your template. To extract the coordinates from a re-centered trajectory, first get the re-centered the trajectory following the standard steps in the EDP section, and then use the extract.inp in "ref-explicit" folder to get the .crd file. The instructions are all in-place within the extract.inp file.
+Two versions of the RDF code are in the "rdf" and "rdf-explicit" subfolders. The two are basically the same, and the only difference between them is the method to separate upper and lower leaflets. In "rdf-explicit", a coordinate file where the membrane is centered at 0 in the z direction is needed as input (usually step5_assembly.psf from CHARMM-GUI meets this requirement). In the worst case, a coordinate file extracted from a re-centered trajectory can be used. In the other one ("rdf"), the user have to find the last atom index of the upper leaflet from the psf file. My experience is that the "rdf" template is good for membrane containing only phospholipids built from CHARMM-GUI. If you have any lipid built by patching two residues (for example, the glycolipids), then you should use the "rdf-explicit" template. To extract the coordinates from a re-centered trajectory (if you really need to), first get the re-centered trajectory following the standard steps in the EDP section, and then use the extract.inp in "ref-explicit" folder to get the .crd file. The instructions are all in-place within the extract.inp file.
 To run the rdf analysis between any two selections, update the corresponding parts in the run_rdf.py file and submit the job by "sbatch run_rdf.sh". Please refer to http://mdtraj.org/1.9.3/atom_selection.html for correct selection word format.
@@ Line 48: / Line 59: @@
 Note this RDF code also allows self-RDF (just make the two atom selection words the same).
 (3) The clustering code
@@ Line 91: / Line 104: @@
 * use lower case for head group name but upper case for atom name.
-b. In **run_cluster.py**, the pair specific Dcut can be implemented by the clfix keyword argument when creating a ClusterLip instance (see the code and you will know what I mean). Other changes should follow the instructions in the code.
+b. In **run_cluster.py**, the pair specific Dcut can be implemented by the **clfix** argument when creating a ClusterLip instance (see the code and you will know what I mean). The minimum number of lipids in a cluster can be modified using the **min_lipids** argument.
-Submit the clustering analysis by "sbatch run_cluster.sh"
+c. Include ALL lipids in your system.
+d. Other changes should follow the instructions in the code.
+(4) The util.py
+a. Make corresponding edits in util.py (could copy from run_cluster.py) before submitting job
+b. Submit the clustering analysis by "sbatch run_cluster.sh"
 **Check the results:**
@@ Line 105: / Line 128: @@
 * in topsize.png / botsize.png, you can find the distribution of the cluster size.
-c. If you want to have visualization of the code, please use the visualize.py file. Usually it takes about 10 seconds to generate cluster pictures for one frame, so that I won't use slurm. However, if you want to create pictures for more than 10 frames, please consider using slurm. Again, the instructions is within the code.
+(5) Visualization
+If you want to visualize the clusters, please use the visualize.py file.
+a. activate conda environment
+b. Make corresponding edits in visualize.py
+c. define colors for different headgroups at the beginning of clustering.py.
+<code>
+def color_(head):
+    # used for cluster visualization
+    if head == 'pe':
+        return 'yellow'
+    elif head == 'pc':
+        return 'yellow'
+    elif head == 'pg':
+        return 'yellow'
+    elif head == 'pi':
+        return 'yellow'
+    elif head == 'ps':
+        return 'yellow'
+    elif head == 'sito':
+        return 'b'
+    elif head == 'cer':
+        return 'r'
+</code>
+* Could also use //"else:"// command to minimize repetitive line entry for colors of lipids not of our interest
+d. You can define what's included in the cluster plot by editing the //top_res_info// and //bot_res_info// in visualize.py. (Note: the code structure is identical to that in //run_cluster.py//, but you do not have to include all lipids in the system and just the lipids of interest.)
+e. Usually it takes about 10 seconds to generate two cluster pictures (upper&lower leaflets) for one frame, so that I won't use slurm. So, just type python visualize.py and wait for a while (remember to activate the conda env you created). However, if you want to create pictures for more than 10 frames, please consider using slurm. Again, the instructions are within the code.