计算重复采样中字符串出现的次数

Question

我正在尝试完成一个计算遗传学作业问题。当我已知的核苷酸频率（A、T、C、G）都为 0.25 并且我对 1000 个字符长的字符串进行 1000 次采样时，我需要计算二核苷酸“AA”出现的次数。最终，确定“AA”二核苷酸的频率。

SeqLength <- 1000
Nreps <- 1000
x <- c("A", "T", "C", "G")
pi <- c(rep(0.25, 4))
counter <- 0

for (i in 1:Nreps){
    X <- sample(x, SeqLength, replace=TRUE, pi)
    for (j in 1:X - 1){
        if (j == "AA"){
            counter <- counter + 1
        }
    }
}
counter / Nreps

我收到以下错误和警告。它们都有意义，我只是不确定如何解决它们。

Error in 1:X : NA/NaN argument
In addition: Warning messages:
1: In 1:X : numerical expression has 1000 elements: only the first used
2: In 1:X : NAs introduced by coercion

我似乎无法将我脑海中想做的事情转化为功能代码。非常感谢任何帮助/指导/课程。谢谢！！

Answer 1

这是一种方法。您可以使用replicate 执行采样Nreps 次数。对于每个采样的序列，您检查每个位置和其后的位置的核苷酸是否为"A"，并将结果逻辑相加得到计数。

SeqLength <- 1000
Nreps <- 1000
x <- c("A", "T", "C", "G")
pi <- c(rep(0.25, 4))
counter <- 0

set.seed(42)
seqs = replicate(n = Nreps,
                 expr = sample(x = x, size = SeqLength, replace = TRUE, prob = pi),
                 simplify = FALSE)
counts = sapply(seqs, function(s){
    sum((s[-1] == "A") & (s[-length(s)] == "A"))
})
mean(counts)
#> [1] 62.321

^{由reprex package (v0.3.0) 于 2019 年 9 月 13 日创建}