一种具有歧义的热编码(核苷酸序列)

One hot encoding with ambiguitiy (nucleotide sequences)

核苷酸序列(或 DNA 序列)通常由 4 个碱基组成:ATGC,这是一种非常好的、简单且有效的编码方式,用于机器学习目的。

sequence = AAATGCC
ohe_sequence = [[1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 0, 1]]

但是,当您考虑到 RNA 序列和测序仪中有时会出现的错误时,会添加字母 UYRWSKMDVHBXN...当对其进行单热编码时,您最终会得到一个包含 17 行的矩阵,其中最后一行13行一般都是0。

这是非常低效的,并且没有赋予这些额外的(不明确的)字母所具有的生物学意义。

例如:

所以我制作了一本代表这种生物学意义的字典

nucleotide_dict = {'A': 'A', 'T':'T', 'U':'T', 'G':'G', 'C':'C', 'Y':['C', 'T'], 'R':['A', 'G'], 'W':['A', 'T'], 'S':['G', 'C'], 'K':['T', 'G'], 'M':['C', 'A'], 'D':['A', 'T', 'G'], 'V':['A', 'G', 'C'], 'H':['A', 'T', 'C'], 'B':['T', 'G', 'C'], 'X':['A', 'T', 'G', 'C'], 'N':['A', 'T', 'G', 'C']}

但我似乎无法弄清楚如何制作一个有效的单热编码脚本(或者是否有一种方法可以使用 scikit 学习模块来做到这一点),它利用这本字典来获得结果是这样的:

sequence = ANTUYCC
ohe_sequence = [[1, 0, 0, 0], [1, 1, 1, 1], [0, 1, 0, 0], [0, 1, 0, 0], [0, 1, 0, 1], [0, 0, 0, 1], [0, 0, 0, 1]]

# or even better:
ohe_sequence = [[1, 0, 0, 0], [0.25, 0.25, 0.25, 0.25], [0, 1, 0, 0], [0, 1, 0, 0], [0, 0.5, 0, 0.5], [0, 0, 0, 1], [0, 0, 0, 1]]

我喜欢后一种方法,因为它更接近真实的含义: Y 并不表示 C T,而是 C T。如果没有进一步的信息可用,假设概率相等(即权重)似乎是合理的。当然,这种 non-standard 编码需要通过损失函数的选择来反映。

回答您的问题:您可以预先计算从字母到编码的映射,然后创建一个 encode 函数,该函数将 sequence 作为输入并 returns 编码序列作为 (len(sequence), 4)np.array 如下:

import numpy as np

nucleotide_dict = {'A':'A', 'T':'T', 'U':'T', 'G':'G', 'C':'C', 'Y':['C', 'T'], 'R':['A', 'G'], 'W':['A', 'T'], 'S':['G', 'C'], 'K':['T', 'G'], 'M':['C', 'A'], 'D':['A', 'T', 'G'], 'V':['A', 'G', 'C'], 'H':['A', 'T', 'C'], 'B':['T', 'G', 'C'], 'X':['A', 'T', 'G', 'C'], 'N':['A', 'T', 'G', 'C']}
index_mapper = {'A': 0, 'T': 1, 'G': 2, 'C': 3}
mapper_dict = dict()
for k, v in nucleotide_dict.items():
    encoding = np.zeros(4)
    p = 1 / len(v)
    encoding[[index_mapper[i] for i in v]] = p
    mapper_dict[k] = encoding

def encode(sequence):
    return np.array([mapper_dict[s] for s in sequence])

这似乎产生了预期的结果,并且可能有点高效。

一个例子:

print(encode('AYSDX'))

打印以下内容:

array([[1.        , 0.        , 0.        , 0.        ],
       [0.        , 0.5       , 0.        , 0.5       ],
       [0.        , 0.        , 0.5       , 0.5       ],
       [0.33333333, 0.33333333, 0.33333333, 0.        ],
       [0.25      , 0.25      , 0.25      , 0.25      ]])

这很有趣!我认为您可以使用具有适当值的字典来做到这一点。我添加了 scikit-learn class 因为你提到你正在使用它。请参阅下面的 transform

from sklearn.base import BaseEstimator, TransformerMixin
import numpy as np

nucleotide_dict = {
    "A": [1, 0, 0, 0],
    "G": [0, 1, 0, 0],
    "C": [0, 0, 1, 0],
    "T": [0, 0, 0, 1],
    "U": [0, 0, 0, 1],
    "Y": [0, 0, 1, 1],
    "R": [1, 1, 0, 0],
    "W": [1, 0, 0, 1],
    "S": [0, 1, 1, 0],
    "K": [0, 1, 0, 1],
    "M": [1, 0, 1, 0],
    "D": [1, 1, 0, 1],
    "V": [1, 1, 1, 0],
    "H": [1, 0, 1, 1],
    "B": [0, 1, 1, 1],
    "X": [1, 1, 1, 1],
    "N": [1, 1, 1, 1],
    "-": [0, 0, 0, 0],
}


class NucleotideEncoder(BaseEstimator, TransformerMixin):
    def __init__(self, norm=True):
        self.norm = norm

    def fit(self, X, y=None):
        return self

    def transform(self, X, y=None):
        f1 = lambda a: list(a)
        f2 = lambda g: nucleotide_dict[g]
        f3 = lambda c: list(map(f2, f1(c[0])))
        f4 = lambda t: np.array(f3(t)) / np.sum(np.array(f3(t)), axis=1)[:, np.newaxis]
        f = f3
        if self.norm:
            f = f4
        return np.apply_along_axis(f, 1, X)


samples = np.array([["AAATGCC"], ["ANTUYCC"]])
print(NucleotideEncoder().fit_transform(samples))

输出:

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

 [[1.   0.   0.   0.  ]
  [0.25 0.25 0.25 0.25]
  [0.   0.   0.   1.  ]
  [0.   0.   0.   1.  ]
  [0.   0.   0.5  0.5 ]
  [0.   0.   1.   0.  ]
  [0.   0.   1.   0.  ]]]

另一个允许不同长度序列的版本:

import numpy as np

nucleotide_dict = {
    "A": [1, 0, 0, 0],
    "G": [0, 1, 0, 0],
    "C": [0, 0, 1, 0],
    "T": [0, 0, 0, 1],
    "U": [0, 0, 0, 1],
    "Y": [0, 0, 1, 1],
    "R": [1, 1, 0, 0],
    "W": [1, 0, 0, 1],
    "S": [0, 1, 1, 0],
    "K": [0, 1, 0, 1],
    "M": [1, 0, 1, 0],
    "D": [1, 1, 0, 1],
    "V": [1, 1, 1, 0],
    "H": [1, 0, 1, 1],
    "B": [0, 1, 1, 1],
    "X": [1, 1, 1, 1],
    "N": [1, 1, 1, 1],
    "-": [0, 0, 0, 0],
}

norm = True
samples = np.array(["AAATGCC", "ANTUYCC", "".join(list(nucleotide_dict.keys()))[:-1]])


def nucleotide_encode(samples, norm=True):
    m = map(list, samples)
    f1 = lambda x: np.array(list(map(nucleotide_dict.get, x)))
    f = f1
    if norm:
        f = lambda x: np.nan_to_num(f1(x) / np.sum(f1(x), axis=1)[:, np.newaxis])
    return list(map(f, m))


for i, j in zip(samples, nucleotide_encode(samples, norm=norm)):
    print(i)
    print(j)

产量:

AAATGCC
[[1. 0. 0. 0.]
 [1. 0. 0. 0.]
 [1. 0. 0. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 1. 0.]]
ANTUYCC
[[1.   0.   0.   0.  ]
 [0.25 0.25 0.25 0.25]
 [0.   0.   0.   1.  ]
 [0.   0.   0.   1.  ]
 [0.   0.   0.5  0.5 ]
 [0.   0.   1.   0.  ]
 [0.   0.   1.   0.  ]]
AGCTUYRWSKMDVHBXN
[[1.         0.         0.         0.        ]
 [0.         1.         0.         0.        ]
 [0.         0.         1.         0.        ]
 [0.         0.         0.         1.        ]
 [0.         0.         0.         1.        ]
 [0.         0.         0.5        0.5       ]
 [0.5        0.5        0.         0.        ]
 [0.5        0.         0.         0.5       ]
 [0.         0.5        0.5        0.        ]
 [0.         0.5        0.         0.5       ]
 [0.5        0.         0.5        0.        ]
 [0.33333333 0.33333333 0.         0.33333333]
 [0.33333333 0.33333333 0.33333333 0.        ]
 [0.33333333 0.         0.33333333 0.33333333]
 [0.         0.33333333 0.33333333 0.33333333]
 [0.25       0.25       0.25       0.25      ]
 [0.25       0.25       0.25       0.25      ]]