R语言如何拆分每行长度不均的单列核苷酸序列为25个单字符列？

直接按顺序运行代码即可，所有步骤均加了注释，替换对应文件路径、文件名就能直接使用。

前置准备 #

先把存储序列的csv文件放到好找的文件夹里，记住文件夹路径；Windows系统注意路径里的默认反斜杠\要手动替换成正斜杠/，否则会识别报错。

步骤1：读取原始数据 #

用基础R函数读取即可，无需额外安装第三方包：

# 把引号内内容替换为你存放文件的文件夹实际路径
setwd("C:/your/sequence/folder/path")
# 把引号内内容替换为你自己的csv文件名
seq_data <- read.csv("your_sequence_file.csv", stringsAsFactors = FALSE)
# 运行该行查看前6行数据，确认存在Sequence列、序列读取正常
head(seq_data)

如果输出结果中能看到Sequence列、列内容为你存储的核苷酸字符串，即可进入下一步。

步骤2：拆分序列为25列（自动截断超长部分） #

代码逻辑为先保留所有序列的前25个字符，超出部分直接丢弃，再将每个位置的核苷酸拆分为单独列：

# 截取所有序列的前25位
trunc_seq <- substr(seq_data$Sequence, start = 1, stop = 25)
# 逐行拆分单个字符，合并为25列的数据框
split_result <- as.data.frame(
  t(
    sapply(trunc_seq, function(x) strsplit(x, "")[[1]], USE.NAMES = FALSE)
  ),
  stringsAsFactors = FALSE
)
# 给新列命名，Pos1对应序列第1位，Pos25对应序列第25位
colnames(split_result) <- paste0("Pos", 1:25)
# 需要保留原序列列就运行这行，不需要可直接跳过
final_data <- cbind(seq_data, split_result)
# 检查拆分结果，确认每个Pos列内仅存储单个核苷酸字符
head(final_data)

如果习惯用tidyverse风格的代码，也可以用以下实现方式，首次使用需要先安装对应包：

# 首次运行才需要执行安装，安装过可直接注释掉该行
install.packages("tidyverse")
library(tidyverse)

final_data <- seq_data %>%
  mutate(trunc_seq = str_sub(Sequence, 1, 25)) %>%
  separate(
    col = trunc_seq,
    into = paste0("Pos", 1:25),
    sep = 1:24,
    remove = FALSE
  )

步骤3：统计每个位置的核苷酸出现频率 #

拆分完成后运行对应代码，即可得到每个位置各核苷酸的计数、占比频率，支持直接导出为csv文件：

基础R版本 #

freq_list <- list()
for (pos_col in colnames(split_result)) {
  nt_count <- table(split_result[[pos_col]])
  nt_freq <- prop.table(nt_count)
  freq_list[[pos_col]] <- data.frame(
    position = pos_col,
    nt = names(nt_count),
    count = as.numeric(nt_count),
    freq = as.numeric(nt_freq)
  )
}
freq_result <- do.call(rbind, freq_list)
# 查看统计结果前6行
head(freq_result)
# 导出结果到csv，文件会存储在之前setwd设置的文件夹内
write.csv(freq_result, "nt_position_freq.csv", row.names = FALSE)

tidyverse版本 #

freq_result <- final_data %>%
  select(starts_with("Pos")) %>%
  pivot_longer(cols = everything(), names_to = "position", values_to = "nt") %>%
  count(position, nt, name = "count") %>%
  group_by(position) %>%
  mutate(freq = count / sum(count)) %>%
  ungroup()
# 导出结果
write.csv(freq_result, "nt_position_freq.csv", row.names = FALSE)

常见报错排查：

• 如果提示找不到Sequence对象，检查原始csv的列名拼写、大小写是否和代码内写法一致
• 如果拆分后出现空值/NA，先运行table(nchar(seq_data$Sequence))查看序列长度分布，过滤掉长度不足25的序列再重新运行：seq_data <- seq_data[nchar(seq_data$Sequence)>=25, ]
• 如果读文件提示路径不存在，检查setwd内的路径是否正确，反斜杠是否全部替换为正斜杠

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