calling bioperl functions from matlab -凯发k8网页登录

this example shows the interoperability between matlab® and bioperl - passing arguments from matlab to perl scripts and pulling blast search data back to matlab.

note: perl and the bioperl modules must be installed to run the perl scripts in this example. since version 1.4, bioperl modules have a warnings.pm dependency requiring at least version 5.6 of perl. if you have difficulty running the perl scripts, make sure your perl5lib environment variable includes the path to your bioperl installation or try running from the bioperl installation directory. see the links at and for current release files and complete installation instructions.

introduction

gleevec™ (sti571 or imatinib mesylate) was the first approved drug to specifically turn off the signal of a known cancer-causing protein. initially approved to treat chronic myelogenous leukemia (cml), it is also effective for treatment of gastrointestinal stromal tumors (gist).

research has identified several gene targets for gleevec including: proto-oncogene tyrosine-protein kinase abl1 (np_009297), proto-oncogene tyrosine-protein kinase kit (np_000213), and platelet-derived growth factor receptor alpha precursor (np_006197).

target_abl1 = 'np_009297';
target_kit = 'np_000213';
target_pdgfra = 'np_006197';

accessing sequence information

you can load the sequence information for these proteins from local genpept text files using genpeptread.

abl1_seq = getfield(genpeptread('abl1_gp.txt'), 'sequence');
kit_seq = getfield(genpeptread('kit_gp.txt'), 'sequence');
pdgfra_seq = getfield(genpeptread('pdgfra_gp.txt'), 'sequence');

alternatively, you can obtain protein information directly from the online genpept database maintained by the national center for biotechnology information (ncbi).

run these commands to download data from ncbi:

% abl1_seq = getgenpept(target_abl1, 'sequenceonly', true);
% kit_seq = getgenpept(target_kit, 'sequenceonly', true);
% pdgfra_seq = getgenpept(target_pdgfra, 'sequenceonly', true);

the matlab whos command gives information about the size of these sequences.

whos abl1_seq
whos kit_seq
whos pdgfra_seq

  name          size              bytes  class    attributes
  abl1_seq      1x1149             2298  char               
  name         size             bytes  class    attributes
  kit_seq      1x976             1952  char               
  name            size              bytes  class    attributes
  pdgfra_seq      1x1089             2178  char

calling perl programs from matlab

from matlab, you can harness existing bioperl modules to run a blast search on these sequences. mw_blast.pl is a perl program based on the remoteblast bioperl module. it reads sequences from fasta files, so start by creating a fasta file for each sequence.

fastawrite('abl1.fa', 'abl1 proto-oncogene tyrosine-protein kinase (np_009297)', abl1_seq);
fastawrite('kit.fa', 'kit proto-oncogene tyrosine-protein kinase (np_000213)', kit_seq);
fastawrite('pdgfra.fa', 'pdgfra alpha precursor (np_006197)', pdgfra_seq);

warning: abl1.fa already exists. the data will be appended to the file. 
warning: kit.fa already exists. the data will be appended to the file. 
warning: pdgfra.fa already exists. the data will be appended to the file.

blast searches can take a long time to return results, and the perl program mw_blast includes a repeating sleep state to await the report. sample results have been included with this example, but if you want to try running the blast search with the three sequences, uncomment the following commands. mw_blast.pl will save the blast results in three files on your disk, abl1.out, kit.out and pdgfra.out. the process can take 15 minutes or more.

% try
%     perl('mw_blast.pl','blastp','pdb','1e-10','abl1.fa','kit.fa','pdgfra.fa');
% catch
%     error(message('bioinfo:bioperldemo:perlerror'))
% end

here is the perl code for mw_blast:

type mw_blast.pl

#!/usr/bin/perl -w
use bio::tools::run::remoteblast;
use strict;
use 5.006;
# a sample blast program  based on the remoteblast.pm bioperl module. takes
# parameters for the blast search program, the database, and the expectation
# or e-value (defaults: blastp, pdb, 1e-10), followed by a list of fasta files
# containing sequences to search.
# 凯发官网入口首页 copyright 2003-2004 the mathworks, inc.
# retrieve arguments and set parameters
my $prog = shift @argv;
my $db   = shift @argv;
my $e_val= shift @argv;
my @params = ('-prog' => $prog,
	      '-data' => $db,
	      '-expect' => $e_val,
	      '-readmethod' => 'searchio' );
# create a remote blast factory	          
my $factory = bio::tools::run::remoteblast->new(@params);
# change a parameter in remoteblast
$bio::tools::run::remoteblast::header{'entrez_query'} = 'homo sapiens [orgn]';
# remove a parameter from remoteblast
delete $bio::tools::run::remoteblast::header{'filter'};
# submit each file
while ( defined($argv[0])) {
    my $fa_file = shift @argv;
    my $str = bio::seqio->new(-file=>$fa_file, '-format' => 'fasta' );    
    my $r = $factory->submit_blast($fa_file);
    # wait for the reply and save the output file
    while ( my @rids = $factory->each_rid ) {
	foreach my $rid ( @rids ) {
	    my $rc = $factory->retrieve_blast($rid);
	    if( !ref($rc) ) {
    		if( $rc < 0 ) {
        	    $factory->remove_rid($rid);
            }
            sleep 5;
	    } else {
            my $result = $rc->next_result();
            my $filename = $result->query_name()."\.out";
            $factory->save_output($filename);
            $factory->remove_rid($rid);
            }
        }
    }
}

the next step is to parse the output reports and find scores >= 100. you can then identify hits found by more than one protein for further research, possibly identifying new targets for drug therapy.

try
   protein_list = perl('mw_parse.pl', which('abl1.out'), which('kit.out'), which('pdgfra.out'))
catch
    error(message('bioinfo:bioperldemo:perlerror'))
end

protein_list =
    '
     /home/data/abl1.out
     1opl, 2584, 0.0, chain a, structural basis for the auto-inhibition of c-abl...
     1fmk, 923, 1e-100, crystal structure of human tyrosine-protein kinase c-src p...
     1qcf, 919, 1e-100, chain a, crystal structure of hck in complex with a src fa...
     1ksw, 916, 1e-100, chain a, structure of human c-src tyrosine kinase (thr338g...
     1ad5, 883, 6e-96, chain a, src family kinase hck-amp-pnp complex pdb|1ad5|b ...
     2abl, 866, 5e-94, sh3-sh2 domain fragment of human bcr-abl tyrosine kinase
     3lck, 666, 9e-71, the kinase domain of human lymphocyte kinase (lck), activa...
     1qpe, 666, 9e-71, chain a, structural analysis of the lymphocyte-specific ki...
     1qpd, 656, 1e-69, chain a, structural analysis of the lymphocyte-specific ki...
     1k2p, 620, 2e-65, chain a, crystal structure of bruton's tyrosine kinase dom...
     1byg, 592, 3e-62, chain a, kinase domain of human c-terminal src kinase (csk...
     1m7n, 561, 1e-58, chain a, crystal structure of unactivated apo insulin-like...
     1jqh, 560, 2e-58, chain a, igf-1 receptor kinase domain pdb|1jqh|b chain b, ...
     1p4o, 560, 2e-58, chain a, structure of apo unactivated igf-1r kinase domain...
     1k3a, 553, 1e-57, chain a, structure of the insulin-like growth factor 1 rec...
     1gjo, 550, 2e-57, chain a, the fgfr2 tyrosine kinase domain
     1fvr, 540, 3e-56, chain a, tie2 kinase domain pdb|1fvr|b chain b, tie2 kinas...
     1ab2, 528, 9e-55, proto-oncogene tyrosine kinase (e.c.2.7.1.112) (src homolo...
     1irk, 525, 2e-54, insulin receptor (tyrosine kinase domain) mutant with cys ...
     1i44, 523, 3e-54, chain a, crystallographic studies of an activation loop mu...
     1ir3, 522, 4e-54, chain a, phosphorylated insulin receptor tyrosine kinase i...
     1fgk, 522, 4e-54, chain a, crystal structure of the tyrosine kinase domain o...
     1p14, 521, 6e-54, chain a, crystal structure of a catalytic-loop mutant of t...
     1m14, 496, 4e-51, chain a, tyrosine kinase domain from epidermal growth fact...
     1pkg, 496, 4e-51, chain a, structure of a c-kit kinase product complex pdb|1...
     1vr2, 463, 3e-47, chain a, human vascular endothelial growth factor receptor...
     1ju5, 330, 8e-32, chain c, ternary complex of an crk sh2 domain, crk-derived...
     1bbz, 317, 3e-30, chain a, crystal structure of the abl-sh3 domain complexed...
     1awo, 303, 1e-28, the solution nmr structure of abl sh3 and its relationship...
     1bbz, 303, 1e-28, chain e, crystal structure of the abl-sh3 domain complexed...
     1g83, 287, 8e-27, chain a, crystal structure of fyn sh3-sh2 pdb|1g83|b chain...
     1lck, 270, 7e-25, chain a, sh3-sh2 domain fragment of human p56-lck tyrosine...
     1muo, 233, 1e-20, chain a, crystal structure of aurora-2, an oncogenic serin...
     1gri, 232, 2e-20, chain a, grb2 pdb|1gri|b chain b, grb2
     1a9u, 220, 4e-19, the complex structure of the map kinase p38sb203580 pdb|1b...
     1bmk, 213, 3e-18, chain a, the complex structure of the map kinase p38sb2186...
     1ian, 209, 8e-18, human p38 map kinase inhibitor complex
     1gz8, 208, 1e-17, chain a, human cyclin dependent kinase 2 complexed with th...
     1ove, 208, 1e-17, chain a, the structure of p38 alpha in complex with a dihy...
     1oit, 207, 1e-17, chain a, imidazopyridines: a potent and selective class of...
     1b38, 206, 2e-17, chain a, human cyclin-dependent kinase 2 pdb|1b39|a chain ...
     1ogu, 206, 2e-17, chain a, structure of human thr160-phospho cdk2cyclin a co...
     1e9h, 206, 2e-17, chain a, thr 160 phosphorylated cdk2 - human cyclin a3 com...
     1jst, 206, 2e-17, chain a, phosphorylated cyclin-dependent kinase-2 bound to...
     1wfc, 206, 2e-17, structure of apo, unphosphorylated, p38 mitogen activated ...
     1qmz, 206, 2e-17, chain a, phosphorylated cdk2-cyclyin a-substrate peptide c...
     1di8, 206, 2e-17, chain a, the structure of cyclin-dependent kinase 2 (cdk2)...
     1h1p, 206, 2e-17, chain a, structure of human thr160-phospho cdk2cyclin a co...
     1di9, 205, 2e-17, chain a, the structure of p38 mitogen-activated protein ki...
     1h4l, 202, 5e-17, chain a, structure and regulation of the cdk5-p25(nck5a) c...
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1h01|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1oir|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1gii|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1csy|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1f3m|c)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1a81|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1h1w|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1b6c|b)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1ig1|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1jkk|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1jow|b)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1bi8|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1o6k|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1gzk|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1gzn|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1o6l|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1bhf|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1lcj|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1pme|)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1cm8|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1a1a|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|3hck|)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1aot|f)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1pmq|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1lkk|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1jnk|)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1shd|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1lkl|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1bm2|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1bmb|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1cwd|l)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1bhh|b)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1ia8|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1fbz|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     --------------------- warning ---------------------
     msg: no hsps for this minimal hit (pdb|1ijr|a)
     if using ncbi blast, check bits() instead
     ---------------------------------------------------
     
     /home/data/kit.out
     1pkg, 974, 1e-106, chain a, structure of a c-kit kinase product complex pdb|1...
     1vr2, 805, 6e-87, chain a, human vascular endothelial growth factor receptor...
     1gjo, 730, 3e-78, chain a, the fgfr2 tyrosine kinase domain
     1fgk, 700, 8e-75, chain a, crystal structure of the tyrosine kinase domain o...
     1opl, 410, 4e-41, chain a, structural basis for the auto-inhibition of c-abl...
     1fvr, 405, 1e-40, chain a, tie2 kinase domain pdb|1fvr|b chain b, tie2 kinas...
     1m7n, 383, 5e-38, chain a, crystal structure of unactivated apo insulin-like...
     1p4o, 383, 5e-38, chain a, structure of apo unactivated igf-1r kinase domain...
     1jqh, 381, 8e-38, chain a, igf-1 receptor kinase domain pdb|1jqh|b chain b, ...
     1qcf, 377, 2e-37, chain a, crystal structure of hck in complex with a src fa...
     1k3a, 371, 1e-36, chain a, structure of the insulin-like growth factor 1 rec...
     1i44, 368, 3e-36, chain a, crystallographic studies of an activation loop mu...
     1irk, 367, 3e-36, insulin receptor (tyrosine kinase domain) mutant with cys ...
     1p14, 361, 2e-35, chain a, crystal structure of a catalytic-loop mutant of t...
     1ir3, 361, 2e-35, chain a, phosphorylated insulin receptor tyrosine kinase i...
     3lck, 354, 1e-34, the kinase domain of human lymphocyte kinase (lck), activa...
     1qpe, 354, 1e-34, chain a, structural analysis of the lymphocyte-specific ki...
     1qpd, 354, 1e-34, chain a, structural analysis of the lymphocyte-specific ki...
     1ad5, 348, 6e-34, chain a, src family kinase hck-amp-pnp complex pdb|1ad5|b ...
     1ksw, 344, 2e-33, chain a, structure of human c-src tyrosine kinase (thr338g...
     1fmk, 344, 2e-33, crystal structure of human tyrosine-protein kinase c-src p...
     1byg, 342, 3e-33, chain a, kinase domain of human c-terminal src kinase (csk...
     1m14, 335, 2e-32, chain a, tyrosine kinase domain from epidermal growth fact...
     1k2p, 294, 1e-27, chain a, crystal structure of bruton's tyrosine kinase dom...
     1h4l, 167, 5e-13, chain a, structure and regulation of the cdk5-p25(nck5a) c...
     1pme, 158, 6e-12, structure of penta mutant human erk2 map kinase complexed ...
     1f3m, 156, 1e-11, chain c, crystal structure of human serinethreonine kinase...
     
     /home/data/pdgfra.out
     1pkg, 625, 5e-66, chain a, structure of a c-kit kinase product complex pdb|1...
     1vr2, 550, 2e-57, chain a, human vascular endothelial growth factor receptor...
     1fgi, 500, 1e-51, chain a, crystal structure of the tyrosine kinase domain o...
     1gjo, 492, 1e-50, chain a, the fgfr2 tyrosine kinase domain
     1fvr, 419, 4e-42, chain a, tie2 kinase domain pdb|1fvr|b chain b, tie2 kinas...
     1qcf, 380, 1e-37, chain a, crystal structure of hck in complex with a src fa...
     1qpe, 364, 9e-36, chain a, structural analysis of the lymphocyte-specific ki...
     1qpd, 364, 9e-36, chain a, structural analysis of the lymphocyte-specific ki...
     3lck, 360, 2e-35, the kinase domain of human lymphocyte kinase (lck), activa...
     1opl, 358, 4e-35, chain a, structural basis for the auto-inhibition of c-abl...
     1fmk, 354, 1e-34, crystal structure of human tyrosine-protein kinase c-src p...
     1ksw, 353, 2e-34, chain a, structure of human c-src tyrosine kinase (thr338g...
     1ad5, 353, 2e-34, chain a, src family kinase hck-amp-pnp complex pdb|1ad5|b ...
     1byg, 352, 2e-34, chain a, kinase domain of human c-terminal src kinase (csk...
     1i44, 351, 3e-34, chain a, crystallographic studies of an activation loop mu...
     1irk, 350, 4e-34, insulin receptor (tyrosine kinase domain) mutant with cys ...
     1m7n, 349, 5e-34, chain a, crystal structure of unactivated apo insulin-like...
     1jqh, 349, 5e-34, chain a, igf-1 receptor kinase domain pdb|1jqh|b chain b, ...
     1p4o, 349, 5e-34, chain a, structure of apo unactivated igf-1r kinase domain...
     1p14, 344, 2e-33, chain a, crystal structure of a catalytic-loop mutant of t...
     1ir3, 343, 2e-33, chain a, phosphorylated insulin receptor tyrosine kinase i...
     1k3a, 338, 9e-33, chain a, structure of the insulin-like growth factor 1 rec...
     1m14, 332, 4e-32, chain a, tyrosine kinase domain from epidermal growth fact...
     1k2p, 315, 4e-30, chain a, crystal structure of bruton's tyrosine kinase dom...
     1pme, 167, 6e-13, structure of penta mutant human erk2 map kinase complexed ...
     1jow, 155, 1e-11, chain b, crystal structure of a complex of human cdk6 and ...
     1bi8, 155, 1e-11, chain a, mechanism of g1 cyclin dependent kinase inhibitio...
     1f3m, 150, 6e-11, chain c, crystal structure of human serinethreonine kinase...
     '

this is the code for mw_parse:

type mw_parse.pl

#!/usr/bin/perl
use bio::searchio;
use strict;
use 5.006;
# a sample blast parsing program based on the searchio.pm bioperl module. takes
# a list of blast report files and prints a list of the top hits from each
# report based on an arbitrary minimum score.
# 凯发官网入口首页 copyright 2003-2012 the mathworks, inc.
# set a cutoff value for the raw score.
my $min_score = 100;
# take each report name and print information about the top hits.
my $seq_count = 0;
while ( defined($argv[0])) {
    my $breport = shift @argv;
    print "\n$breport\n";
    my $in = new bio::searchio(-format => 'blast', 
                               -file   => $breport);
    my $num_hit = 0;
    my $short_desc;
    while ( my $result = $in->next_result) {
	while ( my $curr_hit = $result->next_hit ) {
	    if ( $curr_hit->raw_score >= $min_score ) {
		if (length($curr_hit->description) >= 60) {
		    $short_desc = substr($curr_hit->description, 0, 58)."...";
		} else {
		    $short_desc = $curr_hit->description;
		}
		print $curr_hit->accession, ", ",
		      $curr_hit->raw_score, ", ",
		      $curr_hit->significance, ", ",
		      $short_desc, "\n";
	    }
	    $num_hit  ;
	}
    }
    $seq_count  ;
}

calling matlab functions within perl programs

if you are running on windows®, it is also possible to call matlab functions from perl. you can launch matlab in an automation server mode by using the /automation switch in the matlab startup command (e.g. d:\applications\matlab7x\bin\matlab.exe /automation).

here's a script to illustrate the process of launching an automation server, calling matlab functions and passing variables between perl and matlab.

type matlab_from_perl.pl

#!/usr/bin/perl -w
use win32::ole;
use win32::ole::variant;
# simple perl script to execute commands in matlab.
# note the name win32::ole is misleading and this actually uses com!
#
# use existing instance if matlab is already running.
eval {$matlabapp = win32::ole->getactiveobject('matlab.application')};
die "matlab not installed" if $@;
unless (defined $matlabapp) {
   $matlabapp = win32::ole->new('matlab.application')
      or die "oops, cannot start matlab";
}
# examples of executing matlab commands - these functions execute in
# matlab and return the status.
@exe_commands = ("irk = pdbread('pdb1irk.ent');",
             "lck = pdbread('pdb3lck.ent');",
             "seqdisp(irk)",
             "seqdisp(lck)",
             "[score, alignment] = swalign(irk, lck,'showscore',1);");
# send the commands to matlab
foreach $exe_command (@exe_commands)
{  $status = &send_to_matlab('execute', $exe_command);
   print "matlab status = ", $status, "\n";
}
sub send_to_matlab
{  my ($call, @command) = @_;
   my $status = 0;
   print "\n>> $call( @command )\n";
   $result = $matlabapp->invoke($call, @command);
   if (defined($result))
   {   unless ($result =~ s/^.\?{3}/error:/)
       {  print "$result\n" unless ($result eq "");
       }
       else
       {  print "$result\n";
          $status = -1;
       }
   }
   return $status;
}
# examples of passing variables between matlab and perl.
#
# matlab supoprts passing character arrays directly with the following syntax:
#
# putchararray([in] bstr name, [in] bstr workspace, [in] bstr string);
# getchararray([in] bstr name, [in] bstr workspace, [out] bstr string);
&send_to_matlab('putchararray', 'centraldogma', 'base', 'dna->rna->protein.');
&send_to_matlab('getchararray', 'centraldogma', 'base');
# numeric arrays can be passed by reference in a safearray using the
# putfullmatrix and getfullmatrix functions.
#
# putfullmatrix([in] bstr name, [in] bstr workspace, [in] bstr data);
# getfullmatrix([in] bstr varname, [in] bstr workspace, [out] bstr retdata);
$mreal = variant(vt_array|vt_r8, 4, 4);
$mimag = variant(vt_array|vt_r8, 4, 4);
$mreal->put([[0,0,0,0], [0,0,0,0], [0,0,0,0], [0,0,0,0]]);
print "\n>> putfullmatrix( 'magicarray', 'base', ",'$mreal, $mimag'," )\n";
$matlabapp->putfullmatrix('magicarray', 'base', $mreal, $mimag);
$matlabapp->execute('magicarray = magic(4)');
$m2real = variant(vt_array|vt_r8|vt_byref,4,4);
$m2imag = variant(vt_array|vt_r8|vt_byref,4,4);
print "\n>> getfullmatrix( 'magicarray', 'base', ",'$m2real, $m2imag'," )\n";
$matlabapp->getfullmatrix('magicarray', 'base', $m2real, $m2imag);
for ($i = 0; $i < 4; $i  ) {
	printf "= = = =\n", $m2real->get($i,0), $m2real->get($i,1),
                                $m2real->get($i,2), $m2real->get($i,3);
}
# additionally, you can use variants to send scalar variables by reference
# to matlab for all data types except sparse arrays and function handles through
# putworkspacedata:
# putworkspacedata([in] bstr name, [in] bstr workspace, [in] bstr data);
#
# results are passed back to perl directly with getvariable:
# hresult = getvariable([in] bstr name, [in] bstr workspace); 
   
# create and initialize a date variant.
$dnadate = variant->new(vt_date|vt_byref, 'feb 28, 1953');
&send_to_matlab('putworkspacedata', 'dnadate', 'base', $dnadate);
&send_to_matlab('execute', 'dnadate');
# create and initialize a new string variant.
$aminostring = variant->new(vt_bstr|vt_byref, 'matlap');
&send_to_matlab('putworkspacedata', 'aminoacids', 'base', $aminostring);
# change the value in matlab
&send_to_matlab('execute', "aminoacids = 'arndcqeghilkmfpstwyv';");
# bring the new value back
$aa = $matlabapp->getvariable('aminoacids', 'base');
printf "amino acid codes: %s\n", $aa;
undef $matlabapp; # close matlab if we opened it

protein analysis tools in bioinformatics toolbox™

matlab offers additional tools for protein analysis and further research with these proteins. for example, to access the sequences and run a full smith-waterman alignment on the tyrosine kinase domain of the human insulin receptor (pdb 1irk) and the kinase domain of the human lymphocyte kinase (pdb 3lck), load the sequence data:

irk = pdbread('pdb1irk.ent');
lck = pdbread('pdb3lck.ent');
% run these commands to bring the data from the internet:
% irk = getpdb('1irk');
% lck = getpdb('3lck');

now perform a local alignment with the smith-waterman algorithm. matlab uses blosum 50 as the default scoring matrix for aa strings with a gap penalty of 8. of course, you can change any of these parameters.

[score, alignment] = swalign(irk, lck, 'showscore', true);

matlab and the bioinformatics toolbox™ offer additional tools for investigating nucleotide and amino acid sequences. for example, pdbdistplot displays the distances between atoms and amino acids in a pdb structure, while ramachandran generates a plot of the torsion angle phi and the torsion angle psi of the protein sequence. the toolbox function proteinplot provides a graphical user interface (gui) to easily import sequences and plot various properties such as hydrophobicity.