Bioperl: probe matching
Jong Park
jong@ebi.ac.uk
Thu, 03 Jun 1999 17:46:18 +0000
This is a multi-part message in MIME format.
--------------E48DD72FE129715BCCFF0D96
Content-Type: text/plain; charset=EUC-KR
Content-Transfer-Encoding: 7bit
>
Hi
What you need is a simple Dynamic Programming (or graph algorithm)
to produce an alignment(and match score)
I am sure some people in this list can give you some
simple example(or full) code for it.
I have a perl routines which can produce some kind
of match score based on DP. It was just to experiment
on other things, but if you are reasonably expert, it might
be useful. Just type the following:
make_dynamic_programming_matrix.pl
Jong
>
> I'm not sure if UnivAln is exactly what I'm looking for. I'm working with
> short probes 18 and 20mers. I'm comparing probes to check to make sure they
> are not to similar. So I will put two probes into UnivAln and get the
> consensus. But how would I get the program to slide the probes against each
> other to get the best match?
>
> Example
>
> TATATATATATATATATA
> ATATATATATATATATAT
> seems to produce no match
>
> what I want to see is
> TATATATATATATATAT
> ATATATATATATATATA
> which would produce 16 matches.
>
> Anyhow if this is confusing I apologize, but would appreciate any help you
> might be able to give.
>
> Thank you
> Clark Nelson
> cjnelso3@students.wisc.edu
>
> --
> ^X
>
>
-- Jong Park : a Biology student -- Tel work: 49-4613
... Science is not a craftsmanship nor a simple clever discovery
of nature. It is not a task for any clever mind with a good
memory and a good understanding on the techniques involved with
a good class mark from previous schools. It requires a maturity
as an educated human personality with proper insights on the
processes and consequences of the work involved to the society
and all humanbeings in the far extent.
It is a scholastic activity with unselfish and unbiased causes..
--------------E48DD72FE129715BCCFF0D96
Content-Type: application/x-perl;
name="make_dynamic_programming_matrix.pl"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="make_dynamic_programming_matrix.pl"
#!/usr/bin/perl
#______________________________________________________________________________
# Title : make_dynamic_programming_matrix.pl
# Usage :
# Function :
# Example :
# Keywords :
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
$seq2='HEAGAWGHEE';
$seq1='PAWHEAE';
%seq=('seq', $seq1, 'seq2', $seq2);
%exchange_matrix=%{&get_homology_matrix("m=Blosum62")};
%two_seq_homol_matrix=%{&make_2_seq_comparison_matrix(\%exchange_matrix, \%seq)};
@temp=&make_dynamic_programming_matrix(\%two_seq_homol_matrix, \%seq);
@DP_matrix=@{$temp[0]};
@path_matrix=@{$temp[1]};
&show_array_matrix(\@DP_matrix);
&show_array_matrix(\@path_matrix);
#______________________________________________________________________________
# Title : make_dynamic_programming_matrix
# Usage :
# Function :
# Example :
# Keywords :
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub make_dynamic_programming_matrix{
my(%DP_matrix, @DP_matrix, %two_seq_homol_matrix, $i, $j, $gap_penalty, @seq_names, @seq1,
@seq2, @path_matrix);
$gap_penalty=10; ## 11 is best for Blosum62
%two_seq_homol_matrix=%{$_[0]};
%seq=%{$_[1]};
@seq_names=keys %seq;
@seq1=(split(//, $seq{$seq_names[0]}));
@seq2=(split(//, $seq{$seq_names[1]}));
for($i=0; $i < @seq1; $i++){
for($j=0; $j < @seq2; $j++){
if($i==0 and $j !=0){
$DP_matrix{$i}{$j}= ($j)*( - $gap_penalty);
$DP_matrix[$i][$j]= ($j)*( - $gap_penalty);
$path_matrix[$i][$j]=[$i, 'x'];
}elsif($j==0 and $i !=0){
$DP_matrix{$i}{$j}= ($i)*( - $gap_penalty);
$DP_matrix[$i][$j]= ($i)*( - $gap_penalty);
$path_matrix[$i][$j]=['x', $j];
}elsif($i==0 and $j ==0){
$DP_matrix{$i}{$j}=0;
$DP_matrix[$i][$j]=0;
$path_matrix[$i][$j]=[0, 0];
}else{
$diagonal_case=$DP_matrix{$i-1}{$j-1} + $two_seq_homol_matrix{$seq1[$i]}{$seq2[$j]};
$seq1_gap_case=$DP_matrix{$i-1}{$j} - $gap_penalty;
$seq2_gap_case=$DP_matrix{$i}{$j-1} - $gap_penalty;
$diagonal_case=$DP_matrix[$i-1][$j-1] + $two_seq_homol_matrix{$seq1[$i]}{$seq2[$j]};
$seq1_gap_case=$DP_matrix[$i-1][$j] - $gap_penalty;
$seq2_gap_case=$DP_matrix[$i][$j-1] - $gap_penalty;
if($diagonal_case >= $seq1_gap_case and $diagonal_case >= $seq2_gap_case){
$DP_matrix{$i}{$j}= $diagonal_case;
$DP_matrix[$i][$j]= $diagonal_case;
$path_matrix[$i][$j]=[$i-1, $j-1];
$path_matrix{$i}{$j}=[$i-1, $j-1];
$new_seq1 .=$seq1[$i];
$new_seq2 .=$seq2[$j];
print "\n \$diagonal_case: $DP_matrix{$i-1}{$j-1} + $two_seq_homol_matrix{$seq1[$i]}{$seq2[$j]}\n";
}elsif($seq1_gap_case > $diagonal_case and $seq1_gap_case >= $seq2_gap_case){
$DP_matrix{$i}{$j}= $seq1_gap_case;
$DP_matrix[$i][$j]= $seq1_gap_case;
$path_matrix[$i][$j]=[$i-1, $j];
$path_matrix{$i}{$j}=[$i-1, $j];
print "\n \$seq1_gap_case: $DP_matrix{$i-1}{$j} - $gap_penalty\n";
$new_seq1 .='-';
$new_seq2 .=$seq2[$j];
}elsif($seq2_gap_case > $seq1_gap_case and $seq2_gap_case > $diagonal_case){
$DP_matrix{$i}{$j}= $seq2_gap_case;
$DP_matrix[$i][$j]= $seq2_gap_case;
$path_matrix[$i][$j]=[$i, $j-1];
$path_matrix{$i}{$j}=[$i, $j-1];
#splice(@seq2, $j, 1, '-');
print "\n \$seq2_gap_case: $DP_matrix{$i}{$j-1} - $gap_penalty\n";
$new_seq1 .=$seq1[$i];
$new_seq2 .='-';
}else{
print "\# ERROR \n\n"; exit;
}
}
}
}
print "\n$new_seq1";
print "\n$new_seq2";
#&show_hash(\%DP_matrix);
@keys=keys %DP_matrix;
@values=values %DP_matrix;
print "\nX |\t@seq2\n";
&show_hash_matrix(\%DP_matrix);
&show_array_matrix(\@DP_matrix);
print "\n@keys\n@values";
return(\@DP_matrix, \@path_matrix);
}
#______________________________________________________________________________
# Title : make_2_seq_comparison_matrix
# Usage :
# Function :
# Example :
# Keywords : make_protein_seq_comparison_matrix, make_matrix,
# make_sequence_comparison_matrix
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub make_2_seq_comparison_matrix{
my(%exchange_matrix, %seq, %comparison_matrix, @seq1, @seq2, @seq_names,
$i, $j, $k);
%exchange_matrix=%{$_[0]};
%seq=%{$_[1]};
#&show_hash(\%exchange_matrix);
#&show_hash(\%seq);
@seq_names=keys %seq;
@seq1=split(//, $seq{$seq_names[0]});
@seq2=split(//, $seq{$seq_names[1]});
for($i=0; $i< @seq1; $i++){
for($j=0; $j<@seq2; $j++){
$comparison_matrix{$seq1[$i]}{$seq2[$j]}=$exchange_matrix{$seq1[$i]}{$seq2[$j]};
print "\n$comparison_matrix{$seq1[$i]}{$seq2[$j]}";
}
}
return(\%comparison_matrix);
}
#______________________________________________________________________________
# Title : get_homology_matrix
# Usage :
# Function :
# Example :
# Keywords :
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub get_homology_matrix{
my(%matrix, $matrix_type);
if($_[0]=~/Blosum62/i){ $matrix_type='Blosum62' };
if($matrix_type eq 'Blosum62'){
print "\n \$matrix_type is $matrix_type\n";
$blosum_matrix_string="
# Matrix made by matblas from blosum62.iij
# * column uses minimum score
# BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
# Blocks Database = /data/blocks_5.0/blocks.dat
# Cluster Percentage: >= 62
# Entropy = 0.6979, Expected = -0.5209
# A R N D C Q E G H I L K M F P S T W Y V B Z X *
# A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0 -2 -1 0 -4
# R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3 -1 0 -1 -4
# N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3 3 0 -1 -4
# D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3 4 1 -1 -4
# C 0 -3 -3 -3 9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
# Q -1 1 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2 0 3 -1 -4
# E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4
# G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3 -1 -2 -1 -4
# H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3 0 0 -1 -4
# I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3 -3 -3 -1 -4
# L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1 -4 -3 -1 -4
# K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2 0 1 -1 -4
# M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1 -3 -1 -1 -4
# F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6 -4 -2 -2 1 3 -1 -3 -3 -1 -4
# P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7 -1 -1 -4 -3 -2 -2 -1 -2 -4
# S 1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2 0 0 0 -4
# T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0 -1 -1 0 -4
# W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11 2 -3 -4 -3 -2 -4
# Y -2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2 2 7 -1 -3 -2 -1 -4
# V 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4 -3 -2 -1 -4
# B -2 -1 3 4 -3 0 1 -1 0 -3 -4 0 -3 -3 -2 0 -1 -4 -3 -3 4 1 -1 -4
# Z -1 0 0 1 -3 3 4 -2 0 -3 -3 1 -1 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4
# X 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 0 0 -2 -1 -1 -1 -1 -1 -4
# * -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 1
";
%matrix=%{&read_matrix(\$blosum_matrix_string)};
}
return(\%matrix);
}
#________________________________________________________________________________
# Title : read_matrix
# Usage : %matrix=%{&read_matrix(\$string)};
# Function : Makes similarrity matrix hash(reflexive, so it has AT as well as TA)
# %matrix looks like this: $matrix{X}{Y}= 4
# Example :
# Keywords : get_2D_aa_matrix, read_seq_matrix
# Options :
# $reflexive_combi=r by r -r # both direction AC=-1, CA=-1
# Category :
# Version : 1.0
#--------------------------------------------------------------------------------
sub read_matrix{
my ($reflexive_combi, %matrix, %out_matrix, @matrix_lines,
@residue, $residue_num, $first_residue, $second_residue,
@matrix_val, $paired_residues, @sorted, $i, $j, $k, $l);
$reflexive_combi='r'; ## both direction AC=-1, CA=-1
@matrix_lines=split(/\n/, ${$_[0]});
if(@matrix_lines < 1){ print "\n Error, \@matrix_lines is empty in read_matrix sub\n"; exit }
for($j=0; $j < @matrix_lines; $j++){
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Skipping all non-matrix contents
#_________________________________________________________
if($matrix_lines[$j]=~/^ *# .+\S\S\S/ or $matrix_lines[$j]=~/^ *# *\W/){ next }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# matching top line: ' A R N D C Q E G H I ..' line
#___________________________________________________
if(@residue=$matrix_lines[$j]=~/ ([\*\w])/g){
$residue_num=@residue;
if($residue_num < 20){ next }
#print "\n@residue $residue_num\n,,,,,,,\n";
$j++;
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# matching ' A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -'
# ' R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -'
#___________________________________________________
for($k=0; $k<@residue; $k++){
$first_residue=$residue[$k];
#print "$first_residue\n";
#print $matrix_lines[$j] if $verbose;
if(@matrix_val=$matrix_lines[$j]=~/ *([\-\d]+)/g){
#print "\n-> @matrix_val\n";
for($l=0; $l< @matrix_val; $l++){
$second_residue=$residue[$l];
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Taking only non-reflexive ones
#________________________________________________
if($reflexive_combi){
$paired_residues=join('', $first_residue, $second_residue);
}else{
@sorted = sort($first_residue, $second_residue);
$paired_residues=join('', @sorted);
}
$matrix{$paired_residues}=$matrix_val[$l] unless $matrix{$paired_residues};
#print "\n$paired_residues $matrix{$paired_residues}";
$out_matrix{$first_residue}{$second_residue}=$matrix_val[$l];
#print "\n$out_matrix{$first_residue}{$second_residue}";
}
$j++;
}
}
}
}
if($string_pair_hash_out==1){ ## out looks : $matrix_lines{'XY'} =4
return(\%matrix);
}else{ ## out looks : $matrix_lines{'X'}{'Y'}=4
print "\n read_matrix returns a matrix form of \$matrix\{A\}\{C\} \n";
return(\%out_matrix);
}
}
#________________________________________________________________________
# Title : show_hash
# Usage : &show_hash(\@input_array);
# Function : for debugging purpose. Shows any array elem line by line.
# the line is 60 elements long (uses recursion)
# Example : Output: item1
# Output: item2
# Output: item3
# Warning : There is a global variable: $show_hash_option
# It tries to detect any given sting which is joined by ','
# Keywords :
# Options : -s or -S or s or S for spaced output. Eg)
# seq1 1 1 1 1 1 1 1 1 1 1 1 1
#
# instead of
# seq1 111111111111
#
# -h or -H or h or H for horizontal line of '---------...'
#
# Returns :
# Argument :
# Category :
# Version : 1.9
#--------------------------------------------------------------------
sub show_hash{
my($k, $i, $t, @in2, $in, $LEN, %TEM, $out_in_html_form ); ## You should not put $show_hash_option
my(@in)=@_; ## and $horizontal_line in my !!!
my($KL)=18; # default keys string length;
my($VL)=80; # default values string length;
my($GAP)=2; # default space between keys and values
my($horizontal_line, $show_hash_optionXX, $Hash_counter, @line);
## This is to get the option of 'space' to make spaced output.
for($t=0; $t < @in; $t++){
if($in[$t] =~/^[-]+[sS][pace]*$/){
$show_hash_optionXX = 1; splice(@in, $t, 1);
}elsif(${in[$t]} =~/^[-]+[sS][pace]*$/){
$show_hash_optionXX = 1; splice(@in, $t, 1);
}elsif($in[$t] =~/^[-]+[hH][rR]*$/){
$horizontal_line = 1; splice(@in, $t, 1);
}elsif($in[$t] =~/html/i){
$out_in_html_form="HTML"; splice(@in, $t, 1);
}elsif(${in[$t]} =~/html/i){
$out_in_html_form='HTML'; splice(@in, $t, 1);
}
}
######## Main loop #################
if($horizontal_line ==1){ ## This puts the delimiter '--------------( )'
$Hash_counter ++;
print "\n","-"x78,"(${Hash_counter}th hash)", "\n";
}
for($k=0; $k < @in; $k++){
if(ref($in[$k]) eq 'ARRAY'){ ### When the hashes were given in array ref.
&show_hash(@{$in[$k]}, $show_hash_optionXX, $horizontal_line, $out_in_html_form);
print "\n";
}elsif(ref($in[$k]) eq 'HASH'){ ### recursion
&show_hash(%{$in[$k]}, $out_in_html_form, $show_hash_optionXX, $horizontal_line);
print "\n";
}elsif(ref($in[$k+1]) eq 'HASH'){ ### recursion
&show_hash(%{$in[$k+1]}, $horizontal_line, $out_in_html_form, $show_hash_optionXX); print "\n";
}elsif(ref($in[$k]) eq 'SCALAR'){
if($out_in_html_form){ print ${$_[$k]}, "<br>\n";
}else{ print ${$_[$k]}, "\n"; }
}elsif( !ref($in[$k]) ){
if( !ref($in[$k+1]) && defined($in[$k+1]) ){
if($show_hash_optionXX){ #### space option checking.
%TEM = @in;
$LEN = ${&max_elem_string_array_show_hash(keys %TEM)};
if($LEN > $KL){ $KL = $LEN + $GAP +2};
if($out_in_html_form){
printf ("<br>%-${KL}s ", $in[$k]); $k++;
printf ("%-${VL}s<br>\n","@line");
}else{
printf ("%-${KL}s ", $in[$k]); $k++;
printf ("%-${VL}s\n","@line");
}
}else{ ### If not option is set, just write
%TEM = @in;
$LEN = ${&max_elem_string_array_show_hash( keys %TEM)};
if($LEN > $KL){ $KL = $LEN + $GAP +2};
if($out_in_html_form){
printf ("<br>%-${KL}s ", $in[$k]); $k++; # print $in[$k], "\t"; $k++;
printf ("%-${VL}s<br>\n",$in[$k]); # print $in[$k], "\n";
}else{
printf ("%-${KL}s ", $in[$k]); $k++; # print $in[$k], "\t"; $k++;
printf ("%-${VL}s\n",$in[$k]); # print $in[$k], "\n";
}
}
}
#________________________________________________________
# Title : max_elem_string_array_show_hash
# Keywords : largest string length of array
# Function : gets the largest string length of element of any array of numbers.
# Usage : ($out1, $out2)=@{&max_elem_array(\@array1, \@array2)};
# ($out1) =${&max_elem_array(\@array1) };
# Argument : numerical arrays
# Version : 1.1
#-------------------------------------------------------
sub max_elem_string_array_show_hash{
my(@input, $i, $max_elem); @input = @{$_[0]} || @_;
for($i=0; $i< @input ; $i++){
$max_elem = length($input[0]);
if (length($input[$i]) > $max_elem){
$max_elem = length($input[$i]);
}
}
return(\$max_elem);
}
#####################################insert_gaps_in_seq_hash
}
}
}
#________________________________________________________________________
# Title : show_hash_matrix
# Usage : &show_hash_matrix(\%input_hash);
# Function :
#
# Example : &show_hash_matrix(\%DP_matrix);
# 0 0 -0 -0 -0 -0 -0 -0 -0 -0 -0
# 1 -0 -1 4 0 4 -3 0 -2 -1 -1
# 2 -0 -3 -4 2 -3 15 4 -2 -5 -4
# 3 -0 0 -5 -6 0 4 13 12 1 -5
# 4 -0 5 -1 -7 -7 -3 2 13 17 6
# 5 -0 -1 9 -1 -3 -10 -3 2 12 16
# 6 -0 5 -2 7 -2 -6 -12 -3 7 17
#
# Warning : There is a global variable: $show_matrix_option
# Keywords : show_matrix, show_matrix_hash
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.1
#--------------------------------------------------------------------
sub show_hash_matrix{
my($i, $j, %hash, @keys1, @keys2 ); ## You should not put $show_matrix_option
%hash=%{$_[0]};
@keys1=keys %hash;
print "\n";
for($i=0; $i< @keys1; $i++){
print "$keys1[$i] |\t";
@keys2=keys %{$hash{$keys1[$i]}};
for($j=0; $j< @keys2; $j++){
if(ref($hash{$keys[$i]}{$keys2[$j]}) eq 'ARRAY'){
print "@{$hash{$keys[$i]}{$keys2[$j]}}\t";
}elsif(ref($hash{$keys[$i]}{$keys2[$j]}) eq 'SCALAR'){
print "${$hash{$keys[$i]}{$keys2[$j]}}\t";
}else{
print "$hash{$keys[$i]}{$keys2[$j]}\t";
}
}
print "\n";
}
}
#________________________________________________________________________
# Title : show_array_matrix
# Usage : &show_array_matrix(\@input_array);
# Function :
#
# Example : &show_array_matrix(\@DP_matrix);
# 0 | 0 -0 -0 -0 -0 -0 -0 -0 -0 -0
# 1 | -0 -1 4 0 4 -3 0 -2 -1 -1
# 2 | -0 -3 -4 2 -3 15 4 -2 -5 -4
# 3 | -0 0 -5 -6 0 4 13 12 1 -5
# 4 | -0 5 -1 -7 -7 -3 2 13 17 6
# 5 | -0 -1 9 -1 -3 -10 -3 2 12 16
#
# Warning : There is a global variable: $show_matrix_option
# Keywords : show_matrix, show_matrix_array
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.1
#--------------------------------------------------------------------
sub show_array_matrix{
my($i, $j, @array, @keys1, @keys2 ); ## You should not put $show_matrix_option
@array=@{$_[0]};
print "\n";
for($i=0; $i< @array; $i++){
print "$i |\t";
$dim_size=@{$array[0]};
for($j=0; $j< $dim_size; $j++){
if(ref($array[$i][$j]) eq 'ARRAY'){
print "@{$array[$i][$j]}\t";
}elsif(ref($array[$i][$j]) eq 'SCALAR'){
print "$array[$i][$j]\t";
}else{
print "$array[$i][$j]\t";
}
}
print "\n";
}
}
#________________________________________________________________________
# Title : show_array_path_matrix
# Usage : &show_array_path_matrix(\@input_array);
# Function :
#
# Example : &show_array_path_matrix(\@DP_matrix);
# 0 | 0 -0 -0 -0 -0 -0 -0 -0 -0 -0
# 1 | -0 -1 4 0 4 -3 0 -2 -1 -1
# 2 | -0 -3 -4 2 -3 15 4 -2 -5 -4
# 3 | -0 0 -5 -6 0 4 13 12 1 -5
# 4 | -0 5 -1 -7 -7 -3 2 13 17 6
# 5 | -0 -1 9 -1 -3 -10 -3 2 12 16
#
# Warning : There is a global variable: $show_matrix_option
# Keywords : show_matrix, show_matrix_array
# Options :
# Returns :
# Argument :
# Category :
# Version : 1.0
#--------------------------------------------------------------------
sub show_array_path_matrix{
my($i, $j, @array, @keys1, @keys2 ); ## You should not put $show_matrix_option
@array=@{$_[0]};
print "\n";
for($i=0; $i< @array; $i++){
print "$i |\t";
$dim_size=@{$array[0]};
for($j=0; $j< $dim_size; $j++){
print "@{$array[$i][$j]}\t";
}
print "\n";
}
}
__END__
@seq_names=keys %seq;
$seq1_len=@seq1=(split(//, $seq{$seq_names[0]}));
$seq2_len=@seq2=(split(//, $seq{$seq_names[1]}));
$seq1_len=1;
$seq2_len=1;
%final_seq=%{&align_2_sequences(\@DP_matrix, \%seq, \$seq1_len, \$seq2_len, \$gap_penalty)};
#______________________________________________________________________________
# Title : align_2_sequences
# Usage :
# Function :
# Example :
# Keywords :
# Options :
# Author : jong@salt2.med.harvard.edu,
# Category :
# Version : 1.0
#------------------------------------------------------------------------------
sub align_2_sequences{
@DP_matrix=@{$_[0]};
@seq_names=keys %seq;
@seq1=(split(//, $seq{$seq_names[0]}));
@seq2=(split(//, $seq{$seq_names[1]}));
$seq1_len=${$_[2]};
$seq2_len=${$_[3]};
$gap_penalty=${$_[4]};
if($DP_matrix[$seq1_len][$seq2_len]== $DP_matrix[$seq1_len -1][$seq2_len] - $gap_penalty){
%final_seq=%{&align_2_sequences(\@DP_matrix, \%seq, \$seq1_len, \$seq2_len, \$gap_penalty)};
}
--------------E48DD72FE129715BCCFF0D96--
=========== Bioperl Project Mailing List Message Footer =======
Project URL: http://bio.perl.org/
For info about how to (un)subscribe, where messages are archived, etc:
http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
====================================================================