1 问题的提出
超几何分布检验和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是最出名的,但不是唯一的!
但是对于人和鼠的基因来说,人的基因多数是大写字母表示,而鼠是同源基因多数是首字母大写,其余字母小写。但仅仅首字母大写并不能完全实现人鼠基因转换。可以找到一一对应的基因列表,进行转换。因此本文使用其他方法进行人鼠基因的转换
希望得到:
全部的50个基因集差异情况
下面的表格里面的表头分别是:
人基因集的基因数量
小鼠基因集的基因数量
两个基因集overlap数量
人特有的基因数量
小鼠特有的基因数量:
2 答案
使用msgdb数据库提取hallmarker通路基因集,就可以解决了
rm(list=ls())
library(stringr)
library(msigdbr)
library(dplyr)
library(purrr)
hallmark.mouse=msigdbr(species = "Mus musculus",category = "H") #提取基因信息
hallmark.human=msigdbr(species = "Homo sapiens",category = "H")
v1=hallmark.human%>%group_by(gs_name)%>%summarise(v1=n())
v2=hallmark.mouse%>%group_by(gs_name)%>%summarise(v2=n())
v3=hallmark.human%>%filter(hallmark.human$human_gene_symbol%in%str_to_upper(hallmark.mouse$gene_symbol))%>%group_by(gs_name)%>%summarise(v3=n())
v4=hallmark.human%>%filter(!str_to_title(hallmark.human$human_gene_symbol)%in%hallmark.mouse$gene_symbol)%>%group_by(gs_name)%>%summarise(v4=n())
v5=hallmark.mouse%>%filter(!str_to_upper(hallmark.mouse$gene_symbol)%in%hallmark.human$gene_symbol)%>%group_by(gs_name)%>%summarise(v5=n())
res=reduce(list(v1,v2,v3,v4,v5),.f=full_join,by="gs_name") #reduce批量整合
res=data.frame(res)
res[is.na(res)]<-0 #转换NA为0
res
网友(董理聪)投稿
应曾老师的邀请,这里发邮件说说生成全排列的算法。也只是基本算法的实现而已。
递归函数
如果一个函数在内部调用自身,就是递归函数。
通常来说,一门程序语言的教程在讲到函数的时候,也都会提及到 迭代 与 递归。
最经典的使用递归函数的例子是生成 Fibonacci 数列。
def Fibonacci(n):
if n == 0:
return 0
if n in [1, 2]:
return 1
return Fibonacci(n - 1) + Fibonacci(n - 2)
于是生成全排列的方法也可以使用递归函数。
def permutations(iterable):
if len(iterable) == 1:
yield iterable
# 每次抽出一个元素,其他元素再产生全排列
for i, element in enumerate(iterable):
for p in permutations(iterable[:i] + iterable[i + 1:]):
yield [element, ] + p
permutations(list(range(1, 9)))
这个网友强推python ,一直在使用 Jupyter Notebook;希望我们以后的教程,比如 rm(ls = ls()) 的代码没有必要放在推文里;
参考了学徒的答案:将每次递归的结果保存下来可以使用函数传递参数的方法实现。
perm_wapper=function(strLst){
result<-list()
idx<-1
perm=function(x,temp=c()){
if (length(x)==1) {
result[[idx]]<-c(temp,x)
idx<-idx+1
}else{
for (i in x){
ind=which(i==x)
perm(x[-ind],c(temp,i))
}
}
}
perm(strLst)
return(result)
}
res=perm_wapper(c("a","b","c"))
[[1]]
[1] "a" "b" "c"
[[2]]
[1] "a" "c" "b"
[[3]]
[1] "b" "a" "c"
[[4]]
[1] "b" "c" "a"
[[5]]
[1] "c" "a" "b"
[[6]]
[1] "c" "b" "a"