量化跨膜序列中密码子的频率 - 应用功能？

Question

我正在尝试查看某些蛋白质跨膜结构域内的密码子使用情况。

为此，我有 TM 域的序列，我想在这些序列中搜索某些密码子出现的频率（频率）。

理想情况下，我想将新列添加到现有数据框中，其中包含每个基因每个密码子的计数。比如这个假设数据：

Gene ID	TM_domain_Seq	GGA
ENSG00000003989	TGGAGCCTCGCTC	1
ENSG00000003989	TGGAGCCTCGCTC	1
ENSG00000003989	TGGAGCCTCGCTC	1
ENSG00000003989	TGGAGCCTCGCTC	1
ENSG00000003989	TGGAGCCTCGCTC	1

我尝试了以下方法 - 创建一个函数来计算特定密码子出现的频率，并将其应用于每个 TM 序列。 我遇到的问题是如何将每个密码子的新列添加到我的数据框中，以及如何将密码子频率放入其中。

我试过 for 循环，但它们花费的时间太长

amino_search <- function(seq) {
  
  count <- str_count(seq, pattern = codons)
  return(count)
}

codon_search <- function(TMseq) {
  
 High_cor$Newcol <- unlist(lapply(TMseq, amino_search))
}

如有任何帮助，我们将不胜感激。谢谢！

Answer 1

创建可能组合的向量，然后使用 str_count:

comb <- expand.grid(replicate(3, c("A", "T", "G", "C"), simplify = FALSE)) |>
  apply(MARGIN = 1, FUN = paste, collapse = "")
  #apply(X = _, 1, FUN = paste, collapse = "") #with the new placeholder

df[, comb] <- t(sapply(df$TM_domain_Seq, stringr::str_count, comb))

如果您只需要 in-frame 个密码子，一种方法是每三个字符添加一个 space：

gsub('(.{3})', '\1 ', df$TM_domain_Seq[1])
#[1] "TGG AGC CTC GCT C"

df[, comb] <- t(sapply(gsub('(.{3})', '\1 ', df$TM_domain_Seq), stringr::str_count, comb))

输出

# A tibble: 5 × 66
  Gene_ID TM_domain_Seq   AAA   CAC   GGA   TAA   GAA   CAA   ATA   TTA   GTA   CTA   AGA   TGA   CGA   ACA   TCA
  <chr>   <chr>         <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
1 ENSG00… TGGAGCCTCGCTC     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
2 ENSG00… TGGAGCCTCGCTC     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
3 ENSG00… TGGAGCCTCGCTC     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
4 ENSG00… TGGAGCCTCGCTC     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
5 ENSG00… TGGAGCCTCGCTC     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
# … with 49 more variables: GCA <int>, CCA <int>, AAT <int>, TAT <int>, GAT <int>, CAT <int>, ATT <int>,
#   TTT <int>, GTT <int>, CTT <int>, AGT <int>, TGT <int>, GGT <int>, CGT <int>, ACT <int>, TCT <int>,
#   GCT <int>, CCT <int>, AAG <int>, TAG <int>, GAG <int>, CAG <int>, ATG <int>, TTG <int>, GTG <int>,
#   CTG <int>, AGG <int>, TGG <int>, GGG <int>, CGG <int>, ACG <int>, TCG <int>, GCG <int>, CCG <int>,
#   AAC <int>, TAC <int>, GAC <int>, ATC <int>, TTC <int>, GTC <int>, CTC <int>, AGC <int>, TGC <int>,
#   GGC <int>, CGC <int>, ACC <int>, TCC <int>, GCC <int>, CCC <int>

Answer 2

将问题拆分为 sub-problems，分别解决它们，然后组合解决方案。

第一个子问题是：如何获得给定 (in-frame) 序列的密码子频率？答案是要么使用 pre-made 解决方案（例如 Bioconductor 的 Biostrings:: trinucleotideFrequency(…, steps = 3L)），要么像下面这样快速而肮脏的东西：

codon_frequencies = function (seq) {
    # Take care of incomplete codon at end.
    len = nchar(seq) - (nchar(seq) %% 3L)
    start = seq(1L, len, by = 3L)
    substring(seq, start, start + 2L) |> table()
}

试一试：

codon_frequencies('TGGAGCCTCGCTC')
#
# AGC CTC GCT TGG
#   1   1   1   1

… 顺便说一句，你的序列有碎片密码子是故意的吗？如果是这样，您确定它们总是以完整的密码子开始吗？

好的。下一步是为 table 中的每个基因 ID 调用此函数并收集结果。在这一点上，计数 table 可以转换为整洁的数据框这一事实对我们有所帮助：

data.frame(codon_frequencies('TGGAGCCTCGCTC'))
#   Var1 Freq
# 1  AGC    1
# 2  CTC    1
# 3  GCT    1
# 4  TGG    1

为了我们的目的，这是一种方便的格式，因为它使 table 操作更容易（尤其是在使用整洁的数据格式时，我在下面使用包 'dplyr'、' tidyr' 和 'purrr'):

df |>
    group_by(`Gene ID`) |>
    summarize(map_dfr(TM_domain_Seq, ~ data.frame(codon_frequencies(.x))))
# # A tibble: 20 × 3
# # Groups:   Gene ID [1]
#    `Gene ID`       Var1   Freq
#    <chr>           <fct> <int>
#  1 ENSG00000003989 AGC       1
#  2 ENSG00000003989 CTC       1
#  3 ENSG00000003989 GCT       1
#  4 ENSG00000003989 TGG       1
# …

至此我们可能可以收工了：这是一种方便使用的格式。但是，如果您愿意，也可以将数据转换为宽格式：

    … |>
    pivot_wider(
        id_cols = `Gene ID`,
        names_from = Var1,
        values_from = Freq,
        values_fill = 0L # Otherwise missing codons will be `NA`
    )
# # A tibble: 5 × 11
# # Groups:   Gene ID [5]
#   `Gene ID`         AGC   CTC   GCA   TGG   TGA   GCG   AGT   GAT   TAC   GCT
#   <chr>           <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
# 1 ENSG00000003981     1     1     1     1     0     0     0     0     0     0
# 2 ENSG00000003982     1     1     0     1     1     0     0     0     0     0
# 3 ENSG00000003983     1     1     0     1     0     1     0     0     0     0
# 4 ENSG00000003984     0     0     0     0     0     0     1     1     1     0
# 5 ENSG00000003989     1     1     0     1     0     0     0     0     0     1

（这是使用了一些不同的玩具数据。）

最后，如果您想要所有密码子的列，即使是那些不存在于您的数据中的密码子，您可以对 codon_frequencies 函数进行小的修改：

all_codons = c('A', 'C', 'G', 'T') %>% expand.grid(., ., .) |> apply(1L, paste, collapse = '')

codon_frequencies = function (seq, all = FALSE) {
    # Take care of incomplete codon at end.
    len = nchar(seq) - (nchar(seq) %% 3L)
    start = seq(1L, len, by = 3L)
    codons = substring(seq, start, start + 2L)
    table(if (all) factor(codons, levels = all_codons) else codons)
}

然后在上面的代码中调用为codon_frequencies(.x, all = TRUE)。 pivot_wider 不再需要 values_fill = 0L 参数。

综合起来：

df |>
    group_by(`Gene ID`) |>
    summarize(
        map_dfr(TM_domain_Seq, ~ data.frame(codon_frequencies(.x, all = TRUE))),
        .groups = 'drop'
    ) |>
    pivot_wider(
        id_cols = `Gene ID`,
        names_from = Var1,
        values_from = Freq
    )

量化跨膜序列中密码子的频率 - 应用功能？

Quantifying frequency of codons in a transmembrane sequence - apply function?

r

bioinformatics

apply

dna-sequence

lapply