Vol.2 基于小鼠的基因集数数据库资源
“像清风追着流云 望你时忽远忽近”
About Author
Yifan Fu
Undergraduate, BUAA, Beijing
- Major : Bioinformatics
- Recent focus : single-cell transcriptomics, metagenomics, multi-omics interaction
- Contact : fan@buaa.edu.cn
0x00 写在前面
又开了一个系列新坑, 这一系列的文章主要面向生物信息学, 尝试着从R和统计开始, 自己边学边写, 把遇到的一些问题展现出来.
0x01 问题的提出
超几何分布检验和GSEA的差异
通常拿到了上下调差异基因列表,然后说的GO/KEGG数据库注释,指的是超几何分布检验。
但是如果我们并不是首先自定义阈值,确定上下调差异基因列表,而是根据某个指标(比如logFC)把全部的基因排序,再去进行GO/KEGG数据库注释,一般来说就是GSEA分析啦。
但是数据库不仅仅是GO/KEGG哦
MSigDB(Molecular Signatures Database)数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列(Collection),每个系列分别是:
H: hallmark gene sets (癌症)特征基因集合,共50组,最常用;
C1: positional gene sets 位置基因集合,根据染色体位置,共326个,用的很少;
C2: curated gene sets:(专家)校验基因集合,基于通路、文献等:
C3: motif gene sets:模式基因集合,主要包括microRNA和转录因子靶基因两部分
C4: computational gene sets:计算基因集合,通过挖掘癌症相关芯片数据定义的基因集合;
C5: GO gene sets:Gene Ontology 基因本体论,包括BP(生物学过程biological process,细胞原件cellular component和分子功能molecular function三部分)
C6: oncogenic signatures:癌症特征基因集合,大部分来源于NCBI GEO 发表芯片数据
C7: immunologic signatures: 免疫相关基因集合。
可以看到,GO/KEGG是最出名的,但不是唯一的!
在 wehi的MSigDB发现了两份清单,下载后希望进行数据库的相互匹配与检索
http://bioinf.wehi.edu.au/software/MSigDB/human_H_v5p2.rdata
http://bioinf.wehi.edu.au/software/MSigDB/mouse_H_v5p2.rdata
但是对于人和鼠的基因来说,人的基因多数是大写字母表示,而鼠是同源基因多数是首字母大写,其余字母小写。但仅仅首字母大写并不能完全实现人鼠基因转换。可以找到一一对应的基因列表,进行转换。因此本文使用其他方法进行人鼠基因的转换
希望得到:
全部的50个基因集差异情况
下面的表格里面的表头分别是:
人基因集的基因数量
小鼠基因集的基因数量
两个基因集overlap数量
人特有的基因数量
小鼠特有的基因数量:
0x02 解决方法
目前有的方法主要是大小写转换, 正则表达式匹配等, 但存在如TP53(human) 对应 Trp53(mouse)的情况. 为了解决这一问题, 还是要回到公共数据库进行检索匹配. 使用R包biomaRt进行.借用了Vol.1写的工具,进行了一下改进。
load("human_H_v5p2.rdata") #Hs.H
load("mouse_H_v5p2.rdata") #Mm.H
if(!require(biomaRt))
BiocManager::install("biomaRt")
library(biomaRt)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
h.gene <- lapply(Hs.H, function(x){
bitr(x, #将rdata中的ID转换为gene Symbol
fromType = "ENTREZID",
toType = "SYMBOL",
OrgDb = org.Hs.eg.db)[,2]
})
m.gene <- lapply(Mm.H, function(x){
bitr(x,
fromType = "ENTREZID",
toType = "SYMBOL",
OrgDb = org.Mm.eg.db)[,2]
})
identical(names(h.gene),names(m.gene))#gene分类名相同
human=biomaRt::useMart("ensembl","hsapiens_gene_ensembl")
mouse=biomaRt::useMart("ensembl","mmusculus_gene_ensembl")
Ms2Hs <- function(gene){
require(biomaRt)
geneTrans <- getLDS(
values = gene,mart = mouse,
attributes = "mgi_symbol",filters = "mgi_symbol",
attributesL = "hgnc_symbol",
martL=human)
hs.gene <- as.vector(geneTrans$HGNC.symbol)
names(hs.gene) <- geneTrans$MGI.symbol
return(hs.gene)}
m2h.gene <- lapply(m.gene, function(x){
Ms2Hs(x)
})
result <- data.frame("name"=names(h.gene))
result$m <- as.numeric(lapply(m2h.gene, function(x){length(x)}))
result$h <- as.numeric(lapply(h.gene, function(x){length(x)}))
for( i in 1:50){
result$m.op[i] <- as.numeric((table(m2h.gene[[i]]%in%h.gene[[i]]))["TRUE"])
result$h.op[i] <- as.numeric((table(h.gene[[i]]%in%m2h.gene[[i]]))["TRUE"])
result$m.diff[i] <- as.numeric((table(m2h.gene[[i]] %in% h.gene[[i]]))["FALSE"])
result$h.diff[i] <- as.numeric((table(h.gene[[i]] %in% m2h.gene[[i]]))["FALSE"])
}
#tips:数学加和不等的原因
m2h.gene[[1]][duplicated(m2h.gene[[1]])]
#会出现多个基因对应人类同一个基因的情况,对于这种情况计数会导致差异
0x03 小结
好困, 人还在新主楼,好饿, 今天练胸了, 火锅真香。