热图-标点

给热图表上标点symbol
代码如下

转换基因id
用clusterProfiler包 55.31% of input gene IDs are fail to map...
#用biomaRt包
library('curl')
library(org.Hs.eg.db)
listMarts()
plant<-useMart("ensembl")   
listDatasets(plant)
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
listFilters(mart)
hg_symbols<- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name", "start_position","end_position", "band"), filters= 'ensembl_gene_id', values = DEG$ENSEMBL, mart = mart)
names(DEG)[names(DEG) == "ENSEMBL"] <- "ensembl_gene_id"

DEG_symbolid <- na.omit(DEG_symbolid)
names(DEG)[names(DEG) == "ENSEMBL"] <- "ensembl_gene_id"
DEG_symbolid <- merge(DEG, hg_symbols, by="ensembl_gene_id")

input：

image.png

ggplot画火山图

library(ggplot2)
p <- ggplot(data = DEG_symbolid,
            aes(x = log2FoldChange,
                y = -log10(pvalue))) +
  geom_point(alpha=0.4, size=3.5,
             aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(0.01),lty=4,col="black",lwd=0.8) +
  theme_bw()
DEG_symbolid$change = ifelse(DEG_symbolid$pvalue < 0.000001 &
                               abs(log2(DEG_symbolid$log2FoldChange)) >= 1,
                             ifelse(log2(DEG_symbolid$log2FoldChange)> 1 ,'Up','Down'),
                             'Stable'))
DEG_symbolid$change = ifelse(DEG_symbolid$pvalue < 0.000001 & abs(log2(DEG_symbolid$log2FoldChange)) >= 1,
                     ifelse(log2(DEG_symbolid$log2FoldChange)> 1 ,'Up','Down'),
                     'Stable')

for_label <- DEG_symbolid %>%
  filter(abs(log2FoldChange) >3 & -log10(pvalue)> -log10(0.000001))
p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = hgnc_symbol),
    data = for_label,
    color="black"
  )

input： image.png

和文章里面的点，对不上号？？

image.png

查一下文章里面标出的点FOS基因

"FOS" %in% DEG_symbolid$hgnc_symbol
[1] TRUE
which(DEG_symbolid == "FOS", arr.ind = TRUE)
row col
12789 9455   8
DEG_symbolid[9455,]

>ensembl_gene_id log2FoldChange     lfcSE      stat    pvalue      padj baseMean hgnc_symbol chromosome_name start_position
12789 ENSG00000170345     0.04420872 0.5116253 0.0864084 0.9311418 0.9998393 51.98664         FOS              14       75278826
      end_position  band change
12789     75282230 q24.3 Stable

查出来的pvalue是，0.9311418，应该是不会标出来的呀？

不知道是哪里出问题的了？

①、在转换基因的时候，转换错了。但是百度了一下，ENSG00000170345 ，确实是对应FOS基因的。

②、那就可能是在做差异分析的时候做错了。