适用背景
之前写了一篇跨物种同源基因转换的文章，里面用到biomaRt这个包，但是这个包只能对2个物种进行同源基因转换。如果有多个物种，使用起来就比较麻烦，因此本文将介绍一步转换多个物种同源基因的函数。
useEnsembl和useMart的区别
上一篇文章提到，用useMart时有时候会出现网络连接错误的问题，可以指定host解决，另外根据一篇Bioconductor文章，建议使用useEnsembl替代useMart，或许可以解决网络问题。
查看帮助文档发现，这两者除了名称没有太大差别，如果要深究需要看源代码了。因此，本文将使用useEnsembl而不用useMart，可以把下面脚本的useEnsembl都替换成useMart使用。
- useEnsembl Connects to the selected BioMart database and dataset hosted by Ensembl
- useMart Connects to the selected BioMart database and dataset
关于biomaRt未收录的非模式物种
因为有些非模式物种没有被收纳进biomaRt包，因此在找其它同源基因库时发现一个HGNC网站。
HGNC
这个网站有个同源基因转换工具，虽然物种不是很多，但列出了不少同源基因库的信息，感兴趣可以逐个了解一下。
HCOP
2个物种的同源基因转换函数get2sp_orth_gene
- ref，参考物种的数据集名称，例如小鼠的是mmusculus_gene_ensembl
- que，查询物种的数据集名称
- gene.use，指定基因集，默认为空，使用查询物种的所有基因
- a.ref，参考物种输出的属性类型向量，默认会输出c("ensembl_gene_id","external_gene_name")
- f.ref，参考物种输入的属性类型，默认为基因名"external_gene_name"
- a.que，查询物种输出的属性类型向量，默认输出基因名c("external_gene_name")
- host， 使用的数据库网站版本，默认为"https://dec2021.archive.ensembl.org/"
如果需要查询某个物种的数据集名称以及可用的属性类型，可以查看跨物种同源基因转换。
函数代码如下：
get2sp_orth_gene <- function(ref=NULL,que=NULL,gene.use=NULL,
a.ref=NULL,f.ref="external_gene_name",
a.que=NULL,uniqueRows=T,host = "https://dec2021.archive.ensembl.org/"){
a.ref <- unique(c("ensembl_gene_id","external_gene_name",a.ref))
a.que <- unique(c("external_gene_name",a.que))
mref <- useEnsembl("ensembl",dataset = ref, host = host)
mque <- useEnsembl("ensembl",dataset = que, host = host)
if (is.null(gene.use)) {
g1 <- getBM(attributes = c("external_gene_name"),mart = useEnsembl("ensembl",dataset = que, host = host),uniqueRows = T)
g2 <- g1$external_gene_name
}else{
g2 <- gene.use
}
gene.out <- getLDS(attributes = a.ref, filters = f.ref, values = g2, mart = mref,
attributesL = a.que, martL = mque, uniqueRows = uniqueRows)
nref <- length(a.ref)
nque <- length(a.que)
recol <- colnames(gene.out)
recol[1:2] <- c('ref.gene_ID','ref.gene_name')
recol[nref+1] <- 'que.gene_name'
colnames(gene.out) <- recol
gene.out$ref.specice <- ref
gene.out$que.specice <- que
return(gene.out)
}
- 使用示例
这个函数返回的是一个数据框，第1列ref.gene_ID是参数物种的基因ID，第2列ref.gene_name是参考物种的基因名，第3列que.gene_name是询问物种的基因名，第4列是ref.specice参考物种的数据集名称，第5列是参考物种的数据集名称que.specice。
> g1 <- get2sp_orth_gene(ref="mmusculus_gene_ensembl",que="mfascicularis_gene_ensembl")
> head(g1)
ref.gene_ID ref.gene_name que.gene_name ref.specice
1 ENSMUSG00000038473 Nos1ap NOS1AP mmusculus_gene_ensembl
2 ENSMUSG00000036885 Arhgef26 ARHGEF26 mmusculus_gene_ensembl
3 ENSMUSG00000033285 Wdr3 WDR3 mmusculus_gene_ensembl
4 ENSMUSG00000074480 Mex3a MEX3A mmusculus_gene_ensembl
5 ENSMUSG00000062127 Cttnbp2nl CTTNBP2NL mmusculus_gene_ensembl
6 ENSMUSG00000041355 Ssr2 SSR2 mmusculus_gene_ensembl
que.specice
1 mfascicularis_gene_ensembl
2 mfascicularis_gene_ensembl
3 mfascicularis_gene_ensembl
4 mfascicularis_gene_ensembl
5 mfascicularis_gene_ensembl
6 mfascicularis_gene_ensembl
多个个物种的同源基因转换函数getmtsp_orth_gene
getmtsp_orth_gene调用了前面的get2sp_orth_gene函数，可以一步实现多个物种的同源基因转换，大部分参数与get2sp_orth_gene一样，只有que和gene.use不一样，que需要传入查询物种数据集的向量，例如que=c('mmusculus_gene_ensembl','mfascicularis_gene_ensembl','tguttata_gene_ensembl')，而gene.use如果不为空则需输入基因集列表，例如gene.use= list(c("Gad1","Sst"), c("GAD1","SST"))，基因集列表要一一对应查询物种，如果为空则会使用所有基因。
函数代码如下：
getmtsp_orth_gene <- function(ref=NULL,que=NULL,gene.use=NULL,
a.ref=NULL, f.ref="external_gene_name",
a.que=NULL, uniqueRows=T,host = "https://dec2021.archive.ensembl.org/"){
df2 <- data.frame()
for (i in 1:length(que)) {
if (is.null(gene.use)) {
g1 <- getBM(attributes = c("external_gene_name"),mart = useEnsembl("ensembl",dataset = que[i], host = host),uniqueRows = T)
g2 <- g1$external_gene_name
}else{
g2 <- gene.use[[i]]
}
df1 <- get2sp_orth_gene(ref=ref,que=que[i],gene.use=g2,
a.ref=a.ref,f.ref=f.ref,a.que=a.que,uniqueRows=uniqueRows,host=host)
if (i==1) {
df2 <- df1
}else{
df2 <- merge(df2,df1,all=T)
}
}
return(df2)
}
- 使用示例
函数返回的也是一个数据框，输出结果与get2sp_orth_gene的类似，如果指定了a.que，则会根据a.que增加相应的列。
> g3 <- getmtsp_orth_gene(ref='hsapiens_gene_ensembl', que= c('mmusculus_gene_ensembl','mfascicularis_gene_ensembl','tguttata_gene_ensembl'),
a.que= c("hsapiens_homolog_associated_gene_name","hsapiens_homolog_orthology_type"))
> head(g3)
ref.gene_ID ref.gene_name que.gene_name Human.gene.name Human.homology.type
1 ENSG00000000003 TSPAN6 Tspan6 TSPAN6 ortholog_one2one
2 ENSG00000000003 TSPAN6 TSPAN6 TSPAN6 ortholog_one2one
3 ENSG00000000003 TSPAN6 TSPAN6 TSPAN6 ortholog_one2one
4 ENSG00000000005 TNMD Tnmd TNMD ortholog_one2one
5 ENSG00000000005 TNMD TNMD TNMD ortholog_one2one
6 ENSG00000000005 TNMD TNMD TNMD ortholog_one2one
ref.specice que.specice
1 hsapiens_gene_ensembl mmusculus_gene_ensembl
2 hsapiens_gene_ensembl mfascicularis_gene_ensembl
3 hsapiens_gene_ensembl tguttata_gene_ensembl
4 hsapiens_gene_ensembl mmusculus_gene_ensembl
5 hsapiens_gene_ensembl mfascicularis_gene_ensembl
6 hsapiens_gene_ensembl tguttata_gene_ensembl
-
补充
如果想找一对一同源的同源基因，可以对参数a.que赋值，以人为例，hsapiens_homolog_associated_gene_name表示对应的人同源基因，hsapiens_homolog_orthology_type表示同源类型，ortholog_one2one表示的就是一对一同源，可以根据生成的列进行筛选。
human.name <-c("hsapiens_homolog_associated_gene_name","hsapiens_homolog_orthology_type")
mouse.name <- c("mmusculus_homolog_associated_gene_name","mmusculus_homolog_orthology_type")
网友评论