Perl script for ATM

#!/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";
}