R语言编写卡住了，想要得到火山图

怎么修改，写不下去了，目的是想要分析数据得到火山图

rm(list = ls())
rm(list = ls())
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("GEOquery")

library("GEOquery")
gse<-getGEO('GSE147602',destdir='.',
            AnnotGPL= T,
            getGPL= T)
gse[[1]]
exp<-exprs(gse[[1]])
cli<-pData(gse[[1]])

group<-c(rep("cancer",10),rep("control",10),rep("cancer",10),rep("control",10))
group

exp<-as.data.frame(exp)
exp$ID<-rownames(exp)
GPL<-fData(gse[[1]])

library("stringr")
probe2gene<- GPL[,c(1,9)]
probe2gene$ID<-rownames(probe2gene)
probe2gene$symbol=trimws(str_split(probe2gene$gene_assignment,'//',simplify = T)[,50])
plot(table(table(probe2gene$symbol)),xlim=c(1,50))
head(probe2gene)
eSet<-probe2gene[,c(1,3)]

exp_symbol<-merge(exp,eSet,by="ID")
table(duplicated(exp_symbol$symbol))

library("limma")
exp_unique<-avereps(exp_symbol[,-c(1,ncol(exp_symbol))],ID=exp_symbol$symbol)
table(duplicated(rownames(exp_unique)))


boxplot(exp_unique)
options(stringsAsFactors = F)
group<-c(rep("cancer",10),rep("control",10),rep("cancer",10),rep("control",10))
group

library("limma")
design<-(model.matrix(~0+factor(group)))
colnames(design)<-levels(factor(group))
row.names(design)<-colnames(exp_unique)
exp_unique1<-voom(exp_unique,design,normalize="quantile") 

contrast.matrix<-makeContrasts(cancer-control,levels=design)



fit<-lmFit(exp_unique,design)
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
options(digits = 4)
DEG<-topTable(fit2,coef = 1,n=Inf)
DEG$group<-ifelse(DEG$P.Value>0.05,"no_change",
                  ifelse(DEG$logFC>1,"up",
                         ifelse(DEG$logFC< -1,"down","no_change")))
table(DEG$group)
DEG$gene<-rownames(DEG)
write.table(DEG, file = "DEG163654.txt", row.names = T, quote = T) # 


BiocManager::install("ggplot2")
BiocManager::install("ggpubr")
BiocManager::install("ggrepel")
BiocManager::install("pheatmap")



library(pheatmap)
library(ggplot2)
library(ggpubr)

data<-data.frame(gene=exp_unique["Timp1",],group=group)
ggplot(data=data,aes(x=group,y=gene,fill=group))+geom_boxplot()+
  geom_point(position = "jitter")+stat_compare_means()
火山图
library(ggrepel)
DEG$P<- -log10(DEG$P.Value)
ggplot(data=DEG,aes(x=logFC,y=P,color=group))+geom_point()
ggplot(data=DEG,aes(x=logFC,y=P,color=group))+
  geom_point(data=DEG[DEG$P.Value<0.01&abs(DEG$logFC)>1,],size=3)+
  geom_point(data=DEG[DEG$P.Value>0.01|abs(DEG$logFC)<1,],size=1)+
  scale_color_manual(values = c("limegreen","grey64","red"))+
  xlab("log2Foldchange")+ylab("-log10(p.value)")+
  geom_text_repel(data=DEG[DEG$P.Value<0.001&abs(DEG$logFC)>1.5,],
                  aes(label=gene),
                  size=3,color="black")

library("pheatmap")

X<-DEG$logFC
names(X)<-rownames(DEG)
X[1:8]
upgene<-names(tail(sort(X),8))
downgene<-names(head(sort(X),8))
top20<-c(upgene,downgene)

group<-data.frame(group=group)
rownames(group)<-colnames(exp_unique)

pheatmap(exp_unique[top20,],
         show_rownames =T,
         cluster_rows = F,
         cluster_cols = F,
         show_colnames = F,
         scale = "row",
         annotation_col = group)

pheatmap(exp_unique[top20,],
         show_rownames =T,
         cluster_rows = F,
         cluster_cols = F,
         show_colnames = F,
         scale = "row",
         annotation_col = group,
         colors=c("grey92","skyblue","gold1","orange"))

save(contrast.matrix,DEG,design,eSet,fit,fit2,cli,exp,GPL,exp_unique,exp_symbol,probe2gene,group,gse,file="GSE147602.Rdata")

这是包安装发生错误，viridisLite这个包是麻烦，用回ggplot2画吧

不知道你这个问题是否已经解决, 如果还没有解决的话:

• 你可以参考下这篇文章：R 数据可视化 02 | 火山图
• 除此之外, 这篇博客: R语言绘制表达矩阵的热图与火山图中的 火山图 部分也许能够解决你的问题, 你可以仔细阅读以下内容或者直接跳转源博客中阅读:
  

  仅需要pval 和log2FoldChange

  v_circ$color <- ifelse(v_circ$pval <0.05 & abs(v_circ$log2FoldChange)>= 1,ifelse(v_circ$log2FoldChange > 1,'red','blue'),'gray')
  color <- c(red = "red",gray = "gray",blue = "blue")
  
  circ_vo <- ggplot(v_circ, aes(log2FoldChange, -log10(pval), col = color)) +
    geom_point() +
    theme_bw() +
    scale_color_manual(values = color) +
    labs(x="log2 (fold change)",y="-log10 (p-value)") +
    geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
    geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
    theme(legend.position = "none",
          panel.grid=element_blank(),
          axis.title = element_text(size = 16),
          axis.text = element_text(size = 14))
  
  circ_vo + ggtitle('circRNA')
  circ_vo
  

• 以下回答来自chatgpt:

  根据参考资料中提到的，绘制火山图需要数据类型为data.frame，包含基因名称，PValue、FD等数据，并且需要使用ggplot2包进行绘图。下面给出一个简单的步骤，更改plot的每个参数，以获得更美观的火山图。

  1.导入数据

  首先需要将数据导入R中，可以使用read.table()函数进行导入。假设数据文件名为data.txt，数据中包含三列:基因名称，PValue和FD，可以使用以下代码导入数据：

  data <- read.table("data.txt", header = T, sep = "\t")
  

  2.绘制基础火山图

  使用ggplot2包进行绘图，需要先调用该包。使用aes()函数指定x轴和y轴的变量，geom_point()指定绘制点的形状、大小和颜色等信息。具体代码如下：

  library(ggplot2)
  
  ggplot(data, aes(x = log2(FD), y = -log10(PValue))) +
    geom_point(shape = 8, size = 1, color = "black")
  

  3.按照调控模式上色

  可以按照调控模式（如上调、下调和正常）将点按照不同颜色进行上色。具体代码如下：

  ggplot(data, aes(x = log2(FD), y = -log10(PValue))) +
    geom_point(aes(color = regulate), shape = 8, size = 1) +
    scale_color_manual(values = c("red", "black", "blue"), 
                       labels = c("Up", "Normal", "Down"))
  

  其中color参数指定按照什么变量进行上色，scale_color_manual()是手动指定颜色和标签。

  4.调整坐标轴和标签

  可以使用theme()和labs()函数进行调整。其中，theme()函数可以设定与图片背景，文字有关的参数。具体代码如下：

  ggplot(data, aes(x = log2(FD), y = -log10(PValue))) +
    geom_point(aes(color = regulate), shape = 8, size = 1) +
    scale_color_manual(values = c("red", "black", "blue"), 
                       labels = c("Up", "Normal", "Down")) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.border = element_rect(colour = "black"),
          axis.line = element_line(colour = "black"),
          axis.text = element_text(size = 14),
          axis.title = element_text(size = 14),
          legend.position = "right",
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 16)) +
    labs(x = expression(Log[2]*" FC"), y = expression("-log"[10]*" P Value"))
  

  5.添加参考线和标题

  可以使用geom_hline()和geom_vline()函数添加水平和垂直参考线，使用ggtitle()函数添加标题。具体代码如下：

  ggplot(data, aes(x = log2(FD), y = -log10(PValue))) +
    geom_point(aes(color = regulate), shape = 8, size = 1) +
    scale_color_manual(values = c("red", "black", "blue"), 
                       labels = c("Up", "Normal", "Down")) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.border = element_rect(colour = "black"),
          axis.line = element_line(colour = "black"),
          axis.text = element_text(size = 14),
          axis.title = element_text(size = 14),
          legend.position = "right",
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 16)) +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    geom_vline(xintercept = c(-2, 2), linetype = "dashed") +
    labs(x = expression(Log[2]*" FC"), y = expression("-log"[10]*" P Value")) +
    ggtitle("Volcano Plot")
  

  6.最终代码

  综上所述，火山图的绘制代码如下：

  library(ggplot2)
  
  data <- read.table("data.txt", header = T, sep = "\t")
  
  ggplot(data, aes(x = log2(FD), y = -log10(PValue))) +
    geom_point(aes(color = regulate), shape = 8, size = 1) +
    scale_color_manual(values = c("red", "black", "blue"), 
                       labels = c("Up", "Normal", "Down")) +
    theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.border = element_rect(colour = "black"),
          axis.line = element_line(colour = "black"),
          axis.text = element_text(size = 14),
          axis.title = element_text(size = 14),
          legend.position = "right",
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 16)) +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    geom_vline(xintercept = c(-2, 2), linetype = "dashed") +
    labs(x = expression(Log[2]*" FC"), y = expression("-log"[10]*" P Value")) +
    ggtitle("Volcano Plot")
  

  注：以上代码仅供参考，实际应用需要根据数据情况进行调整。

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^