在每個(gè)cluster注釋好細(xì)胞類型后背稼，需要找出每個(gè)細(xì)胞種類的marker基因并進(jìn)行可視化

讀取數(shù)據(jù)宫纬、創(chuàng)建seurat對(duì)象并進(jìn)行錨點(diǎn)整合

##系統(tǒng)報(bào)錯(cuò)改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())

##改變工作路徑强法，將相關(guān)結(jié)果儲(chǔ)存至對(duì)應(yīng)文件夾
setwd("C:/Users/86269/Desktop/shun.C/single_cell")

#加載所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)

#讀取文件,創(chuàng)建Seurat對(duì)象
dir=c("CAF1/","CAF2/","dapi1/","dapi2/")
names(dir)=c("CAF1","CAF2","dapi1","dapi2")

scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  scRNAlist[[i]] <- CreateSeuratObject(counts,min.cells=3,min.features=300)
}

#歸一化處理與尋找高變基因
for(i in 1:length(scRNAlist)){
  scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]])
  scRNAlist[[i]] <- FindVariableFeatures(scRNAlist[[i]],selection.method='vst',nfeatures=3000)
}

#進(jìn)行錨點(diǎn)整合
library(future)
scRNA.anchor <- FindIntegrationAnchors(object.list = scRNAlist,anchor.features = 3000)
scRNA1 <- IntegrateData(anchorset = scRNA.anchor)

DefaultAssay(scRNA1) <- "integrated"

降維聚類 進(jìn)行SingleR注釋

#尋找高變基因万俗，降維聚類
scRNA1 <- ScaleData(scRNA1)
scRNA1 <- FindVariableFeatures(scRNA1,Selection.method='vst',nFeatures=3000) 
scRNA1 <- RunPCA(scRNA1,npcs=20,verbose=T)
scRNA1 <- FindNeighbors(scRNA1,reduction='pca',dims=1:7) 
scRNA1 <- FindClusters(scRNA1,resolution=0.4) 
scRNA1 <- RunUMAP(scRNA1,resolution='pca',dims=1:7)
scRNA1<- RunTSNE(scRNA1,dims=1:7)
DimPlot(scRNA1,reduction = 'umap')
#運(yùn)用SingleR進(jìn)行注釋
library(SingleR)
#準(zhǔn)備參考數(shù)據(jù)與注釋數(shù)據(jù)
library(celldex)
refdata <- MouseRNAseqData() 
testdata <- GetAssayData(scRNA1,slot="data")


clusters <- scRNA1@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata,ref = refdata,labels = refdata$label.main,
                    method='cluster',clusters = clusters,
                    assay.type.test = "logcounts",assay.type.ref = 'logcounts')
celltype <- data.frame(ClusterID=rownames(cellpred),celltype=cellpred$labels,stringsAsFactors = T)



scRNA1@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
  scRNA1@meta.data[which(scRNA1@meta.data$seurat_clusters==celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]
}



DimPlot(scRNA1,reduction = 'umap',group.by = 'celltype',label=T)
save(scRNA1,file="scRNA1.Rdata")

image.png

該圖展示了5種細(xì)胞種類的聚類情況

展示marker

1.heatmap

ROC方式尋找marker基因

#更改分組信息
Idents(scRNA1)='celltype'
#找到不同細(xì)胞類型的marker
degs <- FindAllMarkers(scRNA1,logfc.threshold = 0.5,
                       test.use = 'roc',
                       return.thresh = 0.25,
                       min.pct = 0.3,
                       only.pos = T)


#對(duì)marker進(jìn)行篩選
degs_sig <- degs %>% filter(pct.1>0.3& power>0.25) %>% filter(cluster != 'others') %>% arrange(cluster,-power)

#篩選出top50
 degs_top50 <- degs_sig %>% group_by(cluster) %>% top_n(50,power) %>% top_n(50,avg_diff) %>% arrange(cluster,-power)
 
 #取平均值,對(duì)每個(gè)基因在每種細(xì)胞內(nèi)的表達(dá)量取平均值，例如A基因在巨噬細(xì)胞內(nèi)的表達(dá)平均值與上皮細(xì)胞內(nèi)的表達(dá)平均值
 avgData <- scRNA1@assays$RNA@data[degs_top50$gene,] %>% apply(1,function(x){
   tapply(x,scRNA1$celltype,mean)
 }) %>% t

 #對(duì)最大值與最小值進(jìn)行scale與限定,最大值與最小值偏離過大會(huì)對(duì)畫圖效果造成影響
 phData <- MinMax(scale(avgData),-2,2)
 
 #可視化
 library(pheatmap)
 rownames(phData)<- 1:nrow(phData)
 phres <- pheatmap(
   phData,
   color = colorRampPalette(c("darkblue","white","red3"))(99),
   scale='row',
   cluster_rows=F,
   cluster_cols=T,#不按行聚類饮怯，按列聚類
   clustering_method='complete',
   show_rownames=F,
   annotation_row = data.frame(cluster=degs_top50$cluster)
 )

image.png

Wilcoxon方法尋找marker基因

Idents(scRNA1) ="celltype"
 marker = FindAllMarkers(scRNA1,logfc.threshold = 0.5)
 top10 <- marker %>% group_by(cluster) %>% top_n(n=10,wt = avg_log2FC)
 
 DimPlot(scRNA1,reduction="tsne",group.by = 'celltype',pt.size = 1)
 DoHeatmap(scRNA1,features = top10$gene,label = F,slot="data")

image.png

heatmap

可以觀察到熱圖各個(gè)列的寬度差別很大闰歪，即各個(gè)種類的細(xì)胞數(shù)量差別很大
接下來每種細(xì)胞隨機(jī)抽取150個(gè)再進(jìn)行熱圖繪制

scRNA1 <- subset(scRNA1,downsample=150)
DoHeatmap(scRNA1,features = top10$gene,label = F,slot="data")

heatmap

自定義方法畫熱圖

#將barcode加入metadata列中
colanno <- scRNA1@meta.data
#按照細(xì)胞類型排序
colanno=colanno%>%arrange(celltype)

#將細(xì)胞類型與barcode一一對(duì)應(yīng)
colanno$celltype =factor(colanno$celltype,levels=unique(colanno$celltype))
colanno1 <- colanno[,6]
colanno1<- as.data.frame(colanno1)
rownames(colanno1)=rownames(colanno)
colnames(colanno1)<-'celltype'

heatmap

2.violin plot

library(reshape2)
#每個(gè)細(xì)胞種類挑選兩個(gè)marker基因
top2 <- marker %>% group_by(cluster) %>% top_n(n=2,wt = avg_log2FC)
#提取10個(gè)marker基因的表達(dá)矩陣
Vln.df <- as.data.frame(scRNA1[["RNA"]]@data[top2$gene[1:10],])
Vln.df$gene=rownames(Vln.df)
Vln.df=melt(Vln.df,id='gene')
colnames(Vln.df)[c(2:3)]=c('barcode','exp')

#加入metadata中的信息
anno=scRNA1@meta.data
anno$barcode=rownames(anno)
Vln.df = inner_join(Vln.df,anno,by='barcode')

#可視化
Vln.df %>% ggplot(aes(celltype,exp))+geom_violin(aes(fill=celltype),scale="width")+
  facet_grid(Vln.df$gene~.,scales = "free_y")+
  scale_fill_brewer(palette = "Set3",direction = 1)+
  scale_x_discrete("")+scale_y_continuous("")+
  theme_bw()+
  theme(
    axis.text.x.bottom = element_text(angle=45,hjust=1,vjust=1),
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),#移除網(wǎng)格線
    legend.position = "none"
  )

violin plot

如果想將顏色調(diào)整為與基因的表達(dá)量相關(guān)，則代碼如下

#將細(xì)胞類型與基因連接組成新的列
Vln.df$celltype_gene=paste(Vln.df$celltype,Vln.df$gene,sep = "_")
#計(jì)算每個(gè)marker基因在每種細(xì)胞內(nèi)的平均表達(dá)量
stat.df=as.data.frame(Vln.df%>%  dplyr::group_by(celltype,gene)%>%dplyr::summarize(mean=mean(exp)))
stat.df$celltype_gene=paste(stat.df$celltype,stat.df$gene,sep = "_")
stat.df=stat.df[,c("mean","celltype_gene")]
#加上其他信息
Vln.df=inner_join(Vln.df,stat.df,by="celltype_gene")
#控制平均數(shù)范圍防止出現(xiàn)極值
Vln.df$mean=ifelse(Vln.df$mean > 3, 3, Vln.df$mean)

Vln.df%>%ggplot(aes(celltype,exp))+geom_violin(aes(fill=mean),scale = "width")+
  facet_grid(Vln.df$gene~.,scales = "free_y")+
  scale_fill_gradient(limits=c(0,3),low = "lightgrey",high = "red")+
  scale_x_discrete("")+scale_y_continuous("",expand = c(0.02,0))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
    axis.text.x.bottom = element_text(angle = 45,hjust = 1,vjust = 1)
  )

violin plot

3.bubble plot

#提取10個(gè)marker基因的表達(dá)矩陣
bubble.df=as.matrix(scRNA1[['RNA']]@data[top2$gene[1:10],])
bubble.df = t(bubble.df)

bubble.df

#添加上其他信息
bubble.df=as.data.frame(scale(bubble.df))
bubble.df $CB=rownames(bubble.df)
bubble.df=merge(bubble.df,scRNA1@meta.data[,c(1,6)],by.x="CB",by.y=0)
bubble.df$CB=NULL

bubble.df

celltype_v=c()
gene_v=c()
mean_v=c()
ratio_v=c()
#兩層循環(huán)蓖墅，第一層針對(duì)celltype库倘，第二層針對(duì)該celltype內(nèi)的marker基因
for(i in unique(bubble.df$celltype)){
  bubble.df_small=bubble.df %>% filter(celltype==i)
  for(j in top2$gene[1:10]){
    exp_mean=mean(bubble.df_small[,j])
    exp_ratio=sum(bubble.df_small[,j] > min(bubble.df_small[,j]))/length(bubble.df_small[,j])
    celltype_v=append(celltype_v,i)
    gene_v=append(gene_v,j)
    mean_v=append(mean_v,exp_mean)
    ratio_v=append(ratio_v,exp_ratio )
    }

  plotdf=data.frame(
  celltype=celltype_v,
  gene=gene_v,
  exp=mean_v,
  ratio=ratio_v
)
 #可視化
library(RColorBrewer)
plotdf$celltype=factor(plotdf$celltype,levels=sort(unique(plotdf$celltype)))
plotdf$gene=factor(plotdf$gene,levels=rev(as.character(top10$gene[1:10])))
#限定基因表達(dá)量最高值
plotdf$exp=ifelse(plotdf$exp>3,3,plotdf$exp)

plotdf %>% ggplot(aes(celltype,gene,size=ratio,color=exp))+geom_point()+
  scale_x_discrete("")+scale_y_discrete("")+
  scale_color_gradientn(colours = rev(c('#FFD92F','#FEE391',brewer.pal(11,"Spectral")[7:11])))+
  scale_size_continuous(limits=c(0,1))+theme_classic()+
  theme(
    axis.text.x.bottom = element_text(hjust=1,vjust=1,angle=90)
  )

image.png

4.feature plot

以基因Rbp1為例，畫TSNE降維結(jié)果的聚類圖

#提取Rbp1基因在每個(gè)細(xì)胞中的表達(dá)量
mat1 = as.data.frame(scRNA1[['RNA']]@data['Rbp1',])
colnames(mat1)='exp'
#提取降維結(jié)果并與Rbp1表達(dá)量合并
mat2=Embeddings(scRNA1,'tsne')
mat3=merge(mat2,mat1,by='row.names')
#可視化
mat3 %>% ggplot(aes(tSNE_1,tSNE_2))+geom_point(aes(color=exp))+
  scale_color_gradient(low="grey",high="purple")+theme_bw()+
  theme(
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
    axis.ticks=element_blank(),axis.text = element_blank(),
    legend.position = 'none',
    plot.title=element_text(hjust=0.5,size=14)
  )+scale_x_continuous("")+scale_y_continuous("")

feature plot

點(diǎn)的顏色越深论矾，基因表達(dá)量越高

5.H形圖

#加載數(shù)據(jù)并進(jìn)行聚類與注釋
##系統(tǒng)報(bào)錯(cuò)改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())

library(dplyr)
library(Seurat)
library(patchwork)
library(reshape2)
library(RColorBrewer)
library(ggplot2)
library(ggrepel)
library(magrittr)
library(data.table)
library(tidyverse)

setwd("C:/Users/86269/Desktop/shun.C/single_cell")
load("scRNA_harmony.Rdata")
#進(jìn)行SingleR注釋
library(SingleR)
refdata <- HumanPrimaryCellAtlasData()
testdata <- GetAssayData(scRNA_harmony, slot="data")
clusters <- scRNA_harmony@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main, 
                    method = "cluster", clusters = clusters, 
                    assay.type.test = "logcounts", assay.type.ref = "logcounts")
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = FALSE)
scRNA_harmony@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}

DimPlot(scRNA_harmony,group.by = "celltype")

#得到細(xì)胞種類
sc.combined=scRNA_harmony
table(sc.combined@meta.data$celltype)

#找出差異基因
type=c("Chondrocytes","Endothelial_cells","Fibroblasts",
       "Macrophage","Monocyte","T_cells",
       "Tissue_stem_cells")
library(future)
r.deg=data.frame()
for (i in 1:7) {
  Idents(sc.combined)="celltype"
  #每個(gè)細(xì)胞種類在sample11相對(duì)于sample3的差異基因
  deg=FindMarkers(sc.combined,ident.1 = "sample_11",ident.2 = "sample_3",
                  group.by = "orig.ident",subset.ident =type[i]   )
  #在deg表中加入細(xì)胞類型
  deg$celltype=type[i]
  #在deg表中加入cluster數(shù)
  deg$unm=i-1
  #合并到一個(gè)表
  r.deg=rbind(deg,r.deg) 
}

r.deg

篩選差異基因

#作初步篩選教翩，篩選出差異明顯的基因
r.deg <- subset(r.deg, p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)
#找出上調(diào)和下調(diào)基因
r.deg$threshold <- as.factor(ifelse(r.deg$avg_log2FC > 0 , 'Up', 'Down'))
#找出顯著變化與不顯著變化基因
r.deg$adj_p_signi <- as.factor(ifelse(r.deg$p_val_adj < 0.01 , 'Highly', 'Lowly'))
#合并得到基因的調(diào)整情況
r.deg$thr_signi <- paste0(r.deg$threshold, "_", r.deg$adj_p_signi)
r.deg$unm %<>% as.vector(.) %>% as.numeric(.)

r.deg

挑選log2FC為top5的基因進(jìn)行展示

#挑選上調(diào)基因(每個(gè)cluster上調(diào)最明顯的5個(gè)基因)
top_up_label <- r.deg %>% 
  subset(., threshold%in%"Up") %>% 
  group_by(unm) %>% 
  top_n(n = 5, wt = avg_log2FC) %>% 
  as.data.frame()
#挑選下調(diào)基因（每個(gè)cluster下調(diào)最明顯的5個(gè)基因）
top_down_label <- r.deg %>% 
  subset(., threshold %in% "Down") %>% 
  group_by(unm) %>% 
  top_n(n = -5, wt = avg_log2FC) %>% 
  as.data.frame()
#合并得到差異明顯的基因
top_label <- rbind(top_up_label,top_down_label)
top_label$thr_signi %<>% factor(., levels = c("Up_Highly","Down_Highly","Up_Lowly","Down_Lowly"))

畫圖準(zhǔn)備

### 準(zhǔn)備繪制暗灰色背景所需數(shù)據(jù) 
#灰色背景的高度
background_position <- r.deg %>% group_by(unm) %>%summarise(Min = min(avg_log2FC) - 0.2, Max=max(avg_log2FC) + 0.2) %>%as.data.frame()
#灰色背景的寬度
background_position$unm %<>% as.vector(.) %>% as.numeric(.)
background_position$start <- background_position$unm - 0.4
background_position$end <- background_position$unm + 0.4

### 準(zhǔn)備繪制中間區(qū)域cluster彩色bar所需數(shù)據(jù)
cluster_bar_position <- background_position
cluster_bar_position$start <- cluster_bar_position$unm - 0.5
cluster_bar_position$end <- cluster_bar_position$unm + 0.5
cluster_bar_position$unm %<>% 
  factor(., levels = c(0:max(as.vector(.))))

## 設(shè)置填充顏色
cols_thr_signi <- c("Up_Highly" = "#d7301f",
                    "Down_Highly" = "#225ea8",
                    "Up_Lowly" = "black",
                    "Down_Lowly" = "black")
cols_cluster <- c("0" = "#35978f",
                  "1" = "#8dd3c7",
                  "2" = "#ffffb3",
                  "3" = "#bebada",
                  "4" = "#fb8072",
                  "5" = "#80b1d3",
                  "6" = "#fdb462",
                  "7" = "#b3de69")
p= ggplot() +
  geom_rect(data = background_position, aes(xmin = start, xmax = end, ymin = Min,ymax = Max),
            fill = "#525252", alpha = 0.1) + ###添加灰色背景色
  geom_jitter(data = r.deg, aes(x =unm, y = avg_log2FC, colour = thr_signi),
              size = 1,position = position_jitter(seed = 1)) +
  scale_color_manual(values = cols_thr_signi)+#改變點(diǎn)的顏色
  scale_x_continuous(limits = c(-0.5, max(r.deg$unm) + 0.5),
                     breaks = seq(0, max(r.deg$unm), 1),
                     label = seq(0, max(r.deg$unm),1)) + #修改坐標(biāo)軸顯示刻度

  # 根據(jù)top_label標(biāo)注基因名
  geom_text_repel(data = top_label, aes(x =unm, y = avg_log2FC, label = gene),
                  position = position_jitter(seed = 1), show.legend = F, size = 2.5,
                  box.padding = unit(0, "lines")) +
#繪制彩色bar
  geom_rect(data = cluster_bar_position, aes(xmin = start, xmax = end, ymin = -0.4,ymax = 0.4, fill = unm), color = "black", alpha = 1, show.legend = F) +
 scale_fill_manual(values = cols_cluster) +
  labs(x = "Cluster", y = "average log2FC") +
  theme_bw()

plot1 <- p + theme(panel.grid.minor = element_blank(), ##去除網(wǎng)格線
                   panel.grid.major = element_blank(),
                   axis.text.y = element_text(colour = 'black', size = 14),
                   axis.text.x = element_text(colour = 'black', size = 14, vjust = 67), #調(diào)整x軸坐標(biāo),vjust的值按照最終結(jié)果稍加調(diào)整
                   panel.border = element_blank(), ## 去掉坐標(biāo)軸
                   axis.ticks.x = element_blank(), ## 去掉的坐標(biāo)刻度線
                   axis.line.y = element_line(colour = "black")) #添加y軸坐標(biāo)軸

plot1

H形圖

6.FoldChange比較圖

加載package與數(shù)據(jù)

##系統(tǒng)報(bào)錯(cuò)改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())

###加載所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)
library(ggrepel)
library(reshape2)
setwd("C:/Users/86269/Desktop/shun.C/single_cell")
load("scRNA_harmony.Rdata")
library(ggplot2)
library(cowplot)

提取數(shù)據(jù)并尋找差異基因

#取cluster數(shù)為1的細(xì)胞
t.cells <- subset(scRNA_harmony, idents = "1")
#按sample分組
Idents(t.cells) <- "orig.ident"
t.cells=subset(t.cells,ident=c("sample_11","sample_3"))
#放縮
avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))#log1p(x)表示log2(x+1)
avg.t.cells$gene <- rownames(avg.t.cells)
#得出兩個(gè)sample在cluster1的差異表達(dá)基因
deg <- FindMarkers(t.cells,ident.1 = "sample_11",ident.2 = "sample_3", logfc.threshold = 2,group.by = "orig.ident")

抽取一部分差異表達(dá)基因在圖中展示

#隨機(jī)抽取15個(gè)基因
genes.to.label = sample(rownames(deg),15)
p1 <- ggplot(avg.t.cells, aes(sample_11,sample_3)) + geom_point() + ggtitle("T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p1 

p1

p1中基因label的指向不明確，下面要將被label的基因和其他基因區(qū)分開

degdf=avg.t.cells
degdf.l=degdf[genes.to.label,]
#提取label基因之外的點(diǎn)
degdf.p=degdf[-match(genes.to.label,rownames(degdf)),]
#增加一列提示基因是否被label
degdf.p$Sig="NO_DIFF"
degdf.l$Sig="DIFF"
degdf=rbind(degdf.p,degdf.l)

avg.t.cells

degdf

接下來進(jìn)行可視化

#用不同顏色將標(biāo)記的基因與未標(biāo)記基因區(qū)分開
p=ggplot(degdf, aes(sample_11,sample_3)) +
  geom_point( aes(color=Sig)) +
  
  scale_color_manual(values=c("red","grey"))+
  ggtitle("Stromal Cells")+
  theme(plot.title = element_text(size =18,hjust = 0.5, face = "bold")) +
  
  theme(axis.title.x =element_text(size=14), axis.title.y=element_text(size=14)) + 
  
  theme_bw()

p

#去除邊框贪壳，加深坐標(biāo)軸
p2 <- p+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))

#對(duì)選擇的基因進(jìn)行標(biāo)記
data_selected <- degdf[genes.to.label,]
p2 <- p2 + geom_label_repel(data=data_selected,
                     aes(label=rownames(data_selected)))

p2

還可以選擇展示某個(gè)特定的基因,以基因TIMP1為例

p +
  geom_point(size = 3, shape =2, data = data_selected[c("TIMP1"),],colour="green") +
  ggrepel::geom_label_repel(
    aes(label =genes.to.label),size =2.5,
    data = data_selected,
    color="black")+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))

展示TIMP1

7.火山圖

尋找差異基因

cd4.naive=subset(scRNA_harmony,idents ="1")
degdf <- FindMarkers(cd4.naive,ident.1 = "sample_11",ident.2 = "sample_3", logfc.threshold = 0.01,group.by = "orig.ident")
degdf$symbol <- rownames(degdf)

degdf

設(shè)定閾值并篩選差異表達(dá)明顯的基因

logFC_t=0
P.Value_t = 1e-28
#log2FC>0且padj<p則是上調(diào)基因,log2FC<0且padj<p則是下調(diào)基因饱亿，padj>p則是差異表達(dá)不明顯的基因
degdf$change = ifelse(degdf$p_val_adj < P.Value_t & degdf$avg_log2FC < 0,"down",
                      ifelse(degdf$p_val_adj < P.Value_t & degdf$avg_log2FC > 0,"up","stable"))

degdf

可視化

p=ggplot(degdf, aes(avg_log2FC,  -log10(p_val_adj))) +
  geom_point(alpha=0.4, size=2.8, aes(color=change)) +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("green", "grey","red"))+
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))

火山圖 越靠近兩側(cè)的基因差異表達(dá)越明顯

#展示感興趣的基因
data_selected <- degdf["TIMP1",]
p + geom_label_repel(data=data_selected,
                     aes(label=rownames(data_selected)))

展示TIMP1

#更改點(diǎn)的形狀
p +
  geom_point(size = 4, shape =2, data = data_selected,colour="black") +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = data_selected,
    color="green")