如何使用 R 或 python 为 DNA 序列生成一种热编码

Question

我想为一系列 DNA 序列生成一个热编码矩阵。我尝试通过以下链接解决我的问题： 如何为 DNA 序列生成一种热编码？但有些解决方案仅针对单个 DNA 序列给出，而不是针对一系列 DNA 序列。

例如

def one_hot_encode(seq):
    mapping = dict(zip("ACGT", range(4)))    
    seq2 = [mapping[i] for i in seq]
    return np.eye(4)[seq2]

one_hot_encode("AACGT")

在上面给出的代码中，如果我运行one_hot_encode("AACGT","GGTAC","CGTAC")它将失败，我也想生成矩阵作为输出。

目前，我正在 R 中工作，下面是我在 R 数据框中的 DNA 序列（单列文件）

ACTTTA
TTGATG
CTTACG
GTACGT

预期产出

1   0   0   0   0   1   0   0   0   0   0   1   0   0   0   1   0   0   0   1   1   0   0   0
0   0   0   1   0   0   0   1   0   0   1   0   1   0   0   0   0   0   0   1   0   0   1   0
0   1   0   0   0   0   0   1   0   0   0   1   1   0   0   0   0   1   0   0   0   0   1   0
0   0   1   0   0   0   0   1   1   0   0   0   0   1   0   0   0   0   1   0   0   0   0   1

在 R 中可以做到这一点吗？

Answer 1

Here\xe2\x80\x99s 是一个使用基本 R 的解决方案，可生成解决方案的转置。也就是说，它为每个单独的字符创建一个 one-hot 列向量，并将这些列连接起来（即，无论字符串的数量如何，总是有四行）。

\n

（有关生成与问题相同格式的解决方案，请参阅底部。）

\n

sequences = c(\'ACTTTA\', \'TTGATG\', \'GATTACA\')\n

\n

strsplit(sequences, \'\') |>\n  lapply(match, table = c(\'A\', \'C\', \'G\', \'T\')) |>\n  lapply(\\(x) {\n    m = diag(0L, nrow = 4L, ncol = length(x))\n    m[cbind(x, seq_along(x))] = 1L\n    m\n  }) |>\n  do.call(\'cbind\', args = _)\n

\n

或者，另一种选择（更短但不一定更具可读性）：

\n

strsplit(sequences, \'\') |>\n  lapply(match, table = c(\'A\', \'C\', \'G\', \'T\')) |>\n  unlist() %>%\n  {diag(1L, 4L)[, .]}\n

\n

不幸的是，这需要 \xe2\x80\x98magrittr\xe2\x80\x99 管道，但我们可以通过将最后一行的表达式抽象为函数来更改它（并使其再次更具可读性）：

\n

num_to_one_hot = function (x, bits) {\n  diag(1L, bits)[, x]\n}\n

\n

这是一个很好的通用函数，可以对数字向量进行一次性编码。此外，我们可以unlist在第一步之后，这使我们能够避免lapply. 这样我们的 DNA 编码代码就变成了：

\n

strsplit(sequences, \'\') |>\n  unlist() |>\n  match(c(\'A\', \'C\', \'G\', \'T\')) |>\n  num_to_one_hot(bits = 4L)\n

\n

\xe2\x80\xa6 是迄今为止所有替代方案中最不言自明（也是最易读的！）的。它也完全矢量化，并且不使用lapply或类似的，因此它也更高效。

\n\n

为了完整起见，这里\xe2\x80\x99s是一个解决方案，通过将先前算法生成的矩阵从列优先方向转换为行优先方向，产生与问题中所要求的相同的输出：

\n

strsplit(sequences, \'\') |>\n  unlist() |>\n  match(c(\'A\', \'C\', \'G\', \'T\')) |>\n  num_to_one_hot(4L) |>\n  matrix(nrow = length(sequences), byrow = TRUE)\n

\n

num_to_one_hot我对当前发布的所有解决方案进行了基准测试，使用（Konrad2及以下）的解决方案Konrad3是最快的：

\n

> bench::mark(Ma\xc3\xabl(), Paul(), Konrad1(), Konrad2(), Konrad3(), GKi1(), GKi2(), Thomas(), check = FALSE)\n# # A tibble: 8 \xc3\x97 13\n#   expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory     time       gc      \n#   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>     <list>     <list>  \n# 1 Ma\xc3\xabl()        776\xc2\xb5s 818.25\xc2\xb5s     1154.    1.36KB    13.3    522     6    452.2ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 2 Paul()       1.07ms   1.17ms      743.    4.38KB     2.04   365     1    491.2ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 3 Konrad1()    21.7\xc2\xb5s   23.7\xc2\xb5s    39401.      432B    15.8   9996     4    253.7ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 4 Konrad2()     4.8\xc2\xb5s    5.7\xc2\xb5s   148619.    1.69KB    14.9   9999     1     67.3ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 5 Konrad3()     6.6\xc2\xb5s    7.9\xc2\xb5s   108540.    2.11KB    10.9   9999     1     92.1ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 6 GKi1()        9.2\xc2\xb5s   10.8\xc2\xb5s    83596.      960B    16.7   9998     2    119.6ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 7 GKi2()      258.3\xc2\xb5s    278\xc2\xb5s     3442.      960B     0     1721     0      500ms <NULL> <Rprofmem> <bench_tm> <tibble>\n# 8 Thomas()       10\xc2\xb5s   11.6\xc2\xb5s    77604.    2.39KB     7.76  9999     1    128.8ms <NULL> <Rprofmem> <bench_tm> <tibble>\n

\n