接上回的故事，公司给了新一批的数据，但是批次效应比较重，需要去批次。去批次之后，fpkm出现了负值，但是counts在很多分析里不能直接用，所以需要counts转fpkm。目前网上关于人类基因组的已经比较多了，但是小鼠的我没咋找到，特别是biomart这个包我用的不熟练，绕了很多弯路，有点挠头。于是有了今天的文章，如何获取小鼠的基因长度。
rm(list=ls())
一键清空工作空间。
listMarts()
这个函数能够返回bioMRT可以连接的数据库列表，就是看看哪些数据库是能用的
mart = useMart('ensembl')
head(listDatasets(mart))
dataset <- listDatasets(mart)
dataset$description[grep("mouse",dataset$description)]
通过useMart函数得到一个mart对象，然后看里面有哪些数据集。这里注意区分数据库和数据集的概念。想想GEO数据库。如果找人类的基因集，这里有现成的GRCh38.p14，需要注意首字母是大写：
> dataset[grep("Human",dataset$description),]
dataset description version
80 hsapiens_gene_ensembl Human genes (GRCh38.p14) GRCh38.p14
然后可以选择自己想要的基因集了。
bmart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "mmusculus_gene_ensembl",
host = "https://www.ensembl.org")
这里把小鼠为对象的mart对象取出来。
Attributes <- listAttributes(bmart)
Attributes[grep("symbol",Attributes$name),]
attributes = c(
"ensembl_gene_id",
#"uniprot_gn_symbol",
"start_position",
"end_position"
)
filters = "ensembl_gene_id"
当然也可以搜其他类型的数据，比如
> Attributes[grep("entrez",Attributes$name),]
name description
58 entrezgene_trans_name EntrezGene transcript name ID
77 entrezgene_description NCBI gene (formerly Entrezgene) description
78 entrezgene_accession NCBI gene (formerly Entrezgene) accession
79 entrezgene_id NCBI gene (formerly Entrezgene) ID
这里要设置检索数据库的格式和参数之类的。下面会用到。
feature_info <- getBM(attributes = attributes,
filters = filters,
values = "",
mart = bmart)
这里的getBM函数是bioMRT查询的主要函数，attributes是属性，就是你需要返回的东西有哪些；filters是过滤器，就类似百度的高级检索，你需要哪些东西满足你的条件这样子；values就是你需要什么样的东西。通过这个函数，你能得到filter这个属性里符合values值的样本/观测。给定一组过滤器和相应的值，它从连接到的BioMart数据库中检索用户指定的属性。所以filters的长度和values的长度一定是相等的。
eff_length <- abs(feature_info$end_position - feature_info$start_position)
names(eff_length) <- feature_info$ensembl_gene_id
write.csv(eff_length, "./gene_length_mmu.csv", row.names = TRUE)
然后这样一减就得到基因长度了！
总结一下心得：
今天看Y叔的推送，谈到这么一个问题，说“访问内部数据和面向对象封装的概念是背道而驰的”。作为一个学临床的同学，特别是对数据的理解不深的时候，其实是不太能理解的。今天重新复习了biomart的操作之后，意识到问题所在。对于一个内部数据较少的对象（比如一个简单的ggplot2对象），访问内部数据可能是OK的。但是对于一个大型的数据，不可能一头扎进去找自己需要的数据。一个封装好的对象，肯定有各种函数能够进行操作或者是访问，需要的就是怎么合理利用好这些函数。比如今天用到的getBM()
