如何用dplyr管道分割字符串并统计每行氨基酸频率

我来帮你用tidyverse处理这个问题，咱们可以分几种方式实现，按需选择就行～

首先先确认下你的原始数据：

library(tidyverse)

dat <- structure(list(fasta_header = c(">seq1", ">seq2"), sequence = c("MPSRGTRPE", "VSSKYTFWNF")), .Names = c("fasta_header", "sequence"), row.names = c(NA, -2L), class = c("tbl_df", "tbl", "data.frame"))

dat
#> # A tibble: 2 x 2
#>   fasta_header sequence
#>   <chr>        <chr>
#> 1 >seq1        MPSRGTRPE
#> 2 >seq2        VSSKYTFWNF

方法1：统计氨基酸出现次数（宽格式） #

这个方式会把每个氨基酸作为单独的列，显示对应序列中的出现次数，没出现的氨基酸自动填0：

aa_counts <- dat %>%
  # 把序列拆分成单个氨基酸的列表
  mutate(aa = str_split(sequence, "")) %>%
  # 展开列表，让每个氨基酸单独占一行
  unnest_longer(aa) %>%
  # 按fasta标题和氨基酸分组，统计出现次数
  group_by(fasta_header, aa) %>%
  count(name = "count") %>%
  # 转成宽格式，缺失的氨基酸填充0
  pivot_wider(names_from = aa, values_from = count, values_fill = 0) %>%
  # 合并回原始数据，保留原序列列
  left_join(dat, ., by = "fasta_header")

# 查看结果
aa_counts

运行后会得到这样的输出：

#> # A tibble: 2 × 12
#>   fasta_header sequence     M     P     S     R     G     T     E     V     K     Y     F     W     N
#>   <chr>        <chr>    <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
#> 1 >seq1        MPSRGTRPE     1     1     1     2     1     1     1     0     0     0     0     0     0
#> 2 >seq2        VSSKYTFWNF    0     0     2     0     0     1     0     1     1     1     2     1     1

方法2：计算氨基酸频率（宽格式） #

如果你需要的是频率（出现次数/序列总长度），可以在上面的基础上调整：

aa_frequencies <- dat %>%
  # 先计算序列长度，方便后续算频率
  mutate(seq_length = str_length(sequence),
         aa = str_split(sequence, "")) %>%
  unnest_longer(aa) %>%
  group_by(fasta_header, aa) %>%
  count(name = "count") %>%
  # 计算频率：次数除以序列总长度
  mutate(frequency = count / first(seq_length)) %>%
  select(-count) %>%
  pivot_wider(names_from = aa, values_from = frequency, values_fill = 0) %>%
  left_join(dat, ., by = "fasta_header")

# 查看结果
aa_frequencies

方法3：长格式结果（适合后续分析） #

要是你更喜欢长格式（每行对应一个氨基酸的统计结果），可以跳过转宽格式的步骤：

long_format <- dat %>%
  mutate(aa = str_split(sequence, "")) %>%
  unnest_longer(aa) %>%
  group_by(fasta_header, sequence, aa) %>%
  count(name = "count") %>%
  mutate(frequency = count / str_length(sequence))

# 查看结果
long_format

这个结果会是这样的：

#> # A tibble: 14 × 5
#> # Groups:   fasta_header, sequence, aa [14]
#>    fasta_header sequence     aa    count frequency
#>    <chr>        <chr>        <chr> <int>     <dbl>
#>  1 >seq1        MPSRGTRPE    E         1     0.111
#>  2 >seq1        MPSRGTRPE    G         1     0.111
#>  3 >seq1        MPSRGTRPE    M         1     0.111
#>  4 >seq1        MPSRGTRPE    P         1     0.111
#>  5 >seq1        MPSRGTRPE    R         2     0.222
#>  6 >seq1        MPSRGTRPE    S         1     0.111
#>  7 >seq1        MPSRGTRPE    T         1     0.111
#>  8 >seq2        VSSKYTFWNF   F         2     0.2  
#>  9 >seq2        VSSKYTFWNF   K         1     0.1  
#> 10 >seq2        VSSKYTFWNF   N         1     0.1  
#> 11 >seq2        VSSKYTFWNF   S         2     0.2  
#> 12 >seq2        VSSKYTFWNF   T         1     0.1  
#> 13 >seq2        VSSKYTFWNF   V         1     0.1  
#> 14 >seq2        VSSKYTFWNF   W         1     0.1  

这样就能灵活满足你的需求啦～

内容的提问来源于stack exchange，提问作者littleworth