getwd()

#下载方式有很多，挑自己能下载的办法下载好几个包就好，我们用自带的数据就好

install.packages("SeuratData")
devtools::install_github('satijalab/seurat-data')
library(SeuratData)
library(Seurat)
library(airway)
library(magrittr)
#加载火山图的包
install.packages('BiocManager')
#install.packages("EnhancedVolcano")
devtools::install_github('kevinblighe/EnhancedVolcano')
library(EnhancedVolcano)
#控制数据大小(读取时间）)
getOption('timeout')
options(otimeout=10000)

#载入PBMC3k数据，已经是一个seurat结构
#InstallData("pbmc3k")
data(pbmc3k)
pbmc3k
str(pbmc3k)
head(pbmc3k,5)#有一个注释群在

#final是处理过的pbmc3k数据
sce <- pbmc3k.final
#展示数据
sce
head(sce,5)#这个东西被降维，计算线粒体基因等一系列操作了
table(Idents(sce))#展示一下ident
pbmc3k@meta.data$seurat_annotations #展示一下注释

#寻找两个群之间的高变基因
deg = FindMarkers(sce,ident.1 = 'NK', ident.2 = 'B')
head(deg[order(deg$p_val),])#根据P值排序

这个包的内容非常丰富，有时间再多整理出来吧
#lab=rownames 设置标签（gene名字）
#xlim = c(-17, 13)、ylim=c(1,6)：设置x、y轴的展示区间
#pCutoff = 0.001、FCcutoff = 2：自定义阈值线划分的范围，也就是那两个虚线。
#cutoffLineCol  自定义阈值线的颜色
#cutoffLineType  自定义阈值线的类型
#cutoffLineWidth    自定义阈值线的宽度
#pointSize 点的大小默认是2
#shape=19
#subtitle 是title下面那个title
#caption 是右下角的title
#titleLabSize,subtitleLabSize,captionLabSize三个标题的大小，默认是14
EnhancedVolcano(deg,
                lab = rownames(deg),
                #xlim = c(-5,7),
                #ylim = c(0,100),
                #pCutoff = 0.01,
                #FCcutoff = 4,
                x = 'avg_log2FC',
                y = 'p_val_adj',
                #subtitle = F,
                #caption = "what"，
                #cutoffLineCol = "red",
                #titleLabSize = 14,
                #shape=16,
                #pointSize=5,
                title = "NK vs B",)

deg$p_val_adj
deg$avg_log2FC

#用ggplo去画图，同一组数据
deg$threshold<-as.factor(ifelse(deg$avg_log2FC >= 2,'Up',ifelse(deg$p_val_adj<0.05 & deg$avg_log2FC <= -2,'Down','Not')))
deg<-data.frame(deg)
ggplot(data=deg, aes(x=avg_log2FC, y=-log10(p_val_adj), colour=threshold, fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(alpha=0.6,size=2) +
  #xlim(c(-6, 6)) +
  #ylim(c(0, 300)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-2,2),lty=1,col=c("blue","red"),lwd=0.6)+
  geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)+
  theme(legend.position="right",
        panel.grid=element_blank(),
        legend.title = element_blank(),
        legend.text= element_text(face="bold", color="black",family = "Times", size=8),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(color="black", size=12),
        axis.text.y = element_text(color="black", size=12),
        axis.title.x = element_text(face="bold", color="black", size=12),
        axis.title.y = element_text(face="bold",color="black", size=12))+
  labs(x="log2 (Fold Change)",y="-log10 (p-value)",title="NK vs B")

 

 ggplot封装的东西比较少，还是没那么方便