如何为DNA序列生成独热编码？

Question

我想为一组DNA序列生成独热编码。例如，序列ACGTCCA可以以下面的方式表示。但是下面的代码将以水平方式生成独热编码，而我更喜欢以垂直形式呈现。有人可以帮助我吗？

ACGTCCA 
1000001 - A
0100110 - C 
0010000 - G
0001000 - T

示例代码：

from sklearn.preprocessing import OneHotEncoder
import itertools

# two example sequences
seqs = ["ACGTCCA","CGGATTG"]


# split sequences to tokens
tokens_seqs = [seq.split("\\") for seq in seqs]

# convert list of of token-lists to one flat list of tokens
# and then create a dictionary that maps word to id of word,
# like {A: 1, B: 2} here
all_tokens = itertools.chain.from_iterable(tokens_seqs)
word_to_id = {token: idx for idx, token in enumerate(set(all_tokens))}

# convert token lists to token-id lists, e.g. [[1, 2], [2, 2]] here
token_ids = [[word_to_id[token] for token in tokens_seq] for tokens_seq in tokens_seqs]

# convert list of token-id lists to one-hot representation
vec = OneHotEncoder(n_values=len(word_to_id))
X = vec.fit_transform(token_ids)

print X.toarray()

然而，这段代码给我输出：

[[ 0.  1.]
 [ 1.  0.]]

预期输出:

[[1. 0. 0. 0. 0. 0. 1. 0. 1. 0. 0. 1. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 1. 0. 0. 0. 0. 1. 1. 0.]]

- Xiong89

Answer 1

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")

## Output: 
array([[1., 0., 0., 0.],
   [1., 0., 0., 0.],
   [0., 1., 0., 0.],
   [0., 0., 1., 0.],
   [0., 0., 0., 1.]])

- DrIDK

Answer 2

我建议用稍微手动一些的方式来完成：

import numpy as np

seqs = ["ACGTCCA","CGGATTG"]

CHARS = 'ACGT'
CHARS_COUNT = len(CHARS)

maxlen = max(map(len, seqs))
res = np.zeros((len(seqs), CHARS_COUNT * maxlen), dtype=np.uint8)

for si, seq in enumerate(seqs):
    seqlen = len(seq)
    arr = np.chararray((seqlen,), buffer=seq)
    for ii, char in enumerate(CHARS):
        res[si][ii*seqlen:(ii+1)*seqlen][arr == char] = 1

print res

这将会给您想要的结果：

[[1 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0]
 [0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1 1 0]]

- John Zwinck

---

谢谢。这个解决方案完美地适用于单个碱基A、T、C、G。我想知道是否可能将其用于二核苷酸，其有16种可能的组合，如AA、AT、AC、AG等。如果它出现为TGA，则TG将被赋值为1，GA也将被赋值为1。 - Xiong89

@Xiong89：我相信你可以做到，但如果我要为你编写所有的代码，我可能会要求按小时计费！我认为你需要用这个主要的构建块来处理di case：arr = np.chararray((3, 2), buffer='asdfds') 然后 np.all(np.chararray((2,), buffer='as') == arr, axis=1) - 这将给你一个数组 [True, False, False]。也就是说，你需要将所有16个两个字符组合作为chararray构建，然后与之匹配。试试看吧。 :) - John Zwinck

Answer 3

from keras.utils import to_categorical

def one_hot_encoding(seq):
    mp = dict(zip('ACGT', range(4)))
    seq_2_number = [mp[nucleotide] for nucleotide in seq]
    return to_categorical(seq_2_number, num_classes=4, dtype='int32').flatten()

- Rajan saha Raju