自作のモジュール（kerasy.ML.HMM.HMM）を使って実装をします。
In [1]:
from kerasy.ML.HMM import HMM
In [2]:
model = HMM()
In [3]:
# 以下のようなパラメタファイルを読み込む。
! cat params.json
In [4]:
model.load_params("params.json")
In [5]:
# インスタンスに読み込まれた各種パラメータ
model.params()
In [6]:
# なお、この時以下のような辞書を同時に作成しています。
print(model.base2int)
print(model.int2base)
In [7]:
# また、パラメータはnumpyの行列の形で読み込まれています。
print(type(model.pi))
print(type(model.A))
print(type(model.phi))
In [8]:
# 以下のようにすれば、学習したパラメタを同じ形式で出力することができます。
model.save_params("out.json")
Viterbi algorithm¶
HMMのViterbiアルゴリズムを実装して、尤もらしいサンプル配列の隠れ状態列を推定します。
In [9]:
from kerasy.utils.bio import alignStr, readMonoSeq, readMultiSeq
# ファイルの中身
cat seq1.fasta
```
sample-RNA
gagaguccuauacaaacuccaaaacacugagaccauacaaguaaaaccagucgagaaaau
agucaacagcacccccguugcguauccugcggagaaugccucgcuaucuguccucacgau
uggucuaaccgcucggcucaggcgugugggccugaaauccgggcagcaaacuacgguaag
uuuucgcguaucaaaaauacauugaugaauguucgcuauuagccgggucgacguuuugau
ggugacuaggagcgaaagugauuuuuuugugagcggugucauagcaggggaucaucugag
gugaacuauggacgggucagucgcccucauuggguuuguuacucagauacgugacacacg
uaaggucgcacggcaguagugauccacgagaaucggcacucuuacgagcaagucuauagc
gacguggcuugcuauuaagacaguaauggacucggacagcuaguacgcgacuaaucauuu
ucgcacgccaucacacuuucaacccaggagcuuacaguguuccaggggccaggcuugggg
cagccucuguggaagugcgagggcugcgcuaguguacauuagccgacccucagacgugaa
auaaagagaugcugcugguugcacagugagcaacguucacucggcaacccgucggauacc
```
In [10]:
# ファイルの読み込み
seq1 = readMonoSeq('seq1.fasta')
print(f"Input Sequence:\n{alignStr(seq1[0])}")
In [11]:
# この時、文字列をリストに格納しています。
print(type(seq1))
In [12]:
pred1 = model.predict(sequences=seq1, asai=True)
In [13]:
# 変数asaiは、隠れ状態がZero-based numberingかどうかを決めるものであり、特に意味はありません。（課題の形式に対応したものがasai=True.）
print(f"Hidden Indexes:\n{alignStr(pred1[0])}")
In [14]:
# もちろんパラメタの更新は行われていません。
model.params()
Baum-Welch algorithm¶
HMMのBaum-Welchアルゴリズムを実装して、モデルを学習させてます。
In [15]:
model = HMM()
model.load_params("params2.json")
model.params() # 文字が "acgu"→"ACGT"に変わっただけです。
cat seq2.txt
```
>Vitis_vinifera_chr14.trna15-AlaAGC (28776087-28776159) Ala (AGC) 73 bp Sc: 68.12
GGGGATGTAGCTCAGATGGTAGAGCGCTCGCTTAGCATGCGAGAGGTACGGGGATCGATA
CCCCGCATCTCCA
>Vitis_vinifera_chr13.trna5-ArgCCG (2691241-2691313) Arg (CCG) 73 bp Sc: 83.85
GATCGCATAGCGGAGTGGATATCGCGTTAGACTCCGAATCTAAAGGTCGTGGGTTCGAAT
CCCACTGCGATCG
>Vitis_vinifera_chr3.trna6-GlnCTG (5807449-5807520) Gln (CTG) 72 bp Sc: 71.93
GGTCCCATGGTCTAGTGGTCAGGACATTGGACTCTGAATCCAGTAACCCGAGTTCAAATC
TCGGTGGGACCT
>Vitis_vinifera_chr14.trna19-LeuAAG (25640120-25640040) Leu (AAG) 81 bp Sc: 62.71
GTGGAGATGGCCGAGTTGGTCTAAGGCGCCAGATTAAGGTTCTTGTCCGAAAGGGCGTGG
GTTCAAATCCCACTCTCCACA
>Vitis_vinifera_chr18.trna9-MetCAT (2896034-2896118) Met (CAT) 85 bp Sc: 60.66
GGGGTGGTGGCGCAGTTGGCTAGCGCGTAGGTCTCATAGCTTCCGAGTAATCCTGAGGTC
GAGAGTTCGAGCCTCTCTCACCCCA
>Vitis_vinifera_chr3.trna30-ProTGG (7498415-7498344) Pro (TGG) 72 bp Sc: 68.98
GGGCATTTGGTCTAGTGGTATGATTCTCGCTTTGGGTGCGAGAGGTCCCGAGTTCGATTC
TCGGAATGACCC
>Vitis_vinifera_chr5.trna13-GluTTC (21745120-21745048) Glu (TTC) 73 bp Sc: 36.67
AGTCCCATAGCCTAGTTGGTCGAGCACAAGGTTTTCAATCTTGTGGTCGTGGGTTCGAGC
CCTATTGGTGGTT
>Vitis_vinifera_chr5.trna3-TyrGTA (7422575-7422666) Tyr (GTA) 92 bp Sc: 33.68
CCGATCTTAACTCAATTGGTAGAGCAGAGGACCGTAGTGGGGTTGACCCAATAGATATCC
TTAGGTCGTTAGTTTGAATCCAACAGGTCTAA
>Vitis_vinifera_chr1.trna8-ValCAC (11396083-11396156) Val (CAC) 74 bp Sc: 66.86
GTCTAGGTGGTATAGTTGGTTATCACGCTAGTCTCACACACTAGAGGTCCCTAGTTCGAA
CCCAGGCTCAGATA
>Vitis_vinifera_chr14.trna2-SerTGA (4021381-4021462) Ser (TGA) 82 bp Sc: 80.82
GTCGATATGTCCGAGTGGTTAAGGAGACAGACTTGAAATCTGTTGGGCTTCGCCCGCGCA
GGTTCGAACCCTGCTGTCGACG
```
In [16]:
seq2 = readMultiSeq('seq2.txt')
In [17]:
for i, seq in enumerate(seq2):
print(f"No.{i:>02}\n{alignStr(seq, width=45)}")
In [18]:
# 学習させる。
model.fit(sequences=seq2, epochs=2000, verbose=False)
In [19]:
import matplotlib.pyplot as plt
In [20]:
def plotHistory(model):
""" 学習の過程をプロットする関数 """
plt.plot(model.history)
plt.title(f"Mean Log Likelihood Transition (eoochs={model.epoch})")
plt.xlabel("epochs")
plt.ylabel("Mean Log Likelihood")
plt.show()
In [21]:
plotHistory(model)
In [22]:
model.params()
※ 隠れ状態ごとにかなり偏った遷移を行なう。また、遷移行列の対角成分がそこまで大きい値では無い。
→ 特定の構造が、ある程度の長さのゆとりを持って繰り返されていると予想される。
In [23]:
pred2 = model.predict(sequences=seq2, asai=True)
In [24]:
for i,pred in enumerate(pred2):
print(f"No.{i:>02} Hidden Indexes:\n{alignStr(pred, width=45)}")
補足：上記プログラムにバグがないかの確認¶
Baum-Welchアルゴリズムが正しく実装できていたか、極端な配列を学習させることで確認します。
In [25]:
model = HMM()
model.load_params("params3.json")
model.params()
In [26]:
seq_extreme = readMonoSeq('seq3.txt')
In [27]:
# かなり極端な配列の例で実装してみる。
print(f"Input Sequence:\n{alignStr(seq_extreme[0])}")
上の配列は、以下の構造を持っていることがわかります。
| index | type |
|---|---|
| 1~150 | pyrimidine base |
| 151~450 | both |
| 451~600 | purine base |
In [28]:
model.fit(sequences=seq_extreme, epochs=100, verbose=False)
In [29]:
plotHistory(model)
In [30]:
model.params()
In [31]:
pred_extreme = model.predict(sequences=seq_extreme)
In [32]:
print(f"Hidden Indexes:\n{alignStr(pred_extreme[0])}")
うまく傾向をつかんでいることから、学習アルゴリズムに大きな誤りはなさそうです。
In [ ]: