This Protocol is divided into following parts.
A. Run setup
B. Analysis of results
C. Energy analysis
D. Filtering best runs from multiple kappas
E. Filtering best runs from XY scaling and Z scaling
F. Helper Scripts

For help with syntax/usage:
Simply run the script without any arguments,
script_path is the path to top level scripts folder
eg. bash script_path/calculation/make_cafeinput_monomer.sh

Interactive login/ dependency:
sinteractive --gres=lscratch:500 --mem=20g
module load python/3.5
module load R/4.0

A) Run setup

#### Protocol for Running Cafemol
First of all, make a project directory and cd to this directory
If you are running for the first time, you have to generate the CG model, so proceed to step 1.
If you already have the CG model of pdb and the ninfo file, proceed to step 2.


1. All atom runs to create CG model and native info file

•       Warning: Please refine you structure in all atom using Xplor, before CG conversion
•	For monomer, copy the script_path/templates/monomer_ini folder to your project directory and cd to monome_ini
•       For dimmer, copy the script_path/templates/dimer_ini folder to your project directory and cd to dimer_ini
•	Specify the necessary file path at the top of the input file
•	No afm data needed for this step
•	The PDB file (all atom) works as is. No further processing needed
        Only keep the last model from the output PDB file, which will be the input CG pdb file for afm run.
        In case of dimer, introduce a TER word between two chain in the PDB for chain termination (cafemol is fussy about this).
•       Run cafemol:  path/cafemol.exe input.inp

2. AFM fitting run
•       Login to interactive node: sinteractive --gres=lscratch:500 --mem=10g
•       Prepare the input AFM file based on afm surface file from MountainSpip

        For monomer, Run: bash script_path/calculation/make_cafeinput_monomer.sh afm_file cg_pdb_file ninfo_file
        For dimer, Run: bash script_path/calculation/make_cafeinput_dimer.sh afm_file cg_pdb_file ninfo_file

        Note: If you are using the simulated data generated using cafemol
        For monomer, Run: bash script_path/calculation/make_cafeinput_monomer_simu.sh afm_file cg_pdb_file ninfo_file
        For dimer, Run: bash script_path/calculation/make_cafeinput_dimer_simu.sh afm_file cg_pdb_file ninfo_file
        Also, please do not remove the "x y z" header line in the simulated data file. No padding of afm image is done.

        afm_file -  is the afm data file from MountainSPIP, in nanometer units
        cg_pdb_file -  is the CG model of PDB obtained from allatom to cg calculation
        ninfo_file - is native info file from allatom to cg conversion calculation, edit ninfo_file if needed
•	Edit the input camemol file "input_cafe.inp" for md_steps, total atoms (look for GROUP(1)) and accordingly kappa, x/y pixel size and pixel numbers n_x, n_y.
        If you want to run in sumulated anneal mode, change i_run_mode to 3 and edit the section "annealing".
        Please note: You do not have to set up any path or raname any PDB or ninfo file. This is all automated now.
•	Edit the qsub script and submit the Batch Job.
        You can quickly check if cafemol run works in interactive mode, before submitting batch job.
•	After the analysis is complete, check the output folder.

3. 
   A) Scale the afm image data file and add gaussian noise to z column.
      module load R/4.0
      Rscript path script_path/analysis/scale.R [a,b,c] [p,q] afm_image.txt
      Usage:
       Rscript scale.R [a,b,c] [p,q] afm_image.txt
       image.txt    afm image file with headers x, y, z
       [a,b,c]      scale factors that multiply x, y, z values resp., scaling is applied before adding noise
       [p,q]        min and max limits of gaussian noise, p & q have to be same numbers with opposite signs
                    eg. [-10,10] will add gaussian distributed noise around mean 0 and limits -10 to +10, sd ~ 10/3 
                    noise is added to all z values (signal as well as background) 
                    Two final files ending in A.txt (angstoms units) and nm.txt (nanometer units) are produced

   B) Scale the afm image file and add uniform distribution noise within certain percent of the z value
   script path: script_path/analysis/addNoise_uniform.py
   Usage: pyXplor addNoise_uniform.py -noisePercent A -xscale P -yscale Q -zscale afm_file.txt
          eg. pyXplor addNoise_uniform.py -noisePercent 0.10 -xscale 2 -yscale 2 -zscale 2 cmpt1_5A_resample.txt
          The afm_file has headers x, y, z.
          Noise is only added if z values greater than 0.          

4. Edit the ninfo file to scale the force coefficients and add tertiary contacts
   (works for both monomer and dimer)
   set python environment: module load python/3.5
   i) To only scale the force coefficients:
      python script_path/calculation/edit_ninfo_full.py scale_factor ninfo_file

   ii) To scale coefficients and add tertiary restraints
      python script_path/calculation/edit_ninfo_full.py scale_factor ninfo_file restraint_file cg_pdbfile

   The scale_factor is given as [x,x,x,x,x,x,x] representing [bond, angle, dihedral_1, dihedral_2,long_range, basepair, basestack]   "x" respresents the multiplication factor, eg. [2,2,2,2,2,2,2]
   Templates for restraint file is given in the directory script_path/calculation

   iii) Separately scale the coefficients for duplex and non duplex regions
      module load python/3.5
      python script_path/calculation/scale_coef_duplex_ninfo.py scale_duplex scale_nonduplex ninfo_file basepair_file cg_pdbfile

      Both scale_duplex, and scale_nonduplex have the form [x,x,x,x,x,x,x] representing [bond, angle, dihedral_1, dihedral_2, long_range, basepair, basestack] and "x" respresents the multiplication factor.
      If all of the residues participating in bond, angle, long range, basepair and basestack belong to duplex, the duplex scaling factor is chosen. For dihedral, if both second and third resiude belong to the duplex, the duplex scaling factor is applied.

   iv) To remove tertiary contacts between a group of residues and rest of the molecule
      module load python/3.5
      python script_path/calculation/remove_contacts_ninfo.py  "A:1-20,B:30-40" ninfo_file cg_pdbfile

      In above example we choose chain A with residues 1 thr 20, and chain B with resiudes 30 thr 40, the script will remove all long range restraints between each of these groups and rest of the molecule.

B) Analysis of results
•       Login to interactive node: sinteractive --gres=lscratch:500 --mem=10g


   1) Plot AFMcc, rmsf, residue residue correlation, and generate related text files.
      Calculation of PDB files that score top N percent in AFMcc, and generate their pseudo images

   Monomer Run: bash script_path/analysis2/analyze_all_monomer_short.sh cafe_infile atom_type top_percent
   Dimer Run: bash script_path/analysis2/analyze_all_dimer_short.sh cafe_infile atom_type top_percent

Usage:
      bash analyze_all_monomer_short.sh/analyze_all_monomer_short.sh cafe_infile atom_type top_percent

      cafe_infile     cafemol input file used for calculation
      atom_type       choose one (sugar, phosphate, base, all)
      top_percent     interger number for extracting PDBs that score top N percent in afmCC

      A new folder top_N%_pdbs is created, which contains all results
      The cafe_infile should contain full path to the output folder containing dcd and timestep file, ninfo file


   2) Find rmsd of all PDBs in the dcd file compared to any given PDB. The atom fields in PDB format should match the DCD format.
   module load R/4.0
   Rscript script_path/analysis/get_rmsd.R pdb_file dcd_file atom_type
Usage:
       Rscript get_rmsd.R pdb_file dcd_file atom_type
       pdb_file     reference PDB file compare against
       dcd_file     trajectory file
       atom_type    choose one (sugar, phosphate, base, all)

   3) Core calculation and PCA analysis of PDB cohorts
      bio3d package was use to identify the stable core residues and carry out downstream PCA analysis.
      Script Path:
      script_path/analysis/core_calc_bio3d.R
      script_path/analysis/analyze_PCA_bio3d.R 



C) Energy analysis
   1) Plot Energy and output various energies as text files
      Script path: script_path/analysis2/plot_energy.R
      Usage:
             Rscript plot_energy.R timestep_file top_percent frames_cutoff

             timestep_file    timestep file from cafemol output
             top_percent      integer number for the top percent of afmcc structures
             frames_cutoff    number of first frames to remove for plotting

             Out files are en_all.txt (all frames), en_topNpercent.txt (top frames) and en_plots.pdf
             total energy vs afmcc plot: red dots showing top X percent pdbs with lowest afmcc
             normCC x Etot vs afmcc plot: red dots showing top 20 pdbs with lowest  normCC x Etot

   2) Plot mean and standard deviation of various energy terms of the top X percent afmcc strutures for each kappa values
   Rscript path: script_path/analysis/get_stats.R
Usage:
       Rscript get_stats.R output1 output2 output3  [k1,k2,k3,..] top_percent

       output1 output2 output3 etc. are the folder names corresonding to kappa values
       [k1,k2,k3,..]    k1, k2, k3 etc are the kappa values
       top_percent      integer number for the top percent of afmcc structures

       The number of folders and kappa values should match
       Each folder contains the timestep file and another folder named top_X%_pdbs (contains the rmsd file)
       The final output is a pdf file with summary statistics figures for energy terms such as correlation,
       total energy, local energy, etc against kappa

   3) Extract PDBs and plot various energy type based on correlation-energy (ce) score
        ce score is the product of normalized afmcc and normalized total Energy
   Rscript path: script_path/analysis/get_stats_CE.R

Usage:
      Rscript get_stats_CE.R input_folder top_percent pdbfile rmsd_file

      input_folder    this folder should contain the timestep and dcd file
      top_percent     integer representing the top percent of PDBs based on ce score
      pdbfile         pdbfile is a template PDB file for extracting coordinates from dcd
      rmsd_file       path to rmsd_file


      Folder top_Xpercent_ce_pdbs and plot top_Xpercent_ce_score_plots.pdf are generated
      ce_score and afmcc against RMSDs are plotted

    4) Calculate mean and sd of various energy terms for each kappa
        for the top X percent structures with lowest afmCC
        Rscript path: script_path/analysis2/get_stats_3D_short.R
Usage:
      Rscript get_stats_3D_short.R kappa_file top_percent
      kappa_file    file with list of kappa values
      top_percent   top X percent structures with lowest afmCC

      Also generate a plot and table of correlation between afmCC and (normalzed afmCC x total Energy) for each kappa
      Output files: 3Dstats_topXpercent.txt, correlation.txt, correlation.pdf

    5) Create a table with afmCC and (normalzed afmCC x total Energy) for each kappa
       Rscript path: script_path/analysis2/get_energyTally_3D.R
       Usage:
              Rscript get_energyTally_3D.R kappa_file frequency_factor
              kappa_file         file with list of kappa values
              frequency_factor   interger to print every Nth value from the kappa list

              Output files: energy_tally_3D.txt

D) Running multiple kappa runs and analysis to get the best kappa combinations

   1.  Setting up runs for multiple kappa values
   i) Generate copies of a template folder to run cafemol runs at various kappa values
   and submit all the runs.
   script path: script_path/calculation/setupNsubmit.sh
   Copy the file setupNsubmit.sh to your top directory containing the template folder.
   Edit the file bash file starting line "sbatch" to edit the batch submission parameters.
   Usage:
    bash setupNsubmit.sh template_directory kappa.txt

    template_directory    template folder that has a complete set up for the cafemol run
    kappa.txt             file with column of kappa values

    example of kappa.txt is given at script_path/calculation/
    please remove any empty lines at the end of file in kappa.txt

   2. Calculate top structures based on kappa filtering and frame filtering
      Kappa filtering is done based on linear regression and outlier removal by Cooks Distance
      get_gbStats_v1.sh does kappa filter sorting based on stack + hbond, get_gbStats_v2.sh uses stack + hbond + repulsion
      Frame filtering is based on sorting by Etot*CC^N or EstackEhbond*CC^N or EstackEhbondErepulsion*CC^N
      script path: script_path/analysis/get_gbStats_v1.sh
      Usage:

      Rscript get_gbStats_v1.sh kappa_file N pdb_cafemol pdb_rmsd atom_type F_cutoff M_cutoff N_kappa N_frame
      kappa_file     file with list of kappa values
      N              power to which afmCC is raised while multipying with Energy Terms
      pdb_cafemol    pdb file listed in the cafemol_input file
      pdb_rmsd       reference pdb file for calculating the rmsd against
      atom_type      atom type from the list (sugar, phosphate, base, all) for rmsd calculation
      F_cutoff       integer value representing the # of beginning frames to cutoff for frame filtering
      M_cutoff       integer value representing the number of frames to include for calculation summary statistics for each kappa
      N_kappa        integer value representing the number of top kappa runs to write out
      N_frame        integer value representing the number of top PDB frames to write out under each kappa
      Calculates mean and std dev summary statistics based on lowest 10 energy terms sorted by Etot - Afmfit
      Output files are in the folder top_statsStHb_PowN, frame filtering results are under subdirectory topPDB_Etot, topPDB_StckHb and topPDB_StckHbRepul
      And also under the output directory of each kappa/scan3D directory (Escore_ccPowN)
   
E)  Runing multiple runs with varying X,Y pixel scaling and Z scaling to find the best scaling paramters
    
   1. Setting up runs for multiple X,Y pixel scaling, Z scaling and kappa scaling.
   Generate copies of a template folder to run cafemol runs at various XY and Z scale factors and kappa values
   and submit all the runs

   Usage: 
      bash scan3D.sh '[x1,x2,dx] [z1,z2,dz] [k1,k2,dk] c templt_folder'

      x1,x2,dx              x1,x2 lower and upperbound for scale factor for x and y direction, dx is the increment interval
                            both x and y pixel coordinates are scaled by the same amount
      z1,z2,dz              z1, z2 are lower and upperbound for scale factor for z direction, dz is the increment interval
      k1,k2,dk              k1, k2 are lower and upperbound for kappa values, dk is the kappa increment interval
      c                     cutoff value for z, z values are first scaled and any scaled value of z below cutoff is set to 0
      templt_folder         template folder that has a complete set up for the cafemol run

      This script also generates a file called scan3Dlist.txt that contain all combinatios of xy and z scale factors and kappa values

   2. Calculate top structures based on filtering XY and Z Scaling Combinations for specific kappa
      Filtering XY and Z Scaling Combinations is based on sorting of mean afmfit energy and its sd
      Frame filtering is based on sorting by Etot*CC^N or EstackEhbond*CC^N or EstackEhbondErepulsion*CC^N
      Usage:
      script path: script_path/analysis/get_gbStats_scan2D.sh

       Rscript get_gbStats_scan2D.sh scan3Dlist.txt N pdb_cafemol pdb_rmsd atom_type F_cutoff M_cutoff K_input N_frame
       scan3Dlist.txt file with list of 3D scan directories values
       N              power to which afmCC is raised while multipying with Energy Terms
       pdb_cafemol    pdb file listed in the cafemol_input file
       pdb_rmsd       reference pdb file for calculating the rmsd against
       atom_type      atom type from the list (sugar, phosphate, base, all) for rmsd calculation
       F_cutoff       integer value representing the # of beginning frames to cutoff for frame filtering
       M_cutoff       integer value representing the number of frames to include for calculation summary statistics for each kappa
       K_input        input kappa value upon which the filtering of XY and Z pixel combinations is carried
       N_comb        integer value representing the number of top XY and Z pixel combinations to write out
       N_frame        integer value representing the number of top PDB frames to write out under each kappa
       Calculates mean and std dev summary statistics based on lowest 10 energy terms sorted by Etot - Afmfit
       Output files are in the folder top_statsStHb11_PowN, frame filtering results are under subdirectory topPDB_Etot, topPDB_StckHb and topPDB_StckHbRepul
       And also under the output directory of each kappa/scan3D directory (Escore_ccPowN, FrameFilter_ccPowN

       The energy and rmsd values can be used to generate contour plots by using the follwing script.


       Generate the contour plots for afmfit energy and RMSD.
       i)Generate a contour plot based on mean afmfit energy + scaled SD
       script path: script_path/analysis
       Usage:
       Rscript contour_plot_afmfit.R summary_stats_file cutoff_val 
       summary_stats_file    file containing the sumarry_statistics information containing means energies and SD
       cutoff_val            value abouve which the contour plot is saturated in color 
       A pdf file called contour_cutoff_val is generated 
       The summary_stats_file is summary_stats_sort.txt generated by the script get_gbStats_scan2D.sh 

       ii)Generate a contour plot based on RMSD
       Usage:
       Rscript contour_plot_afmfit.R rmsd_file cutoff_val 
       rmsd_file     file containing the RMSD
       cutoff_val    value above which the contour plot is saturated in color 
       A pdf file called contour_cutoff_val is generated 
       The rmsd_file is topPDB_StckHbRepul.txt generated by the script get_gbStats_scan2D.sh 

     
F) Helper scripts


1. Generate 3D open structure from given secondary structure information

   Script path: script_path/rna2d3d/convert2d3d.sh
Usage:
      bash convert2d3d.sh input.bpl

      input.bpl      input file with the sequence and base pairs

      see the template input.bpl at script_path/rna2d3d/input.bpl
      the nulceotide sequence should have a maximum of 50 spaces in each line
      the output pdb file has extension is _open.pdb

2. Use Xplot-NIH to refine the 3D open structure from rna2d3d

   Make a directory called "refine" and copy the files refine.py, qsub1.sh, qsub2.sh from script_path/rna2d3d/refine
   to this directory.
   provide the files in current directory: input_open.pdb and input.bpl, you can edit their names in file qsub1.sh.
   As the name suggest input_open.pdb is the 3D open structure from rna2d3d
   input.bpl is the base pair list file used in rna2d3d before.
   Run the sbatch command at the bottom of the file qsub1.sh.
   The final refined file generated is refine_final.pdb

3. Glob to all atom conversion and refinement
   Script path: script_path/cg2aa/
   For Momomer:
   Carry out the glob to all atom conversion and xplor refinement

   Usage:
         bash master_aa_monomer.sh input.bpl cgRNA_list_file

         input.pdb         base pair list file
         cgRNA_list_file   file containging full path of a list of cg PDBs
                           example at script_path/cg2aa/cgPDB_list.txt
                           glob to allatom conversion and Xplor refinement is iterated for each PDB
                           The refined PDB files have the format refine_**.pdb, will produce 5 PDBs per refinement
  For Dimer:
  Carry out the glob to all atom conversion and xplor refinement

  Usage:
        bash master_aa_dimer.sh input.bpl cgRNA_list_file

        input.pdb         base pair list file
        cgRNA_list_file   file containging full path of a list of cg PDBs
                          example at script_path/cg2aa/cgPDB_list.txt
                          glob to allatom conversion and Xplor refinement is iterated for each PDB
                          The refined PDB files have the format refine_**.pdb, will produce 5 PDBs per refinement

4. Fix covalent geometry and minimize the RNA structure
   Regularizes bond, angle and improper terms, not allowing significant motion of the P atoms
   Usage: xplor script_path/analysis/regularize.py FILE.pdb > regularize.out
   where you specify FILE.pdb. The regularized structure is written to FILE.pdb.out

   Refinement of group of PDBs:
   Script path: script_path/analysis/master_refine.sh
   Carry out all atom refinement using Xplor
   regularizes bond, angle and improper terms, not allowing significant motion of the P atoms

   Usage: 
      bash master_refine.sh pdb_list

      pdb_list          list of all atom PDBs with full path
                        example at script_path/analysis/aaPDBlist.txt
                        The refined PDB files have extension _ref.pdb

5. Calculate the AFM image for individual coarse grained PDB file
   Script path: script_path/analysis/
   Process a single PDB to produce the afm image data
   Usage:
         ## Create the AFM image from PDB
         bash script_path/analysis/pdb2image_monomer_unit.sh pdb_file cafemol_input_file (Monomer)
         bash script_path/analysis/pdb2image_dimer_unit.sh pdb_file cafemol_input_file (Dimer)

         pdb_file              coarse grained pdb file recentered in the afm field
         cafemol_input_file    input file for cafemol

6. Calculate the AFM image for individual coarse grained PDB file based on user defined image size
   Script path: script_path/analysis2/
Usage:
      bash pdb2image_monomer_setField.sh pixel_size max_dim pdb_file cafemol_inputfile (Monomer)

      pixel_size            pixel size along x and y (integer or decimal)
      max_dim               maximum dimension along x and y (integer or decimal)
      pdb_file              coarse grained pdb file
      cafemol_input_file    input file for cafemol

7. Format CG model for cafemol run
   Usage:
       Rscript format_pdb_cg.R input.pdb 

       input.pdb           PDB file to format
       Input PDB file is formated in place, i.e. modified pdb has the same name
8. Rotate the given PDB by x,y,z axes 
   script path:   script_path/calculation/rotate_PDB.R
   Usage:
       Rscript rotate_PDB.R [a,b,c] pdb_file
       [a,b,c]     a,b,c are the angles in degrees(+ve or -ve) to rotate about x,y and z axes resp.
       pdb_file    the pdb file to rotate
       The output pdb file is given extension of _rot 
       The order of rotation is X, Y and Z axes resp. 
       [0,0,0] denotes no rotation 

9. Subtract a constant number from X,Y, or Z columns in a PDB file 
   script path:   script_path/analysis/edit_pdb.R
   Usage: 
      Rscript edit_pdb.R pdb_file axes Num
 
      pdb_file   input pdb file 
      axes       either X, Y or Z 
      Num        constant number to subtract
      Example: Rscript edit_pdb.R pdb_file Z 20 

