如何在R语言中将4×16三核苷酸计数表转为4×4汇总表

问题描述 #

现有如下三核苷酸计数数据：

aaa  aac  aag  aat  aca  acc  acg  act  aga  agc  agg  agt  ata  atc  atg  att  caa  cac  cag  cat  cca  ccc  ccg  cct 
1501  895  660  989  759  858  658  524  667  878  608  621  897  983  901 1030  669  713 1381  708  610  514 1308  542 
 cga  cgc  cgg  cgt  cta  ctc  ctg  ctt  gaa  gac  gag  gat  gca  gcc  gcg  gct  gga  ggc  ggg  ggt  gta  gtc  gtg  gtt 
 516  888 1398  664  196  410 1371  523 1071  624  519 1118 1045  881 1007  835  894 1037  734  915  618  660  862  863 
 taa  tac  tag  tat  tca  tcc  tcg  tct  tga  tgc  tgg  tgt  tta  ttc  ttg  ttt 
 805  567  214  996 1057  721  493  599 1255  965  840  803  871  816  729 1364 

已通过以下R代码生成4×16的矩阵R：

triple <- count(morseq, wordsize=3)
triple
R <- matrix(triple,4,16,byrow=T)
colnames(R) <- c("aa","ac","ag","at","ca","cc","cg","ct","ga","gc","gg","gt","ta","tc","tg","tt") ## 给列命名
row.names(R) <- c("a","c","g","t")

需要将该矩阵简化为4×4矩阵：行代表三核苷酸的第一个碱基，列代表第三个碱基，对4×16矩阵中首尾碱基匹配的对应值求和，预期格式如下：

a   c    g    t   
a   x   x    x    x  
c   x   x    x    x  
g   x   x    x    x  
t   x   x    x    x 

解决方案 #

方法1：基础R实现 #

通过提取列名的第三个碱基，按行分组求和：

# 提取每个列名的第三个碱基（列名的第二个字符）
third_bases <- substr(colnames(R), 2, 2)

# 初始化4×4结果矩阵，设置行名和列名
result_matrix <- matrix(0, nrow = 4, ncol = 4, 
                        dimnames = list(c("a", "c", "g", "t"), c("a", "c", "g", "t")))

# 遍历每一行，按第三个碱基分组求和并赋值
for (row_name in rownames(result_matrix)) {
  row_data <- R[row_name, ]
  # 按第三个碱基分组计算总和
  sum_values <- tapply(row_data, third_bases, sum)
  # 将结果填入对应位置
  result_matrix[row_name, names(sum_values)] <- sum_values
}

# 查看最终结果
result_matrix

方法2：tidyverse风格实现 #

如果熟悉dplyr和tidyr，可以用更简洁的代码：

library(dplyr)
library(tidyr)
library(tibble)

# 将矩阵转换为长格式，分组求和后转回宽格式矩阵
result_matrix <- R %>% 
  as_tibble(rownames = "first_base") %>%
  pivot_longer(-first_base, names_to = "last_two", values_to = "count") %>%
  mutate(third_base = substr(last_two, 2, 2)) %>%
  group_by(first_base, third_base) %>%
  summarise(total = sum(count), .groups = "drop") %>%
  pivot_wider(names_from = third_base, values_from = total) %>%
  column_to_rownames("first_base") %>%
  as.matrix()

# 查看结果
result_matrix

两种方法都能得到符合要求的4×4矩阵，其中每个位置的值是对应第一个碱基和第三个碱基的所有三核苷酸计数之和。

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