Perl script for ATM

<pre>#!/usr/bin/perl -w

use strict;

doMain();

sub doMain
{
my $baseDirPath = shift @ARGV || "";
my $pirAlignmentFile = shift @ARGV || "";
my $targetProteinCode = shift @ARGV || "";
my $pdbFiles_ref = \@ARGV || ();

print "----\n";
print "\nChecking command-line parameters ...\n";
printUsageAndExit() unless ( $baseDirPath && -d $baseDirPath );
$baseDirPath =~ s/\/$//; # remove the ending slash if exists

printUsageAndExit() unless ( $pirAlignmentFile && -f "$baseDirPath/$pirAlignmentFile" && $pirAlignmentFile =~ /\.pir$/ );
printUsageAndExit() unless ( $targetProteinCode && length($targetProteinCode) == 4 );
printUsageAndExit() unless ( @$pdbFiles_ref > 0 );

foreach my $pdbFile ( @$pdbFiles_ref )
{
printUsageAndExit() unless ( -f "$baseDirPath/$pdbFile" && $pdbFile =~ /\.pdb$/ );
}

print "Making modeller input files ...\n";
makeModellerInput( $baseDirPath, $pirAlignmentFile, $targetProteinCode, $pdbFiles_ref );
}

sub printUsageAndExit
{
print STDERR "Error!\n\n";
print STDERR "Usage: perl -w $0 A1_baseDirectoryPath A2_pirAlignmentFile A3_targetProteinCode A4_pdbFiles\n\n";
print STDERR " A1_baseDirectoryPath: directory path containing all the PIR alignment file and pdb files\n\n";
print STDERR " A2_pirAlignmentFile: PIR-format output file (with the extension, '.pir') of multiple sequence alignment including a four-letter code in each header line (e.g. cdk2HumanHomologs)\n\n";
print STDERR " A3_targetProteinCode: the same four-letter code of your target protein as that in the alignment file (e.g. trgt or mine)\n\n";
print STDERR " A4_pdbFiles: pdb files (with the extension, '.pdb') separated by blank of which file name consists of the same four-letter code as that in the alignment file (e.g. 1erk.pdb 1p38.pdb)\n\n";
print STDERR "----\n";
exit;
}

sub makeModellerInput
{
my ( $baseDirPath, $pirAlignmentFile, $targetProteinCode, $pdbFiles_ref ) = @_;
print "\tchanging the directory to $baseDirPath ...\n";
chdir $baseDirPath or die "Dir Change Error, $baseDirPath: $!";

my ( $outputPrefix ) = $pirAlignmentFile =~ /^(\S+)\.\S+$/;

# modeller input files
my $alignmentFile = $outputPrefix. ".ali";
my $topScriptFile = $outputPrefix. ".top";
my %pdbCode2AtomFile = get_pdbCode2File( $pdbFiles_ref, "atm" );

my %pdbCode2PdbFile = get_pdbCode2File( $pdbFiles_ref, "pdb" );
my %pdbCode2PirAlignment = get_pdbCode2PirAlignment( $pirAlignmentFile, $targetProteinCode, \%pdbCode2PdbFile );

# parse the target protein's alignment and information first
my $alignmentFileContent = getTargetProteinAlignment( $targetProteinCode, $pdbCode2PirAlignment{$targetProteinCode} );

# parse alignments and information of known proteins and make atom files of them
foreach my $pdbCode ( keys %pdbCode2PdbFile )
{
my $atomFile = $pdbCode2AtomFile{ $pdbCode };
my $pdbFile = $pdbCode2PdbFile{ $pdbCode };
my $pirAlignment = $pdbCode2PirAlignment{ $pdbCode };

# make the atom files
print "\tmaking the atom file, $atomFile ...\n";
makeAtomFile( $pdbFile, $atomFile );

$alignmentFileContent .= get_eachAlignmentContent( $pdbCode, $pirAlignment, $pdbFile );
}

# make the alignment file
print "\tmaking the alignment file, $alignmentFile ...\n";
makeAlignmentFile( $alignmentFile, $alignmentFileContent );

# make the top script file
print "\tmaking the top script file, $topScriptFile ...\n";
makeTopScriptFile( $topScriptFile, $alignmentFile, $targetProteinCode, \%pdbCode2AtomFile );

print "Finished!\n\n";
print "----\n";
}

sub get_pdbCode2File
{
my ( $pdbFiles_ref, $fileExtension ) = @_;
my %pdbCode2File = ();

foreach my $pdbFile ( @$pdbFiles_ref )
{
my ( $pdbCode ) = $pdbFile =~ /(\w{4})\.pdb$/; # extract only the four-letters
my $file = $pdbCode. ".". $fileExtension;
$pdbCode2File{ $pdbCode } = $file;
}

return %pdbCode2File;
}

sub get_pdbCode2PirAlignment
{
my ( $pirAlignmentFile, $targetProteinCode, $pdbCode2PdbFile_ref ) = @_;
my %pdbCode2PirAlignment = ();

my @pdbCodes = ();
push @pdbCodes, $targetProteinCode;
push @pdbCodes, (sort keys %$pdbCode2PdbFile_ref );

$/ = "*\n"; # change the reading-block separator
open FI, $pirAlignmentFile or die "File Open Error, $pirAlignmentFile: $!\n";

while ( my $entry = <fi> )
{
chomp $entry;
next if $entry =~ /^\s*$/;

my @lines = split /\n/, $entry;

my $pdbCode = "";
foreach my $code ( @pdbCodes )
{
if ( $lines[0] =~ /$code/ )
{
$pdbCode = $code;
last;
}
}

if ( $pdbCode eq "" )
{
print STDERR "PDB code mistmatch with pdb files on this entry!\n";
print STDERR "pdbCode: $pdbCode\n";
print STDERR "entry: $entry\n";
}

$pdbCode2PirAlignment{ $pdbCode } = "$entry*\n";
print "\nPDB code To PIR alignment found: $pdbCode =>\n$entry\n";
}

close FI;
$/ = "\n"; # restore the reading-block separator

return %pdbCode2PirAlignment;
}

sub getTargetProteinAlignment
{
my ( $targetProteinCode, $targetProteinAlignment ) = @_;
my $targetProteinAlignmentContent = "";

my @lines = split /\n/, $targetProteinAlignment;
shift @lines; # remove header line
shift @lines; # remove blank line
my $sequenceAlignmentOnly = join "\n", @lines;

$targetProteinAlignmentContent .= ">P1;$targetProteinCode\n";
$targetProteinAlignmentContent .= "sequence: : : : : : : : :\n";
$targetProteinAlignmentContent .= "$sequenceAlignmentOnly\n";

return $targetProteinAlignmentContent;
}

sub makeAtomFile
{
my ( $pdbFile, $atomFile ) = @_;
 !system("cp $pdbFile $atomFile") or die "File Copy Error: $!";
print "\tcopying $pdbFile to $atomFile ...\n";
}

sub get_eachAlignmentContent
{
my ( $pdbCode, $pirAlignment, $pdbFile ) = @_;
my $eachAlignmentContent = "";

my @lines = split /\n/, $pirAlignment;
shift @lines; # remove header line
shift @lines; # remove blank line
pop @lines; # remove ending *
my $sequenceAlignmentOnly = join "\n", @lines;
$sequenceAlignmentOnly =~ s/[Xx]/-/g;

my ( $structureInfoLine, $adjustedAlignmentSequence )
= get_structureInfoAndAdjustedSequence( $pdbCode, $sequenceAlignmentOnly, $pdbFile );
$eachAlignmentContent .= ">P1;$pdbCode\n";
$eachAlignmentContent .= $structureInfoLine;
$eachAlignmentContent .= "$adjustedAlignmentSequence\n*\n";

return $eachAlignmentContent;
}

sub get_structureInfoAndAdjustedSequence
{
my ( $pdbCode, $sequenceAlignmentOnly, $pdbFile ) = @_;
my ( $structureInfoLine, $adjustedAlignmentSequence ) = ( "", "" );

my ( $startResNum, $startChain, $endResNum, $endChain ) = ("", "", "", "");
my ( $structureXorN, $proteinName, $source, $resolution, $rFactor ) = ("", "", "", "", "");
my $allResiduesInPdb = "";
my $includedResNum = 0;
my $isChainExist = 1; # true


open FI, $pdbFile or die "File Open Error, $pdbFile: $!\n";
while ( my $line = <fi> )
{
chomp $line;

if ($line =~ /^COMPND.*MOLECULE:\s*(\S.*\S);\s*$/)
{
$proteinName = $1;
}
elsif ($line =~ /^COMPND.*CHAIN:\s+NULL;\s*$/)
{
$isChainExist = 0; # false
}
elsif ($line =~ /ORGANISM_SCIENTIFIC:\s*(\S.*\S);\s*$/)
{
$source = $1;
}
elsif ($line =~ /^EXPDTA\s+(\S.*\S)\s*$/)
{
my $method = $1;
$structureXorN =
( $method =~ /X-RAY/ ) ? "structureX" :
( $method =~ /NMR/ ) ? "structureN" :
"structure";
}
elsif ($line =~ /^REMARK.*RESOLUTION.\s+([\d\.]+)\s+ANGSTROM.*$/)
{
$resolution = $1;
}
elsif ($line =~ /^REMARK\s+\d+\s+R VALUE\s+\(WORKING SET\)\s*:\s*([\d\.\-]+)\s*$/)
{
$rFactor = $1;
}
elsif ($line =~ /^ATOM/)
{
my @values = split /\s+/, $line;
my $threeLetterAminoAcidCode = $values[3];
my $oneLetterAminoAcidCode = convertThreeToOneLetterAminoAcidCode($threeLetterAminoAcidCode);

if ( $startResNum eq "" && $isChainExist )
{
( $startChain, $startResNum ) = ( $values[4], $values[5] );
}
elsif ( $startResNum eq "" )
{
$startResNum = $values[4];
}

if ( $isChainExist )
{
( $endChain, $endResNum ) = ( $values[4], $values[5] );
}
else
{
$endResNum = $values[4];
}

if ( $endResNum ne $includedResNum )
{
$allResiduesInPdb .= $oneLetterAminoAcidCode;
$includedResNum = $endResNum;
}
}
elsif ($line =~ /^TER/) # read only the atom list ( if multi-chain protein, only the first chain is read )
{
last;
}
}
close FI;

$structureInfoLine = "$structureXorN:$pdbCode:$startResNum :$startChain :$endResNum :$endChain ".
":$proteinName :$source :$resolution :$rFactor \n";
$adjustedAlignmentSequence = get_adjustedAlignmentSequence( $sequenceAlignmentOnly, $allResiduesInPdb );

return ( $structureInfoLine, $adjustedAlignmentSequence );
}

sub convertThreeToOneLetterAminoAcidCode
{
my ( $threeLetterAminoAcidCode ) = @_;
return
( $threeLetterAminoAcidCode =~ /gly/i ) ? "G" :
( $threeLetterAminoAcidCode =~ /pro/i ) ? "P" :
( $threeLetterAminoAcidCode =~ /ala/i ) ? "A" :
( $threeLetterAminoAcidCode =~ /val/i ) ? "V" :
( $threeLetterAminoAcidCode =~ /leu/i ) ? "L" :
( $threeLetterAminoAcidCode =~ /ile/i ) ? "I" :
( $threeLetterAminoAcidCode =~ /met/i ) ? "M" :
( $threeLetterAminoAcidCode =~ /cys/i ) ? "C" :
( $threeLetterAminoAcidCode =~ /phe/i ) ? "F" :
( $threeLetterAminoAcidCode =~ /tyr/i ) ? "Y" :
( $threeLetterAminoAcidCode =~ /trp/i ) ? "W" :
( $threeLetterAminoAcidCode =~ /his/i ) ? "H" :
( $threeLetterAminoAcidCode =~ /lys/i ) ? "K" :
( $threeLetterAminoAcidCode =~ /arg/i ) ? "R" :
( $threeLetterAminoAcidCode =~ /gln/i ) ? "Q" :
( $threeLetterAminoAcidCode =~ /asn/i ) ? "N" :
( $threeLetterAminoAcidCode =~ /glu/i ) ? "E" :
( $threeLetterAminoAcidCode =~ /asp/i ) ? "D" :
( $threeLetterAminoAcidCode =~ /ser/i ) ? "S" :
( $threeLetterAminoAcidCode =~ /thr/i ) ? "T" :
"-";
}

sub get_adjustedAlignmentSequence
{
my ( $sequenceAlignmentOnly, $allResiduesInPdb ) = @_;
my $adjustedAlignmentSequence = "";

my @lines = split /\n/, $sequenceAlignmentOnly;
my $charactersPerLine = length $lines[0];
$sequenceAlignmentOnly =~ s/\n//g; # remove new line characters

# substitute front and end residues without atom coordinates to '-'s
my $startPos = get_startPos( $sequenceAlignmentOnly, $allResiduesInPdb );
my $endPos = get_endPos( $sequenceAlignmentOnly, $allResiduesInPdb );
$adjustedAlignmentSequence = "-" x ($startPos - 1).
substr( $sequenceAlignmentOnly, $startPos - 1, $endPos - $startPos + 1 ). # substr is zero-based
"-" x ( length($sequenceAlignmentOnly) - $endPos );

# substitute loop residues without atom coordinates to '-'s
my @loopPositions = get_loopPositions( $adjustedAlignmentSequence, $allResiduesInPdb );
foreach my $loopPosition ( @loopPositions )
{
my ( $loopStartPos, $loopEndPos ) = split /-/, $loopPosition; # e.g. loop residues 45-56
$adjustedAlignmentSequence =
substr( $adjustedAlignmentSequence, 0, $loopStartPos - 1 ).
"-" x ( $loopEndPos - $loopStartPos + 1).
substr( $adjustedAlignmentSequence, $loopEndPos, length($adjustedAlignmentSequence) - $loopEndPos );
}

my $dashRemovedAdjustedAlignmentSequence = $adjustedAlignmentSequence;
$dashRemovedAdjustedAlignmentSequence =~ s/-//g;

if ( $dashRemovedAdjustedAlignmentSequence ne $allResiduesInPdb )
{
print STDERR "\nError!!!\n";
print STDERR "Number of residues in the PIR-alignment file and the PDB file mismatches !\n";
print STDERR "Substitue residues in the PIR-alignment file without atom coordinates to dashes('-'s) manually!\n\n";
print STDERR "dashRemovedAdjustedAlignmentSequence: $dashRemovedAdjustedAlignmentSequence\n";
print STDERR "allResiduesInPdb: $allResiduesInPdb\n";
}

return formatFastaSequence($adjustedAlignmentSequence, $charactersPerLine);
}

sub get_startPos
{
my ( $sequenceAlignmentOnly, $allResiduesInPdb ) = @_;

# Hypothesis: if ten residues match in the front end, we can determine the start position in the alignment
# matching with the residues in the PDB file
my $firstTenAminoAcidsOfAllResiduesInPdb = substr( $allResiduesInPdb, 0, 10 );

my $startPos = "";
my @allAminoAcidsInTheAlignmentSequence = split //, $sequenceAlignmentOnly;
my $pos = 0;

foreach my $aa ( @allAminoAcidsInTheAlignmentSequence )
{
$pos++;
next if $aa eq "-";

my $frontRemovedSeq = substr( $sequenceAlignmentOnly, $pos - 1, length($sequenceAlignmentOnly) - $pos + 1 );
$frontRemovedSeq =~ s/-//g; # remove all dashes

if ( $frontRemovedSeq =~ /^$firstTenAminoAcidsOfAllResiduesInPdb/i )
{
$startPos = $pos;
print "\tstart position or non-loop start position determined: $pos\n";
last;
}
}

if ( $startPos eq "" )
{
print STDERR "\nError!\n";
print STDERR "Start position mismatches!\n";
print STDERR "Check the start amino acids of the alignment!\n\n";
print STDERR "sequenceAlignment: $sequenceAlignmentOnly\n";
print STDERR "residuesInPdb: $allResiduesInPdb\n\n";
}

return $startPos;
}

sub get_endPos
{
my ( $sequenceAlignmentOnly, $allResiduesInPdb ) = @_;

# Hypothesis: if ten residues match in the rear end, we can determine the end position in the alignment
# matching with the residues in the PDB file
my $reversedLastTenAminoAcidsOfAllResiduesInPdb =
reverseMyString(substr( $allResiduesInPdb, length($allResiduesInPdb) - 10, 10 ));

my $endPos = "";
my @allAminoAcidsInTheAlignmentSequence = split //, $sequenceAlignmentOnly;
my $pos = @allAminoAcidsInTheAlignmentSequence + 1;

foreach my $aa ( reverse @allAminoAcidsInTheAlignmentSequence )
{
$pos--;
next if $aa eq "-";

my $endRemovedReverseSeq = reverseMyString( substr( $sequenceAlignmentOnly, 0, $pos ) );
$endRemovedReverseSeq =~ s/-//g; # remove all dashes

if ( $endRemovedReverseSeq =~ /^$reversedLastTenAminoAcidsOfAllResiduesInPdb/i )
{
$endPos = $pos;
print "\tend position determined: $pos\n";
last;
}
}

if ( $endPos eq "" )
{
print STDERR "\nError!\n";
print STDERR "End position mismatches between the alignment file and the pdb file!\n";
print STDERR "Check the end amino acids of the alignment!\n\n";
print STDERR "sequenceAlignmentOnly: $sequenceAlignmentOnly\n";
print STDERR "allResiduesInPdb: $allResiduesInPdb\n\n";
}

return $endPos;
}

sub reverseMyString
{
my ( $string ) = @_;
my @allChars = split //, $string;

return join("", reverse @allChars );
}

sub get_loopPositions
{
my ( $adjustedAlignmentSequence, $allResiduesInPdb ) = @_;
my @loopPositions = ();

my $dashRemovedAdjustedSequence = $adjustedAlignmentSequence;
$dashRemovedAdjustedSequence =~ s/-//g;

if ( $dashRemovedAdjustedSequence eq $allResiduesInPdb )
{
return @loopPositions;
}

my @allAminoAcidsInTheAlignmentSequence = split //, $adjustedAlignmentSequence;
my @allAminoAcidsInThePdbFile = split //, $allResiduesInPdb;
my ( $alignmentResiduePos, $pdbFileResiduePos ) = ( 0, 0 );

while ( $alignmentResiduePos <= $#allAminoAcidsInTheAlignmentSequence ) # '$#' means the zero-based last entry number
{
$alignmentResiduePos++;
next if $allAminoAcidsInTheAlignmentSequence[$alignmentResiduePos-1] eq "-";
$pdbFileResiduePos++;

if ( $allAminoAcidsInTheAlignmentSequence[$alignmentResiduePos-1] ne
$allAminoAcidsInThePdbFile[$pdbFileResiduePos - 1] )
{
my $frontRemovedAlignmentSeq =
substr( $adjustedAlignmentSequence,
$alignmentResiduePos - 1,
length($adjustedAlignmentSequence) - $alignmentResiduePos + 1 );
my $frontRemovedPdbResidueSeq =
substr( $allResiduesInPdb,
$pdbFileResiduePos - 1,
length($allResiduesInPdb) - $pdbFileResiduePos + 1 );

my $loopLength = get_startPos($frontRemovedAlignmentSeq, $frontRemovedPdbResidueSeq) - 1;
my ( $loopStartPos, $loopEndPos ) = ( $alignmentResiduePos, $alignmentResiduePos + $loopLength - 1);

if ( $loopLength ne "" && $loopStartPos <= $loopEndPos )
{
push @loopPositions, "$loopStartPos-$loopEndPos";
print "\tloop pisition found: $loopStartPos-$loopEndPos, length=$loopLength\n";

$alignmentResiduePos += $loopLength - 1;
$pdbFileResiduePos--;
}
}
}

return @loopPositions;
}

sub formatFastaSequence
{
my ($seq, $charactersPerLine) = @_;
$seq =~ s/(.{$charactersPerLine})/$1\n/g;
return $seq;
}

sub makeAlignmentFile
{
my ( $alignmentFile, $alignmentFileContent ) = @_;

open FO, ">$alignmentFile" or die "File Write Error, $alignmentFile: $!";
print FO $alignmentFileContent;
close FO;
}

sub makeTopScriptFile
{
my ( $topScriptFile, $alignmentFile, $targetProteinCode, $pdbCode2AtomFile_ref ) = @_;

open FO, ">$topScriptFile" or die "File Write Error, $topScriptFile: $!";
print FO printTopScriptInvariableContentFront();

# print variable contents of the top script file
print FO "SET ALNFILE = '$alignmentFile' # alignment filename\n";
print FO "SET KNOWNS = ";
foreach my $pdbCode ( sort keys %$pdbCode2AtomFile_ref )
{
print FO "'$pdbCode' ";
}
print FO " # codes of the templates\n";
print FO "SET SEQUENCE = '$targetProteinCode' # code of the target\n";

print FO printTopScriptInvariableContentEnd();
close FO;
}

sub printTopScriptInvariableContentFront
{
return "".
"# Homology modelling by the MODELLER TOP routine 'model'.\n".
"\n".
"INCLUDE # Include the predefined TOP routines\n".
"\n".
"SET OUTPUT_CONTROL = 1 1 1 1 1 # uncomment to produce a large log file\n".
"\n";
}

sub printTopScriptInvariableContentEnd
{
return "".
"SET ATOM_FILES_DIRECTORY = './' # directories for input atom files\n".
"SET STARTING_MODEL= 1 # index of the first model\n".
"SET ENDING_MODEL = 1 # index of the last model\n".
" # (determines how many models to calculate)\n".
"CALL ROUTINE = 'model' # do homology modelling\n".
"\n";
}

</fi></fi></pre>