不使用工具包读取多序列FastA文件的优化方案咨询

问题背景 #

用户有如下格式的FASTA多序列文件：

>gi|1079586|gb|AAA82053.1| polymerase [Human endogenous retrovirus K]
IQKTSGRWRMLTDLRAVNAVIQPKGHLQPGLPSPAMIPKDWPLIIIDLKDCFFTIPLEEQDFEKFAFTIP
AINNKEPATRF
>gi|138362618|gb|ABO76712.1| polymerase [Kodoko virus]
SERESNSESLSKALSLTKCLSAALKNLCFYSEESPTSYTSVGPDSGRLKFALSYKEQVGGNRELYIGDLR
TKMFTRLVEDYFESFTSFFHGSCLNDEKEFENAILSMTVNVRQGYLSYSMDHSKWGPM
>gi|387147|gb|AAA37560.1| polymerase [Mus musculus]
PSLQAHLQALQAVQREVWKPLAAAYQDQQDQPVIPHPFRVGDTVWVRRHQTKNLEPRWKGPYTVLLTTPT
ALKVDGIAAWIHAAHVKAATTPPAGTASGPTWKVQRSQNPLKIRLTRGAP
>gi|74474929|dbj|BAE44450.1| polymerase [Hepatitis B virus]
MPLSYQHFRKLLLLDDEAGPLEEELPRLADEGLNRRVAEDLNLGNLNVSIPWTHKVGNFTGLYSSTVPVF
NPEWQTPSFPNIHLQEDIINRCQQYVGPLTVNEKRRLKLIMPARFYPNLTKYLPLDKGIKPYYPEHAVNH
YFKTRHYLHTLWKAGILYKRETTRSASFCGSPYSWEQELQHGRLVFQTSTRHGDESFCSQSSGILSRSPV
GPCVRSQLKQSRLGLQPQQGSLARGKSGRSGSIRARVHPTTRRSFGVEPSGSGHIDNSASSASSCLHQSA
VRKTAYSHLSTSKRQSSSGHAMELHHIPPSSARSQSEGPIFSCWWLQFRNSKPCSDYCLTHIVNLLEDWG
PCTEHGEHNIRIPRTPARVTGGVFLVDKNPHNTTESRLVVDFSQFSRGSTHVSWPKFAVPNLQSLTNLLS
SNLSWLSLDVSAAFYHIPLHPAAMPHLLVGSSGLPRYVARLSSTSRNINYQHGTMQDLHDSCSRNLYVSL
LLLYKTFGRKLHLYSHPIILGFRKIPMGVGLSPFLLGQFTSAICSVVRRAFPHCLAFSYMDDVVLGAKSV
QHLESLFTSITNFLLSLGIHLNPNKTKRWGYSLNFMGYVIGSWGTLPQEHIVLKLKQCFRKLPVNRPIDW
KVCQRIVGLLGFAAPFTQCGYPALMPLYACIQAKQAFTFSPTYKAFLCQQYLHLYPVGRQRSGLCQVFAD
ATPTGWGLAIGHRRMRGTFVAPLPIHTAELLAACFARSRSGAKLIGTDNSVVLSRKYTSFPWLLGCAANW
ILRGTSFVYVPSALNPADDPSRGRLGLYRPLLHLPFRPTTGRTSLYAVSPSVPSHLPDRVHFASPLHVAW
RPP

需求为：

• 不依赖第三方工具包读取文件
• 验证文件是否符合标准FASTA格式
• 将内容存储为可单独获取头部（header）和序列（sequence）的结构

之前的尝试存在问题：

1. 单序列处理代码会将所有序列合并，无法区分多条目
2. 用strsplit拆分的方式会丢失分隔符，产生无效空元素，且无法为列表元素命名

解决方案 #

1. 标准FASTA格式验证规则 #

解析前先验证文件格式，确保符合以下规则：

• 文件非空，且至少包含一条序列
• 所有header行必须以>开头，且不能为空
• header行之后的序列行仅包含合法字符（氨基酸序列允许A-Z/a-z，核酸序列允许A/T/C/G/U/a/t/c/g/u，可根据需求调整）
• 不能出现连续的header行（即两个>行之间必须有序列内容）

2. 解析与存储实现代码 #

以下是完整的R代码，包含格式验证和多序列解析功能，最终将数据存储为命名列表，每个列表元素是包含header和sequence的子列表：

# 读取文件（本地路径或URL）
con <- "URL或本地文件路径"
lines <- readLines(con)
# 过滤空行
lines <- lines[nchar(lines) > 0]

# --------------------------
# 格式验证
# --------------------------
# 检查是否有header行
header_indices <- grep("^>", lines)
if (length(header_indices) == 0) {
  stop("无效FASTA文件：未找到以>开头的头部行")
}
# 检查是否有连续的header行
if (any(diff(header_indices) == 1)) {
  stop("无效FASTA文件：存在连续的头部行，缺少对应序列")
}
# 检查文件是否以header开头
if (header_indices[1] != 1) {
  stop("无效FASTA文件：第一条非空行必须是头部行")
}
# 检查序列行是否包含非法字符（以氨基酸为例，可根据需求修改正则）
sequence_indices <- setdiff(1:length(lines), header_indices)
invalid_sequence <- grepl("[^A-Za-z]", lines[sequence_indices])
if (any(invalid_sequence)) {
  invalid_lines <- paste(sequence_indices[invalid_sequence], collapse = ", ")
  stop(paste("无效FASTA文件：序列行包含非法字符，行号：", invalid_lines))
}

# --------------------------
# 解析多序列
# --------------------------
fasta_data <- list()
# 遍历每个header，获取对应的序列
for (i in seq_along(header_indices)) {
  current_header <- lines[header_indices[i]]
  # 确定当前序列的结束位置
  if (i == length(header_indices)) {
    seq_lines <- lines[(header_indices[i]+1):length(lines)]
  } else {
    seq_lines <- lines[(header_indices[i]+1):(header_indices[i+1]-1)]
  }
  # 合并序列行
  current_sequence <- paste(seq_lines, collapse = "")
  # 提取header的标识部分作为列表名称（比如gi|xxx部分）
  name <- sub("^>(.*?)\\|.*", "\\1", current_header)
  # 添加到结果列表
  fasta_data[[name]] <- list(
    header = current_header,
    sequence = current_sequence
  )
}

# 示例：获取第一条序列的header和sequence
fasta_data[["gi"]]$header
fasta_data[["gi"]]$sequence

代码说明 #

• 格式验证：通过检查header行的位置、是否连续、序列字符合法性等，确保输入符合标准FASTA格式
• 解析逻辑：通过定位所有header行的索引，分割每个条目对应的序列行，合并后存储为子列表
• 命名规则：提取header中的gi标识作为列表元素名称，也可根据需求修改为其他标识（比如gb编号）
• 访问方式：通过fasta_data[[名称]]$header和fasta_data[[名称]]$sequence单独获取每个条目的头部和序列

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