Difference between revisions of "MutEn.pl"
Line 1: | Line 1: | ||
NOTE: Do not forget to substitute "2FBT" with your own PDB code. |
NOTE: Do not forget to substitute "2FBT" with your own PDB code. |
||
The mutant structures must be in the same folder as the script and named PDB_SOMETHING.pdb (where "PDB" must be the PDB code for your protein and "SOMETHING" can be anything to differ each mutation from the rest). |
The mutant structures must be in the same folder as the script and named PDB_SOMETHING.pdb (where "PDB" must be the PDB code for your protein and "SOMETHING" can be anything to differ each mutation from the rest). |
||
+ | Do not forget to save the script with a perl extension (example: script.pl) |
||
+ | To get this script, the easiest way is to click the "edit" button and copy everything below this sentence. |
||
− | |||
− | To copy this script, click the "edit" button and copy everything below this sentence. |
||
Revision as of 19:03, 3 July 2011
NOTE: Do not forget to substitute "2FBT" with your own PDB code. The mutant structures must be in the same folder as the script and named PDB_SOMETHING.pdb (where "PDB" must be the PDB code for your protein and "SOMETHING" can be anything to differ each mutation from the rest). Do not forget to save the script with a perl extension (example: script.pl)
To get this script, the easiest way is to click the "edit" button and copy everything below this sentence.
#!/usr/bin/perl -w
my $PDB='2FBT'; # Insert PDB code for protein my @fileslist = <${PDB}_*>; # Create an array with all the mutant structures.
foreach my $file ( @fileslist ) {
print "--> FILE ${file}\n\n\n";
print "Creating ID from file name\n\n";
# Create ID from file name
my @fullname = split(/.pdb/, $file); # Splits the file name in two parts. my $id = $fullname[0]; # Identifier of file name is the first part (before ".pdb").
print "Creating sequence file as an input for SCWRL\n\n";
# Creates sequence file as an input for SCWRL:
# Extract sequence from PDB file: `repairPDB ${id}.pdb -seq >> ${id}.fasta`; # Creates sequence file.
# It is now necessary to convert the sequence to lower case: # Start open(SEQ, "<${id}.fasta"); # Opens the sequence file. open(SEQLC, ">${id}_LC.fasta"); # Opens a new file to store the lower case sequence. my $sequence=<SEQ>; # Stores the sequence in a string. my $sequencelc = lc($sequence); # Converts the sequence to lower case. print SEQLC $sequencelc; # Prints the lower case sequence string to file. close SEQ; close SEQLC; # End
print "Runing SCWRL to assert the position of side-chains\n\n";
# Runs SCWRL to assert the position of side-chains:
`scwrl -s ${id}_LC.fasta -i ${id}.pdb -o ${id}_Scwrl.pdb`;
print "Remove Solvent molecules and Hydrogen atoms to prepare the structure for energy minimization\n\n";
# Remove Solvent molecules and Hydrogen atoms to prepare the structure for energy minimization:
`repairPDB ${id}_Scwrl.pdb -nosol -noh >> ${id}_NoSolNoH.pdb`; # Removes Water molecules and Hydrogen atoms from structure.
print "Runing minimization with Gromacs\n";
# Run minimization with Gromacs:
print "\tRuning PDB2GMX\n";
# PDB2GMX
`pdb2gmx -f ${id}_NoSolNoH.pdb -o ${id}.gro -p ${id}.top -water tip3p -ff amber03`;
print "\tCreating MDP file\n";
# Create MDP file:
open(MDP,">Minimize.mdp"); print MDP "title = PBSA minimization in vacuum\n"; print MDP "cpp = /usr/bin/cpp\n"; print MDP "define = -DFLEXIBLE -DPOSRES\n"; print MDP "implicit_solvent = GBSA\n"; print MDP "integrator = steep\n"; print MDP "emtol = 1.0\n"; print MDP "nsteps = 500\n"; print MDP "nstenergy = 1\n"; print MDP "energygrps = System\n"; print MDP "ns_type = grid\n"; print MDP "coulombtype = cut-off\n"; print MDP "rcoulomb = 1.0\n"; print MDP "rvdw = 1.0\n"; print MDP "constraints = none\n"; print MDP "pbc = no"; close MDP;
print "\tRuning GROMPP\n";
# GROMPP
`grompp -v -f Minimize.mdp -c ${id}.gro -p ${id}.top -o ${id}.tpr`;
print "\tRuning MDRUN\n";
# MDRUN
`mdrun -v -deffnm ${id}`;
print "\tAnalyzing\n\n\n\n";
# Analyze
`g_energy -f ${id}.edr -o ${id}_Energy.xvg`;
} </nowiki>