前言

Hello小伙伴們大家好，我是生信技能樹的托更小能手”我才不吃蛋黃“哈哈哈哈。最近在科室匯報(bào)了單細(xì)胞的學(xué)習(xí)成果,發(fā)現(xiàn)長時(shí)間沒有學(xué)習(xí)新知識，匯報(bào)的效果并未達(dá)到預(yù)期。要想給沒有生信基礎(chǔ)的聽眾把單細(xì)胞講明白，確實(shí)是需要一定的水平。只能怪自己學(xué)藝不精。

前段時(shí)間由于熬夜過猛，開始頻繁的出現(xiàn)胸痛，有一天早上上班時(shí)間持續(xù)胸痛了4小時(shí)，嚇得趕緊去查了心肌酶譜和冠脈CTA,萬幸檢查結(jié)果正常（回家被“蛋黃都給你吃”教育了一頓）。于是休息了（不熬夜）1周，終于緩過來了。外科大夫可太辛苦了，生信人也是經(jīng)常需要熬夜爆肝。為了保持競爭力，近半年一直持續(xù)扮演外科大夫+生信人的角色，buff可算是疊滿了。身體累垮一次之后，算是敲響了警鐘。怎么說呢，學(xué)會高效的工作、學(xué)習(xí)和休息也是一種必不可少的能力。2025年工作得繼續(xù)，學(xué)習(xí)也得繼續(xù)。

今天是肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)系列第五期。第四期我們繪制餅圖、堆積柱狀圖、箱線圖、氣泡圖等，比較不同分組之間細(xì)胞比例差異（見：肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)(四)：細(xì)胞比例可視化）。本期，我們將在單細(xì)胞水平進(jìn)行組間差異分析，并繪制多種圖形可視化差異基因。由于篇幅過長，推文分上下兩期，干貨滿滿，歡迎持續(xù)追更。

1.背景介紹

差異基因分析（Differential Gene Expression Analysis）是轉(zhuǎn)錄組研究中一種常見的分析方法，用于比較不同實(shí)驗(yàn)條件、時(shí)間點(diǎn)或處理組之間基因的表達(dá)差異。其目的是識別在不同條件下顯著變化的基因（稱為差異表達(dá)基因，Differentially Expressed Genes，DEGs），從而揭示與生物過程或表型相關(guān)的分子機(jī)制。

Bulk轉(zhuǎn)錄組差異分析常用的R包為DESeq2（適用于RNA-Seq數(shù)據(jù)的負(fù)二項(xiàng)分布建模）、EdgeR（適用于計(jì)數(shù)數(shù)據(jù)的廣義線性模型）和limma（用于微陣列和RNA-Seq數(shù)據(jù)的線性模型分析）。單細(xì)胞常用的是Seurat包內(nèi)置FindMarkers函數(shù)。

單細(xì)胞測序技術(shù)通過在單個(gè)細(xì)胞水平上對基因表達(dá)進(jìn)行定量分析，能夠揭示細(xì)胞群體內(nèi)部的異質(zhì)性和復(fù)雜性。進(jìn)行差異基因分析是單細(xì)胞測序研究中的關(guān)鍵步驟，其重要性主要體現(xiàn)在以下幾個(gè)方面：

揭示細(xì)胞異質(zhì)性：單細(xì)胞測序可以檢測到細(xì)胞群體中不同細(xì)胞的基因表達(dá)差異，這些差異可能與細(xì)胞的類型、狀態(tài)或功能有關(guān)。

識別關(guān)鍵生物標(biāo)志物：差異表達(dá)的基因可能與特定的生物學(xué)過程或疾病狀態(tài)相關(guān)，因此可以作為生物標(biāo)志物或治療靶點(diǎn)。

理解細(xì)胞發(fā)育和分化：通過比較不同發(fā)育階段或不同細(xì)胞分化路徑的細(xì)胞，差異基因分析有助于揭示細(xì)胞命運(yùn)決定的分子機(jī)制。

發(fā)現(xiàn)新的細(xì)胞類型：在復(fù)雜的組織或疾病微環(huán)境中，單細(xì)胞測序可以發(fā)現(xiàn)新的細(xì)胞亞群，這些亞群可能在以往的bulk測序中被忽視。

研究疾病機(jī)理：在疾病研究中，差異基因分析有助于理解疾病發(fā)生發(fā)展的分子機(jī)制，特別是在腫瘤微環(huán)境、自身免疫疾病等領(lǐng)域。

驗(yàn)證和擴(kuò)展bulk測序結(jié)果：單細(xì)胞測序可以與bulk測序結(jié)果相結(jié)合，提供更精細(xì)的基因表達(dá)模式，驗(yàn)證bulk測序發(fā)現(xiàn)的差異表達(dá)基因，并揭示更多細(xì)節(jié)。

時(shí)間序列分析：在動態(tài)過程中，如胚胎發(fā)育、組織再生或疾病進(jìn)程中，差異基因分析可以揭示隨時(shí)間變化的基因表達(dá)模式。

減少技術(shù)噪聲的影響：單細(xì)胞數(shù)據(jù)可能受到技術(shù)噪聲的影響，如dropout現(xiàn)象，差異基因分析有助于識別和排除這些噪聲，提高分析結(jié)果的可靠性。

提供更精確的生物學(xué)解釋：通過比較不同條件下的細(xì)胞，如治療前后或不同環(huán)境條件下，差異基因分析可以提供更精確的生物學(xué)解釋。

指導(dǎo)后續(xù)實(shí)驗(yàn)設(shè)計(jì)：差異基因分析的結(jié)果可以指導(dǎo)后續(xù)的實(shí)驗(yàn)設(shè)計(jì)，如功能驗(yàn)證實(shí)驗(yàn)、基因編輯實(shí)驗(yàn)等。

綜上所述，差異基因分析在單細(xì)胞測序中扮演著至關(guān)重要的角色，它不僅能夠提供細(xì)胞水平上的基因表達(dá)信息，還能夠?yàn)樯飳W(xué)研究提供新的視角和深入的見解。

2.差異分析

首先加載R包，加載細(xì)胞注釋后的數(shù)據(jù)：

rm(list=ls())
source('scRNA_scripts/lib.R')
dir.create("6-DEG")
setwd('6-DEG/')
sce.all=readRDS( "../3-Celltype/sce_celltype.rds")
sce.all

使用FindMarkers函數(shù)使用FindMarkers函數(shù)對肺癌三個(gè)組別IAC, MIA, AIS進(jìn)行差異分析:

library(tidyverse)
library(tinyarray)
scRNA = sce.all
head(scRNA@meta.data)
Idents(scRNA) = scRNA$sample
ct <- c('IAC', 'MIA', 'AIS')
if (!file.exists('markers_list.Rdata')) {
  markers_list <- list()
  for (i in 1:(length(ct)-1)) {
    for (j in (i+1):length(ct)) {
      markers <- FindMarkers(scRNA, group.by = "sample",
                             logfc.threshold = 0.1,
                             ident.1 = ct[i],
                             ident.2 = ct[j])
      markers_list[[paste0(ct[i], "_vs_", ct[j])]] <- markers
    }
  }
  save(markers_list,file = 'markers_list.Rdata')
} else {
  load('markers_list.Rdata')
}

查看差異分析的結(jié)果：在這里，我們設(shè)置的差異表達(dá)基因的閾值是p_val_adj < 0.01，Up(avg_log2FC>1)，Down(avg_log2FC< -1)。

length(markers_list)
lapply(markers_list,nrow)

all_markers_sig = lapply(markers_list, function(x){
  markers_sig <- subset(x, p_val_adj < 0.01)
})

marker_stat = as.data.frame(lapply(all_markers_sig,function(x){
  # x=all_markers_sig[[1]]
  Up = sum(x$avg_log2FC>1)
  Down = sum(x$avg_log2FC< -1)
  Total = Up+Down
  return(c(Up, Down, Total))
}))

查看差異基因數(shù)量：

rownames(marker_stat) = c("Up","Down","Total")
marker_stat

三組差異分析的結(jié)果分別標(biāo)記上UP和DOWN：

library(gridExtra)
IAC_vs_MIA = all_markers_sig[[1]]
IAC_vs_AIS = all_markers_sig[[2]]
MIA_vs_AIS = all_markers_sig[[3]]
##每一組加上group上下調(diào)
deg <- function(dat){
  dat$group = NA
  dat$group[(dat$avg_log2FC < -1) & (dat$p_val_adj < 0.01)] <- "down"
  dat$group[(dat$avg_log2FC > 1) & (dat$p_val_adj < 0.01)] <- "up"
  table(dat$group )
  dat = na.omit(dat)
  dat
}

查看上下調(diào)基因的數(shù)量：

P001 = deg(IAC_vs_MIA)
table(P001$group)
P002 = deg(IAC_vs_AIS)
table(P002$group)
P003 = deg(MIA_vs_AIS)
table(P003$group)

3.差異分析可視化

veen圖：

####veen圖####
###上調(diào)
IAC_vs_MIA_up = rownames(P001[P001$group == 'up',]) 
IAC_vs_AIS_up = rownames(P002[P002$group == 'up',]) 
MIA_vs_AIS_up = rownames(P003[P003$group == 'up',]) 
inter_upgene <- intersect(intersect(IAC_vs_MIA_up, IAC_vs_AIS_up), MIA_vs_AIS_up)
#三元#
#BiocManager::install("VennDetail")
library(VennDetail)
library(VennDiagram) 
#IAC_vs_MIA_up,IAC_vs_AIS_up,MIA_vs_AIS_up
venn <- venndetail(list(IAC_vs_MIA_up = IAC_vs_MIA_up, IAC_vs_AIS_up = IAC_vs_AIS_up, MIA_vs_AIS_up= MIA_vs_AIS_up))
detail(venn) 

venn.diagram函數(shù)及代碼注釋：

venn.diagram(x=list(IAC_vs_MIA_up,IAC_vs_AIS_up,MIA_vs_AIS_up),
             scaled = F, # 根據(jù)比例顯示大小
             alpha= 0.5, #透明度
             lwd=1,lty=1,col=c('#FFFFCC','#CCFFFF',"#FFCCCC"), #圓圈線條粗細(xì)、形狀、顏色；1 實(shí)線, 2 虛線, blank無線條
             label.col ='black' , # 數(shù)字顏色abel.col=c('#FFFFCC','#CCFFFF',......)根據(jù)不同顏色顯示數(shù)值顏色
             cex = 2, # 數(shù)字大小
             fontface = "bold",  # 字體粗細(xì)；加粗bold
             fill=c('#FFFFCC','#CCFFFF',"#FFCCCC"), # 填充色 配色https://www.58pic.com/
             category.names = c("IAC_vs_MIA_up", "IAC_vs_AIS_up","MIA_vs_AIS_up") , #標(biāo)簽名
             cat.dist = 0.07, # 標(biāo)簽距離圓圈的遠(yuǎn)近
             cat.pos = c(-30, -330, -180), # 標(biāo)簽相對于圓圈的角度cat.pos = c(-10, 10, 135)
             cat.cex = 1, #標(biāo)簽字體大小
             cat.fontface = "bold",  # 標(biāo)簽字體加粗
             cat.col='black' ,   #cat.col=c('#FFFFCC','#CCFFFF',.....)根據(jù)相應(yīng)顏色改變標(biāo)簽顏色
             cat.default.pos = "outer",  # 標(biāo)簽位置, outer內(nèi);text 外
             output=TRUE,
             filename='veenup.png',# 文件保存
             imagetype="png",  # 類型（tiff png svg）
             resolution = 400,  # 分辨率
             compression = "lzw",# 壓縮算法
             height = 2100 ,   # 高度
             width = 2300
             
)

韋恩餅圖：

plot(venn, type = "vennpie")

韋恩條形圖:

dplot(venn, order = TRUE, textsize = 4)

upset圖:

plot(venn, type = "upset")

下調(diào)基因的韋恩圖和上調(diào)基因一致：

IAC_vs_MIA_down = rownames(P001[P001$group == 'down',]) 
IAC_vs_AIS_down = rownames(P002[P002$group == 'down',]) 
MIA_vs_AIS_down = rownames(P003[P003$group == 'down',]) 

inter_downgene <- intersect(intersect(IAC_vs_MIA_down, IAC_vs_AIS_down), MIA_vs_AIS_down)

venn <- venndetail(list(IAC_vs_MIA_down = IAC_vs_MIA_down, IAC_vs_AIS_down = IAC_vs_AIS_down, MIA_vs_AIS_down= MIA_vs_AIS_down))
detail(venn) 

# 韋恩圖
venn.diagram(x=list(IAC_vs_MIA_down,IAC_vs_AIS_down,MIA_vs_AIS_down),
             scaled = F, # 根據(jù)比例顯示大小
             alpha= 0.5, #透明度
             lwd=1,lty=1,col=c('#FFFFCC','#CCFFFF',"#FFCCCC"), #圓圈線條粗細(xì)、形狀、顏色；1 實(shí)線, 2 虛線, blank無線條
             label.col ='black' , # 數(shù)字顏色abel.col=c('#FFFFCC','#CCFFFF',......)根據(jù)不同顏色顯示數(shù)值顏色
             cex = 2, # 數(shù)字大小
             fontface = "bold",  # 字體粗細(xì)；加粗bold
             fill=c('#FFFFCC','#CCFFFF',"#FFCCCC"), # 填充色 配色https://www.58pic.com/
             category.names = c("IAC_vs_MIA_down", "IAC_vs_AIS_down","MIA_vs_AIS_down") , #標(biāo)簽名
             cat.dist = 0.07, # 標(biāo)簽距離圓圈的遠(yuǎn)近
             cat.pos = c(-30, -330, -180), # 標(biāo)簽相對于圓圈的角度cat.pos = c(-10, 10, 135)
             cat.cex = 1, #標(biāo)簽字體大小
             cat.fontface = "bold",  # 標(biāo)簽字體加粗
             cat.col='black' ,   #cat.col=c('#FFFFCC','#CCFFFF',.....)根據(jù)相應(yīng)顏色改變標(biāo)簽顏色
             cat.default.pos = "outer",  # 標(biāo)簽位置, outer內(nèi);text 外
             output=TRUE,
             filename='veendown.png',# 文件保存
             imagetype="png",  # 類型（tiff png svg）
             resolution = 400,  # 分辨率
             compression = "lzw",# 壓縮算法
             height = 2100 ,   # 高度
             width = 2300
             
)
###結(jié)果和文章中不一致，沒關(guān)系繼續(xù)往下

# 韋恩餅圖
plot(venn, type = "vennpie")

# 韋恩條形圖
dplot(venn, order = TRUE, textsize = 4)

# downset圖
plot(venn, type = "upset")

3D PCA:

繪制PCA圖，可以展示樣本間的差異。

數(shù)據(jù)準(zhǔn)備：

##3D PCA
set.seed(98765)
##PCA showing the DEGs among the three groups. The distance between dots represents the difference between groups
table(scRNA$patient)
meta = scRNA@meta.data
meta$sample_patient = paste(meta$sample,meta$patient,sep = '_')
table(meta$sample_patient)
#Idents(scRNA)
sce = scRNA
colnames(sce)
rownames(sce)
sce@meta.data = meta
avg <- AverageExpression(object = sce, group.by = "sample_patient")
a = avg$RNA
b = as.data.frame(a)

提取出差異基因：

DEGgene = c(rownames(IAC_vs_AIS),rownames(IAC_vs_MIA),rownames(MIA_vs_AIS))
##有一些重復(fù)的差異基因你
DEG = b[rownames(b) %in% DEGgene,]
str(DEG)
dat = as.data.frame(t(DEG))
pca.res <- prcomp(dat, scale. = T, center = T)
pca.res
tmp <- as.data.frame(pca.res$x)
head(tmp)

加載R包，繪制PCA圖：

library(scatterplot3d)
rownames(tmp)
col = ifelse(str_detect(rownames(tmp),"AIS"),'purple',
             ifelse(str_detect(rownames(tmp),"IAC"),'orange','green'))
scatterplot3d(tmp[,1:3],color=col,
              pch = 16,angle=30,
              box=T,type="p",
              lty.hide=2,lty.grid = 2)
legend("topleft",c('AIS','IAC','MIA'),
       fill=c('purple','orange','green'),box.col=NA,cex=0.7)

3D PCA：

二由老師：文章中PCA顯示MIA和IAC在轉(zhuǎn)錄組水平上非常接近；我這里顯示的是IAC和AIS更為接近，倒是和前面的韋恩圖對應(yīng)上了；但是這更能說明IAC是最終狀態(tài)啊，期間的MIA是中間狀態(tài)，感覺也是可以解釋的。

結(jié)語

本期，我們在單細(xì)胞水平進(jìn)行組間差異分析，并繪制了差異基因韋恩圖。下一期，我們講繪制兩組差異分析共顯示dotplot圖及三組差異分析火山圖。干貨滿滿，歡迎大家持續(xù)追更，謝謝！