The HMM package

my_state <- c("AT-rich", "GC-rich")
at_tran_prob <- c(0.7, 0.3)
gc_tran_prob <- c(0.1, 0.9)
tran_matrix <- matrix(c(at_tran_prob, gc_tran_prob), 2, 2, byrow = TRUE)
rownames(tran_matrix) <- my_state
colnames(tran_matrix) <- my_state

my_nuc <- c("A", "C", "G", "T")
at_emi_prob <- c(0.39, 0.1, 0.1, 0.41)
gc_emi_prob <- c(0.1, 0.41, 0.39, 0.1)
emi_matrix <- matrix(c(at_emi_prob, gc_emi_prob), 2, 4, byrow = TRUE)
rownames(emi_matrix) <- my_state
colnames(emi_matrix) <- my_nuc

generatehmmseq <- function(tran_matrix, emi_matrix, init_prob, seq_len){
  new_seq <- character()
  new_seq_state <- character()
  first_state <- sample(my_state, 1, prob=init_prob)
  my_prob <- emi_matrix[first_state,]
  first_nuc <- sample(my_nuc, 1, prob=my_prob)
  new_seq[1] <- first_nuc
  new_seq_state[1] <- first_state
  
  for (i in 2:seq_len){
    prev_state <- new_seq_state[i-1]
    state_probs <- tran_matrix[prev_state,]
    cur_state <- sample(my_state, 1, prob=state_probs)
    my_prob <- emi_matrix[cur_state,]
    nuc <- sample(my_nuc, 1, prob=my_prob)
    new_seq[i] <- nuc
    new_seq_state[i] <- cur_state
  }
  
  my_result <- list(sequence = new_seq, state = new_seq_state)
}

set.seed(31)
init_prob <- c(0.5, 0.5)
my_seq <- generatehmmseq(tran_matrix, emi_matrix, init_prob, 30)
my_seq

library(HMM)
my_hmm <- initHMM(my_state, my_nuc, transProbs = tran_matrix, emissionProbs = emi_matrix)

my_most_prob_state <- viterbi(hmm = my_hmm, observation = my_seq$sequence)
data.frame(nuc = my_seq$sequence, real = my_seq$state, inferred = my_most_prob_state, match = my_seq$state == my_most_prob_state)