写在前面：经常做转录组分析，就是把差异基因搞个火山图和Venn图看各组差异基因的共有和特有情况。看见有个比较好的选择，能直观比较各种处理带来的影响，如下：
image.png
来自Nature plants的一篇文章
Ref:
https://github.com/YulongNiu/MPIPZ_microbe-host_homeostasis
https://www.nature.com/articles/s41477-021-00920-2
这个图很牛逼啊，表示的信息量很全，值得学习。去扒作者的代码，复现出了大部分
所需文件：
总的基因丰度表，即各个基因在每个样品中的丰度
image.png
各个样品的基因差异表达情况
举一个例子，这是Deseq2算出来的：
image.png
1. 导入数据并做一些处理
setwd("C:/Users/yjk/Desktop/cluster_DEGs")
library("dplyr")
library("ComplexHeatmap")
library("tibble")
library("RColorBrewer")
library("ggplot2")
my_theme()
all_genes <- read.table("all_genes.txt", header = TRUE) # fpkm
# DESeq2获得的差异表达基因(DEGs), |log2foldchange| > 1, padj < 0.05
HF12N_Rs <- read.table("HF12N_Rs.txt", header = TRUE, sep = "\t")
HF12N_S <- read.table("HF12N_S.txt", header = TRUE, sep = "\t")
HF12Rs_N_S <- read.table("HF12Rs_N_S.txt", header = TRUE, sep = "\t")
HF12S_Rs <- read.table("HF12S_Rs.txt", header = TRUE, sep = "\t")
HG64N_Rs <- read.table("HG64N_Rs.txt", header = TRUE, sep = "\t")
HG64N_S <- read.table("HG64N_S.txt", header = TRUE, sep = "\t")
HG64Rs_N_S <- read.table("HG64Rs_N_S.txt", header = TRUE, sep = "\t")
HG64S_Rs <- read.table("HG64S_Rs.txt", header = TRUE, sep = "\t")
# 合并差异表达和基因丰度
all <- HF12N_Rs %>%
full_join(HF12N_S) %>%
full_join(HF12Rs_N_S) %>%
full_join(HF12S_Rs) %>%
full_join(HG64N_Rs) %>%
full_join(HG64N_S) %>%
full_join(HG64Rs_N_S) %>%
full_join(HG64S_Rs) %>%
left_join(all_genes)
all[is.na(all)] <- "No"
## mean func
meanGp <- function(v) {
require('magrittr')
res <- v %>%
split(rep(1 : 8, each = 3)) %>%
sapply(mean, na.rm = TRUE)
return(res)
}
all_for_cluster <- select(all, -contains("vs")) # 只选择原始fpkm数据
rownames(all_for_cluster) <- all_for_cluster$id
all_for_cluster <- all_for_cluster[-1]
## sample name
sampleN <- c("HG64NRs","HG64SRs","HG64N","HG64S",
"HF12NRs","HF12SRs","HF12N","HF12S")
meanCount <- all_for_cluster %>%
apply(1, meanGp) %>%
t
colnames(meanCount) <- sampleN
## scale
scaleCount <- meanCount %>%
t %>%
scale %>%
t
scaleCount %<>% .[complete.cases(.), ]
2. 然后要先计算多少个cluster合适：
# set.seed(123) 聚类结果一致
set.seed(123)
## choose cluster num
## 1. sum of squared error
wss <- (nrow(scaleCount) - 1) * sum(apply(scaleCount, 2, var))
for (i in 2:20) {
wss[i] <- sum(kmeans(scaleCount,
centers=i,
algorithm = 'MacQueen')$withinss)
}
ggplot(tibble(k = 1:20, wss = wss), aes(k, wss)) +
geom_point(colour = '#D55E00', size = 3) +
geom_line(linetype = 'dashed') +
xlab('Number of clusters') +
ylab('Sum of squared error')
# ggsave("Sum_of_squared_error.pdf", height = 3, width = 4)
## 2. Akaike information criterion
kmeansAIC = function(fit){
m = ncol(fit$centers)
n = length(fit$cluster)
k = nrow(fit$centers)
D = fit$tot.withinss
return(D + 2*m*k)
}
aic <- numeric(20)
for (i in 1:20) {
fit <- kmeans(x = scaleCount, centers = i, algorithm = 'MacQueen')
aic[i] <- kmeansAIC(fit)
}
ggplot(tibble(k = 1:20, aic = aic), aes(k, wss)) +
geom_point(colour = '#009E73', size = 3) +
geom_line(linetype = 'dashed') +
xlab('Number of clusters') +
ylab('Akaike information criterion')
# ggsave("Akaike_information_criterion.pdf", height = 3, width = 4)
# choose cluster 10
withinss: Vector of within-cluster sum of squares, one component per cluster.
我们上面计算的是withinss，即cluster内的误差平方和，同一个cluster趋势越一致则越小。所以cluster越多则这个值越小，这和我们认知一致。但是，cluster太多了信息很杂，太少了就成了生拉硬扯成cluster了
image.png
可以看出选则10比较合适
利用另外一个参数：赤池信息量准则（Akaike information criterion，简称AIC）来衡量
AIC介绍，https://www.biaodianfu.com/aic-bic.html
理论上来讲，增加自由参数的数目可以提高拟合的优良性，但是过多，模型过于复杂容易造成过拟合现象。因此，AIC不仅要提高模型拟合度（极大似然），而且引入了惩罚项，使模型参数尽可能少，有助于降低过拟合的可能性。
image.png
可以看出，两种方法，同样结果。选择10
3. 然后聚类：
kclust10 <- kmeans(scaleCount, centers = 10, algorithm = "MacQueen", nstart = 1000, iter.max = 20)
cl <- as.data.frame(kclust10$cluster) %>%
rownames_to_column("id") %>%
set_colnames(c("id","cl"))
heat_all <- all %>% left_join(cl)
# 把所有DEGs的cluster信息保存，为后面富集分析所用
degs_cl <- heat_all %>%
select(c("id","cl"))
write.table(degs_cl, "./enrichment/degs_cl.txt", sep = "\t", quote = FALSE)
scaleC <- heat_all %>%
select(-contains("vs")) %>%
select(-c("id","cl")) %>%
t %>%
scale %>%
t %>%
as_tibble %>%
bind_cols(heat_all %>% select(id, cl))
# 排序
scaleC <- scaleC[,c("HG64S1", "HG64S2", "HG64S3",
"HG64N1", "HG64N2", "HG64N3",
"HG64SRs1", "HG64SRs2", "HG64SRs3",
"HG64NRs1", "HG64NRs2", "HG64NRs3",
"HF12S1", "HF12S2", "HF12S3",
"HF12N1", "HF12N2", "HF12N3",
"HF12SRs1", "HF12SRs2", "HF12SRs3",
"HF12NRs1", "HF12NRs2", "HF12NRs3",
"id", "cl")]
## col annotation
NatSoil <- HeatmapAnnotation(NatSoil = c(rep(rep(c("No", "Yes", "No", "Yes"), each = 3),2)),
col = list(NatSoil = c("Yes" = "grey", "No" = "white")),
gp = gpar(col = "black"))
Rs <- HeatmapAnnotation(Rs = c(rep(rep(c("No", "Yes"), each = 6),2)),
col = list(Rs = c("Yes" = "grey", "No" = "white")),
gp = gpar(col = "black"))
## DEG annotation
deg <- heat_all %>% select(matches("vs"))
# order
deg <- deg[,c("HG64N_vs_HG64S", "HF12N_vs_HF12S", "HG64NRs_vs_HG64SRs", "HF12RsN_vs_HF12RsS",
"HG64NRs_vs_HG64N", "HF12NRs_vs_HF12N", "HG64SRs_vs_HG64S", "HF12SRs_vs_HF12S")]
Heatmap(matrix = scaleC[1:24],
name = 'Scaled Counts',
row_split = scaleC$cl,
row_gap = unit(2, "mm"),
column_order = 1:24,
# column_split = rep(c("HG64S","HG64N","HG64SRs","HG64NRs",
# "HF12S","HF12N","HF12SRs","HF12NRs"), each = 3),
column_split = factor(rep(c("HG64","HF12"), each = 12), levels = c("HG64","HF12")),
show_column_names = FALSE,
col = colorRampPalette(rev(brewer.pal(n = 10, name = 'Spectral'))[c(-3,-8)])(100),
top_annotation = c(NatSoil, Rs),
row_title_gp = gpar(fontsize = 10),
use_raster = FALSE) +
Heatmap(deg,
col = c('Down' = '#00bbf9', 'No' = 'white', 'Up' = '#f94144'),
column_names_gp = gpar(fontsize = 6),
heatmap_legend_param = list(title = 'DEGs'),
cluster_columns = FALSE,
column_names_rot = 45,
use_raster = FALSE)
image.png
4. 然后每个cluster有什么功能呢？富集分析
#####################
setwd("C:/Users/yjk/Desktop/cluster_DEGs/enrichment")
library("clusterProfiler")
library("magrittr")
library("tidyverse")
library("RColorBrewer")
library("AnnotationHub")
my_theme()
# > packageVersion("AnnotationHub")
# [1] ‘3.0.2’
# packageVersion("clusterProfiler")
# [1] ‘4.0.5’
hub <- AnnotationHub()
# snapshotDate(): 2021-05-18
query(hub, "solanum")
# AnnotationHub with 9 records
# # snapshotDate(): 2021-05-18
# # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, Inparanoid8, WikiPathways
# # $species: Solanum lycopersicum, Solanum tuberosum, Solanum pennellii, Solanum pennelli, Sola...
# # $rdataclass: OrgDb, Inparanoid8Db, Tibble
# # additional mcols(): taxonomyid, genome, description, coordinate_1_based, maintainer,
# # rdatadateadded, preparerclass, tags, rdatapath, sourceurl, sourcetype
# # retrieve records with, e.g., 'object[["AH10593"]]'
#
# title
# AH10593 | hom.Solanum_lycopersicum.inp8.sqlite
# AH10606 | hom.Solanum_tuberosum.inp8.sqlite
# AH91815 | wikipathways_Solanum_lycopersicum_metabolites.rda
# AH94101 | org.Solanum_pennellii.eg.sqlite
# AH94102 | org.Solanum_pennelli.eg.sqlite
# AH94160 | org.Solanum_tuberosum.eg.sqlite
# AH94246 | org.Solanum_esculentum.eg.sqlite
# AH94247 | org.Solanum_lycopersicum.eg.sqlite
# AH94248 | org.Solanum_lycopersicum_var._humboldtii.eg.sqlite
sly.db <- hub[["AH94247"]]
这里如果遇到提示
hub <- AnnotationHub()
snapshotDate(): 2021-05-18
Warning message:DEPRECATION: As of AnnotationHub (>2.23.2), default caching location has changed.
Problematic cache: C:\Users\yjk\AppData\Local/AnnotationHub/AnnotationHub/Cache
See https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/TroubleshootingTheCache.html#default-caching-location-update
运行下面的命令即可
moveFiles<-function(package){
olddir <- path.expand(rappdirs::user_cache_dir(appname=package))
newdir <- tools::R_user_dir(package, which="cache")
dir.create(path=newdir, recursive=TRUE)
files <- list.files(olddir, full.names =TRUE)
moveres <- vapply(files,
FUN=function(fl){
filename = basename(fl)
newname = file.path(newdir, filename)
file.rename(fl, newname)
},
FUN.VALUE = logical(1))
if(all(moveres)) unlink(olddir, recursive=TRUE)
}
下面就可以进行GO和KEGG富集分析了
kmeansRes <- read.table("degs_cl.txt")
prefix <- 'kmeans10'
savepath <- "C:/Users/yjk/Desktop/cluster_DEGs/enrichment"
for (i in kmeansRes$cl %>% unique) {
## BP
goBP <- enrichGO(gene = kmeansRes %>% filter(cl == i) %>% .$id,
OrgDb = sly.db,
keyType= 'SYMBOL',
ont = 'BP',
pAdjustMethod = 'BH',
pvalueCutoff = 0.05,
qvalueCutoff = 0.1)
goBPSim <- clusterProfiler::simplify(goBP,
cutoff = 0.5,
by = 'p.adjust',
select_fun = min)
## check and plot
write.table(as.data.frame(goBPSim),
paste0(prefix, '_cl', i, '_cp_BP.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
## KEGG
kk2 <- enrichKEGG(gene = kmeansRes %>% filter(cl == i) %>% .$id %>%
bitr(.,"SYMBOL", "ENTREZID", sly.db) %>% dplyr::select("ENTREZID") %>%
unlist(),
organism = 'sly',
pvalueCutoff = 0.05)
write.table(as.data.frame(kk2),
paste0(prefix, '_cl', i, '_cp_KEGG.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
}
kall <- lapply(kmeansRes$cl %>% unique, function(x) {
eachG <- kmeansRes %>% filter(cl == x) %>% .$id
return(eachG)
}) %>%
set_names(kmeansRes$cl %>% unique %>% paste0('cl', .))
kallGOBP <- compareCluster(geneCluster = kall,
fun = 'enrichGO',
OrgDb = sly.db,
keyType= 'SYMBOL',
ont = 'BP',
pAdjustMethod = 'BH',
pvalueCutoff=0.01,
qvalueCutoff=0.1)
kallGOBPSim <- clusterProfiler::simplify(kallGOBP,
cutoff = 0.9,
by = 'p.adjust',
select_fun = min)
dotplot(kallGOBPSim, showCategory = 10) +
ggtitle("Biological process")
# ggsave("kallGOBPSim.pdf", width = 6.4, height = 5.4)
kallGOCC <- compareCluster(geneCluster = kall,
fun = 'enrichGO',
OrgDb = sly.db,
keyType= 'SYMBOL',
ont = "CC",
pAdjustMethod = 'BH',
pvalueCutoff=0.01,
qvalueCutoff=0.1)
kallGOCCSim <- clusterProfiler::simplify(kallGOCC,
cutoff = 0.9,
by = 'p.adjust',
select_fun = min)
dotplot(kallGOCCSim, showCategory = 20) +
ggtitle("Cellular component")
# ggsave("kallGOCCSim.pdf", width = 6.4, height = 4)
kallGOMF <- compareCluster(geneCluster = kall,
fun = 'enrichGO',
OrgDb = sly.db,
keyType= 'SYMBOL',
ont = "MF",
pAdjustMethod = 'BH',
pvalueCutoff=0.01,
qvalueCutoff=0.1)
kallGOMFSim <- clusterProfiler::simplify(kallGOMF,
cutoff = 0.9,
by = 'p.adjust',
select_fun = min)
dotplot(kallGOMFSim, showCategory = 10) +
ggtitle("Molecular function")
ggsave("kallGOMFSim.pdf", width = 6.9, height = 4)
kallGOBP %>%
as.data.frame %>%
write.table('kmeans10_GOBP.txt', quote = FALSE, sep = "\t")
emapplot(kallGOBP,
showCategory = 5,
pie='count',
pie_scale=1,
layout='kk')
kallKEGG_input <- lapply(kmeansRes$cl %>% unique, function(x) {
eachG <- kmeansRes %>% filter(cl == x) %>%
.$id %>%
bitr(.,"SYMBOL", "ENTREZID", sly.db) %>%
dplyr::select("ENTREZID") %>%
unlist()
return(eachG)
}) %>%
set_names(kmeansRes$cl %>% unique %>% paste0('cl', .))
kallKEGG <- compareCluster(geneCluster = kallKEGG_input,
fun = 'enrichKEGG',
organism = "sly",
pvalueCutoff = 0.05)
dotplot(kallKEGG)
# ggsave('kmeans10_KEGGALL.pdf', width = 8, height = 4)
kallKEGG %>%
as.data.frame %>%
write.table('kmeans10_KEGG.txt', quote = FALSE, sep = "\t")
列一个例子，GO富集的生物学过程
image.png
到此结束
网友评论