DESCRIPTION
 
 This script plots the density profile along the z axis of groups of particles present
 around transmembrane proteins. Density profiles are broken down by the size of the
 TM clusters.

 A file containing the charged particles can also be supplied to calculate the density
 of charges.
 
 density calculation
 -------------------
 Density is calculated based on a cylinder centered on the protein cluster center of 
 geometry. All the particles present in this cylinder are binned into cylinder slices
 based on their distance (along z) to the center of the lipid bilayer. The cylinder
 radius is controlled via --slices_radius and the slices thickness via --slices_thick.

 (a) particles selection
   You can specify the particles for which to plot the density by supplying a file
   via the --particles option. Each line of this file should follow the following
   format (without quotation marks):
    -> 'group,label,colour,MDAnalysis selection string'
 
   where 'group' is used to normalise the densities of several particles. By default
   each particle type is normalised with respect to itself and the following densities
   are used:
    -> protein,protein,#262626,protein
    -> CHOL,CHOL,#bdbdbd,resname CHOL and name ROH
    -> POPC,POPC,#41ab5d,resname POPC and name PO4
    -> POPE,POPE,#6a51a3,resname POPE and name PO4
    -> POPS,POPS,#cc4c02,resname POPS and name PO4
    -> water,water,#1d91c0,resname W
    -> Na+,Na+,#7bccc4,name NA+
    -> Cl-,Cl-,#fa9fb5,name CL-
  
   Note that the only acceptable protein selection is the whole "protein" string and
   that it should be labeled as "protein" and in a group of its own also labeled
   "protein". If you're interested in particular sub-selections of proteins, see (b)
   below.
   
 (b) protein details: residue types
   Density profiles can be broken down further for peptides to show the density
   distribution amongst different residue types. You can group the residues as you
   wish by supplying a file via the --residues options. Each line of this file
   should follow the following format:
   -> 'label,colour,resname1,resname2,...'
  
   By default the following residue types are used:
    -> basic,#253494,ARG,LYS
    -> polar,#a1d99b,SER,THR,ASN,GLN,HIS		
    -> negative,#99d8c9,ASP, GLU
    -> hydrophobic,#993404,VAL,ILE,LEU,MET,PRO,CYS,PHE,TYR,TRP
    -> backbone,#969696,ALA,GLY

   WARNING: you need to have defined a particle type called "protein" in (a) in order
   to show residue details.
   If you do not want to show residue details, just use: '--residues no'.

 (c) charge density
   You can specify which particles to take into account for the calculation of the total
   charge density by supplying a file via the --charges option. Each line of this file
   should follow the format (without quotation marks):
    -> 'group_name,colour,charge_name,charge,MDAnalysis selection string for charge'

   Note that the MDAnalysis selection string should not contain any commas.

   The total charge for each group will be plotted on the charge density profile. The
   group colour must be specified for each charge.
  
   By default charges are defined using the Martini 2.1 settings:
    -> solvent,#52A3CC,Na+,1,name NA+
    -> solvent,#52A3CC,CL-,-1,name CL-
    -> lipids,#b2182b,-1,phosphate,name PO4
    -> lipids,#b2182b,1,amine_choline,name NH3 or name NC3
    -> protein,#053061,Lys,1,resname LYS and name SC2
    -> protein,#053061,Arg,1,resname ARG and name SC2
    -> protein,#053061,Asp,-1,resname ASP and name SC1
    -> protein,#053061,Glu,-1,resname GLU and name SC1
  
   Another default set of charges can be used by specifiying --charges 2.2P :
    -> solvent,#52A3CC,Na+,1,name NA+
    -> solvent,#52A3CC,CL-,-1,name CL-
    -> solvent,#52A3CC,WP,0.46,name WP
    -> solvent,#52A3CC,WM,-0.46,name WM
    -> lipids,#b2182b,-1,phosphate,name PO4
    -> lipids,#b2182b,1,amine_choline,name NH3 or name NC3    
    -> protein,#053061,Lys,1,resname LYS and name SCP
    -> protein,#053061,Arg,1,resname ARG and name SCP
    -> protein,#053061,Asp,-1,resname ASP and name SCN
    -> protein,#053061,Glu,-1,resname GLU and name SCN
    -> protein,#053061,Asn_p,0.46,resname ASN and name SCP
    -> protein,#053061,Asn_n,-0.46,resname ASN and name SCN
    -> protein,#053061,Gln_p,0.46,resname GLN and name SCP
    -> protein,#053061,Gln_n,-0.46,resname GLN and name SCN
    -> protein,#053061,Thr_p,0.31,resname THR and name SCP
    -> protein,#053061,Thr_n,-0.31,resname THR and name SCN
    -> protein,#053061,Ser_p,0.4,resname SER and name SCP
    -> protein,#053061,Ser_n,-0.4,resname SER and name SCN

   By default proteins termini are considered uncapped and +1 and -1 are added to the
   first and last backbone beads ("name BB") respectively. If this is not what you want
   just use the option --capped.
   
   Note that for systems containing different protein species automatic addition of
   charges at the termini has not yet been implemented: you should use --capped and
   provide a file via --charges specifying all charged particles in the system, including
   termini beads.

   If you do not want to calculate charge density, just use: '--charges no'

 (d) colour definition
   Colours can be specified using single letter code (rgbcmykw), hex code  or the name of
   a colour map (see the matplotlib website for a list of the available colour maps).
   In case a colour map is used, its name must be specified as the only colour.

 detection of transmembrane protein clusters
 -------------------------------------------
 Two clustering algorithms can be used to identify protein clusters.

 (a) Connectivity based (relies on networkX module):
   A protein is considered in a cluster if it is within a distance less than --nx_cutoff
   from another protein. This means that a single protein can act as a connector between
   two otherwise disconnected protein clusters.
   This algorithm can be ran using either the minimum distante between proteins (default, 
   --algorithm 'min') or the distance between their center of geometry (--algorithm 'cog').
   The 'min' option scales as the square of the number of proteins and can thus be very
   slow for large systems.

 (b) Density based (relies on the sklearn module and its implementation of DBSCAN):
   A protein is considered in a cluster if is surrounded by at least --db_neighbours other
   proteins within a radius of --db_radius.
   This density based approach is usually less suited to the detection of protein
   clusters but as a general rule the more compact the clusters, the smaller --db_radius
   the higher --db_neighbours can be - for details on this algorithm see its online
   documentation.
   This algorithm is selected by setting the --algorithm option to 'density'.

 The identified protein clusters are considered to be transmembrane only if the closest
 lipid headgroup neighbours to the cluster particles are all within the same leaflet.
 In addition to the sizes identified, size groups can be defined - see note 7.


REQUIREMENTS

 The following python modules are needed :
  - MDAnalysis
  - matplotlib
  - numpy
  - scipy
  - networkX (if option --algorithm is set to 'min' or 'cog')
  - sklearn (if option --algorithm is set to 'density')


NOTES

 1. The density is calculated with respect to the z axis, not the bilayer normal. So the
    more your system deforms the noiser the less meaningful the results get.

 2. Identification of the bilayer leaflets can be controlled via 3 options:
    (a) beads
     By default, the particles taken into account to define leaflet are:e
     -> name PO4 or name PO3 or name B1A
   
     Note that only lipids which contain one of the beads mentioned in the selection string
     will be taken into account. If you wish to specify your own selection string (e.g. to
     choose different beads or add a bead not in the default list in order to take into
     account a particular lipid specie) you can do so by supplying a file via the --beads
     option. This file should contain a single line that can be passed as the argument
     to MDAnalysis selectAtoms() routine and should not contain any quotation marks, e.g.:
      -> name PO4 or name PO3 or name B1A or name AM1
        
    (b) leaflet finding method
     By default leaflets are identified using the MDAnalysis LeafletFinder routine and the
     the optimum cutoff to identify 2 lipids groups is determined using the optimize_cutoff
     routine.
     This optimisation process can take time in large systems and you can specify your own
     cutoff value to skip this step. For instance to use a 15 Angstrom cutoff value:
      -> '--leaflet 15'
    
     In very large systems (more then ~50,000 phospholipids) LeafletFinder (or rather the
     networkX module that it relies on) can fail. To  avoid this you can choose not to use
     this routine by specifying:
      -> '--leaflet large'
     In this case lipids whose headgroups z value is above the average lipids z value will
     be considered to make up the upper leaflet and those whose headgroups z value is below
     the average will be considered to be in the lower leaflet.
     This means that the bilayer should be as flat as possible in the 1st frame of the xtc
     file supplied in order to get a meaningful outcome. 

	 NOTE: By default the gro file is only used as a topology file and the 1st frame of the
	 xtc is used to identify leaflets. If you wish to use the gro file instead, for instance
	 in the case that the 1st frame of the xtc is not flat, you need to specify the --use_gro
	 flag: be warned that this might take a few minutes longer on large systems.

    (c) flipflopping lipids
     In case lipids flipflop during the trajectory, a file listing them can be supplied
     with the --flipflops option. Each line of this file should follow the format:
      -> 'resname,resid,starting_leaflet,z_bead'
     where starting_leaflet is either 'upper' or 'lower' - e.g. 'POPC,145,lower,PO4'. The
     z_bead is used to track the position of the lipid.
     If flipflopping lipids are not specified they may add significant noise to results as
     they prevent the correct identification of TM clusters.
     Note that, even when specified, flipflopping lipids will be taken into account when
     calculating densities and charges.   

 3. Proteins are detected automatically but you can specify an input file to define your
    own selection with the --proteins option.
    In this case the supplied file should contain on each line a protein selection string
    that can be passed as the argument of the MDAnalysis selectAtoms() routine - for 
    instance 'bynum 1:344'.
   
 4. The densities are calculated for each TM cluster size identified but can also be
    binned into size groups.
    The size groups are defined by supplying a file with --groups, whose lines all
    follow the format:
     -> 'lower_size,upper_size'

    Size groups definition should follow the following rules:
     -to specify an open ended group use 'max', e.g. '3,max'
     -groups should be ordered by increasing size and their boundaries should not overlap
     -boundaries are inclusive so you can specify one size groups with 'size,size,colour'
     -any cluster size not covered will be labeled as 'other'
 
 5. There are 3 possible options to determine the local normal to the bilayer. These are
    controlled with the flags --normal and --normal_d:
    (a) 'z': the bilayer is assumed flat in the x,y plane and the z axis is taken to be the
     normal. Best for systems without significant curvature and local deformations. In this
     case the --normal_d flag is ignored.

    (b) 'cog': in this case neighbourhing particles to current cluster of interest are
     identified in the lower and upper leaflet. The local normal is then considered to be the
     vector going from the cog of the lower ones to the cog of the upper ones. In each leaflet,
     neighbouring particles are the particles selected by --beads which are within --normal_d
     Angstrom of the cog of the protein cluster of interest.

    (c) 'svd': in this case neighbourhing particles to current cluster of interest are
     identified in the lower and upper leaflet as in (b) above. The normal of the best fitting
     plane to these particles is obtained by performing a singular value decomposition of their
     coordinate matrix.

 
USAGE

Option	      Default  	Description                    
-----------------------------------------------------
-f			: structure file [.gro] (required)
-x			: trajectory file [.xtc]
-o			: name of output folder
-b			: beginning time (ns) (the bilayer must exist by then!)
-e			: ending time (ns)	
-t 		10	: process every t-frames
--particles		: definition of particles, see 'DESCRIPTION'
--types			: definition of residue groups, see 'DESCRIPTION'
--charges	2.1	: definition of charged particles, see 'DESCRIPTION' 
--capped		: assumes protein termini are capped
 
Density profile options
-----------------------------------------------------
--range		40 	: distance spanned on either side of the bilayer center
--slices_thick	0.5 	: z thickness of the slices (Angstrom)
--slices_radius	30 	: radius of the slices (Angstrom)
--normal	svd	: local normal to bilayer ('z', 'cog' or 'svd'), see note 5
--normal_d	50	: distance of points to take into account for local normal, see note 5

Lipids identification  
-----------------------------------------------------
--beads			: leaflet identification technique, see note 2(a)
--flipflops		: input file with flipflopping lipids, see note 2(c)
--leaflets	optimise: leaflet identification technique, see note 2(b)
--use_gro		: use gro file instead of xtc, see note 2(b)

Protein clusters identification
-----------------------------------------------------
--groups		: cluster groups definition file, see note 4
--proteins		: protein selection file, (optional, see note 3)
--algorithm	min	: 'cog','min' or 'density', see 'DESCRIPTION'
--nx_cutoff 	6	: networkX cutoff distance for protein-protein contact (Angstrom)
--db_radius 	20	: DBSCAN search radius (Angstrom)
--db_neighbours	3	: DBSCAN minimum number of neighbours within a circle of radius --db_radius	

Other options
-----------------------------------------------------
--version		: show version number and exit
-h, --help		: show this menu and exit