我们一般做测序的时候，如果没有明确要求进行mRNA富集，一般给出的测序数据是包含编码和非编码RNA的。对于不同课题组经常会有不同的研究倾向，需要从数据中将所有非编码RNA提取或剔除出去。这个其实很简单，可以自己下载GFT文件，用Excel VLOOKUP函数进行比对，然后复制出来。咱们今天讲一下如何用R语言实现。
开始之前，大家先了解一下什么是Ensembl ID、Entrez ID和Gene symbol？
Gene ID 也称Entrez ID/EntrezGene ID ，是 NCBI 使用的能够对众多数据库进行联合搜索的搜索引擎, 其对不同的 Gene 进行了编号, 每个 gene 的编号就是 entrez gene id，就是一串数字，比如：TP53 的Gene ID是:7157。因为entrez ID相对稳定, 所以也被其他数据库, 如 KEGG 等采用。不同物种的同一个基因的Gene ID是不同的。NCBI的RefSeq数据库ID，一般是两个大写首字母，加下划线，后面为数字。两个首字母 ”NC”、”NM”、”NP_”分别代表DNA、mRNA、Protein。
Gene Symbol ，是HGNC数据库为基因提供的官方名称，HGNC是人类基因命名委员会(HUGO Gene Nomenclature Committee);人类基因组命名委员会。有专门的数据库：https://www.genenames.org/。主要是按基因的功能起的名字，字母一般为英文全称的缩写，由大写字母和数字组成，如TP53基因的Official Symbol就是由HGNC所提供。需要注意的是，因为基因命名的变化，可能一个基因有多个别名，甚至一些别名与其他基因名字一样，因此使用基因Symbol一定要注意，不建议一开始分析就用symbol。
Ensembl ID，是在Ensembl数据库中对基因的命名，常见的物种前缀：“ENS“表示Homo sapiens (Human)，”ENSMUS“表示Mus musculus (Mouse)，”ENSDAR“表示Danio rerio (Zebrafish)；而常见的序列类型用G、P、T、分别表示gene、protein和transcript。这个和Entrez ID一样比较稳定，甚至优于Entrez，在版本更新之后会在相应ID后面添加一个版本号.**。TCGA下载的数据一般就是用的Ensembl ID。
接下来开始使用BioMart包提取基因注释信息：
# Load the required packages
suppressPackageStartupMessages({
library(biomaRt)
library(tidyverse)
})
# Use biomaRt functions to create the conversion table
biomart_table <- getBM(
attributes = c(
"external_gene_name",
"ensembl_gene_id",
"entrezgene_id",
"gene_biotype"
),
mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
)
#########对于不同物种，请修改dataset，可以使用如下代码查看：
##mart = useMart('ensembl')
##listDatasets(mart)
##比如小鼠：mmusculus_gene_ensembl，斑马鱼：drerio_gene_ensembl
# Replace empty values with NA
biomart_table[biomart_table == ""] <- NA
biomart_table <- biomart_table %>%
dplyr::select(
"Symbol" = external_gene_name,
"Ensembl" = ensembl_gene_id,
"Entrez" = entrezgene_id,
"Biotype" = gene_biotype
) %>%
arrange(Symbol, Ensembl, Entrez, Biotype)
# Save the output for use in the app
saveRDS(biomart_table, "app_data/biomart_table_homo.Rds")
至此，就可以得到注释文件，接下来的任务就是匹配，可以使用Excel，方便的话还是用R：
biomart_table <- readRDS("app_data/biomart_table_homo.Rds")
results <- read.delim("diffexpr-results.txt") #数据第一列为Ensembl ID
results2 <- merge(results,biomart_table[c(2,1,4)],by = "Ensembl")
results2 <- merge(biomart_table[c(2,1,4)],by = "Ensembl",results)
write.csv(results2,"results.csv")
此外，还可以根据这个包进行人鼠同源基因名字转化，使用getLDS函数，它是biomaRt查询的主要功能，连接两个数据集，并从这些链接的biomaRt数据集检索信息。在Ensembl中，这转化为同源映射。
## using the asia mirror
human <- useEnsembl('ensembl',dataset = "hsapiens_gene_ensembl", mirror = "asia")
mouse <- useEnsembl('ensembl',dataset = "mmusculus_gene_ensembl", mirror = "asia")
#比如我有一系列小鼠基因：
mouse.genes <- c("Psme2b","Zfp661","Tsc22d1","Prox2","Eml3","Creb1")
#将其映射到人的基因上
MtoH <- getLDS(attributes = "mgi_symbol", # 要转换符号的属性，这里基因名（第3步是基因名）
filters = "mgi_symbol", #参数过滤
mart = mouse, #需要转换的基因名的种属来源，也就是第2步的mouse
values = mouse.genes, #要转换的基因集
attributesL = "hgnc_symbol", #要同源转换的目标属性，这里还是转为基因名，也可加其他
martL = human, #要同源转换的目标种属，也就是第2步的human
uniqueRows = TRUE)
###Error: biomaRt has encountered an unexpected server error.换成其它镜像也不行
###有人给出了权宜之计，用既往存档的数据,参考https://support.bioconductor.org/p/9144001/
##https://www.ensembl.org/info/website/archives/assembly.html
human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = "https://dec2021.archive.ensembl.org/")
mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl", host = "https://dec2021.archive.ensembl.org/")
##换成2022的archive，又不行了
human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = "https://oct2022.archive.ensembl.org/")
mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl", host = "https://oct2022.archive.ensembl.org/")
##烦躁，
打赏赞(13)微海报分享