Difference between revisions of "MutEn.pl"
(3 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
+ | NOTE: Do not forget to substitute "2FBT" with your own PDB code. |
||
− | #!/usr/bin/perl -w |
||
+ | 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. |
||
− | my $PDB='2FBT'; # Insert PDB code for protein |
||
+ | |||
− | my @fileslist = <${PDB}_*>; # Create an array with all the mutant structures. |
||
+ | |||
+ | |||
+ | #!/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 ) { |
foreach my $file ( @fileslist ) { |
||
Line 9: | Line 17: | ||
print "Creating ID from file name\n\n"; |
print "Creating ID from file name\n\n"; |
||
− | # Create ID from file name |
+ | # Create ID from file name |
my @fullname = split(/.pdb/, $file); # Splits the file name in two parts. |
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"). |
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"; |
print "Creating sequence file as an input for SCWRL\n\n"; |
||
− | # Creates sequence file as an input for SCWRL: |
+ | # Creates sequence file as an input for SCWRL: |
− | + | # Extract sequence from PDB file: |
|
− | `repairPDB ${id}.pdb -seq >> ${id}.fasta`; # Creates sequence file. |
+ | `repairPDB ${id}.pdb -seq >> ${id}.fasta`; # Creates sequence file. |
+ | |||
− | |||
− | + | # It is now necessary to convert the sequence to lower case: |
|
− | # Start |
+ | # Start |
− | open(SEQ, "<${id}.fasta"); # Opens the sequence file. |
+ | open(SEQ, "<${id}.fasta"); # Opens the sequence file. |
− | open(SEQLC, ">${id}_LC.fasta"); # Opens a new file to store the lower case sequence. |
+ | 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 $sequence=<SEQ>; # Stores the sequence in a string. |
− | my $sequencelc = lc($sequence); # Converts the sequence to lower case. |
+ | my $sequencelc = lc($sequence); # Converts the sequence to lower case. |
− | print SEQLC $sequencelc; # Prints the lower case sequence string to file. |
+ | print SEQLC $sequencelc; # Prints the lower case sequence string to file. |
close SEQ; close SEQLC; |
close SEQ; close SEQLC; |
||
− | # End |
+ | # End |
print "Runing SCWRL to assert the position of side-chains\n\n"; |
print "Runing SCWRL to assert the position of side-chains\n\n"; |
||
− | # Runs SCWRL to assert the position of side-chains: |
+ | # Runs SCWRL to assert the position of side-chains: |
`scwrl -s ${id}_LC.fasta -i ${id}.pdb -o ${id}_Scwrl.pdb`; |
`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"; |
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: |
+ | # 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. |
`repairPDB ${id}_Scwrl.pdb -nosol -noh >> ${id}_NoSolNoH.pdb`; # Removes Water molecules and Hydrogen atoms from structure. |
||
print "Runing minimization with Gromacs\n"; |
print "Runing minimization with Gromacs\n"; |
||
− | # Run minimization with Gromacs: |
+ | # Run minimization with Gromacs: |
print "\tRuning PDB2GMX\n"; |
print "\tRuning PDB2GMX\n"; |
||
− | # PDB2GMX |
+ | # PDB2GMX |
`pdb2gmx -f ${id}_NoSolNoH.pdb -o ${id}.gro -p ${id}.top -water tip3p -ff amber03`; |
`pdb2gmx -f ${id}_NoSolNoH.pdb -o ${id}.gro -p ${id}.top -water tip3p -ff amber03`; |
||
print "\tCreating MDP file\n"; |
print "\tCreating MDP file\n"; |
||
− | # Create MDP file: |
+ | # Create MDP file: |
open(MDP,">Minimize.mdp"); |
open(MDP,">Minimize.mdp"); |
||
print MDP "title = PBSA minimization in vacuum\n"; |
print MDP "title = PBSA minimization in vacuum\n"; |
||
Line 64: | Line 72: | ||
print "\tRuning GROMPP\n"; |
print "\tRuning GROMPP\n"; |
||
− | # GROMPP |
+ | # GROMPP |
`grompp -v -f Minimize.mdp -c ${id}.gro -p ${id}.top -o ${id}.tpr`; |
`grompp -v -f Minimize.mdp -c ${id}.gro -p ${id}.top -o ${id}.tpr`; |
||
print "\tRuning MDRUN\n"; |
print "\tRuning MDRUN\n"; |
||
− | # MDRUN |
+ | # MDRUN |
`mdrun -v -deffnm ${id}`; |
`mdrun -v -deffnm ${id}`; |
||
print "\tAnalyzing\n\n\n\n"; |
print "\tAnalyzing\n\n\n\n"; |
||
− | # Analyze |
+ | # Analyze |
− | `g_energy -f ${id}.edr -o ${id}_Energy.xvg`; |
+ | `g_energy -f ${id}.edr -o ${id}_Energy.xvg <<EOF\n1\n2\n12\n0\nEOF`; |
} |
} |
Latest revision as of 19:29, 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 <<EOF\n1\n2\n12\n0\nEOF`;
} </nowiki>