Добрый день, коллеги. Здесь же можно задавать вопросы о качестве

кода в плане восприятия/оптимальности и вообще стилистики? Поскольку задача решена, но нет уверенности, что так делать хорошо и правильно.
Суть в следующем: есть две последовательности ДНК (строки), в которых есть различия на определенных позициях. Последовательности должны быть одной длины и моя задача определить, в каких позициях произошли замены. Однако, поскольку это ДНК, то в зависимости от позиции относительно рамки считывания, будет различаться конечный результат замены. Нуклеотиды считываются по 3, поэтому есть три возможных ситуации: замена в первой позиции, во второй и в третьей относительно старта. В своей реализации использую case_when для проверки остатка от деления на три и выписываю результат.

Вопрос, насколько это адекватно и понятно?
r
library(purrr)
library(stringr)
library(dplyr)

ref_seq <- 'ATGACCCCAATACGCAA'
mut_seq <- 'ATGTCCGCAATGCGCAA'

ref_seq <- unlist(strsplit(ref_seq, ''))
mut_seq <- unlist(strsplit(mut_seq, ''))

res <- map2_chr(ref_seq, mut_seq, .f = ~ifelse(.x == .y, "", paste(.x, .y, sep = "_to_")))
mutation_indeces <- str_which(res, '')
gene_code <- list('UUU' = 'Phe', 'UUC'='Phe', 'UCU'='Ser', 'UCC'='Ser',
'UAU'='Tyr', 'UAC'='Tyr', 'UGU'='Cys', 'UGC'='Cys',
'UUA'='Leu', 'UCA'='Ser', 'UAA'='stop', 'UGA'='stop',
'UUG'='Leu', 'UCG'='Ser', 'UAG'='stop', 'UGG'='Trp',
'CUU'='Leu', 'CUC'='Leu', 'CCU'='Pro', 'CCC'='Pro',
'CAU'='His', 'CAC'='His', 'CGU'='Arg', 'CGC'='Arg',
'CUA'='Leu', 'CUG'='Leu', 'CCA'='Pro', 'CCG'='Pro',
'CAA'='Gln', 'CAG'='Gln', 'CGA'='Arg', 'CGG'='Arg',
'AUU'='Ile', 'AUC'='Ile', 'ACU'='Thr', 'ACC'='Thr',
'AAU'='Asn', 'AAC'='Asn', 'AGU'='Ser', 'AGC'='Ser',
'AUA'='Ile', 'ACA'='Thr', 'AAA'='Lys', 'AGA'='Arg',
'AUG'='Met', 'ACG'='Thr', 'AAG'='Lys', 'AGG'='Arg',
'GUU'='Val', 'GUC'='Val', 'GCU'='Ala', 'GCC'='Ala',
'GAU'='Asp', 'GAC'='Asp', 'GGU'='Gly', 'GGC'='Gly',
'GUA'='Val', 'GUG'='Val', 'GCA'='Ala', 'GCG'='Ala',
'GAA'='Glu', 'GAG'='Glu', 'GGA'='Gly', 'GGG'='Gly')
translate_codon <- function(start_coord, end_coord, string) {
gene_code[[gsub('T', 'U', paste0(string[start_coord:end_coord], collapse = ''))]]
}
map(mutation_indeces, .f = ~case_when(
.x %% 3 == 1 ~ c(res[.x], .x, translate_codon(.x, (.x+2), ref_seq),
translate_codon(.x, (.x+2), mut_seq)),
.x %% 3 == 2 ~ c(res[.x], .x, translate_codon((.x-1), (.x+1), ref_seq),
translate_codon((.x-1), (.x+1), mut_seq)),
.x %% 3 == 0 ~ c(res[.x], .x, translate_codon((.x-2), (.x), ref_seq),
translate_codon((.x-2), (.x), mut_seq))))
#> [[1]]
#> [1] "A_to_T" "4" "Thr" "Ser"
#>
#> [[2]]
#> [1] "C_to_G" "7" "Pro" "Ala"
#>
#> [[3]]
#> [1] "A_to_G" "12" "Ile" "Met"


^{Created on 2022-05-18 by the reprex package (v2.0.1)}

Можно попробовать подключить Biostrings. Там большинство стандартных задач для работы с биологическими последовательностями решаются (та же трансляция). Уж генетический код вручную точно можно не задавать.

один минус — подключение биокондуктора сносит весь устоявшийся набор пакетов и глаголов в небытие. там полностью свой мир.

Это да, и минус не единственный, конечно.

Глянул. Я бы написал примерно так. Атомарные верифицируемые шаги. Исходный вариант хорош для демонстрации владения различными методами. Но там отсутствует валидация и многие шаги тяжело проследить без дебага или знания нюансов функций, таких как "An empty pattern, '', is equivalent to boundary('character')." в str_which, а также избыточности под капотом. library(data.table) library(magrittr) ref_seq <- 'ATGACCCCAATACGCAA' mut_seq <- 'ATGTCCGCAATGCGCAA' gene_code <- list('UUU'='Phe', 'UUC'='Phe', 'UCU'='Ser', 'UCC'='Ser', 'UAU'='Tyr', 'UAC'='Tyr', 'UGU'='Cys', 'UGC'='Cys', 'UUA'='Leu', 'UCA'='Ser', 'UAA'='stop', 'UGA'='stop', 'UUG'='Leu', 'UCG'='Ser', 'UAG'='stop', 'UGG'='Trp', 'CUU'='Leu', 'CUC'='Leu', 'CCU'='Pro', 'CCC'='Pro', 'CAU'='His', 'CAC'='His', 'CGU'='Arg', 'CGC'='Arg', 'CUA'='Leu', 'CUG'='Leu', 'CCA'='Pro', 'CCG'='Pro', 'CAA'='Gln', 'CAG'='Gln', 'CGA'='Arg', 'CGG'='Arg', 'AUU'='Ile', 'AUC'='Ile', 'ACU'='Thr', 'ACC'='Thr', 'AAU'='Asn', 'AAC'='Asn', 'AGU'='Ser', 'AGC'='Ser', 'AUA'='Ile', 'ACA'='Thr', 'AAA'='Lys', 'AGA'='Arg', 'AUG'='Met', 'ACG'='Thr', 'AAG'='Lys', 'AGG'='Arg', 'GUU'='Val', 'GUC'='Val', 'GCU'='Ala', 'GCC'='Ala', 'GAU'='Asp', 'GAC'='Asp', 'GGU'='Gly', 'GGC'='Gly', 'GUA'='Val', 'GUG'='Val', 'GCA'='Ala', 'GCG'='Ala', 'GAA'='Glu', 'GAG'='Glu', 'GGA'='Gly', 'GGG'='Gly') translateCodon <- function(vec){ gene_code[[gsub('T', 'U', paste0(vec, collapse = ''))]] } dt <- data.table(ref = ref_seq, mut = mut_seq) %>% .[, lapply(.SD, function(x){unlist(strsplit(x, ""))})] %>% .[, `:=`(grp = (.I - 1) %/% 3, is_noneq = (ref != mut))] %>% .[, `:=`(n_diff = sum(is_noneq), n = .N, idx = .I), by = grp] %>% # санация неполных триплетов + оставляем только измененные наборы # вопрос про допустимое количество изменений! .[n == 3 & n_diff == 1] # соберем словарик генов по группам dict_dt <- dt[, .(ref_gene = translateCodon(ref), mut_gene = translateCodon(mut)), by = grp] # делаем свертку res_dt <- merge(dt[is_noneq == TRUE], dict_dt, by = "grp") glue::glue_data(res_dt, "{ref}->{mut}, @{idx}, {ref_gene}->{mut_gene}") #> A->T, @4, Thr->Ser #> C->G, @7, Pro->Ala #> A->G, @12, Ile->Met

Спасибо, решение действительно изящное, правда синтаксис дт пока еще не освоен, поэтому читалось сложно, но при применении шагов все достаточно понятно. Решение конечно отправлю свое, но на заметку взяла такой дататейбл подход