指定通路绘制gsea图热图和火山图
前面在 所有的肿瘤都有恶性增殖的特性吗，我们发现了绝大部分癌症都有MKI67和TOP2A这样的细胞增殖通路相关基因的高表达，最后的gsea分析结果里面展示的通路包括：
2.4 Replication and repair
03030 DNA replication
03410 Base excision repair
03420 Nucleotide excision repair
03430 Mismatch repair
03440 Homologous recombination
03450 Non-homologous end-joining
03460 Fanconi anemia pathway
4.2 Cell growth and death
04110 Cell cycle
但是我们直接是对gsea分析结果的最终es值在可视化，所以是行是通路，列是癌症的，数值是gsea的es打分的矩阵。对初学者来说， 跳过了大量细节，所以跟这个教程会比较吃力，有粉丝就提问了希望可以对这些通路在在具体的癌症里面细化展示，比如绘制gsea图，热图和火山图。
gsea分析大家都没有问题：
nrDEG=DEG_limma_voom[,c(1,4)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
head(nrDEG)
library(org.Hs.eg.db)
library(clusterProfiler)
gene <- bitr(rownames(nrDEG), fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
gene$logfc <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
head(gene)
geneList=gene$logfc
names(geneList)=gene$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)
R.utils::setOption( "clusterProfiler.download.method",'auto' )
library(clusterProfiler)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
tmp=kk@result
pro='test_gsea'
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
对上面的结果，进行查看统计学显著的上下调通路。
# 这里如果找不到显著下调的通路，可以选择调整阈值，或者其它。
kk_gse=kk
down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.5,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.5,];up_kegg$group=1
save(up_kegg,kk,file = 'up_kegg.by.gsea.Rdata')
首先批量针对每个通路绘制gsea图：
pro='gsea'
lapply(1:nrow(up_kegg), function(i){
gseaplot2(kk,up_kegg$ID[i],
title=up_kegg$Description[i],pvalue_table = T)
ggsave(paste0(pro,'_up_kegg_',
gsub('/','-',up_kegg$Description[i]),
'.pdf'))
})
然后 批量针对每个通路绘制热图，需要提取每个通路里面的基因列表：
lapply(1:nrow(up_kegg), function(i){
cg = strsplit(up_kegg$core_enrichment[i],'/')[[1]]
dat=log2(edgeR::cpm(exprSet)+1)
dat[1:4,1:4]
table(group_list)
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = T,
annotation_col=ac)
pheatmap(n,show_colnames =F,show_rownames = T,
annotation_col=ac,
filename = paste0(pro,'_up_kegg_pheatmap_',
gsub('/','-',up_kegg$Description[i]),
'.pdf'))
})
然后 批量针对每个通路绘制火山图，把每个通路里面的基因列表标记在火山图里面，这个时候仍然是分成两步走，首先绘制一个火山图 (不同的包做差异分析得到的矩阵列名不一样，下面是DEseq2的结果举例哦 )：
## for volcano
logFC_cutoff <- with(DEG_DEseq2,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) );logFC_cutoff
logFC_cutoff
logFC_t=3
DEG_DEseq2$change = as.factor(ifelse(DEG_DEseq2$pvalue < 0.05 & abs(DEG_DEseq2$log2FoldChange) > logFC_t,
ifelse(DEG_DEseq2$log2FoldChange > logFC_t ,'up','down'),'not')
)
table(DEG_DEseq2$change)
library(ggplot2)
library(ggstatsplot)
library(ggsci)
p_v = ggplot(data=DEG_DEseq2,
aes(x=log2FoldChange, y=-log10(pvalue),
color=change)) +
geom_point(alpha=0.4, size=1.75) +
xlab("logFC") +
ylab("-log10 pvalue") +
theme_ggstatsplot()+
geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) +
theme(plot.title = element_text(size=12,hjust = 0.5),
panel.grid = element_line(colour = 'white'))+
scale_color_manual(values=c("#34bfb5", "#828586","#ff6633")
)
dev.off()
p_v
ggsave(p_v,filename = 'step4_deg_volcano.png',height = 6,width = 5)
然后循环感兴趣的每个 通路，去叠加到上面的火山图里面 ：
lapply(1:nrow(up_kegg), function(i){
cg = strsplit(up_kegg$core_enrichment[i],'/')[[1]]
cgd = DEG_DEseq2[cg,]
head(cgd)
cgd$gene = rownames(cgd)
p_v + ggrepel::geom_text_repel(data=cgd,
aes(x=log2FoldChange,
y=-log10(padj),
label= gene) ,
min.segment.length = 0,
max.overlaps=30
)
ggsave(paste0(pro,'_up_kegg_volcano_',
gsub('/','-',up_kegg$Description[i]),
'.pdf'))
})
我们来展示一下前面在 所有的肿瘤都有恶性增殖的特性吗提到的细胞周期通路看看。代码如下所示：
pro='BRCA'
brca_gsea = gsea_list[[pro]]
brca_gsea = brca_gsea [match(cg,brca_gsea$Description),]
brca_gsea
brca_deg = deg_list[[pro]]
load('Rdata/TCGA-BRCA.htseq_counts.Rdata')
group_list = ifelse(
as.numeric(substring(colnames(pd_mat),14,15)) < 10,
'tumor','normal'
)
table(group_list)
lapply(1:nrow(brca_gsea), function(i){
cg = strsplit(brca_gsea$core_enrichment[i],'/')[[1]]
dat=log2(edgeR::cpm(pd_mat)+1)
dat[1:4,1:4]
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = T,
annotation_col=ac)
pheatmap(n,show_colnames =F,show_rownames = T,
annotation_col=ac,
filename = paste0(pro,'_up_kegg_pheatmap_',
gsub('/','-',brca_gsea$Description[i]),
'.pdf'))
})
可以看到，虽然一千多个肿瘤样品跟一百多个正常样品进行差异分析，确实细胞增殖相关的基因是统计学显著的高表达，但也有异质性，并不是所有的肿瘤个体都是恶性增殖，其实我们推翻了前面的结论。我们前面在 所有的肿瘤都有恶性增殖的特性吗，得到的结论其实是绝大部分肿瘤从整体上来说，恶性增殖都是大概率事件：
肿瘤不同病人基因异质性
好比我们说北京人都很有钱，并不是说每个北京人都是富人，只不过是北京这个地区相比其它城市来说，整体上富裕的人比较多而已。我们说北京的高考比较容易，也不是说每个人都能上清华北大，其内部也需要竞争。同理，肿瘤确实是有一个很显著的特征就是恶性增殖，但是并不是每个肿瘤类型的每个肿瘤样品都是如此。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2022-06-26，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
腾讯云开发者
