如何使用 R 进行采样，从而最大化相关性-解网

问：

我有以下参考序列：

reference_seq <- "KPAACQHRQDKWKNSHWNRFKAYFVVIKKK"

使用此函数可以简单地计算氨基酸序列的组成：

calculate_aa_content <- function (x) 
{

  AADict <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", 
              "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
  AAC <- summary(factor(strsplit(x, split = "")[[1]], levels = AADict), 
                 maxsum = 21)/nchar(x)
  AAC
}

我可以得到参考序列的 AA 组成：

> calculate_aa_content(reference_seq)
         A          R          N          D          C          E 
0.10000000 0.06666667 0.06666667 0.03333333 0.03333333 0.00000000 
         Q          G          H          I          L          K 
0.06666667 0.00000000 0.06666667 0.03333333 0.00000000 0.23333333 
         M          F          P          S          T          W 
0.00000000 0.06666667 0.03333333 0.03333333 0.00000000 0.06666667 
         Y          V 
0.03333333 0.06666667

然后我有种子序列：

seed          <- "FKDHKHIDVKDRHRTRHLAK??????????"

我想做的是从 10 个氨基酸载体中采样 10aa，由 20 个问号（？）表示：

 AADict <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", 
                  "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")

因此，从种子与参考序列中最终构建的序列的 Pearson 相关性最大化（如果不可能获得绝对最大值，则允许一定次数的迭代）。

例如，种子采样的样本之一可以是：

> sample_1 <- "FKDHKHIDVKDRHRTRHLAKHHHHHHHHHH"
> cor(calculate_aa_content(reference_seq), calculate_aa_content(sample_1))
[1] 0.2955238

但相关性很低。我想做的是在一定次数的迭代后找到从种子构造的字符串与参考序列的最大相关性。

如果与电流最大值的差值在某个阈值（例如 0.01）内，则具有提前停止的附加功能将不胜感激。

请注意，最终序列不应与参考序列相同。

如何使用 R 有效地实现这一点？

R 字符串序列相关采样

calculate_aa_content <- function (x) {
  AADict <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", 
              "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
  AAC <- summary(factor(strsplit(x, split = "")[[1]], levels = AADict), 
                 maxsum = 21)/nchar(x)
  AAC
}

fun <- function(x, ref, dict = AADict) {
  i <- regexpr("\\?", x)
  n <- nchar(x) - i + 1L
  new <- sample(dict, n, TRUE)
  new <- paste(new, collapse = "")
  substring(x, i) <- new
  list(
    new = x,
    cor = cor(calculate_aa_content(ref), calculate_aa_content(x))
  )
}

reference_seq <- "KPAACQHRQDKWKNSHWNRFKAYFVVIKKK"
seed <- "FKDHKHIDVKDRHRTRHLAK??????????"

AADict <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", 
            "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")

set.seed(2022)
R <- 1e4
res <- replicate(R, fun(seed, reference_seq), simplify = FALSE)
cor_vec <- sapply(res, `[[`, 'cor')
res[which.max(cor_vec)]
#> [[1]]
#> [[1]]$new
#> [1] "FKDHKHIDVKDRHRTRHLAKKKWYKEKWFN"
#> 
#> [[1]]$cor
#> [1] 0.7902037

^{创建于 2022-10-21 使用 reprex v2.0.2}

如果有多个值等于最大值，请改用

res[cor_vec == max(cor_vec)]

为什么你的代码不能采用这种模式？它产生这个作为它的新序列。在该示例中，只需对????进行采样(4 ?)在中间。我怎样才能修改你的代码，以便它可以泛化为具有任何位置的种子。其他示例包括或"FKDHKHIDVKDRHRTHLAK????RTRHLAK""FKDHKHIDVKDRHRTHLAKWKNKYCFNAWF"?????HIDVKDRHRTHLAKKKKRTRHLAKKDRH?THLAKKKKRT?HLAK

0赞 littleworth 10/22/2022

我专门为那个 stackoverflow.com/questions/74162714/ 提出了一个新问题......

1赞 Rui Barradas 10/23/2022

@littleworth在那里回答。

上一个：使用 Biopython 保存一个 MODIED 的 FASTA 文件

下一个：如何判断每个子序列是否为递增阶 python

如何使用 R 进行采样，从而最大化相关性

How to maximize correlation by sampling using R

评论

评论