Skip to content

Commit

Permalink
added custom atom selection option (selection text) to lipid contacts…
Browse files Browse the repository at this point in the history 
… and updated the user manual
  • Loading branch information
@bernhardtna
bernhardtna committed Dec 1, 2022
1 parent 1004484 commit 1421d52
Show file tree
Hide file tree
Showing 6 changed files with 142 additions and 8 deletions.
Binary file modified mosaics_user_manual.docx
View file
Binary file not shown.
Binary file modified mosaics_user_manual.pdf
View file
Binary file not shown.
3 changes: 3 additions & 0 deletions src/MosAT/program_variables/pv_lipid_contacts.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ struct program_variables
string leaflet_finder_param_name; //Name of the leaflet finder param file
string protein_finder_param_name; //Name of the protein finder param file
string solvent_finder_param_name; //Name of the solvent finder param file
string selection_text_file_name; //Name of the selection card with the selection text
int b_sf_param; //Tells if the user included a solvent types parameter file
int b_pf_param; //Tells if the user included a protein types parameter file
int b_lf_param; //Tells if the user included a lipid types parameter file
Expand All @@ -43,6 +44,7 @@ struct program_variables
int b_stdev; //Compute the z-coord standard deviation?
int b_clean; //Delete single frame rmsd files after computing average?
int b_mol; //Count molecules instead of contacts?
int b_sel_text; //Refine the selection of protein atoms?
double APS; //This is the area of a grid square
double radius; //Radius of the atom
double cutoff; //Percentage of average rho used for excluding data
Expand Down Expand Up @@ -86,6 +88,7 @@ void initialize_program_variables(program_variables *p)
p->b_lf_pdb = 0;
p->b_lf_param = 0;
p->b_mol = 0;
p->b_sel_text = 0;

//here we set the program description
p->program_description = "Lipid Contacts is a program designed for quantifying contacts between lipids and other molecules. The program can measure lipid-lipid, lipid-protein, and lipid-solvent contacts or a mixture of the three. The number of contacts formed is projected onto the XY plane.";
Expand Down
3 changes: 3 additions & 0 deletions src/MosAT/program_variables/pv_lipid_contacts_3d.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ struct program_variables
string leaflet_finder_param_name; //Name of the leaflet finder param file
string protein_finder_param_name; //Name of the protein finder param file
string solvent_finder_param_name; //Name of the solvent finder param file
string selection_text_file_name; //Name of the selection card with the selection text
int b_sf_param; //Tells if the user included a solvent types parameter file
int b_pf_param; //Tells if the user included a protein types parameter file
int b_lf_param; //Tells if the user included a lipid types parameter file
Expand All @@ -42,6 +43,7 @@ struct program_variables
int b_stdev; //Compute the z-coord standard deviation?
int b_clean; //Delete single frame rmsd files after computing average?
int b_mol; //Count molecules instead of contacts?
int b_sel_text; //Refine the selection of protein atoms?
double APS; //This is the area of a grid square
double radius; //Radius of the atom
double cutoff; //Percentage of average rho used for excluding data
Expand Down Expand Up @@ -88,6 +90,7 @@ void initialize_program_variables(program_variables *p)
p->b_lf_param = 0;
p->b_mol = 0;
p->ex_val = 0.0;
p->b_sel_text = 0;

//here we set the program description
p->program_description = "Lipid Contacts 3d is a program designed for quantifying contacts between lipids and other molecules. The program can measure lipid-lipid, lipid-protein, and lipid-solvent contacts or a mixture of the three. The number of contacts formed is projected onto a 3d lattice.";
Expand Down
72 changes: 68 additions & 4 deletions src/MosAT/programs/lipid_contacts.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,11 @@ using namespace std;
#include "headers/switch.h" //This defines a switch (on, off)
#include "headers/file_reader.h" //This has basic routines for reading text files
#include "headers/vector_mpi.h" //This has routines for collecting vector data
#include "headers/mosat_routines.h" //This is where most of the functions called in main are located
#include "headers/mosat_routines.h" //This is where most of the functions called in main are located
#include "headers/file_naming.h" //This has routines for added tags to an existing file name
#include "headers/file_naming_mpi.h" //This has routines for added tags to an existing file name (mpi)
#include "headers/command_line_args_mpi.h" //This has routines for adding command line arguments
#include "MosAT/program_variables/pv_lipid_contacts.h" //This has the variables specific to the analysis program
#include "MosAT/program_variables/pv_lipid_contacts.h" //This has the variables specific to the analysis program
#include "headers/array.h" //This has routines used for working with arrays
#include "headers/performance.h" //This has a class for logging performance data
#include "headers/index.h" //This has a class for working with index files
Expand All @@ -35,13 +35,14 @@ using namespace std;
#include "headers/protein.h" //This has routines used for working with protein data
#include "headers/force_serial.h" //This has routines used for forcing the code to run on a single mpi process
#include "headers/param.h" //This has routines used for reading complex parameter data
#include "headers/atom_select.h" //This has routines used for making atom selections using a selection text

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// //
// This function computes the number of lipid contacts and adds it to the grid //
// //
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &param,Grid &lcont)
void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &param,Grid &lcont,iv1d &custom_sel)
{
int i = 0; //standard variable used in loops
int j = 0; //standard variable used in loops
Expand Down Expand Up @@ -189,6 +190,42 @@ void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &p
}
}
}

//count contacts with a custom selection
if(p.b_sel_text == 1)
{
for(m=0; m<traj.atoms(); m++) //loop over system atoms
{
if(custom_sel[m] == 1)
{
if(b_counted[m] == 0 || p.b_mol == 0) //check that molecule has not been added already
{
double dif_x = traj.r[k][0] - traj.r[m][0];
double dif_y = traj.r[k][1] - traj.r[m][1];
double dif_z = traj.r[k][2] - traj.r[m][2];

distance = sqrt(dif_x*dif_x + dif_y*dif_y + dif_z*dif_z);

if(distance < p.contact_cutoff)
{
contacts++;

//flag that residue has been counted
if(p.b_mol == 1)
{
int min_1 = traj.res_start[traj.res_nr[m]-1];
int max_1 = traj.res_end[traj.res_nr[m]-1];

for(n=min_1; n<=max_1; n++) //loop over current residue
{
b_counted[n] = 1;
}
}
}
}
}
}
}
}
}
}
Expand Down Expand Up @@ -332,6 +369,7 @@ int main(int argc, const char * argv[])
add_argument_mpi_i(argc,argv,"-stdev", &p.b_stdev, "Compute the STDEV and STEM? (0:no 1:yes)", s.world_rank, nullptr, 0);
add_argument_mpi_i(argc,argv,"-clean", &p.b_clean, "Remove single frame files? (0:no 1:yes)", s.world_rank, nullptr, 0);
add_argument_mpi_i(argc,argv,"-mol", &p.b_mol, "Count residues instead of contacts? (0:no 1:yes)", s.world_rank, nullptr, 0);
add_argument_mpi_s(argc,argv,"-sel", p.selection_text_file_name, "Input file with the atom selection text (sel)", s.world_rank, &p.b_sel_text,0);
conclude_input_arguments_mpi(argc,argv,s.world_rank,s.program_name);

//create a trajectory
Expand Down Expand Up @@ -377,6 +415,10 @@ int main(int argc, const char * argv[])
{
check_extension_mpi(s.world_rank,"-sf_prm",p.solvent_finder_param_name,".prm");
}
if(p.b_sel_text == 1)
{
check_extension_mpi(s.world_rank,"-sel",p.selection_text_file_name,".sel");
}

//create parameter files
Param param;
Expand Down Expand Up @@ -424,6 +466,28 @@ int main(int argc, const char * argv[])

//print info about the grid
lcont.print_dim();

//create object to hold protein atom refinement
iv1d custom_sel(traj.atoms(),0);

if(p.b_sel_text == 1)
{
//create a object to hold an atom selection
Selection this_sel;

//select the atoms
this_sel.get_selection(traj,p.selection_text_file_name);

//generate pdb file name for highlighting the selection
string pdb_filename = chop_and_add_tag(p.selection_text_file_name,".pdb");

//highlight the selected atoms
this_sel.highlight_sel(traj,pdb_filename);

//refine the selection
custom_sel = this_sel.tag_atoms(traj);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////

//print info about the worlk load distribution
Expand All @@ -440,7 +504,7 @@ int main(int argc, const char * argv[])

traj.do_fit();

lip_cont(traj,s,p,param,lcont);
lip_cont(traj,s,p,param,lcont,custom_sel);

traj.set_beta_lf();

Expand Down
72 changes: 68 additions & 4 deletions src/MosAT/programs/lipid_contacts_3d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,11 @@ using namespace std;
#include "headers/switch.h" //This defines a switch (on, off)
#include "headers/file_reader.h" //This has basic routines for reading text files
#include "headers/vector_mpi.h" //This has routines for collecting vector data
#include "headers/mosat_routines.h" //This is where most of the functions called in main are located
#include "headers/mosat_routines.h" //This is where most of the functions called in main are located
#include "headers/file_naming.h" //This has routines for added tags to an existing file name
#include "headers/file_naming_mpi.h" //This has routines for added tags to an existing file name (mpi)
#include "headers/command_line_args_mpi.h" //This has routines for adding command line arguments
#include "MosAT/program_variables/pv_lipid_contacts_3d.h" //This has the variables specific to the analysis program
#include "MosAT/program_variables/pv_lipid_contacts_3d.h" //This has the variables specific to the analysis program
#include "headers/array.h" //This has routines used for working with arrays
#include "headers/performance.h" //This has a class for logging performance data
#include "headers/index.h" //This has a class for working with index files
Expand All @@ -36,13 +36,14 @@ using namespace std;
#include "headers/protein.h" //This has routines used for working with protein data
#include "headers/force_serial.h" //This has routines used for forcing the code to run on a single mpi process
#include "headers/param.h" //This has routines used for reading complex parameter data
#include "headers/atom_select.h" //This has routines used for making atom selections using a selection text

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// //
// This function computes the number of lipid contacts and adds it to the grid //
// //
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &param,Grid_3d &lcont)
void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &param,Grid_3d &lcont,iv1d &custom_sel)
{
int i = 0; //standard variable used in loops
int j = 0; //standard variable used in loops
Expand Down Expand Up @@ -191,6 +192,42 @@ void lip_cont(Trajectory &traj,system_variables &s,program_variables &p,Param &p
}
}
}

//count contacts with a custom selection
if(p.b_sel_text == 1)
{
for(m=0; m<traj.atoms(); m++) //loop over system atoms
{
if(custom_sel[m] == 1)
{
if(b_counted[m] == 0 || p.b_mol == 0) //check that molecule has not been added already
{
double dif_x = traj.r[k][0] - traj.r[m][0];
double dif_y = traj.r[k][1] - traj.r[m][1];
double dif_z = traj.r[k][2] - traj.r[m][2];

distance = sqrt(dif_x*dif_x + dif_y*dif_y + dif_z*dif_z);

if(distance < p.contact_cutoff)
{
contacts++;

//flag that residue has been counted
if(p.b_mol == 1)
{
int min_1 = traj.res_start[traj.res_nr[m]-1];
int max_1 = traj.res_end[traj.res_nr[m]-1];

for(n=min_1; n<=max_1; n++) //loop over current residue
{
b_counted[n] = 1;
}
}
}
}
}
}
}
}
}
}
Expand Down Expand Up @@ -337,6 +374,7 @@ int main(int argc, const char * argv[])
add_argument_mpi_i(argc,argv,"-clean", &p.b_clean, "Remove single frame files? (0:no 1:yes)", s.world_rank, nullptr, 0);
add_argument_mpi_i(argc,argv,"-mol", &p.b_mol, "Count residues instead of contacts? (0:no 1:yes)", s.world_rank, nullptr, 0);
add_argument_mpi_d(argc,argv,"-ex_val", &p.ex_val, "Set excluded lattice points to this value", s.world_rank, nullptr, 0);
add_argument_mpi_s(argc,argv,"-sel", p.selection_text_file_name, "Input file with the atom selection text (sel)", s.world_rank, &p.b_sel_text,0);
conclude_input_arguments_mpi(argc,argv,s.world_rank,s.program_name);

//create a trajectory
Expand Down Expand Up @@ -382,6 +420,10 @@ int main(int argc, const char * argv[])
{
check_extension_mpi(s.world_rank,"-sf_prm",p.solvent_finder_param_name,".prm");
}
if(p.b_sel_text == 1)
{
check_extension_mpi(s.world_rank,"-sel",p.selection_text_file_name,".sel");
}

//create parameter files
Param param;
Expand Down Expand Up @@ -436,6 +478,28 @@ int main(int argc, const char * argv[])

//print info about the grid
lcont.print_dim();

//create object to hold protein atom refinement
iv1d custom_sel(traj.atoms(),0);

if(p.b_sel_text == 1)
{
//create a object to hold an atom selection
Selection this_sel;

//select the atoms
this_sel.get_selection(traj,p.selection_text_file_name);

//generate pdb file name for highlighting the selection
string pdb_filename = chop_and_add_tag(p.selection_text_file_name,".pdb");

//highlight the selected atoms
this_sel.highlight_sel(traj,pdb_filename);

//refine the selection
custom_sel = this_sel.tag_atoms(traj);
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////

//print info about the worlk load distribution
Expand All @@ -452,7 +516,7 @@ int main(int argc, const char * argv[])

traj.do_fit();

lip_cont(traj,s,p,param,lcont);
lip_cont(traj,s,p,param,lcont,custom_sel);

traj.set_beta_lf();

Expand Down

0 comments on commit 1421d52

Please sign in to comment.