使用Biostrings包vmatchPattern匹配密码子DNAStringSetList报错如何解决

vmatchPattern函数要求输入的pattern参数必须是单个字符串或XString对象，你从codons中取出的codons[[i]]是DNAStringSet类型，其中像精氨酸、亮氨酸等氨基酸对应多个密码子pattern，属于长度大于1的DNAStringSet，不符合参数要求，因此触发报错。另外你原代码中的酪氨酸和色氨酸的密码子序列写反了，需要先修正避免后续结果错误。

前置修正 #

首先修正密码子定义笔误，给列表加名称方便后续结果查看：

# 修正笔误：酪氨酸对应TAY，色氨酸对应TGG
tyrosine <- DNAStringSet("TAY")
tryptophan <- DNAStringSet("TGG")

# 其余密码子定义保留原有写法，重新生成codons列表并命名
codons <- DNAStringSetList(list(alanine, arginine, asparginine, aspartic_acid, asparagine_or_aspartic_acid,
                           cysteine, glutamine, glutamic_acid, glutamine_or_glutamic_acid, glycine,
                           histidine, start, isoleucine, leucine, lysine, methionine, phenylalanine,
                           proline, serine, threonine, tyrosine, tryptophan, valine, stop))
names(codons) <- c("丙氨酸","精氨酸","天冬酰胺","天冬氨酸","天冬酰胺/天冬氨酸","半胱氨酸",
                   "谷氨酰胺","谷氨酸","谷氨酰胺/谷氨酸","甘氨酸","组氨酸","起始密码子",
                   "异亮氨酸","亮氨酸","赖氨酸","甲硫氨酸","苯丙氨酸","脯氨酸","丝氨酸",
                   "苏氨酸","酪氨酸","色氨酸","缬氨酸","终止密码子")

方案1：for循环版本（适配你的使用习惯） #

library(Biostrings)
# 读取新冠参考基因组
filepath = "你的本地fasta文件路径"
covid <- readDNAStringSet(filepath)
# 提取第一条序列作为匹配主体，readDNAStringSet读入的结果为DNAStringSet，取第一个元素才是XString类型
reference_seq <- covid[[1]]

codon_locations <- list()
for (i in seq_along(codons)) {
  # 取出当前氨基酸对应的所有密码子pattern
  current_patterns <- codons[[i]]
  # 逐个匹配所有pattern，再合并结果
  current_matches <- lapply(current_patterns, function(p) vmatchPattern(p, reference_seq, fixed = FALSE))
  codon_locations[[i]] <- current_matches
}
# 给结果列表命名
names(codon_locations) <- names(codons)

注意这里加了fixed = FALSE参数兼容密码子中的简并碱基（比如N、R、Y等）。

方案2：lapply版本（代码更简洁） #

codon_locations <- lapply(codons, function(aa_patterns) {
  lapply(aa_patterns, function(p) vmatchPattern(p, reference_seq, fixed = FALSE))
})

如果需要匹配固定阅读框的密码子，可通过匹配结果的起始位置过滤：比如仅保留start(ranges(match_result)) %% 3 == 1的结果，即可得到正链第一个阅读框的密码子匹配结果。定位SNP位点所属密码子时，直接比对SNP坐标是否落在匹配到的密码子区间内即可。

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