Bioperl Scripts

<br />#!/usr/bin/perl<br /># Extract sequence chains from PDB file<br /><br /><br /># Read in PDB file: Warning - some files are very large!<br />print "Enter Pdb file name,\n"; <br />my $filename =<STDIN>;<br /><br /> chomp $filename;<br /> unless(-e $filename) {<br /> print "$filename does not exist!\n";<br /> exit;<br /> }<br /><br /><br />my @file = get_file_data($filename);<br /><br /># Parse the record types of the PDB file<br />my %recordtypes = parsePDBrecordtypes(@file);<br /><br /># Extract the amino acid sequences of all chains in the protein<br />my @chains = extractSEQRES( $recordtypes{'SEQRES'} );<br /><br /># Translate the 3-character codes to 1-character codes, and print<br />foreach my $chain (@chains) {<br /> print "$chain\n";<br /> print iub3to1($chain), "\n";<br />}<br />print "Press <Enter> to exit";<br /><STDIN>;<br />exit;<br /><br /><br /><br /># parsePDBrecordtypes<br />#<br />#-given an array of a PDB file, return a hash with<br /># keys = record type names<br /># values = scalar containing lines for that record type <br /><br />sub parsePDBrecordtypes {<br /><br /> my @file = @_;<br /><br /> use strict;<br /> use warnings;<br /> <br /> my %recordtypes = ( );<br /> <br /> foreach my $line (@file) {<br /> <br /> # Get the record type name which begins at the<br /> # start of the line and ends at the first space<br /><br /> # The pattern (\S+) is returned and saved in $recordtype<br /> my($recordtype) = ($line =~ /^(\S+)/);<br /> <br /> # .= fails if a key is undefined, so we have to<br /> # test for definition and use either .= or = depending<br /> if(defined $recordtypes{$recordtype} ) {<br /> $recordtypes{$recordtype} .= $line;<br /> }else{<br /> $recordtypes{$recordtype} = $line;<br /> }<br /> }<br /> <br /> return %recordtypes;<br />}<br /><br /># extractSEQRES<br />#<br />#-given an scalar containing SEQRES lines,<br /># return an array containing the chains of the sequence<br /><br />sub extractSEQRES {<br /><br /> use strict;<br /> use warnings;<br /><br /> my($seqres) = @_;<br /><br /> my $lastchain = '';<br /> my $sequence = '';<br /> my @results = ( );<br /> # make array of lines<br /><br /> my @record = split ( /\n/, $seqres);<br /> <br /> foreach my $line (@record) {<br /> # Chain is in column 12, residues start in column 20<br /> my ($thischain) = substr($line, 11, 1);<br /> my($residues) = substr($line, 19, 52); # add space at end<br /> <br /> # Check if a new chain, or continuation of previous chain<br /> if("$lastchain" eq "") {<br /> $sequence = $residues;<br /> }elsif("$thischain" eq "$lastchain") {<br /> $sequence .= $residues;<br /> <br /> # Finish gathering previous chain (unless first record)<br /> }elsif ( $sequence ) {<br /> push(@results, $sequence);<br /> $sequence = $residues;<br /> }<br /> $lastchain = $thischain;<br /> }<br /><br /> # save last chain<br /> push(@results, $sequence);<br /> <br /> return @results;<br />}<br /><br /># iub3to1<br />#<br />#-change string of 3-character IUB amino acid codes (whitespace separated)<br /># into a string of 1-character amino acid codes<br /><br />sub iub3to1 {<br /><br /> my($input) = @_;<br /> <br /> my %three2one = (<br /> 'ALA' => 'A',<br /> 'VAL' => 'V',<br /> 'LEU' => 'L',<br /> 'ILE' => 'I',<br /> 'PRO' => 'P',<br /> 'TRP' => 'W',<br /> 'PHE' => 'F',<br /> 'MET' => 'M',<br /> 'GLY' => 'G',<br /> 'SER' => 'S',<br /> 'THR' => 'T',<br /> 'TYR' => 'Y',<br /> 'CYS' => 'C',<br /> 'ASN' => 'N',<br /> 'GLN' => 'Q',<br /> 'LYS' => 'K',<br /> 'ARG' => 'R',<br /> 'HIS' => 'H',<br /> 'ASP' => 'D',<br /> 'GLU' => 'E',<br /> );<br /><br /> # clean up the input<br /> $input =~ s/\n/ /g;<br /><br /> my $seq = '';<br /> <br /> # This use of split separates on any contiguous whitespace<br /> my @code3 = split(' ', $input);<br /><br /> foreach my $code (@code3) {<br /> # A little error checking<br /> if(not defined $three2one{$code}) {<br /> print "Code $code not defined\n";<br /> next;<br /> }<br /> $seq .= $three2one{$code};<br /> }<br /><br />sub get_file_data {<br /><br /> my($filename) = @_;<br /><br /> use strict;<br /> use warnings;<br /><br /> # Initialize variables<br /> my @filedata = ( );<br /><br /> unless( open(GET_FILE_DATA, $filename) ) {<br /> print STDERR "Cannot open file \"$filename\"\n\n";<br /> exit;<br /> }<br /><br /> @filedata = <GET_FILE_DATA>;<br /><br /> close GET_FILE_DATA;<br /><br /> return @filedata;<br />}<br /><br /><br /><br /> return $seq;<br />}<br /><br /><br /><br />

Links

Extract sequence chains from PDB file