Vol.2 基于小鼠的基因集数数据库资源

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和统计开始, 自己边学边写, 把遇到的一些问题展现出来.

这一系列可能会要求有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 小结

好困, 人还在新主楼，好饿, 今天练胸了, 火锅真香。