差异基因得到的方法有很多,例如DESeq2、EdgeR、Wilcox、SCDE等等,真正有意义的差异基因,即使使用不同的方法,最后结果也不会相差太多
有了差异基因,然后就是去注释
rm(list=ls())
options(stringsAsFactors = F)
library(clusterProfiler)
load(file = 'step3.1-DEG-monocle_summary.Rdata')
de_genes <- subset(de_clusters, qval<0.05)
> length(de_genes$genes)
[1] 4435
entrez_genes <- bitr(de_genes$genes, fromType="SYMBOL",
toType="ENTREZID",
OrgDb="org.Mm.eg.db")
entrez_genes <- entrez_genes[!entrez_genes$ENTREZID %in% "101055843",]
> length(entrez_genes$SYMBOL)
[1] 4281
因为有一些de_genes
的SYMBOL没有对应的ENTREZID,因此看到少了100多个基因。
de_gene_clusters <- de_genes[de_genes$genes %in% entrez_genes$SYMBOL,
c("genes", "cluster")]
# 保持de_gene_clusters$genes的顺序不变,将他的symbol变成entrez ID
de_gene_clusters <- data.frame(
ENTREZID=entrez_genes$ENTREZID[entrez_genes$SYMBOL %in% de_gene_clusters$genes],
cluster=de_gene_clusters$cluster
)
也就是拆分成4个cluster组成的列表
list_de_gene_clusters <- split(de_gene_clusters$ENTREZID,
de_gene_clusters$cluster)
formula_res <- compareCluster(
ENTREZID~cluster,
data=de_gene_clusters,
fun="enrichGO",
OrgDb="org.Mm.eg.db",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05
)
# org.Mm.eg.db更新于2019年4月
# 可视化
pdf('DEG_GO_each_cluster.pdf',width = 11,height = 6)
dotplot(formula_res, showCategory=5)
dev.off()
The simplified version of enriched result is more clear and give us a more comprehensive view of the whole story https://guangchuangyu.github.io/2015/10/use-simplify-to-remove-redundancy-of-enriched-go-terms/
start_time <- Sys.time()
lineage1_ego <- simplify(
formula_res,
cutoff=0.5,
by="p.adjust",
select_fun=min
)
end_time <- Sys.time()
(end_time - start_time)
# Time difference of 1.588762 mins
pdf('DEG_GO_each_cluster_simplified.pdf',width = 11,height = 6)
dotplot(lineage1_ego, showCategory=5)
dev.off()
write.csv(formula_res@compareClusterResult,
file="DEG_GO_each_cluster.csv")
# 简化版本
write.csv(lineage1_ego@compareClusterResult,
file="DEG_GO_each_cluster_simplified.csv")
看第一行:背景基因是20000多个,其中属于这个GO通路的有256个基因;然后C1这个cluster这里总共得到的差异基因是1236个,其中属于这个GO通路的是66个
library(org.Mm.eg.db)
# 结果有3百多万个
go2gene <- toTable(org.Mm.egGO2ALLEGS)
其实看到里面的基因ID存在重复,一个go_id会对应多个同样的gene_id,那么就要去重复,获取unique gene id
uni_gene <- unique(go2gene$gene_id[go2gene$go_id == 'GO:0140014'])
> length(uni_gene)
[1] 256
结果和?的BgRatio
中记录的一致
c1_genes <- list_de_gene_clusters[['C1']]
> length(c1_genes)
[1] 1284
这里的1284比之前的1236多个几十个,说明存在几十个基因没有GO注释
其实就是找c1_genes
和go2gene
的交叉
overlap_genes <- intersect(c1_genes,uni_gene)
> length(overlap_genes)
[1] 66
也就是对应上面GO通路的66个基因
Monocle的结果会用到上一篇的 1.6 用包装好的pheatmap画热图 的代码
rm(list = ls())
options(warn=-1)
options(stringsAsFactors = F)
source("../analysis_functions.R")
# 加载RPKM表达矩阵
load('../female_rpkm.Rdata')
# 加载ROTS 的差异分析结果
load('DEG-ROTS_summary.Rdata')
head(summary_pop1)
head(summary_pop2)
head(summary_pop3)
head(summary_pop4)
# 6个发育时期获取
head(colnames(females))
female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
names(females) <- colnames(females)
table(female_stages)
# 4个cluster获取
cluster <- read.csv('female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)
RPKM.full=females
population_subset<-c(rownames(summary_pop1[summary_pop1$ROTS.statistic<0,])[1:18],
rownames(summary_pop2[summary_pop2$ROTS.statistic<0,])[1:18],
rownames(summary_pop3[summary_pop3$ROTS.statistic<0,])[1:18],
rownames(summary_pop4[summary_pop4$ROTS.statistic<0,])[1:18])
总共得到了72个基因,然后取出
RPKM_heatmap<-RPKM.full[population_subset,]
将小表达矩阵的列名重新按照4个cluster排序,并且log转换
RPKM_heatmap<-RPKM_heatmap[,order(female_clustering)]
RPKM_heatmap<-log2(RPKM_heatmap+1)
每个群赋予不同颜色
popul.col<-sort(female_clustering);table(popul.col)
popul.col<-replace(popul.col, popul.col=='C1',"#1C86EE" )
popul.col<-replace(popul.col, popul.col=='C2',"#00EE00" )
popul.col<-replace(popul.col, popul.col=='C3',"#FF9912" )
popul.col<-replace(popul.col, popul.col=='C4',"#FF3E96" )
画图
library(gplots)
png("ROST-DEG-heatmap.png")
heatmap.2(as.matrix(RPKM_heatmap),
ColSideColors = as.character(popul.col), tracecol = NA,
dendrogram = "none",col=bluered, labCol = FALSE,
scale="none", key = TRUE, symkey = F, symm=F, key.xlab = "",
key.ylab = "", density.info = "density",
key.title = "log2(RPKM+1)", keysize = 1.2, denscol="black", Colv=FALSE)
dev.off()
rm(list=ls())
options(stringsAsFactors = F)
load(file = 'DEG-monocle_summary.Rdata')
load(file = 'DEG-ROTS_summary.Rdata')
# monocle差异基因
mnc_DEG <- subset(de_clusters, qval<0.05)
# ROTS差异基因
head(summary_pop1)
head(summary_pop2)
head(summary_pop3)
head(summary_pop4)
g1=rownames(summary_pop1[summary_pop1$FDR>0.05,])
g2=rownames(mnc_DEG[mnc_DEG$cluster=='C1',])
> length(intersect(g1,g2));length(g1);length(g2)
[1] 194
[1] 7006
[1] 1319
看到ROTS得到了7006个差异基因,monocle得到1319个,它们共有194个。就数字来讲,ROTS挑选的有些”粗糙“,一般差异基因都不会挑选太多
## 对第二群比较
g1=rownames(summary_pop2[summary_pop2$FDR>0.05,])
g2=rownames(mnc_DEG[mnc_DEG$cluster=='C2',])
> length(intersect(g1,g2));length(g1);length(g2)
[1] 620
[1] 7015
[1] 1100
# 第二群,二者的重叠度较高
## 对第三群比较
g1=rownames(summary_pop3[summary_pop3$FDR>0.05,])
g2=rownames(mnc_DEG[mnc_DEG$cluster=='C3',])
> length(intersect(g1,g2));length(g1);length(g2)
[1] 205
[1] 7909
[1] 1129
## 对第四群比较
g1=rownames(summary_pop2[summary_pop2$FDR>0.05,])
g2=rownames(mnc_DEG[mnc_DEG$cluster=='C2',])
> length(intersect(g1,g2));length(g1);length(g2)
[1] 620
[1] 7015
[1] 1100