マルコフ連鎖確率確率行列を構築する方法

Question

私は自分自身でRを学習していますが、markovchainパッケージを使用してRstudioで遷移確率行列を構築しようとするといくつかの問題があります。まず、DNA配列の転移確率を計算しようとしました。

ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT 

しかし、どのように推移確率行列は、このような順序で計算することができ、私はRのインデックスを使用して考えていたが、私は実際にそれらの遷移確率を計算する方法がわかりません。

Rでこれを行う方法はありますか？ 私は行列のそれらの確率の出力はこのようなものでなければならないこと推測しています：DNA

 A T C G 
    A 0.60 0.10 0.10 0.20 
    T 0.10 0.50 0.30 0.10 
    C 0.05 0.20 0.70 0.05 
    G 0.40 0.05 0.05 0.50 

Answer 1

作成

DNA <- "ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT" 

文字によって分割されて文字が

DNA_list <- unlist(strsplit(DNA, split = "")) 

ユニークな要素を取得します。

DNA_unique <- unique(DNA_list) 

空行列を作成

matrix <- matrix(0, ncol = length(DNA_unique), nrow=length(DNA_unique)) 

それを埋める：iと要素I + 1をELTとマトリックスの対応するセル内の1つを追加します。

for (i in 1:(length(DNA_list) - 1)){ 
    index_of_i <- DNA_unique == DNA_list[i] 
    index_of_i_plus_1 <- DNA_unique == DNA_list[i + 1] 
    matrix[index_of_i, index_of_i_plus_1] = matrix[index_of_i, index_of_i_plus_1] + 1 
} 

matrix <- matrix/rowSums(matrix) 

> matrix 
      [,1]  [,2]  [,3]  [,4] 
[1,] 0.3000000 0.2166667 0.2250000 0.2583333 
[2,] 0.3068182 0.2954545 0.2159091 0.1818182 
[3,] 0.2857143 0.2142857 0.2619048 0.2380952 
[4,] 0.3764706 0.2235294 0.1882353 0.2117647 

NBをそれを正規化：あなたは計算する非常に大きなDNAを持っている場合、より高速な方法でそれを実行する方法があるかもしれません。しかし、ここでそれは十分に速く見えます。

Answer 2

markovchainパッケージを使用してこれを支援することができます。まず、あなたのデータ

seq <- "ATTCAACACATCCAGCCACATGCTCCGAGAGGAGGCAGAGGGCCCCCGGAATGATGCTTACCGAGATTCTTGTTTTTATCCTCGTGGTTGTTTAAAAACGAGTTGAAACTGACGGCATGTCGGACTATAAGCTACTTACTCACCATAGACGTGACCATAGGCCCTAAAACGTTACCGAGATATTCACTTCTAATAACAGTTGTCGGCAGAGCCAAAAGGCCGGGTGATAATACTTTAAAAAGGGAGTTGATTGTTGTATCTAATCCTAGAATGTCAAGAGCGACCATAACAAGATAATTCGGCAGAGCCAGAAAGCGTTCAAGGACTAGAACCATACCGAGACGCAAACGTTCAGGTCGAACTCTAATACCGATTAGT" 

その後

library(markovchain) 
base_sequence <- strsplit(seq, "")[[1]] 
mcX <- markovchainFit(base_sequence)$estimate 
mcX 

#   A   C   G   T 
# A 0.3000000 0.2250000 0.2583333 0.2166667 
# C 0.2857143 0.2619048 0.2380952 0.2142857 
# G 0.3764706 0.1882353 0.2117647 0.2235294 
# T 0.3068182 0.2159091 0.1818182 0.2954545