如何从Biostrings的vmatchPattern结果中提取基因名称

我正在开展5'UTR结合分析，需要从中提取基因名称，编写的代码在vmatchPattern步骤之前运行正常：

library(biomaRt)
library(GenomicFeatures)
library(XVector)
library(Biostrings)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(BSgenome.Mmusculus.UCSC.mm10)

fUTR <- fiveUTRsByTranscript(TxDb.Mmusculus.UCSC.mm10.knownGene)

Mmusculus <- BSgenome.Mmusculus.UCSC.mm10
seqlevelsStyle(Mmusculus) <- 'ensembl'
seqlevelsStyle(fUTR) <- 'ensembl'

Seq <- getSeq(Mmusculus, fUTR)
Pbind <- RNAString('UGUGUGAAHAA')

Match <- vmatchPattern(Pbind, unlist2(Seq), max.mismatch = 0, min.mismatch = 0, with.indels = F, fixed = T, algorithm = 'auto')

后续需要提取基因名称生成列表用于Python的RNA-seq分析，但尝试三种方法均失败：

##How to get gene names from the match Pattern
#1
matches <- unlist(Match, recursive = T, use.names = T)
m <- as.matrix(matches)
subseq(genes[rownames(m),], start = m[rownames(m),1], width = 20)

#2
transcripts(TxDb.Mmusculus.UCSC.mm10.knownGene, columns = c('tx_id', 'tx_name', 'gene_id'))

#3
count_index <- countIndex(Match)
wh <- which(count_index > 0)
result_list = list()
for(i in 1: length(wh))  
{
  result_list[[i]] = Views(subject[[wh[i]]], mindex[[wh[i]]])
  
}
names(result_listF) = nm[wh]

问题分析 #

你之前的方法存在几个关键问题：

• 方法1中genes对象未定义，无法直接调用
• 方法2仅提取了转录本与基因的映射关系，但未关联到你的匹配结果
• 方法3中subject、mindex、nm均为未定义变量，导致运行报错

可行解决方案 #

核心思路是：从匹配结果中拿到对应的转录本ID，通过TxDb关联到基因ID，再用biomaRt转换为基因名称（symbol），具体代码如下：

# 1. 从匹配结果中提取有匹配的转录本ID，去重避免重复
matched_tx_ids <- unique(names(unlist(Match)))

# 2. 从TxDb中获取转录本ID对应的基因ID（Entrez ID）
tx2gene <- select(TxDb.Mmusculus.UCSC.mm10.knownGene,
                  keys = matched_tx_ids,
                  columns = "GENEID",
                  keytype = "TXID")

# 3. 用biomaRt将Entrez ID转换为基因名称（external_gene_name即基因symbol）
# 连接小鼠的Ensembl数据库
mart <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
gene_mapping <- getBM(attributes = c("entrezgene_id", "external_gene_name"),
                      filters = "entrezgene_id",
                      values = tx2gene$GENEID,
                      mart = mart)

# 4. 合并转录本、基因ID、基因名称的关联结果
final_table <- merge(tx2gene, gene_mapping, 
                     by.x = "GENEID", by.y = "entrezgene_id")

# 5. 提取唯一的基因名称列表，用于后续Python分析
unique_gene_symbols <- unique(final_table$external_gene_name)

# 可选：将结果保存为csv文件，方便Python读取
write.csv(unique_gene_symbols, "matched_genes.csv", row.names = FALSE)

补充说明 #

• 如果你不需要保留转录本和基因的对应关系，仅需要基因名称，直接取unique_gene_symbols即可
• 若biomaRt连接失败，可尝试指定历史镜像host：useMart("ensembl", dataset = "mmusculus_gene_ensembl", host = "https://dec2023.archive.ensembl.org")
• 保存为csv后，Python可直接用pandas.read_csv()读取

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