一简僧、準備工作

1. 軟件包安裝

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("pathview")
install.packages("msigdbr")
install.packages("ggnewscale")

library(Seurat)
library(tidyverse)
library(clusterProfiler)
library(pathview)
library(enrichplot)
library(msigdbr)

2. 尋找差異表達基因

差異分析參考

讀取數(shù)據(jù)

path_to_pbmc_rds = "sce_final.rds"
pbmc <- readRDS(path_to_pbmc_rds)

篩選差異基因

## 在尋找差異基因之前涎劈，把默認的assay切換為RNA。
DefaultAssay(pbmc) <- 'RNA'

## 定義好你想要在哪一個分群基礎(chǔ)上找差異表達基因
Idents(pbmc) <- 'seurat_clusters' 

## 在不同cluster/或者celltype中找差異表達基因
only.pos = TRUE           # 只找上調(diào)的差異表達基因
logfc.threshold = 0.25    # 差異基因的avg_log2FC必須要大于0.25

# 對所有Cluster依次進行差異表達分析谅海，將所有Cluster分為兩類蹦浦，我與非我
# all.clusters.markers <- FindAllMarkers(pbmc, only.pos = T, logfc.threshold = logfc.threshold)
# head(all.clusters.markers)

# 只對指定的Cluster進行差異表達分析
markers <- FindMarkers(pbmc, ident.1 = 1, ident.2 = 2, only.pos = T)
head(markers)

markers = markers %>% rownames_to_column('gene') %>% filter(p_val_adj < 0.05)
head(markers)

注：

• readRDS() 函數(shù)讀取*.rds對象。.rds與.rdata都是R語言特有的數(shù)據(jù)保存格式蝌诡，前者只能保存單一對象枫吧，此處用于獲取上次教程保存的Seurat對象九杂；后者可以保存某次R運行過程中所有的環(huán)境變量，可以下次重復(fù)使用甥捺，不用重復(fù)從頭執(zhí)行代碼镰禾。

• DefaultAssay() 函數(shù)常見的參數(shù)有"RNA"和"integrated"唱逢。在尋找差異基因之前惶我，把默認的assay切換為RNA，意味著后續(xù)的處理將使用原始Count值盯蝴。設(shè)置為integrated意味著使用整合后的數(shù)據(jù)進行后續(xù)操作听怕，此類數(shù)據(jù)可以認為是經(jīng)過批處理過的尿瞭。

• FindAllMarkers()與FindMarkers()函數(shù)：

  • FindAllMarkers()對所有Cluster依次進行差異表達分析，將所有Cluster分為兩類黑竞，我與非我很魂。比如對Cluster0進行分析時檐涝，其余Cluster1~7被認為是另一個整體。
  • FindMarkers()僅對指定的Cluster進行分析凡纳。

• 差異表達分析理論：

  • FoldChange(FC帝蒿，姑且翻譯為差異倍數(shù))葛超，相同基因在不同樣本組（可以是疾病組與正常組，也可以是這里的兩個Clusters）中平均表達量差異倍數(shù)，是A/B的倍數(shù)關(guān)系胖替。
  • 常見的差異倍數(shù)是將上述平均表達量的倍數(shù)關(guān)系進行取log2（FC）豫缨。
  • 根據(jù)對數(shù)圖像可以發(fā)現(xiàn)好芭，當x趨近于0時，y趨近于-∞招狸，因此看到的最終差異倍數(shù)往往是log2[(A+1)/(B+1)]裙戏。
  • 最終差異表達基因根據(jù)提前設(shè)定好的logFC和p值閾值即可得到厕诡。

差異表達分析結(jié)果解讀：

• 對Cluster1和2進行差異表達分析灵嫌，行名是差異表達基因，列名分別為p值猖凛，log2FC平均值形病，pct.1與pct.2表示相應(yīng)基因在兩個Cluster中的表達比例（即多少細胞表達了該基因），矯正p值（是一種更嚴格的p值量瓜，可替代p值作為篩選條件）绍傲。
• “markers = markers %>% rownames_to_column('gene') %>% filter(p_val_adj < 0.05)”作用依次為將行名作為列內(nèi)容添加到dataframe中耍共，相應(yīng)的數(shù)據(jù)維度增加一列试读，該新增列列名為'gene'钩骇，最后根據(jù)矯正p值小于0.05篩選差異基因。
• “%>%”符號類似于Linux中管道符的作用银亲，將數(shù)據(jù)從上游傳到下游纽匙。

Cluster1和2差異表達分析結(jié)果
差異表達基因閾值過濾

ID 轉(zhuǎn)換

# 加載物種包
# 根據(jù)物種來指定
OrgDb = "org.Hs.eg.db" 

# 基因Symbol轉(zhuǎn)化為基因ENTREZID
gene_convert <- bitr(markers$gene, fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = OrgDb)

# 合并兩個表格
markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL"))
head(markers )

二烛缔、富集分析

1. GO富集分析

注意：進行GO分析時需要使用轉(zhuǎn)換后的基因 ENTREZID

# GO
# posiible value: BP, CC, MF, all
ont = "BP" 

go.results <- enrichGO(markers$ENTREZID, keyType="ENTREZID",
                       ont="BP",OrgD = OrgDb, readable = TRUE)
go.results


dotplot(go.results)
barplot(go.results, showCategory = 10)
emapplot(pairwise_termsim(go.results))

GO富集分析可分為三個層次BP, CC, MF力穗，實際使用enrichGO() 函數(shù)進行富集分析時可通過指定BP, CC, MF,all四個參數(shù)当窗。此處僅使用BP舉例說明崖面。

markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL")) 作用：

• “gene_convert”數(shù)據(jù)框中包含基因Symbol和基因ENTREZID信息，markers數(shù)據(jù)框中包含差異表達基因信息庶香，兩者均含有基因但是列名不同赶掖。通過管道符及左聯(lián)函數(shù)將兩個數(shù)據(jù)框信息合并，見下圖陪白。

基因名轉(zhuǎn)換及合并

三個層次結(jié)果圖片一起畫

參考
org.Hs.eg.db

#GO富集分析
go.results <- enrichGO(gene = markers$ENTREZID, 
               OrgDb = OrgDb , 
               pvalueCutoff =0.05, 
               qvalueCutoff = 0.05,
               ont="all",
               readable =T)
write.table(go.results,file="GO.txt",sep="\t",quote=F,row.names = F)          #保存富集結(jié)果
 
#柱狀圖
# tiff(file="barplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="barplot.pdf",width = 10,height = 8)
barplot(go.results, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
 
#氣泡圖
# tiff(file="dotplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="dotplot.pdf",width = 10,height = 8)
dotplot(go.results,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()

#熱圖
tiff(file="heatplot.tiff",width = 40,height = 20,units ="cm",compression="lzw",bg="white",res=600)
heatplot(go.results,showCategory =20, foldChange=2,color = "pvalue")
dev.off()

柱狀圖：這張圖可以看出在CC,MF,BP中各個功能的terms富集的基因程度。橫坐標是富集在GO term中的基因數(shù)左邊的是GO的功能序厉，右邊是GO屬于什么數(shù)據(jù)庫以及可以看出顏色所代表的含義毕箍，越紅代表越顯著.

氣泡圖：可以到CC,BP,MF下的各個terms富集的基因數(shù)量和表達情況而柑。橫坐標代表基因所占的比例牺堰，右邊可以看出點的大小所代表的含義颅围，點越大睛榄，富集的基因越多，顏色越紅代表富集越顯著。

熱圖：可以看到每個基因在不同的GO的terms的表達程度茎辐，這可以選擇表達程度最好的幾個基因進行分析拖陆。

2. KEGG富集分析

KEGG富集分析過程與GO類似

organism = "hsa" # 對應(yīng)KEGG數(shù)據(jù)庫里的物種縮寫
kegg.results <- enrichKEGG(markers$ENTREZID, organism = organism)
write.table(kegg.results,file="KEGG.txt",sep="\t",quote=F,row.names = F)   # 保存富集結(jié)果
kegg.results <- setReadable(kegg.results, OrgDb = OrgDb, keyType='ENTREZID')

#柱狀圖
#tiff(file="barplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="barplot.pdf",width = 10,height = 7)
barplot(kegg.results, drop = TRUE, showCategory = 30)
dev.off()

#氣泡圖
#tiff(file="dotplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="bubble.pdf",width = 10,height = 7)
dotplot(kegg.results, showCategory = 30)
dev.off()

#熱圖
tiff(file="heatplot.tiff",width = 25,height = 15,units ="cm",compression="lzw",bg="white",res=600)
heatplot(kegg.results, showCategory =20, foldChange=cor)
dev.off()


#通路圖
library("pathview")
keggxls=read.table("KEGG.txt",sep="\t",header=T)
for(i in keggxls$ID){
  pv.out <- pathview(gene.data = cor, pathway.id = i, species = "hsa", out.suffix = "pathview")
}

kegg的通路里有的只是基因的id卻不是基因的名字依啰，所以需要轉(zhuǎn)化基因的id店枣，參考。

可視化其中一個通路：

diff_genes_avg_logFC = markers$avg_log2FC
head(diff_genes_avg_logFC)

names(diff_genes_avg_logFC) = markers$ENTREZID
head(diff_genes_avg_logFC)

pathview(gene.data = diff_genes_avg_logFC, species = organism, pathway.id = kegg.results@result$ID[3])

3. GSEA富集分析

GSEA富集分析與GO/KEGG通路富集分析的區(qū)別：

• 兩者的區(qū)別在于Pathway包含基因間互作信息长豁，Gene set只是一個單純的彼此獨立的集合蕉斜。就像新生第一次組成一個班宅此，彼此不認識爬范，相互沒聯(lián)系青瀑，就是Gene set，經(jīng)過一段時間相處彼此熟絡(luò)起來枝嘶，就是Pathway群扶。

• GSEA使用的背景基因稱為Gene set竞阐，GO/KEGG使用的背景基因稱為Pathway暑劝。
• 另外担猛，GSEA富集分析不用進行差異表達基因的尋找，但是需要將基因按照log2FC進行排序智嚷。

將基因按照log2FC進行排序

# markers for GSEA
## 至少多少比例的細胞表達這個基因盏道，過濾一些只在極少數(shù)細胞中有表達的基因
min.pct = 0.01 
## 過濾掉在兩組中幾乎沒有差異的基因
logfc.threshold = 0.01 
markers.for.gsea <- FindMarkers(pbmc, ident.1 = 1, ident.2 = 2, min.pct = min.pct, logfc.threshold=logfc.threshold)
# GSEA 要求輸入的是一個排好序的列表
Markers_genelist <- markers.for.gsea$avg_log2FC
names(Markers_genelist)= rownames(markers.for.gsea)
Markers_genelist <- sort(Markers_genelist, decreasing = T)

雖然說GSEA分析不用進行差異表達基因的尋找猜嘱，但是為了減少計算量還是提前過濾只在極少數(shù)細胞中表達的基因。

• 通過“min.pct = 0.01”和“l(fā)ogfc.threshold = 0.01”參數(shù)指定弦撩。
• 如果擔心過濾掉有意義的基因此處可以將兩個參數(shù)設(shè)為0益楼，此處進行FindMarkers（）操作主要是想得到基因的log2FC用于排序点晴。

進行GSEA富集分析

GSEA原理

MSigDB數(shù)據(jù)庫可選基因集

導(dǎo)入MSigDB

m_df = msigdbr(species = 'Homo sapiens' , category = "C2")
mf_df= m_df %>% dplyr::select(gs_name,gene_symbol)
colnames(mf_df)<-c("term","gene")

gsea.results <- GSEA(Markers_genelist, TERM2GENE = mf_df)
gsea.results

gseaplot(gsea.results, gsea.results@result$ID[3])

# 保存富集結(jié)果
write_csv(gsea.results %>% data.frame, "gsea_results.csv")

• msigdbr()作用：獲取接下來用于注釋的MSigDB背景數(shù)據(jù)集，可以通過“category”參數(shù)指定數(shù)據(jù)集族跛。

• GSEA（）函數(shù)作用：進行GSEA富集分析锐墙，參數(shù)“Markers_genelist”中包含按log2FC排序好的基因溪北，參數(shù)“TERM2GENE”中包含基因集與基因集中所有的基因信息。

• gseaplot（）函數(shù)作用：結(jié)果圖有兩個部分，上半部分橫坐標為待富集的基因敦锌，縱坐標是log2FC乙墙，可以看到基因是按log2FC降序排序好的生均。

• 下半部分马胧，橫坐標反應(yīng)基因在基因集中命中的情況，如果命中坐標軸上就出現(xiàn)一天豎線蛙粘，縱坐標表示富集得分（ES,Enrich Score）,命中就+某值出牧，沒命中就-某值。

• 結(jié)合上下兩部分來看评抚，我們知道此處觀察Cluster1與2之間的差異表達情況慨代，當log2FC大于0時边翼，差異基因列表中的基因更多出現(xiàn)在C2參考基因集中组底，此時FC比值大于1，即基因在Cluster1中更顯著表達江滨。[圖片上傳中...(image-e53601-1679814952079-0)]

• 可視化“gsea.results@result$ID[3]”該條基因集的富集結(jié)果唬滑。

• 結(jié)果圖解讀：

這個圖分為三個部分

• 第一部分：為基因 Enrichment Score 的折線圖晶密，橫軸為該基因下的每個基因模她，縱軸為對應(yīng)的 Running ES侈净，在折線圖中有個峰值，該峰值就是這個基因集的 Enrichemnt score畜侦，峰值之前的基因就是該基因集下的核心基因旋膳，核心基因就是在單細胞很重要很關(guān)鍵的基因，如果基因在頂部富集咏连，我們可以說祟滴，從總體上看垄懂，該基因集是上調(diào)趨勢，反之桶蛔，如果在底部富集仔雷，則是下調(diào)趨勢舔示；

• 第二部分：也就是hit值惕稻，也就是圖中的黑線出現(xiàn)的位置俺祠，這個黑線可以看出第一部分每次折現(xiàn)點的基因的位置

• 第三部分：這是rank值圖，采用了 Signal2Noise 算法淌铐，我對他的理解是腿准，縱軸它反映了基因集中的每個基因的相對表達量，橫軸可以看出這是這個基因在基因數(shù)據(jù)集中出現(xiàn)的位置加叁。

三它匕、GO/KEGG富集分析與GSEA區(qū)別

1. 概念

GO：是基因本體聯(lián)合會所建立的數(shù)據(jù)庫豫柬，旨在建立一個適用于各種物種的扑浸，對基因和蛋白質(zhì)功能進行限定和描述的喝噪，并能隨著研究不斷深入而更新的語義詞匯標準酝惧。GO 提供了一系列的語義用于描述基因功能的概念/類晚唇，以及這些概念之間的關(guān)系盗似。

KEGG：（京都基因與基因組百科全書）是了解高級功能和生物系統(tǒng)赫舒，從分子水平信息，尤其是大型分子數(shù)據(jù)集生成的基因組測序和其他高通量實驗技術(shù)的實用程序數(shù)據(jù)庫資源并鸵，是國際最常用的生物信息數(shù)據(jù)庫之一园担，以“理解生物系統(tǒng)的高級功能和實用程序資源庫”著稱弯汰。

GSEA：基因集富集分析湖雹，用于確定先驗基因集是否在兩種生物狀態(tài)（例如表型）之間差異顯著摔吏。

2. 區(qū)別

GO/KEGG差異基因的一刀切法——僅關(guān)注少數(shù)幾個顯著上調(diào)或下調(diào)的基因征讲，容易遺漏部分差異表達不顯著卻有重要生物學(xué)意義的基因诗箍，忽略一些基因的生物特性、基因調(diào)控網(wǎng)絡(luò)之間的關(guān)系及基因功能和意義等有價值的信息

GSEA不需要指定明確的差異基因閾值，算法根據(jù)實際整體趨勢分析筷狼。

參考：
https://mp.weixin.qq.com/s/ov-tLoxxK2laBm4NHjUS0Q
https://mp.weixin.qq.com/s/kP7qORDcLNFg4IDoy7m0PQ
https://mp.weixin.qq.com/s/oyj45DLIb-MNAw0LV1-VIA
https://mp.weixin.qq.com/s/pBvhxcTv0crH2VgwJRUN-Q
https://mp.weixin.qq.com/s/RvltqYd4PeECgW6yLoB2Lw
https://mp.weixin.qq.com/s/IowGCKgy-re2tzWDk35lJw
https://mp.weixin.qq.com/s/9jLQFuVATuXcvH-bVvwsRQ
https://mp.weixin.qq.com/s/JhFJJiQL-9Z6uDq4DjsgVA
https://mp.weixin.qq.com/s/QbSgYG_y1wsrMqdzGBRrQg