接上回的故事，公司给了新一批的数据，但是批次效应比较重，需要去批次。去批次之后，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()