Site WWW de Laurent Bloch
Slogan du site

ISSN 2271-3905
Cliquez ici si vous voulez visiter mon autre site, orienté vers des sujets moins techniques.

Pour recevoir (au plus une fois par semaine) les nouveautés de ce site, indiquez ici votre adresse électronique :

ChatGPT m’aide à programmer en R
Article mis en ligne le 10 août 2026

par Laurent Bloch

R est-il un langage de programmation ?

Le système R est un logiciel de calcul statistique, mais s’agit-il d’un langage de programmation ? Je dois dire que je le considérais comme un objet hybride, un logiciel de calcul avec des capacités de programmation, mais pas vraiment Turing-équivalent, comme SAS ou SQL. Mais en fait ces systèmes ont tous évolué pour devenir Turing-équivalents, c’est-à-dire capables de calculer tout algorithme.

Dans mon exploration des langages de programmation j’hésitais quant au choix du suivant après TypeScript/JavaScript. Puis j’ai consulté l’index TIOBE du palmarès de popularité des langages de programmation, et constaté la présence de R au neuvième rang, entre SQL et Rust. Cela m’a convaincu d’approfondir avec R, d’autant plus qu’à l’origine des temps je suis statisticien, et que R est très présent dans le domaine de la bioinformatique, auquel je participe régulièrement.

Débuter la programmation R avec ChatGPT

Bonjour ChatGPT,

Ces jours-ci je me lance dans l’analyse de séquences biologiques avec le langage R.

Voici une exemple de fichier qui contient deux séquences à aligner, sequences-1.txt :

ACGTC
AGTC

Voici le programme de lecture du fichier et d’appel du programme d’analyse :

#!/usr/bin/env Rscript

args <- commandArgs(trailingOnly = TRUE)

seqfile <- args[1]
gap <- as.numeric(args[2])
mismatch <- as.numeric(args[3])
match <- as.numeric(args[4])

source("./Needleman-Wunsch-LB.R")

Sequences <- readLines(seqfile)

NeedlemanWunsch(Sequences[1], Sequences[2], gap, mismatch, match)

invoqué ainsi :

./NW.R sequences-1.txt -1 -1 1

Voici le source de la première version du programme, largement emprunté à un exemple de Julius Kittler. C’est un programme monolithique, qui :

– dans premier temps initialise la matrice M des scores et la matrice D des décalages pour alignement ;
– dans un second temps calcule les scores et les décalages ;
– dans un troisième temps effectue le retour sur trace (backtracking) pour obtenir l’alignement.

# References: https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus

NeedlemanWunsch = function(seq1, seq2, gap, mismatch, match){
  
  # Stop conditions
  stopifnot(gap <= 0) # check if penalty negative
  stopifnot(mismatch <= 0)  # check if penalty negative
  stopifnot(match >= 0)  # check if score positive
  
  # Initialize col and rownames for matrices
  len1 = nchar(seq1); len2 = nchar(seq2) # Save number of chars in each sequence
  seq1 = unlist(strsplit(seq1, split = "")) # convert seq to character vector
  seq2 = unlist(strsplit(seq2, split = "")) # convert seq to character vector

  # Initialize matrix M (for scores)
  M = matrix(0, nrow = len1 + 1, ncol = len2 + 1) # Initialize matrix
  rownames(M) = c("-", seq1) # assign seq chars to matrix names
  colnames(M) = c("-", seq2) # assign seq chars to matrix names
  M[1, ] = cumsum(c(0, rep(gap, len2))) # Fill 1st row with gap penalites
  M[, 1] = cumsum(c(0, rep(gap, len1))) # Fill 1st col with gap penalites
  
  # Initialize matrix D (for directions)
  D = matrix(0, nrow = len1 + 1, ncol = len2 + 1) # Initialize matrix
  rownames(D) = c("-", seq1) # assign seq chars to matrix names
  colnames(D) = c("-", seq2) # assign seq chars to matrix names
  D[1, ] = rep("hor") # Fill 1st row with "hor" for horizontal moves
  D[, 1] = rep("ver") # Fill 1st col with "ver" for vertical moves
  type = c("dia", "hor", "ver") # Lookup vector

  # Compute scores and save moves
  for (i in 2:(len1 + 1)){# for every (initially zero) row
    for (j in 2:(len2 + 1)){# for every (initially zero) col
      hor = M[i, j - 1] + gap # horizontal move = gap for seq1
      ver = M[i - 1, j] + gap # vertical move = gap for seq2 
      dia = ifelse(rownames(M)[i] == colnames(M)[j], # diagonal = ifelse(chars equal, match, mismatch) 
                   M[i - 1, j - 1] + match, 
                   M[i - 1, j - 1] + mismatch)
      M[i, j] = max(dia, hor, ver) # Save current (best) score in M
      D[i, j] = type[which.max(c(dia, hor, ver))] # Save direction of move in D
    }
  }
    
  # Backtracing
  align1 = c(); align2 = c() # Note: length of final alignments is unknown at this point
  
  while(i > 1 && j > 1){
    
    if(D[i, j] == "dia") {
      align1 = c(rownames(M)[i], align1)
      align2 = c(colnames(M)[j], align2)
      j = j - 1; i = i - 1  # update indices
    } else if (D[i, j] == "ver") {
      align1 = c(rownames(M)[i], align1)
      align2 = c("-", align2) # vertical movement = gap for seq2
      i = i - 1 # update indices
    } else if (D[i, j] == "hor") {
      align1 = c("-", align1) # horizontal movement = gap for seq1
      align2 = c(colnames(M)[j], align2) 
      j = j - 1 # update indices
    } 
    
  }
  
  # Prepare output
  return(list(aligned_seqs = matrix(c(align1, align2), byrow = TRUE, nrow = 2),
              score = M[nrow(M), ncol(M)], score_matrix = M, movement_matrix = D))
 
}

# Pour que TestThat fonctionne :
# 1 sudo apt install r-base libuv1-dev
# 2 sudo R
# install.packages("testthat")

# Test case 1: Wiki example (https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ACGTC", "AGTC", gap = -1, mismatch = -1, match = 0)$aligned_seqs
  expected = matrix(c("A", "C", "G", "T", "C", 
                      "A", "-", "G", "T", "C"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 2: Problem 2.9 (M. Borodovsky, S. Ekisheva., 2006, 
# Problems and Solutions in Biological Sequence Analysis, Cambridge University Press)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("GAATTC", "GATTA", gap = -2, mismatch = -1, match = 2)$aligned_seqs
  expected = matrix(c("G", "A", "A", "T", "T", "C",
                      "G", "-", "A", "T", "T", "A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 3: Fig. 5-7 (J. Momand, A. McCurdy., 2017, Concepts in Bioinformatics  
# and Genomics, Oxford Uni- versity Press.)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ADCDNRCKCRWP", "AWCNDRQCLCRP", 
                     gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","D","C","D","N","-","R","-","C","K","C","R","W","P",
                      "A","W","C","-","N","D","R","Q","C","L","C","R","-","P"), 
                    byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 4: Ch. 5, Problem 5 MODIFIED (!) (J. Momand, A. McCurdy., 2017, Concepts in 
# Bioinformatics and Genomics, Oxford Uni- versity Press.)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ATAGC", "ATATGA", gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","T","A","-","G","C" ,
                      "A","T","A","T","G","A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

Ce programme fonctionne convenablement. Je veux maintenant introduire de la modularité, en séparant du programme principal les fonctions d’initialisation des matrices, regroupées dans le sous-programme Matrices-Init.R :

Matrices = function(len1, len2, seq1, seq2, gap) {

# Initialize matrix M (for scores)
  M = matrix(0, nrow = len1 + 1, ncol = len2 + 1) # Initialize matrix
  rownames(M) = c("-", seq1) # assign seq chars to matrix names
  colnames(M) = c("-", seq2) # assign seq chars to matrix names
  M[1, ] = cumsum(c(0, rep(gap, len2))) # Fill 1st row with gap penalites
  M[, 1] = cumsum(c(0, rep(gap, len1))) # Fill 1st col with gap penalites
  
  # Initialize matrix D (for directions)
  D = matrix(0, nrow = len1 + 1, ncol = len2 + 1) # Initialize matrix
  rownames(D) = c("-", seq1) # assign seq chars to matrix names
  colnames(D) = c("-", seq2) # assign seq chars to matrix names
  D[1, ] = rep("hor") # Fill 1st row with "hor" for horizontal moves
  D[, 1] = rep("ver") # Fill 1st col with "ver" for vertical moves

  return(list(score_matrix = M, movement_matrix = D))

}

Le programme précédemment monolithique devient :

# References: https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus

NeedlemanWunsch = function(seq1, seq2, gap, mismatch, match){
  
  # Stop conditions
  stopifnot(gap <= 0) # check if penalty negative
  stopifnot(mismatch <= 0)  # check if penalty negative
  stopifnot(match >= 0)  # check if score positive
  
  # Initialize col and rownames for matrices
  len1 = nchar(seq1); len2 = nchar(seq2) # Save number of chars in each sequence
  seq1 = unlist(strsplit(seq1, split = "")) # convert seq to character vector
  seq2 = unlist(strsplit(seq2, split = "")) # convert seq to character vector

source("Matrices-Init.R")
    resultM <- Matrices(len1, len2, seq1, seq2, gap)
    M <- resultM$score_matrix
    D <- resultM$movement_matrix

  # Compute scores and save moves
    type = c("dia", "hor", "ver") # Lookup vector
  for (i in 2:(len1 + 1)){# for every (initially zero) row
    for (j in 2:(len2 + 1)){# for every (initially zero) col
      hor = M[i, j - 1] + gap # horizontal move = gap for seq1
      ver = M[i - 1, j] + gap # vertical move = gap for seq2 
      dia = ifelse(rownames(M)[i] == colnames(M)[j], # diagonal = ifelse(chars equal, match, mismatch) 
                   M[i - 1, j - 1] + match, 
                   M[i - 1, j - 1] + mismatch)
      M[i, j] = max(dia, hor, ver) # Save current (best) score in M
      D[i, j] = type[which.max(c(dia, hor, ver))] # Save direction of move in D
    }
  } 
    
  # Backtracing
  align1 = c(); align2 = c() # Note: length of final alignments is unknown at this point
  
  while(i > 1 && j > 1){
    
    if(D[i, j] == "dia") {
      align1 = c(rownames(M)[i], align1)
      align2 = c(colnames(M)[j], align2)
      j = j - 1; i = i - 1  # update indices
    } else if (D[i, j] == "ver") {
      align1 = c(rownames(M)[i], align1)
      align2 = c("-", align2) # vertical movement = gap for seq2
      i = i - 1 # update indices
    } else if (D[i, j] == "hor") {
      align1 = c("-", align1) # horizontal movement = gap for seq1
      align2 = c(colnames(M)[j], align2) 
      j = j - 1 # update indices
    } 
    
  }
  
  # Prepare output
  return(list(aligned_seqs = matrix(c(align1, align2), byrow = TRUE, nrow = 2),
              score = M[nrow(M), ncol(M)], score_matrix = M, movement_matrix = D))
 
}

# Pour que TestThat fonctionne :
# 1 sudo apt install r-base libuv1-dev
# 2 sudo R
# install.packages("testthat")

# Test case 1: Wiki example (https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ACGTC", "AGTC", gap = -1, mismatch = -1, match = 0)$aligned_seqs
  expected = matrix(c("A", "C", "G", "T", "C", 
                      "A", "-", "G", "T", "C"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 2: Problem 2.9 (M. Borodovsky, S. Ekisheva., 2006, 
# Problems and Solutions in Biological Sequence Analysis, Cambridge University Press)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("GAATTC", "GATTA", gap = -2, mismatch = -1, match = 2)$aligned_seqs
  expected = matrix(c("G", "A", "A", "T", "T", "C",
                      "G", "-", "A", "T", "T", "A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 3: Fig. 5-7 (J. Momand, A. McCurdy., 2017, Concepts in Bioinformatics  
# and Genomics, Oxford Uni- versity Press.)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ADCDNRCKCRWP", "AWCNDRQCLCRP", 
                     gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","D","C","D","N","-","R","-","C","K","C","R","W","P",
                      "A","W","C","-","N","D","R","Q","C","L","C","R","-","P"), 
                    byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 4: Ch. 5, Problem 5 MODIFIED (!) (J. Momand, A. McCurdy., 2017, Concepts in 
# Bioinformatics and Genomics, Oxford Uni- versity Press.)
testthat::test_that("needlemanreturns correct alignment", {
  solution = NeedlemanWunsch("ATAGC", "ATATGA", gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","T","A","-","G","C" ,
                      "A","T","A","T","G","A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

Jusque là cela fonctionne.

Ensuite je veux isoler la fonction de calcul des scores et des décalages, dans la fonction ScoreMoves :

ScoresMoves = function(M, D, len1, len2, gap, match, mismatch) {

    type = c("dia", "hor", "ver") # Lookup vector

# Compute scores and save moves
  for (i in 2:(len1 + 1)){# for every (initially zero) row
    for (j in 2:(len2 + 1)){# for every (initially zero) col
      hor = M[i, j - 1] + gap # horizontal move = gap for seq1
      ver = M[i - 1, j] + gap # vertical move = gap for seq2 
      dia = ifelse(rownames(M)[i] == colnames(M)[j], # diagonal = ifelse(chars equal, match, mismatch) 
                   M[i - 1, j - 1] + match, 
                   M[i - 1, j - 1] + mismatch)
      M[i, j] = max(dia, hor, ver) # Save current (best) score in M
      D[i, j] = type[which.max(c(dia, hor, ver))] # Save direction of move in D 
    }
  }
  
}

Le programme devient alors :

# References: https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus

NeedlemanWunsch = function(seq1, seq2, gap, mismatch, match){
  
  # Stop conditions
  stopifnot(gap <= 0) # check if penalty negative
  stopifnot(mismatch <= 0)  # check if penalty negative
  stopifnot(match >= 0)  # check if score positive
  
  # Initialize col and rownames for matrices
  len1 = nchar(seq1); len2 = nchar(seq2) # Save number of chars in each sequence
  seq1 = unlist(strsplit(seq1, split = "")) # convert seq to character vector
  seq2 = unlist(strsplit(seq2, split = "")) # convert seq to character vector

source("Matrices-Init.R")
    resultM <- Matrices(len1, len2, seq1, seq2, gap)
    M <- resultM$score_matrix
    D <- resultM$movement_matrix

source("ScoresMoves.R")
    ScoresMoves(M, D, len1, len2, gap, match, mismatch)
    
  # Backtracing
  align1 = c(); align2 = c() # Note: length of final alignments is unknown at this point
  
  while(i > 1 && j > 1){
    
    if(D[i, j] == "dia") {
      align1 = c(rownames(M)[i], align1)
      align2 = c(colnames(M)[j], align2)
      j = j - 1; i = i - 1  # update indices
    } else if (D[i, j] == "ver") {
      align1 = c(rownames(M)[i], align1)
      align2 = c("-", align2) # vertical movement = gap for seq2
      i = i - 1 # update indices
    } else if (D[i, j] == "hor") {
      align1 = c("-", align1) # horizontal movement = gap for seq1
      align2 = c(colnames(M)[j], align2) 
      j = j - 1 # update indices
    } 
    
  }
  
  # Prepare output
  return(list(aligned_seqs = matrix(c(align1, align2), byrow = TRUE, nrow = 2),
              score = M[nrow(M), ncol(M)], score_matrix = M, movement_matrix = D))
 
}

# Pour que TestThat fonctionne :
# 1 sudo apt install r-base libuv1-dev
# 2 sudo R
# install.packages("testthat")

# Test case 1: Wiki example (https://de.wikipedia.org/wiki/Needleman-Wunsch-Algorithmus)
testthat::test_that("NeedlemanWunsch returns correct alignment", {
  solution = NeedlemanWunsch("ACGTC", "AGTC", gap = -1, mismatch = -1, match = 0)$aligned_seqs
  expected = matrix(c("A", "C", "G", "T", "C", 
                      "A", "-", "G", "T", "C"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 2: Problem 2.9 (M. Borodovsky, S. Ekisheva., 2006, 
# Problems and Solutions in Biological Sequence Analysis, Cambridge University Press)
testthat::test_that("NeedlemanWunsch returns correct alignment", {
  solution = NeedlemanWunsch("GAATTC", "GATTA", gap = -2, mismatch = -1, match = 2)$aligned_seqs
  expected = matrix(c("G", "A", "A", "T", "T", "C",
                      "G", "-", "A", "T", "T", "A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 3: Fig. 5-7 (J. Momand, A. McCurdy., 2017, Concepts in Bioinformatics  
# and Genomics, Oxford University Press.)
testthat::test_that("NeedlemanWunsch returns correct alignment", {
  solution = NeedlemanWunsch("ADCDNRCKCRWP", "AWCNDRQCLCRP", 
                     gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","D","C","D","N","-","R","-","C","K","C","R","W","P",
                      "A","W","C","-","N","D","R","Q","C","L","C","R","-","P"), 
                    byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

# Test case 4: Ch. 5, Problem 5 MODIFIED (!) (J. Momand, A. McCurdy., 2017, Concepts in 
# Bioinformatics and Genomics, Oxford University Press.)
testthat::test_that("NeedlemanWunsch returns correct alignment", {
  solution = NeedlemanWunsch("ATAGC", "ATATGA", gap = 0, mismatch = 0, match = 1)$aligned_seqs
  expected = matrix(c("A","T","A","-","G","C" ,
                      "A","T","A","T","G","A"), byrow = TRUE, nrow =2)
  testthat::expect_equal(solution, expected)
})

Et là, pour une raison que je ne réussis pas à élucider, le programme appelé ainsi :

./NW.R sequences-1.txt -1 -1 1

échoue avec les messages suivants :

── Error: NeedlemanWunsch returns correct alignment ────────────────────────────
Error in `NeedlemanWunsch("ACGTC", "AGTC", gap = -1, mismatch = -1, match = 0)`: objet 'i' introuvable
Backtrace:
    ▆
 1. └─global NeedlemanWunsch("ACGTC", "AGTC", gap = -1, mismatch = -1, match = 0)
Error:
! Test failed with 1 failure and 0 successes.
Backtrace:
     ▆
  1. ├─base::source("./Needleman-Wunsch-LB.R")
  2. │ ├─base::withVisible(eval(ei, envir))
  3. │ └─base::eval(ei, envir)
  4. │   └─base::eval(ei, envir)
  5. └─testthat::test_that(...)
  6.   └─testthat:::test_code(code, parent.frame())
  7.     ├─withr (local) `<fn>`()
  8.     └─reporter$end_test(context = reporter$.context, test = test)
  9.       └─cli::cli_abort(...)
 10.         └─rlang::abort(...)
Exécution arrêtée

Je reste perplexe devant cet “objet ’i’ introuvable”. Aurais-tu une idée ?

Merci d’avance !

Réponse :

Bonjour Laurent,

Oui : l’erreur « objet ’i’ introuvable » est parfaitement explicable, et elle révèle en fait deux petits problèmes liés à la modularisation en R.

Le premier est celui qui provoque immédiatement votre erreur.

1. Où est passé i ?

Dans votre programme monolithique, vous aviez :

for (i in 2:(len1 + 1)){
  for (j in 2:(len2 + 1)){
    ...
  }
}

while(i > 1 && j > 1){
    ...
}

En R, après l’exécution des for, les variables i et j existent encore dans le même environnement. Vous pouviez donc réutiliser i et j dans le while.

Mais vous avez maintenant déplacé les boucles dans :

ScoresMoves = function(M, D, len1, len2, gap, match, mismatch) {
    ...
    for (i in 2:(len1 + 1)){
        for (j in 2:(len2 + 1)){
            ...
        }
    }
}

i et j sont alors des variables locales à la fonction ScoresMoves.

Lorsque ScoresMoves() se termine, elles ne sont plus accessibles dans NeedlemanWunsch().

C’est exactement ce que vous dit R :

objet 'i' introuvable

Le while rencontre :

while(i > 1 && j > 1){

et cherche i dans l’environnement de NeedlemanWunsch(). Or cette variable n’y existe plus.

2. Il y a toutefois un deuxième problème

Même si vous ajoutiez simplement :

i = len1 + 1
j = len2 + 1

avant le while, votre programme ne fonctionnerait toujours pas correctement.

Pourquoi ?

Parce que vous faites :

ScoresMoves(M, D, len1, len2, gap, match, mismatch)

mais vous ne récupérez aucun résultat.

Dans ScoresMoves, vous modifiez :

M[i, j] = ...
D[i, j] = ...

mais la bonne pratique en R est de considérer qu’une fonction doit retourner les objets qu’elle a produits ou modifiés, plutôt que de compter sur une modification implicite des objets du programme appelant.

Vous aviez très bien fait cela avec Matrices :

return(list(score_matrix = M, movement_matrix = D))

Il faut faire la même chose avec ScoresMoves.

Je vous proposerais donc cette version pour ScoresMoves.R :

ScoresMoves = function(M, D, len1, len2, gap, match, mismatch) {

  type = c("dia", "hor", "ver")

  # Compute scores and save moves
  for (i in 2:(len1 + 1)) {
    for (j in 2:(len2 + 1)) {

      hor = M[i, j - 1] + gap
      ver = M[i - 1, j] + gap

      dia = ifelse(
        rownames(M)[i] == colnames(M)[j],
        M[i - 1, j - 1] + match,
        M[i - 1, j - 1] + mismatch
      )

      M[i, j] = max(dia, hor, ver)
      D[i, j] = type[which.max(c(dia, hor, ver))]
    }
  }

  return(list(
    score_matrix = M,
    movement_matrix = D
  ))
}

Puis dans NeedlemanWunsch :

source("ScoresMoves.R")

resultS <- ScoresMoves(
  M, D, len1, len2,
  gap, match, mismatch
)

M <- resultS$score_matrix
D <- resultS$movement_matrix

Ainsi, vous avez une architecture assez propre :

NeedlemanWunsch
       │
       ├── Matrices()
       │      └── produit M et D initialisées
       │
       ├── ScoresMoves()
       │      └── produit M et D calculées
       │
       └── Backtracking
              └── produit l'alignement

3. Et il faut initialiser i et j avant le backtracking

C’est un point qui était un peu caché dans votre version monolithique.

À la fin des boucles :

for (i in 2:(len1 + 1)){
    for (j in 2:(len2 + 1)){
        ...
    }
}

vous aviez implicitement :

i = len1 + 1
j = len2 + 1

à cause du fonctionnement de for.

Mais maintenant, ce n’est plus le cas.

Il faut donc écrire explicitement :

i = nrow(M)
j = ncol(M)

avant le retour sur trace (backtracking).

4. Votre programme principal devient donc :

La partie essentielle pourrait être :

NeedlemanWunsch = function(seq1, seq2, gap, mismatch, match){

  # Stop conditions
  stopifnot(gap <= 0)
  stopifnot(mismatch <= 0)
  stopifnot(match >= 0)

  # Sequence lengths
  len1 = nchar(seq1)
  len2 = nchar(seq2)

  # Convert sequences to character vectors
  seq1 = unlist(strsplit(seq1, split = ""))
  seq2 = unlist(strsplit(seq2, split = ""))

  # Initialize matrices
  source("Matrices-Init.R")

  resultM <- Matrices(len1, len2, seq1, seq2, gap)

  M <- resultM$score_matrix
  D <- resultM$movement_matrix

  # Compute scores and moves
  source("ScoresMoves.R")

  resultS <- ScoresMoves(
    M, D, len1, len2,
    gap, match, mismatch
  )

  M <- resultS$score_matrix
  D <- resultS$movement_matrix

  # Backtracing
  align1 = c()
  align2 = c()

  i = nrow(M)
  j = ncol(M)

  while(i > 1 && j > 1){

    if(D[i, j] == "dia") {

      align1 = c(rownames(M)[i], align1)
      align2 = c(colnames(M)[j], align2)

      i = i - 1
      j = j - 1

    } else if(D[i, j] == "ver") {

      align1 = c(rownames(M)[i], align1)
      align2 = c("-", align2)

      i = i - 1

    } else if(D[i, j] == "hor") {

      align1 = c("-", align1)
      align2 = c(colnames(M)[j], align2)

      j = j - 1
    }
  }

  # Prepare output
  return(list(
    aligned_seqs = matrix(
      c(align1, align2),
      byrow = TRUE,
      nrow = 2
    ),
    score = M[nrow(M), ncol(M)],
    score_matrix = M,
    movement_matrix = D
  ))
}

Avec cette modification, vos tests devraient à nouveau passer.

5. Une petite subtilité supplémentaire dans votre while :

Il y a néanmoins un autre détail intéressant dans votre code initial :

while(i > 1 && j > 1)

Cela signifie : tant que les deux séquences ont encore des caractères à traiter.

Mais imaginez qu’au cours du backtracking vous arriviez ici :

    - A C G T C
    A - G T C

et que j soit arrivé à 1 alors que i > 1.

Il reste alors des caractères dans la première séquence qui doivent nécessairement être alignés avec des gaps dans la seconde.

Avec :

while(i > 1 && j > 1)

vous sortez immédiatement de la boucle et abandonnez ces caractères.

Pour un Needleman-Wunsch général, il est donc plus correct de faire :

while(i > 1 || j > 1)

puis de traiter les trois cas.

On peut par exemple écrire :

while(i > 1 || j > 1) {

  if(i > 1 && j > 1 && D[i, j] == "dia") {

    align1 = c(rownames(M)[i], align1)
    align2 = c(colnames(M)[j], align2)

    i = i - 1
    j = j - 1

  } else if(i > 1 && D[i, j] == "ver") {

    align1 = c(rownames(M)[i], align1)
    align2 = c("-", align2)

    i = i - 1

  } else if(j > 1 && D[i, j] == "hor") {

    align1 = c("-", align1)
    align2 = c(colnames(M)[j], align2)

    j = j - 1
  }
}

C’est un peu plus robuste.
En résumé

Votre erreur actuelle vient très précisément de la modularisation :

Avant :

NeedlemanWunsch()
 ├── for(i...)
 ├── for(j...)
 └── while(i...)

i et j appartiennent au même environnement.

Maintenant :

NeedlemanWunsch()
 ├── ScoresMoves()
 │    ├── for(i...)
 │    └── for(j...)
 │
 └── while(i...)

i et j appartiennent à l’environnement local de ScoresMoves() et ne sont donc pas disponibles dans NeedlemanWunsch().

Et, indépendamment de cela, ScoresMoves() doit retourner M et D, puisque c’est elle qui les calcule.

C’est d’ailleurs une bonne illustration d’un principe assez important en R : une fonction gagne généralement à recevoir ses données en arguments et à retourner explicitement son résultat, plutôt que de modifier silencieusement l’état du programme appelant.

Je pense que votre découpage Matrices → ScoresMoves → Backtracking est tout à fait pertinent pour la suite de votre développement.

Et votre exemple est particulièrement instructif : en passant d’un programme monolithique à trois fonctions, vous venez de rencontrer l’un des points fondamentaux de R : les environnements et la portée des variables (scope).

Votre i « disparu » est en réalité très pédagogique : il n’a pas disparu, il appartient à l’environnement local de ScoresMoves(), et cet environnement cesse d’être accessible de la même manière lorsque la fonction retourne.

C’est une différence assez importante avec ce que l’on peut spontanément attendre lorsqu’on vient de langages comme C, Python ou JavaScript.

Bon courage pour la suite de votre plongée dans R et la bio-informatique !