2009-09-18 21 views
1

J'ai besoin d'un peu d'aide avec ce code. Je connais les sections qui devraient être récursives, ou du moins je pense que je le fais, mais je ne suis pas sûr de la façon de l'appliquer. J'essaie de mettre en œuvre un programme de recherche de chemin à partir d'une matrice d'alignement qui va trouver plusieurs routes à la valeur zéro. Par exemple, si vous excutez mon code et insérez CGCA comme première séquence et CACGTAT comme deuxième séquence, et 1, 0 et -1 comme score de correspondance, de discordance et d'écart. Le programme donne au large de la voie que HDHHDD et Aligement commeTechniques de récursivité Perl?

CACGTAT

CGC - A-.

Cependant il y a des chemins et aligments plus possible, alors cela, sauf que je ne sais pas combien. Ce que je veux faire est d'avoir une partie de mon code en boucle sur lui-même et trouver d'autres chemins et des alignements, en utilisant le même code que la première fois, jusqu'à ce qu'il soit à court d'alignements possibles. Le meilleur moyen que j'ai trouvé sur le net pour faire ceci est la récursivité, sauf que personne ne peut expliquer comment le faire. Dans ce cas, il devrait y avoir deux autres chemins et aligmens HDDDHHD et CACGTAT, et C - GCA-, et. HDDDDHH, CACGTAT et --CGCA-. Je ne sais pas comment coder pour effectuer cette tâche.

# Implementation of Needleman and Wunsch Algorithm 

my($seq1, $len1, $seq2, $len2, $data, @matrix, $i, $j, $x, $y, $val1, $val2); 
my($val3, $pathrow, $pathcol, $seq1loc, $seq2loc, $gapscore, $matchscore, $mismatchscore); 

#first obtain the data from the user. 
print "Please enter the first sequence for comaprsion\n"; 
$seq1=<STDIN>; 
chomp $seq1; 

print "Please enter the second sequence for comparsion\n"; 
$seq2=<STDIN>; 
chomp $seq2; 


# adding extra characters so sequences align with matrix 
# saves some calculations later on 
$seq1 = " " . $seq1; 
$seq2 = " " . $seq2; 
$len1 = length($seq1); 
$len2 = length($seq2); 

print "Enter the match score: "; 
$matchscore=<STDIN>; 
chomp $matchscore; 

print "Enter the mismatch score: "; 
$mismatchscore=<STDIN>; 
chomp $mismatchscore; 

print "Enter the gap score: "; 
$gapscore=<STDIN>; 
chomp $gapscore; 

# declare a two dimensional array and initialize to spaces 
# array must contain one extra row and one extra column 
@matrix =(); 
for($i = 0; $i < $len1; $i++){ 
    for($j = 0; $j < $len2; $j++){ 
     $matrix[$i][$j] = ' '; 
    } 
} 

# initialize 1st row and 1st column of matrix 
$matrix[0][0] = 0; 
for ($i = 1; $i < $len1; $i ++){ 
    $matrix[$i][0] = $matrix[$i-1][0] + $gapscore; 
} 
for ($i = 1; $i < $len2; $i ++){ 
    $matrix[0][$i] = $matrix[0][$i-1] + $gapscore; 
} 


# STEP 1: 
# Fill in rest of matrix using the following rules: 
# determine three possible values for matrix[x][y] 
# value 1 = add gap score to matrix[x][y-1] 
# value 2 = add gap score to matrix[x-1][y] 
# value 3 = add match score or mismatch score to 
#   matrix[x-1][y-1] depending on nucleotide 
#   match for position x of $seq1 and position y 
#   of seq2 
# place the largest of the three values in matrix[x][y] 
# 
# Best alignment score appears in matrix[$len1][$len2]. 

for($x = 1; $x < $len1; $x++){ 
    for($y = 1; $y < $len2; $y++){ 
$val1 = $matrix[$x][$y-1] + $gapscore; 
$val2 = $matrix[$x-1][$y] + $gapscore; 
if (substr($seq1, $x, 1) eq substr($seq2, $y, 1)){ 
      $val3 = $matrix[$x-1][$y-1] + $matchscore; 
} 
else{ 
    $val3 = $matrix[$x-1][$y-1] + $mismatchscore; 
} 
if (($val1 >= $val2) && ($val1 >= $val3)){ 
    $matrix[$x][$y] = $val1; 
} 
elsif (($val2 >= $val1) && ($val2 >= $val3)){ 
    $matrix[$x][$y] = $val2; 
} 
else{ 
    $matrix[$x][$y] = $val3; 
} 
    } 
} 

# Display scoring matrix 
print "MATRIX:\n"; 
for($x = 0; $x < $len1; $x++){ 
    for($y = 0; $y < $len2; $y++){ 
print "$matrix[$x][$y] "; 
    } 
    print "\n"; 
} 
print "\n";  

# STEP 2: 
# Begin at matrix[$len1][$len2] and find a path to 
# matrix[0][0]. 
# Build string to hold path pattern by concatenating either 
# 'H' (current cell produced by cell horizontally to left), 
# 'D' (current cell produced by cell on diagonal), 
# 'V' (current cell produced by cell vertically above) 
# ***This is were I need help I need this code to be recursive, so I can find more then one path*** 

$pathrow = $len1-1; 
$pathcol = $len2-1; 

while (($pathrow != 0) || ($pathcol != 0)){ 
    if ($pathrow == 0){ 
     # must be from cell to left 
     $path = $path . 'H'; 
     $pathcol = $pathcol - 1; 
    } 
    elsif ($pathcol == 0){ 
     # must be from cell above 
     $path = $path . 'V'; 
     $pathrow = $pathrow - 1; 
    } 
    # could be from any direction 
    elsif (($matrix[$pathrow][$pathcol-1] + $gapscore) == $matrix[$pathrow][$pathcol]){ 
     # from left 
     $path = $path . 'H'; 
     $pathcol = $pathcol - 1; 
    } 
    elsif (($matrix[$pathrow-1][$pathcol] + $gapscore) == $matrix[$pathrow][$pathcol]){ 
     # from above 
     $path = $path . 'V'; 
     $pathrow = $pathrow - 1; 
    } 
    else{ 
     # must be from diagonal 
     $path = $path . 'D'; 
     $pathrow = $pathrow - 1; 
     $pathcol = $pathcol - 1; 
    } 


} 



print "Path is $path\n"; 

# STEP 3: 
# Determine alignment pattern by reading path string 
# created in step 2. 
# Create two string variables ($alignseq1 and $alignseq2) to hold 
# the sequences for alignment. 
# Reading backwards from $seq1, $seq2 and path string, 
# if string character is 'D', then 
#  concatenate to front of $alignseq1, last char in $seq1 
#  and to the front of $alignseq2, last char in $seq2 
# if string character is 'V', then 
#  concatenate to front of $alignseq1, last char in $seq1 
#  and to the front of $alignseq2, the gap char 
# if string character is 'H', then 
#  concatenate to front of $alignseq1 the gap char 
#  and to the front of $alignseq2, last char in $seq2 
# Continue process until path string has been traversed. 
# Result appears in string $alignseq1 and $seq2 
***#I need this code to be recursive as well.*** 

$seq1loc = $len1-1; 
$seq2loc = $len2-1; 
$pathloc = 0; 
print length($path); 
while ($pathloc < length($path)){ 
    if (substr($path, $pathloc, 1) eq 'D'){ 
     $alignseq1 = substr($seq1, $seq1loc, 1) . $alignseq1; 
     $alignseq2 = substr($seq2, $seq2loc, 1) . $alignseq2; 
     $seq1loc--; 
     $seq2loc--; 
    } 
    elsif (substr($path, $pathloc, 1) eq 'V'){ 
     $alignseq1 = substr($seq1, $seq1loc, 1) . $alignseq1; 
     $alignseq2 = '-' . $alignseq2; 
     $seq1loc--; 
    }  
    else{ # must be an H 
     $alignseq1 = '-' . $alignseq1; 
     $alignseq2 = substr($seq2, $seq2loc, 1) . $alignseq2; 
     $seq2loc--; 
    }  
    $pathloc++; 
} 

print "\nAligned Sequences:\n"; 
print "$alignseq2 \n"; 
print "$alignseq1 \n"; 

# statement may be needed to hold output screen 
print "Press any key to exit program"; 
$x = <STDIN>; 

Si quelqu'un se demande s'il s'agit d'un algorithme needleman-wunsch. Toute aide ici serait grandement sollicitée.

+1

Que diriez-vous d'un sous-programme? –

Répondre

7

Je ne peux pas donner une réponse, parce que je ne comprends pas exactement ce que vous essayez de faire, mais je peux offrir des conseils généraux.

commencer à organiser votre code en sous-routines discrets qui exécutent des tâches strictement définies. En outre, les sous-programmes qui implémentent vos algorithmes centraux ne doivent pas être orientés vers la réception d'une entrée du clavier et la production de sortie à l'écran; ils devraient plutôt recevoir des données en tant qu'arguments et renvoyer leurs résultats. S'il y a un besoin d'entrée d'utilisateur ou de sortie d'écran, ces tâches devraient être dans des sous-routines séparées, non associées à vos algorithmes principaux.

Une première étape (et partielle) de ce chemin est de prendre votre programme entier, de le placer dans une définition de sous-routine, puis d'appeler le sous-programme avec les arguments requis. Au lieu d'imprimer ses résultats, le sous-programme doit les retourner - en particulier, une référence à @matrix ainsi que les valeurs pour $path, $alignseq1, $alignseq2.

sub NW_algo { 
    my ($seq1, $seq2, $matchscore, $mismatchscore, $gapscore) = @_; 

    # The rest of your code here, but with all print 
    # statements and <STDIN> inputs commented out. 

    return \@matrix, $path, $alignseq1, $alignseq2; 
} 

my(@return_values) = NW_algo('CGCA', 'CACGTAT', 1, 0, -1); 

Print_matrix($return_values[0]); 

sub Print_matrix { 
    for my $m (@{$_[0]}){ 
     print join(' ', @$m), "\n"; 
    } 
} 

À ce stade, vous aurez un algorithme qui peut être invoqué par un autre code, rendant plus facile à tester et déboguer votre programme à l'avenir. Par exemple, vous pouvez définir différents ensembles de données d'entrée et exécuter NW_algo() sur chaque ensemble. Alors seulement, il sera possible de penser à la récursivité ou à d'autres techniques.

5

Depuis Needleman-Wunsch est un algorithme dynamique de programmation, la plupart du travail est déjà fait au moment où vous calculez votre matrice DP. Une fois que vous avez votre matrice DP, vous êtes censé faire marche arrière à travers la matrice pour trouver l'alignement optimal. Le problème est un peu comme la géométrie des taxis, sauf que vous pouvez vous déplacer en diagonale. Essentiellement, quand vous avez besoin de faire marche arrière à travers la matrice, au lieu de choisir entre monter, gauche ou diagonalement, vous faites tous les trois en faisant trois appels récursifs, et chacun d'eux s'appelle pour chacun, vers le haut ou vers la gauche, vous atteignez votre point de départ. Le chemin tracé par chaque brin de récursivité attirera chaque alignement.

EDIT: Donc, fondamentalement, vous avez besoin de mettre l'étape 2 dans une sous-procédure (qui prend position et la voie tracée jusqu'à présent), il peut donc s'appeler plus et plus.Bien sûr, après avoir défini la procédure que vous devez faire un appel à pour commencer réellement le processus de traçage:

sub tracePaths { 
    $x = shift; 
    $y = shift; 
    $pathSoFar = shift; # assuming you're storing your path as a string 
    # 
    # ... do some work ... 
    # 
    tracePaths($x - 1, $y, $pathSoFar . something); 
    tracePaths($x, $y - 1, $pathSoFar . somethingelse); 
    tracePaths($x - 1, $y - 1, $pathSoFar . somethingelselse); 
    # 
    # 
    # 
    if(reached the end) return $pathSoFar; 
} 
# launch the recursion 
tracePaths(beginningx, beginningy, ""); 
4

Cela ne parle pas spécifiquement à votre problème, mais vous devriez peut-être vérifier le livre Higher Order Perl. Il explique comment utiliser beaucoup de techniques de niveau supérieur (telles que la récursivité).