<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
-->
<title></title>
<meta name="generator" content="MATLAB">
<meta name="title" content="BIOPERLDEMO: Example of calling Bioperl functions in MATLAB">
<meta name="description" content="This demonstration illustrates interoperability between MATLAB and Bioperl--passing arguments from MATLAB to perl scripts and pulling BLAST search data back to MATLAB."><style>
body {
background-color: white;
}
h1 {
color: #990000;
font-size: x-large;
}
h2 {
color: #990000;
font-size: medium;
}
p.footer {
text-align: right;
font-size: xx-small;
font-weight: lighter;
font-style: italic;
color: gray;
}
pre.codeinput {
margin-left: 30px;
}
span.keyword {color: blue}
span.comment {color: green}
span.string {color: #B20000}
span.untermstring {color: purple}
span.syscmd {color: orange}
pre.showbuttons {
margin-left: 30px;
border: solid black 2px;
padding: 4px;
background: #EBEFF3;
}
pre.codeoutput {
color: gray;
font-style: italic;
}
pre.error {
color: red;
}
</style></head>
<body>
<h2>Contents</h2>
<ul>
<li><a href="#1">BIOPERLDEMO: Example of calling Bioperl functions in MATLAB</a></li>
<li><a href="#2">Introduction</a></li>
<li><a href="#4">Accessing On-line Sequence Information</a></li>
<li><a href="#7">Calling Perl Programs from MATLAB</a></li>
<li><a href="#12">Protein Analysis Tools in the Bioinformatics Toolbox</a></li>
</ul>
<h2>BIOPERLDEMO: Example of calling Bioperl functions in MATLAB<a name="1"></a></h2>
<p>This demonstration illustrates interoperability between MATLAB and Bioperl--passing arguments from MATLAB to perl scripts
and pulling BLAST search data back to MATLAB.
</p><pre class="codeinput"><span class="comment">% PLEASE NOTE: Perl and the Bioperl modules must have been previously</span>
<span class="comment">% installed to run the perl scripts in this demonstration. Depending on</span>
<span class="comment">% your system and network configurations, this may not be a trivial</span>
<span class="comment">% installation.</span>
</pre><h2>Introduction<a name="2"></a></h2>
<p>STI571 (imatinib mesylate) is 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 has also been found effective for treatment of gastrointestinal
stromal tumor (GIST).
</p><pre class="codeinput">web(<span class="string">'http://www.cancer.gov/clinicaltrials/digestpage/gleevec'</span>)
</pre><p>Research has identified three gene targets for STI571: Proto-oncogene tyrosine-protein kinase ABL1 (P00519), Mast/stem cell
growth factor receptor precursor (proto-oncogene tyrosine-protein kinase Kit) (P10721), and Platelet-derived growth factor
receptor alpha precursor (NP_006197).
</p><pre class="codeinput">target_ABL1 = <span class="string">'P00519'</span>;
target_Kit = <span class="string">'P10721'</span>;
target_PDGF = <span class="string">'NP_006197'</span>;
</pre><h2>Accessing On-line Sequence Information<a name="4"></a></h2>
<p>Sequence information for these proteins can be obtained directly from the on-line GenPept database maintained by the National
Center for Biotechnology Information (NCBI).
</p><pre class="codeinput">ABL1_seq = getgenpept(target_ABL1, <span class="string">'SequenceOnly'</span>, true);
Kit_seq = getgenpept(target_Kit, <span class="string">'SequenceOnly'</span>, true);
PDGF_seq = getgenpept(target_PDGF, <span class="string">'SequenceOnly'</span>, true);
</pre><p>Alternatively, you can load the protein information from a GenPept text file using <b>genpeptread</b>.
</p><pre class="codeinput"><span class="comment">% Uncomment the following lines to pull the information from local files:</span>
<span class="comment">% ABL1_seq = getfield(genpeptread('ABL1_seq'), 'Sequence');</span>
<span class="comment">% Kit_seq = getfield(genpeptread('Kit_seq'), 'Sequence');</span>
<span class="comment">% PDGF_seq = getfield(genpeptread('PDGF_seq'), 'Sequence');</span>
</pre><p>The MATLAB <b>whos</b> command gives information about the size of these sequences.
</p><pre class="codeinput">whos ABL1_seq
whos Kit_seq
whos PDGF_seq
</pre><pre class="codeoutput"> Name Size Bytes Class
ABL1_seq 1x1130 2260 char array
Grand total is 1130 elements using 2260 bytes
Name Size Bytes Class
Kit_seq 1x976 1952 char array
Grand total is 976 elements using 1952 bytes
Name Size Bytes Class
PDGF_seq 1x1089 2178 char array
Grand total is 1089 elements using 2178 bytes
</pre><h2>Calling Perl Programs from MATLAB<a name="7"></a></h2>
<p>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.
</p><pre class="codeinput">fastawrite(<span class="string">'ABL1.fa'</span>, <span class="string">'ABL1 Tyrosine-protein kinase (P00519)'</span>, ABL1_seq);
fastawrite(<span class="string">'Kit.fa'</span>, <span class="string">'Kit Tyrosine-protein kinase (P10721)'</span>, Kit_seq);
fastawrite(<span class="string">'PDGF.fa'</span>, <span class="string">'PDGF alpha precursor (NP_006197)'</span>, PDGF_seq);
</pre><p>The actual BLAST searches can take a long time to return results, and the perl prgram MW_BLAST includes a repeating sleep
state to await the report. Sample results have been included with this demo, but if you have an internet connection and want
to try running the BLAST with the three sequences, uncomment this line (this may take 15 minutes or more):
</p><pre class="codeinput"><span class="comment">% perl('MW_BLAST.pl', 'blastp', 'pdb', '1e-10', 'ABL1.fa', 'Kit.fa', 'PDGF.fa');</span>
</pre><p>Here is the perl code for MW_BLAST:</p><pre class="codeinput">type MW_BLAST.pl
</pre><pre class="codeoutput">
#!/usr/bin/perl -w
use Bio::Tools::Run::RemoteBlast;
use strict;
# 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.
# 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' );
# Remote BLAST factory creation
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
# Changing a paramter in RemoteBlast
$Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'Homo sapiens [ORGN]';
# Removing a parameter from RemoteBlast
delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};
# Take each file name and submit it
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);
}
}
}
}
</pre><p>The next step is to parse the output reports, and filter for scores >= 100. You can then identify hits found by more than
one protein for further research--possibly identifying new targets for drug therapy.
</p><pre class="codeinput">protein_list = perl(<span class="string">'MW_parse.pl'</span>, <span class="string">'ABL1.out'</span>, <span class="string">'Kit.out'</span>, <span class="string">'PDGF.out'</span>)
</pre><pre class="codeoutput">
protein_list =
ABL1.out
1OPL, 999, 0.0, Chain A, Structural Basis For The Auto-Inhibition Of C-Abl...
1FMK, 360, 1e-100, Crystal Structure Of Human Tyrosine-Protein Kinase C-Src p...
1QCF, 358, 1e-100, Chain A, Crystal Structure Of Hck In Complex With A Src Fa...
1KSW, 357, 1e-100, Chain A, Structure Of Human C-Src Tyrosine Kinase (Thr338g...
1AD5, 344, 6e-96, Chain A, Src Family Kinase Hck-Amp-Pnp Complex pdb|1AD5|B ...
2ABL, 338, 5e-94, Sh3-Sh2 Domain Fragment Of Human Bcr-Abl Tyrosine Kinase
3LCK, 261, 9e-71, The Kinase Domain Of Human Lymphocyte Kinase (Lck), Activa...
1QPE, 261, 9e-71, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
1QPD, 257, 1e-69, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
1K2P, 243, 2e-65, Chain A, Crystal Structure Of Bruton's Tyrosine Kinase Dom...
1BYG, 232, 3e-62, Chain A, Kinase Domain Of Human C-Terminal Src Kinase (Csk...
1M7N, 220, 1e-58, Chain A, Crystal Structure Of Unactivated Apo Insulin-Like...
1JQH, 220, 2e-58, Chain A, Igf-1 Receptor Kinase Domain pdb|1JQH|B Chain B, ...
1P4O, 220, 2e-58, Chain A, Structure Of Apo Unactivated Igf-1r Kinase Domain...
1K3A, 217, 1e-57, Chain A, Structure Of The Insulin-Like Growth Factor 1 Rec...
1GJO, 216, 2e-57, Chain A, The Fgfr2 Tyrosine Kinase Domain
1FVR, 212, 3e-56, Chain A, Tie2 Kinase Domain pdb|1FVR|B Chain B, Tie2 Kinas...
1AB2, 207, 9e-55, Proto-Oncogene Tyrosine Kinase (E.C.2.7.1.112) (Src Homolo...
1IRK, 206, 2e-54, Insulin Receptor (Tyrosine Kinase Domain) Mutant With Cys ...
1I44, 206, 3e-54, Chain A, Crystallographic Studies Of An Activation Loop Mu...
1IR3, 205, 4e-54, Chain A, Phosphorylated Insulin Receptor Tyrosine Kinase I...
1FGK, 205, 4e-54, Chain A, Crystal Structure Of The Tyrosine Kinase Domain O...
1P14, 205, 6e-54, Chain A, Crystal Structure Of A Catalytic-Loop Mutant Of T...
1M14, 195, 4e-51, Chain A, Tyrosine Kinase Domain From Epidermal Growth Fact...
1PKG, 195, 4e-51, Chain A, Structure Of A C-Kit Kinase Product Complex pdb|1...
1VR2, 182, 3e-47, Chain A, Human Vascular Endothelial Growth Factor Receptor...
1JU5, 131, 8e-32, Chain C, Ternary Complex Of An Crk Sh2 Domain, Crk-Derived...
1BBZ, 126, 3e-30, Chain A, Crystal Structure Of The Abl-Sh3 Domain Complexed...
1AWO, 121, 1e-28, The Solution Nmr Structure Of Abl Sh3 And Its Relationship...
1BBZ, 121, 1e-28, Chain E, Crystal Structure Of The Abl-Sh3 Domain Complexed...
1G83, 115, 8e-27, Chain A, Crystal Structure Of Fyn Sh3-Sh2 pdb|1G83|B Chain...
1LCK, 108, 7e-25, Chain A, Sh3-Sh2 Domain Fragment Of Human P56-Lck Tyrosine...
Kit.out
1PKG, 379, 1e-106, Chain A, Structure Of A C-Kit Kinase Product Complex pdb|1...
1VR2, 314, 6e-87, Chain A, Human Vascular Endothelial Growth Factor Receptor...
1GJO, 285, 3e-78, Chain A, The Fgfr2 Tyrosine Kinase Domain
1FGK, 274, 8e-75, Chain A, Crystal Structure Of The Tyrosine Kinase Domain O...
1OPL, 162, 4e-41, Chain A, Structural Basis For The Auto-Inhibition Of C-Abl...
1FVR, 160, 1e-40, Chain A, Tie2 Kinase Domain pdb|1FVR|B Chain B, Tie2 Kinas...
1M7N, 152, 5e-38, Chain A, Crystal Structure Of Unactivated Apo Insulin-Like...
1P4O, 152, 5e-38, Chain A, Structure Of Apo Unactivated Igf-1r Kinase Domain...
1JQH, 151, 8e-38, Chain A, Igf-1 Receptor Kinase Domain pdb|1JQH|B Chain B, ...
1QCF, 149, 2e-37, Chain A, Crystal Structure Of Hck In Complex With A Src Fa...
1K3A, 147, 1e-36, Chain A, Structure Of The Insulin-Like Growth Factor 1 Rec...
1I44, 146, 3e-36, Chain A, Crystallographic Studies Of An Activation Loop Mu...
1IRK, 145, 3e-36, Insulin Receptor (Tyrosine Kinase Domain) Mutant With Cys ...
1P14, 143, 2e-35, Chain A, Crystal Structure Of A Catalytic-Loop Mutant Of T...
1IR3, 143, 2e-35, Chain A, Phosphorylated Insulin Receptor Tyrosine Kinase I...
3LCK, 140, 1e-34, The Kinase Domain Of Human Lymphocyte Kinase (Lck), Activa...
1QPE, 140, 1e-34, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
1QPD, 140, 1e-34, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
1AD5, 138, 6e-34, Chain A, Src Family Kinase Hck-Amp-Pnp Complex pdb|1AD5|B ...
1KSW, 137, 2e-33, Chain A, Structure Of Human C-Src Tyrosine Kinase (Thr338g...
1FMK, 137, 2e-33, Crystal Structure Of Human Tyrosine-Protein Kinase C-Src p...
1BYG, 136, 3e-33, Chain A, Kinase Domain Of Human C-Terminal Src Kinase (Csk...
1M14, 133, 2e-32, Chain A, Tyrosine Kinase Domain From Epidermal Growth Fact...
1K2P, 117, 1e-27, Chain A, Crystal Structure Of Bruton's Tyrosine Kinase Dom...
PDGF.out
1PKG, 245, 5e-66, Chain A, Structure Of A C-Kit Kinase Product Complex pdb|1...
1VR2, 216, 2e-57, Chain A, Human Vascular Endothelial Growth Factor Receptor...
1FGI, 197, 1e-51, Chain A, Crystal Structure Of The Tyrosine Kinase Domain O...
1GJO, 194, 1e-50, Chain A, The Fgfr2 Tyrosine Kinase Domain
1FVR, 166, 4e-42, Chain A, Tie2 Kinase Domain pdb|1FVR|B Chain B, Tie2 Kinas...
1QCF, 150, 1e-37, Chain A, Crystal Structure Of Hck In Complex With A Src Fa...
1QPE, 144, 9e-36, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
1QPD, 144, 9e-36, Chain A, Structural Analysis Of The Lymphocyte-Specific Ki...
3LCK, 143, 2e-35, The Kinase Domain Of Human Lymphocyte Kinase (Lck), Activa...
1OPL, 142, 4e-35, Chain A, Structural Basis For The Auto-Inhibition Of C-Abl...
1FMK, 140, 1e-34, Crystal Structure Of Human Tyrosine-Protein Kinase C-Src p...
1KSW, 140, 2e-34, Chain A, Structure Of Human C-Src Tyrosine Kinase (Thr338g...
1AD5, 140, 2e-34, Chain A, Src Family Kinase Hck-Amp-Pnp Complex pdb|1AD5|B ...
1BYG, 140, 2e-34, Chain A, Kinase Domain Of Human C-Terminal Src Kinase (Csk...
1I44, 139, 3e-34, Chain A, Crystallographic Studies Of An Activation Loop Mu...
1IRK, 139, 4e-34, Insulin Receptor (Tyrosine Kinase Domain) Mutant With Cys ...
1M7N, 139, 5e-34, Chain A, Crystal Structure Of Unactivated Apo Insulin-Like...
1JQH, 139, 5e-34, Chain A, Igf-1 Receptor Kinase Domain pdb|1JQH|B Chain B, ...
1P4O, 139, 5e-34, Chain A, Structure Of Apo Unactivated Igf-1r Kinase Domain...
1P14, 137, 2e-33, Chain A, Crystal Structure Of A Catalytic-Loop Mutant Of T...
1IR3, 136, 2e-33, Chain A, Phosphorylated Insulin Receptor Tyrosine Kinase I...
1K3A, 134, 9e-33, Chain A, Structure Of The Insulin-Like Growth Factor 1 Rec...
1M14, 132, 4e-32, Chain A, Tyrosine Kinase Domain From Epidermal Growth Fact...
1K2P, 125, 4e-30, Chain A, Crystal Structure Of Bruton's Tyrosine Kinase Dom...
</pre><p>This is the code for MW_parse:</p><pre class="codeinput">type MW_parse.pl
</pre><pre class="codeoutput">
#!/usr/bin/perl -w
use Bio::SearchIO;
use strict;
# 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.
# 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++;
}
</pre><h2>Protein Analysis Tools in the Bioinformatics Toolbox<a name="12"></a></h2>
<p>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), try the following:
</p><pre class="codeinput">insulin = getpdb(<span class="string">'1IRK'</span>);
LCK = getpdb(<span class="string">'3LCK'</span>);
</pre><p>PDB often includes a wildcard character to denote something unusual in a sequence, such as a '?' to represent an autophosphorylated
AA. Before you can align the sequences, these need to be converted to valid codes. Use 'X' as a replacement.
</p><pre class="codeinput">clean_insulin = strrep(insulin.Sequence.Sequence,<span class="string">'?'</span>,<span class="string">'X'</span>)
clean_LCK = strrep(LCK.Sequence.Sequence,<span class="string">'?'</span>,<span class="string">'X'</span>)
</pre><pre class="codeoutput">
clean_insulin =
VFPSSVFVPDEWEVSREKITLLRELGQGSFGMVYEGNARDIIKGEAETRVAVKTVNESASLRERIEFLNEASVMKGFTCHHVVRLLGVVS...
KGQPTLVVMELMAHGDLKSYLRSLRPEAENNPGRPPPTLQEMIQMAAEIADGMAYLNAKKFVHRDLAARNCMVAHDFTVKIGDFGMTRDI...
YETDYYRKGGKGLLPVRWMAPESLKDGVFTTSSDMWSFGVVLWEITSLAEQPYQGLSNEQVLKFVMDGGYLDQPDNCPERVTDLMRMCWQ...
FNPKMRPTFLEIVNLLKDDLHPSFPEVSFFHSEENK
clean_LCK =
KPWWEDEWEVPRETLKLVERLGAGQFGEVWMGYYNGHTKVAVKSLKQGSMSPDAFLAEANLMKQLQHQRLVRLYAVVTQEPIYIITEYME...
NGSLVDFLKTPSGIKLTINKLLDMAAQIAEGMAFIEERNYIHRDLRAANILVSDTLSCKIADFGLARLIEDNEXTAREGAKFPIKWTAPE...
AINYGTFTIKSDVWSFGILLTEIVTHGRIPYPGMTNPEVIQNLERGYRMVRPDNCPEELYQLMRLCWKERPEDRPTFDYLRSVLEDFFTAT
</pre><p>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.
</p><pre class="codeinput">[Score, Alignment] = swalign(clean_insulin, clean_LCK)
</pre><pre class="codeoutput">
Score =
213.3333
Alignment =
DEWEVSREKITLLRELGQGSFGMVYEGNARDIIKGEAETRVAVKTVNESASLRERIEFLNEASVMKGFTCHHVVRLLGVVSKGQPTLVVM...
||||| || : |:::|| |:|| |: | : | :|:||||:::: :|: || ||::|| : :::||| :||:: :| :: ...
DEWEVPRETLKLVERLGAGQFGEVWMG-YYN---G--HTKVAVKSLKQ-GSMSPD-AFLAEANLMKQLQHQRLVRLYAVVTQ-EPIYIIT...
ELMAHGDLKSYLRSLRPEAENNPGRPPPTLQEMIQMAAEIADGMAYLNAKKFVHRDLAARNCMVAHDFTVKIGDFGMTRDIYETDY-YRK...
| | :|:| ::|:: | : | |::::::|||:||:|||::: ::::|||| | | :|: :: ||:|||::| | ::: |:...
EYMENGSLVDFLKT--P----S-G-IKLTINKLLDMAAQIAEGMAFIEERNYIHRDLRAANILVSDTLSCKIADFGLARLIEDNEXTARE...
VLKFVMDGGY-LDQPDNCPERVTDLMRMCWQFNPKMRPTFLEIVNLLKDGGKGLLPVRWMAPESLKDGVFTTSSDMWSFGVVLWEITSLA...
|:: :: || : :||||||:: :|||:||: |: |||| : ::|:||:| :|::| |||::: |:|| :||:||||::| ||:: :...
GAK--FPIKWTAPEAINYGTFTIKSDVWSFGILLTEIVTHGRIPYPGMTNPEVIQ-NLERGYRMVRPDNCPEELYQLMRLCWKERPEDRP...
EQPYQGLSNEQ
: || |::| :
TFDYLRSVLED
</pre><p>MATLAB and the Bioinformatics Toolbox offer additional tools for investigation of nucleotide and amino acid sequences. For
example, <b>pdbdistplot</b> displays the distances between atoms and amino acids in a PDB structure, while <b>ramachandran</b> generates a plot of the torsion angle PHI and the torsion angle PSI of the protein sequence. The Toolbox function <b>proteinplot</b> provides a graphical user interface (GUI) to easily import sequences and plot various properties such as hydrophobicity.
</p>
<p class="footer"><br></p>
<!--
##### SOURCE BEGIN #####
%% BIOPERLDEMO: Example of calling Bioperl functions in MATLAB
% This demonstration illustrates interoperability between MATLAB and
% BioperlREPLACE_WITH_DASH_DASHpassing arguments from MATLAB to perl scripts and pulling BLAST
% search data back to MATLAB.
% PLEASE NOTE: Perl and the Bioperl modules must have been previously
% installed to run the perl scripts in this demonstration. Depending on
% your system and network configurations, this may not be a trivial
% installation.
%% Introduction
% STI571 (imatinib mesylate) is 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 has also been found
% effective for treatment of gastrointestinal stromal tumor (GIST).
web('http://www.cancer.gov/clinicaltrials/digestpage/gleevec')
%%
% Research has identified three gene targets for STI571: Proto-oncogene
% tyrosine-protein kinase ABL1 (P00519), Mast/stem cell growth factor
% receptor precursor (proto-oncogene tyrosine-protein kinase Kit) (P10721),
% and Platelet-derived growth factor receptor alpha precursor (NP_006197).
target_ABL1 = 'P00519';
target_Kit = 'P10721';
target_PDGF = 'NP_006197';
%% Accessing On-line Sequence Information
% Sequence information for these proteins can be obtained directly from the
% on-line GenPept database maintained by the National Center for
% Biotechnology Information (NCBI).
ABL1_seq = getgenpept(target_ABL1, 'SequenceOnly', true);
Kit_seq = getgenpept(target_Kit, 'SequenceOnly', true);
PDGF_seq = getgenpept(target_PDGF, 'SequenceOnly', true);
%%
% Alternatively, you can load the protein information from a GenPept text
% file using *genpeptread*.
% Uncomment the following lines to pull the information from local files:
% ABL1_seq = getfield(genpeptread('ABL1_seq'), 'Sequence');
% Kit_seq = getfield(genpeptread('Kit_seq'), 'Sequence');
% PDGF_seq = getfield(genpeptread('PDGF_seq'), 'Sequence');
%%
% The MATLAB *whos* command gives information about the size of these
% sequences.
whos ABL1_seq
whos Kit_seq
whos PDGF_seq
%% 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 Tyrosine-protein kinase (P00519)', ABL1_seq);
fastawrite('Kit.fa', 'Kit Tyrosine-protein kinase (P10721)', Kit_seq);
fastawrite('PDGF.fa', 'PDGF alpha precursor (NP_006197)', PDGF_seq);
%%
% The actual BLAST searches can take a long time to return results, and the
% perl prgram MW_BLAST includes a repeating sleep state to await the
% report. Sample results have been included with this demo, but if you have
% an internet connection and want to try running the BLAST with the three
% sequences, uncomment this line (this may take 15 minutes or more):
% perl('MW_BLAST.pl', 'blastp', 'pdb', '1e-10', 'ABL1.fa', 'Kit.fa', 'PDGF.fa');
%%
% Here is the perl code for MW_BLAST:
type MW_BLAST.pl
%%
% The next step is to parse the output reports, and filter for scores >=
% 100. You can then identify hits found by more than one protein for
% further researchREPLACE_WITH_DASH_DASHpossibly identifying new targets for drug therapy.
protein_list = perl('MW_parse.pl', 'ABL1.out', 'Kit.out', 'PDGF.out')
%%
% This is the code for MW_parse:
type MW_parse.pl
%% Protein Analysis Tools in the 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), try the following:
insulin = getpdb('1IRK');
LCK = getpdb('3LCK');
%%
% PDB often includes a wildcard character to denote something unusual in a
% sequence, such as a '?' to represent an autophosphorylated AA. Before you
% can align the sequences, these need to be converted to valid codes. Use
% 'X' as a replacement.
clean_insulin = strrep(insulin.Sequence.Sequence,'?','X')
clean_LCK = strrep(LCK.Sequence.Sequence,'?','X')
%%
% 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(clean_insulin, clean_LCK)
%%
% MATLAB and the Bioinformatics Toolbox offer additional tools for
% investigation of 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.
##### SOURCE END #####
-->
</body>
</html>