最近導(dǎo)師讓我跑模型，生物信息方向的，我一個(gè)學(xué)計(jì)算機(jī)的，好多東西都看不明白。現(xiàn)在的方向大致是，用深度學(xué)習(xí)的模型預(yù)測(cè)病毒感染人類(lèi)的風(fēng)險(xiǎn)。

既然是病毒，就需要拿到它的DNA，也就是堿基序列，然后把這些ACGT序列丟進(jìn)模型里面，然后就是預(yù)測(cè)能不能感染人類(lèi)，說(shuō)實(shí)話(huà)，估計(jì)結(jié)果不會(huì)好，現(xiàn)在啥都是transformer，而且我看的這篇論文，我認(rèn)為僅僅從DNA序列大概預(yù)測(cè)不出什么東西。

但是就那樣吧，現(xiàn)在數(shù)據(jù)去哪里下載，需要下載什么樣的數(shù)據(jù)，下載完成后怎么處理我還是一臉懵逼，但是假設(shè)上面都處理好了，然后即使把數(shù)據(jù)丟給模型，跑就完了。

也不是沒(méi)進(jìn)度，目前了解到的是，我應(yīng)該使用一種叫fasta格式的文件，然后把里面的一大串ACGT序列拿出來(lái)，轉(zhuǎn)為模型可以處理的數(shù)據(jù)。然后，以后再說(shuō)。

現(xiàn)在假設(shè)我已經(jīng)有了ACGT的序列，然后把它轉(zhuǎn)為模型可以處理的矩陣。

這里，我隨機(jī)生成長(zhǎng)度為131072的基因序列，為什么是這個(gè)數(shù)字呢，因?yàn)檫@是之前看的 論文里的值，，暫時(shí)按照這個(gè)來(lái)做。

實(shí)現(xiàn):

首先是導(dǎo)入庫(kù)

import numpy as np
import random
import tensorflow as tf
import inspect
from typing import Any, Callable, Dict, Optional, Text, Union, Iterable
import os

然后，定義一個(gè)生成長(zhǎng)度為131072bp的函數(shù)：

#隨機(jī)生成131072的dna序列
length = 131072
def randomSeq(length):
  return ''.join([random.choice('ACGT') for i in range(length)])

這個(gè)函數(shù)的返回結(jié)果是長(zhǎng)度為length的字符串，類(lèi)似ACGTTGC這樣。

然后這種序列模型是沒(méi)辦法處理的，所以需要把它變成矩陣，也就用one-hot編碼。

比如ACGT這個(gè)序列，編碼成：

[ [1,0,0,0],

[0,1,0,0],

[0,0,1,0],

[0,0,0,1] ]

這樣的一個(gè)矩陣，這個(gè)就不細(xì)說(shuō)了，網(wǎng)上很多資料。

然后，我從別人的代碼中抄了一個(gè)函數(shù)，很好用。

#DNA序列轉(zhuǎn)為one-hot編碼，可以直接拿來(lái)用
def one_hot_encode(sequence: str,
                   alphabet: str = 'ACGT',
                   neutral_alphabet: str = 'N',
                   neutral_value: Any = 0,
                   dtype=np.float32) -> np.ndarray:
  """One-hot encode sequence."""
  def to_uint8(string):
    return np.frombuffer(string.encode('ascii'), dtype=np.uint8)
  hash_table = np.zeros((np.iinfo(np.uint8).max, len(alphabet)), dtype=dtype)
  hash_table[to_uint8(alphabet)] = np.eye(len(alphabet), dtype=dtype)
  hash_table[to_uint8(neutral_alphabet)] = neutral_value
  hash_table = hash_table.astype(dtype)
  return hash_table[to_uint8(sequence)]

這是一個(gè)嵌套函數(shù)了，仔細(xì)研究下還是可以理解的，我就不說(shuō)了，會(huì)用就行了。

簡(jiǎn)單講一下參數(shù)的意思：

sequence:字符串類(lèi)型，就是輸入的堿基序列。

alphabet: str = ‘ACGT’ ：詞表，一共只需要這四個(gè)詞

neutral_alphabet: str = ‘N’,

neutral_value: Any = 0,

上面這兩一起用，就是說(shuō)遇到N這個(gè)堿基就會(huì)編碼成[0,0,0,0]的向量。

dtype=np.float32，這個(gè)就是內(nèi)部元素值的類(lèi)型。

簡(jiǎn)單生成一下：

然后輸入序列長(zhǎng)度是131072bp,所以輸入的矩陣就是131072x4的矩陣，現(xiàn)在來(lái)把序列變?yōu)榫仃嚒?/p>

編碼成one-hot矩陣

dnaVec = one_hot_encode(dna)

現(xiàn)在DNA序列已經(jīng)變成了矩陣，接下來(lái)需要把這一條序列，也就是一個(gè)樣本數(shù)據(jù)，變成TensorFlow中的TFRecord文件格式。TFRecord 是 TensorFlow 中的數(shù)據(jù)集存儲(chǔ)格式。當(dāng)我們將數(shù)據(jù)集整理成 TFRecord 格式后，TensorFlow 就可以高效地讀取和處理這些數(shù)據(jù)集，從而幫助我們更高效地進(jìn)行大規(guī)模的模型訓(xùn)練。

關(guān)于tfr文件的處理，我就不在細(xì)說(shuō)了，總之現(xiàn)在我們需要構(gòu)建example。

在此之前，我們需要先這么做：

#給出結(jié)果的tfr文件的路徑
path = '/content/drive/MyDrive/test_Enformer/result.tfr'
#dna的numpy數(shù)組轉(zhuǎn)成字節(jié)流，這樣才能存儲(chǔ)
dnaVec = dnaVec.tobytes()

接下來(lái)就是把這個(gè)字節(jié)流數(shù)據(jù)寫(xiě)入到tfr文件中，這里同時(shí)寫(xiě)入這條數(shù)據(jù)的label中，我的問(wèn)題是給一個(gè)Dna序列，預(yù)測(cè)是或者不是的二分類(lèi)問(wèn)題，所以我同時(shí)把這條dna序列對(duì)應(yīng)的真實(shí)標(biāo)簽也寫(xiě)進(jìn)去，但是我是隨機(jī)從0,1中選擇一個(gè)。

from tensorflow.core.example.feature_pb2 import BytesList
with tf.io.TFRecordWriter(path) as writer:
  feature = {
      #序列使用的是tf.train.BytesList類(lèi)型
      'sequence':tf.train.Feature(bytes_list=tf.train.BytesList(value=[dnaVec])),
      #label是隨機(jī)生成的0，或者1
      'label':tf.train.Feature(int64_list=tf.train.Int64List(value=[np.random.choice([0,1])]))
  }
  example = tf.train.Example(features=tf.train.Features(feature=feature))
  writer.write(example.SerializeToString())

這部分的代碼執(zhí)行結(jié)束后，就已經(jīng)把dna序列以及對(duì)應(yīng)的標(biāo)簽寫(xiě)入了tfr文件中，不過(guò)這個(gè)tfr文件中只有一個(gè)example，你可以寫(xiě)更多個(gè)。

剛剛寫(xiě)入的tfr文件

到這里，相當(dāng)于已經(jīng)把數(shù)據(jù)準(zhǔn)備好了，接下來(lái)就是讀取數(shù)據(jù)。

#從剛才的路徑中加載數(shù)據(jù)集
dataset = tf.data.TFRecordDataset(path)

#定義Feature結(jié)構(gòu)，告訴解碼器每個(gè)Feature的類(lèi)型是什么
feature_description = {"sequence": tf.io.FixedLenFeature((), tf.string),
              "label": tf.io.FixedLenFeature((), tf.int64)}

 #將 TFRecord 文件中的每一個(gè)序列化的 tf.train.Example 解碼
def parse_example(example_string):
  #解析之后得到的example
  example = tf.io.parse_single_example(example_string,feature_description)
  #example['sequence']還是字節(jié)流的形式，重新轉(zhuǎn)為數(shù)字向量
  sequence = tf.io.decode_raw(example['sequence'], tf.float32)
  sequence = tf.reshape(sequence,(length,4))  #形狀需要重塑，不然就是一個(gè)長(zhǎng)向量
  label = tf.cast(example['label'],tf.int64)  #標(biāo)簽對(duì)應(yīng)的類(lèi)型轉(zhuǎn)換
#每一天example解析后返回對(duì)應(yīng)的一個(gè)字典
  return {
      'sequence':sequence,
      'label': label
  }

#把parse_example函數(shù)映射到dataset中的每個(gè)example,
#這里的dataset中只有一個(gè)example
dataset = dataset.map(parse_example)

此時(shí)的dataset是一個(gè)可以遍歷的對(duì)象，內(nèi)部元素可以認(rèn)為是解析完成后的example。

這個(gè)字典有兩個(gè)鍵sequence和lable，對(duì)應(yīng)著序列矩陣和標(biāo)簽值

這就是可以用來(lái)訓(xùn)練的數(shù)據(jù)。

到此這篇關(guān)于Python實(shí)現(xiàn)將DNA序列存儲(chǔ)為tfr文件并讀取流程介紹的文章就介紹到這了,更多相關(guān)Python存儲(chǔ)tfr文件內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論