singleseqgset | 單細(xì)胞RNA-Seq基因集富集分析

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測序分析?(重磅綜述：三萬字長文讀懂單細(xì)胞RNA測序分析的最佳實(shí)踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

singleseqgset是用于單細(xì)胞RNA-seq數(shù)據(jù)的基因集富集分析的軟件包。它使用簡單的基礎(chǔ)統(tǒng)計(jì)量（variance inflated Wilcoxon秩和檢驗(yàn)）來確定不同cluster中感興趣的基因集的富集。

Installation

library(devtools) install_github("arc85/singleseqgset") library(singleseqgset)

Introduction

基因集富集分析是一種無處不在的工具，用于從基因組數(shù)據(jù)集中提取具有生物學(xué)意義的信息。有許多不同類型的工具可用于基因集富集分析，以前的基因組富集分析工具的關(guān)鍵概念是比較組內(nèi)基因與組外基因的差異（例如，各組之間基因表達(dá)的log fold change變化）。但最值得肯定的是Subramanian等人（PNAS 2005）(https://www.ncbi.nlm.nih.gov/pubmed/16199517)的開創(chuàng)性工作-Gene Set Enrichment Analysis（GSEA富集分析 - 界面操作）。

GSEA分析中已知功能的基因集可以從許多地方進(jìn)行獲取，包括Broad的Molecular Signatures Database（MSigDB）和The Gene Ontology Resource，關(guān)鍵是選擇最能回答當(dāng)前生物學(xué)問題的基因集。舉個(gè)栗子，通常您不會(huì)使用Broad的C7（免疫學(xué)基因集）中的基因集來回答有關(guān)神經(jīng)元發(fā)育的問題（除非有充分的理由！）。

singleseqgset研究者在這里開發(fā)的基因集富集方法靈感來自Di Wu和Gordon K Smyth在Nucleic Acids Research 2012(https://www.ncbi.nlm.nih.gov/pubmed/16199517)上的工作，而Wu和Smyth在這個(gè)工作的基礎(chǔ)上又開發(fā)了相關(guān)調(diào)整的MEan RAnk基因集測試（CAMERA）。

CAMERA是一項(xiàng)競爭性基因集富集測試，可控制基因集內(nèi)的基因間相關(guān)性。研究團(tuán)隊(duì)之所以選擇使用CAMERA，是因?yàn)樗灰蕾囉诨颡?dú)立性的假設(shè)（因?yàn)槲覀冎阑蛲ǔ＞哂薪Y(jié)構(gòu)化的共表達(dá)模式），也不依賴于基因標(biāo)記的排列或測試樣品的重采樣。CAMERA通過基于基因間相關(guān)程度生成方差膨脹因子來控制基因間相關(guān)性，并將這種方差膨脹納入Wilcoxon秩和檢驗(yàn)中。

此workflow演示了對模擬數(shù)據(jù)使用singleseqgset的標(biāo)準(zhǔn)用法。我們將執(zhí)行以下步驟：

使用splatter R包模擬表達(dá)式數(shù)據(jù);

使用msigdbr下載感興趣的基因集;

將特定基因集添加到我們的模擬數(shù)據(jù)中（就是讓我們的模擬數(shù)據(jù)可以富含某基因集）;

使用標(biāo)準(zhǔn)的Seurat工作流程（v.2.3.4）處理我們的數(shù)據(jù);

使用singleseqgset進(jìn)行基因集富集分析;

將結(jié)果繪制在熱圖中;

Simulate data using splatter

首先，加載所有必要的程序包，并使用splatter模擬數(shù)據(jù)。

suppressMessages({ library(splatter) library(Seurat) library(msigdbr) library(singleseqgset) library(heatmap3) }) #Splatter是用于模擬單細(xì)胞RNA測序count數(shù)據(jù)的軟件包。 #創(chuàng)建參數(shù)并模擬數(shù)據(jù) sc.params <- newSplatParams(nGenes=1000,batchCells=5000) sim.groups <- splatSimulate(params=sc.params,method="groups",group.prob=c(0.4,0.2,0.3,0.1),de.prob=c(0.20,0.20,0.1,0.3),verbose=F) sim.groups #Check out the SingleCellExperiment object with simulated dataset## class: SingleCellExperiment ## dim: 1000 5000 ## metadata(1): Params ## assays(6): BatchCellMeans BaseCellMeans ... TrueCounts counts ## rownames(1000): Gene1 Gene2 ... Gene999 Gene1000 ## rowData names(8): Gene BaseGeneMean ... DEFacGroup3 DEFacGroup4 ## colnames(5000): Cell1 Cell2 ... Cell4999 Cell5000 ## colData names(4): Cell Batch Group ExpLibSize ## reducedDimNames(0): ## spikeNames(0):colData(sim.groups) #The colData column "Group" contains the groups of simulated cells## DataFrame with 5000 rows and 4 columns ## Cell Batch Group ExpLibSize ## <factor> <character> <factor> <numeric> ## Cell1 Cell1 Batch1 Group1 73324.2901799195 ## Cell2 Cell2 Batch1 Group1 71115.1438815874 ## Cell3 Cell3 Batch1 Group3 89689.371004738 ## Cell4 Cell4 Batch1 Group3 44896.460028516 ## Cell5 Cell5 Batch1 Group3 53956.0668720417 ## ... ... ... ... ... ## Cell4996 Cell4996 Batch1 Group2 79041.7558615305 ## Cell4997 Cell4997 Batch1 Group4 66684.3797711752 ## Cell4998 Cell4998 Batch1 Group3 70407.9613848431 ## Cell4999 Cell4999 Batch1 Group3 64645.3257427214 ## Cell5000 Cell5000 Batch1 Group3 62266.9119469426#我們將從sim.groups對象中提取模擬的計(jì)數(shù)和組別 sim.counts <- assays(sim.groups)$counts groups <- colData(sim.groups)$Group names(groups) <- rownames(colData(sim.groups)) table(groups)## groups ## Group1 Group2 Group3 Group4 ## 1977 992 1536 495

有次我們創(chuàng)建了一些模擬數(shù)據(jù)，這些數(shù)據(jù)由4個(gè)cluster組成，cluster中不同基因的表達(dá)比例不同。

Download gene sets of interest using msigdbr

在這里，我們將使用misigdbr軟件包下載Hallmark Gene Sets。另外，我們可以直接從Broad網(wǎng)站（http://software.broadinstitute.org/gsea/msigdb/index.jsp）下載基因集，然后使用`GSEABase`軟件包將其讀入R。

注意：必須為您的項(xiàng)目選擇適當(dāng)?shù)幕蜃⑨寴邮?#xff08;gene symbols或entrez gene ID；在這里我們只是gene symbols）和物種（在這里我們使用human）。否則，您的基因組中的基因名稱將與您的數(shù)據(jù)中的基因名稱不匹配，并且所有基因組都將獲得“0”富集得分。

h.human <- msigdbr(species="Homo sapiens",category="H")h.names <- unique(h.human$gs_name)h.sets <- vector("list",length=length(h.names)) names(h.sets) <- h.namesfor (i in names(h.sets)) {h.sets[[i]] <- pull(h.human[h.human$gs_name==i,"gene_symbol"]) }

Add gene sets to simulated clusters

我們將分配5個(gè)基因集在每個(gè)cluster中專門表達(dá)（如果不進(jìn)行分配，則可能出現(xiàn)無功能集在cluster中富集的情況。）。隨機(jī)選擇這20個(gè)集合，并從zero-inflated的Poisson分布中抽取每個(gè)基因，成功概率等于50％，λ值為10，這將模擬中等dropout和相對較低的基因表達(dá)量。

sets.to.use <- sample(seq(1,50,1),20,replace=F) sets.and.groups <- data.frame(set=sets.to.use,group=paste("Group",rep(1:4,each=5),sep=""))for (i in 1:nrow(sets.and.groups)) {set.to.use <- sets.and.groups[i,"set"]group.to.use <- sets.and.groups[i,"group"]gene.set.length <- length(h.sets[[set.to.use]])blank.matrix <- matrix(0,ncol=ncol(sim.counts),nrow=gene.set.length)rownames(blank.matrix) <- h.sets[[sets.to.use[i]]]matching <- rownames(blank.matrix) %in% rownames(sim.counts)if (any(matching==TRUE)) {matching.genes <- rownames(blank.matrix)[matching]nonmatching.genes <- setdiff(rownames(blank.matrix),matching.genes)blank.matrix <- blank.matrix[!matching,]sim.counts <- rbind(sim.counts,blank.matrix)} else {sim.counts <- rbind(sim.counts,blank.matrix)matching.genes <- rownames(blank.matrix)}group.cells <- colnames(sim.counts)[groups==group.to.use]for (z in group.cells) {if (any(matching==TRUE)) {sim.counts[matching.genes,z] <- ifelse(rbinom(length(matching.genes),size=1,prob=0.5)>0,0,rpois(length(matching.genes),lambda=10))sim.counts[nonmatching.genes,z] <- ifelse(rbinom(length(nonmatching.genes),size=1,prob=0.5)>0,0,rpois(length(nonmatching.genes),lambda=10))} else {sim.counts[matching.genes,z] <- ifelse(rbinom(length(matching.genes),size=1,prob=0.5)>0,0,rpois(length(matching.genes),lambda=10))}}}#Here, we will check out the sum of expression for the first gene set len.set1 <- length(h.sets[[sets.to.use[[1]]]]) plot(apply(sim.counts[1001:(1000+len.set1),],2,sum),xlim=c(0,5000),xlab="Cells",ylab="Sum of Gene Set 1 Expression")

我們已經(jīng)成功地修改了模擬計(jì)數(shù)矩陣，使其包含特定cluster中感興趣的基因集。例如，我們知道Group1細(xì)胞應(yīng)顯示出以下基因的富集：

HALLMARK_APICAL_SURFACE, HALLMARK_COMPLEMENT, HALLMARK_DNA_REPAIR, HALLMARK_IL2_STAT5_SIGNALING, HALLMARK_UNFOLDED_PROTEIN_RESPONSE。

Seurat workflow on simulated data

#構(gòu)建seurat object ser <- CreateSeuratObject(raw.data=sim.counts)#歸一化 ser <- NormalizeData(ser) ser@var.genes <- rownames(ser@raw.data)[1:1000] ser <- ScaleData(ser,genes.use=ser@var.genes)## Scaling data matrix#我們會(huì)假設(shè)1000個(gè)模擬基因是“可變基因”， #我們將跳過Seurat的FindVariableGenes ser <- RunPCA(ser,pc.genes=ser@var.genes,do.print=FALSE)PCElbowPlot(ser)

#其中top4PCs代表了多數(shù)的差異 #我們將會(huì)選擇top5 PCs;ser <- RunTSNE(ser,dims.use=1:5)#將模擬的dataID號加入Seurat object data.to.add <- colData(sim.groups)$Group names(data.to.add) <- ser@cell.names ser <- AddMetaData(ser,metadata=data.to.add,col.name="real_id") ser <- SetAllIdent(ser,id="real_id")#我們將會(huì)跳過使用Seurat聚類的步驟，因?yàn)槲覀円阎勒嬲木垲怚DsTSNEPlot(ser,group.by="real_id")

Use singleseqgset to perform gene set enrichment analysis

首先，我們將基于logFC計(jì)算基因集富集測試的矩陣。我們選擇使用已針對library大小進(jìn)行校正的log-normalized數(shù)據(jù)，并存儲(chǔ)在Seurat的“@data”slot中。

logfc.data <- logFC(cluster.ids=ser@meta.data$real_id,expr.mat=ser@data) names(logfc.data)## [1] "cluster.cells" "log.fc.cluster"

logFC函數(shù)返回一個(gè)長度為2的列表，其中包含每個(gè)cluster中的細(xì)胞以及cluster之間基因的對數(shù)倍變化。下一步需要此數(shù)據(jù)計(jì)算enrichment scores和p值。

gse.res <- wmw_gsea(expr.mat=ser@data,cluster.cells=logfc.data[[1]],log.fc.cluster=logfc.data[[2]],gene.sets=h.sets)## [1] "Removed 1 rows with all z-scores equal to zero."names(gse.res)## [1] "GSEA_statistics" "GSEA_p_values"

此函數(shù)返回兩個(gè)列表，第一個(gè)包含測試統(tǒng)計(jì)信息“GSEA_statistics”，第二個(gè)包含p值“GSEA_p_values”。我們可以將這些結(jié)果繪制為熱圖，以可視化差異調(diào)節(jié)基因集。

res.stats <- gse.res[["GSEA_statistics"]] res.pvals <- gse.res[["GSEA_p_values"]]res.pvals <- apply(res.pvals,2,p.adjust,method="fdr") #Correct for multiple comparisonsres.stats[order(res.stats[,1],decreasing=TRUE)[1:10],] #Top gene sets enriched by z scores## Group1 Group2 Group3 ## HALLMARK_IL2_STAT5_SIGNALING 10.51698615 -7.55892580 -8.113191 ## HALLMARK_DNA_REPAIR 10.47159479 -7.62264059 -7.326454 ## HALLMARK_COMPLEMENT 10.17982979 -3.68698969 -8.347697 ## HALLMARK_UNFOLDED_PROTEIN_RESPONSE 9.54278550 -8.98321228 -7.275756 ## HALLMARK_APICAL_SURFACE 7.63384635 -5.70116613 0.000000 ## HALLMARK_ANGIOGENESIS 1.32886963 -0.91447563 0.000000 ## HALLMARK_HEDGEHOG_SIGNALING 0.77628519 0.66226233 0.000000 ## HALLMARK_KRAS_SIGNALING_UP 0.59844933 -0.20261236 -2.365948 ## HALLMARK_COAGULATION 0.57787404 9.03118414 -6.988972 ## HALLMARK_INFLAMMATORY_RESPONSE 0.05206778 -0.08661822 -1.260939 ## Group4 ## HALLMARK_IL2_STAT5_SIGNALING -8.4400626 ## HALLMARK_DNA_REPAIR -10.3905166 ## HALLMARK_COMPLEMENT -12.5780500 ## HALLMARK_UNFOLDED_PROTEIN_RESPONSE -5.8693210 ## HALLMARK_APICAL_SURFACE 0.0000000 ## HALLMARK_ANGIOGENESIS -0.6935659 ## HALLMARK_HEDGEHOG_SIGNALING 0.3350711 ## HALLMARK_KRAS_SIGNALING_UP -1.1822606 ## HALLMARK_COAGULATION -4.3269207 ## HALLMARK_INFLAMMATORY_RESPONSE -3.3347607res.pvals[order(res.stats[,1],decreasing=TRUE)[1:10],] #Top gene sets by p values## Group1 Group2 Group3 ## HALLMARK_IL2_STAT5_SIGNALING 2.858190e-24 2.489271e-13 3.451510e-15 ## HALLMARK_DNA_REPAIR 2.858190e-24 1.739767e-13 1.286641e-12 ## HALLMARK_COMPLEMENT 3.985134e-23 6.949503e-04 6.821492e-16 ## HALLMARK_UNFOLDED_PROTEIN_RESPONSE 1.362655e-20 2.577150e-18 1.687982e-12 ## HALLMARK_APICAL_SURFACE 1.594960e-13 5.830539e-08 1.000000e+00 ## HALLMARK_ANGIOGENESIS 3.003553e-01 5.697704e-01 1.000000e+00 ## HALLMARK_HEDGEHOG_SIGNALING 5.642487e-01 7.318339e-01 1.000000e+00 ## HALLMARK_KRAS_SIGNALING_UP 6.589077e-01 1.000000e+00 4.406075e-02 ## HALLMARK_COAGULATION 6.589077e-01 2.080345e-18 1.130707e-11 ## HALLMARK_INFLAMMATORY_RESPONSE 1.000000e+00 1.000000e+00 3.631134e-01 ## Group4 ## HALLMARK_IL2_STAT5_SIGNALING 1.942653e-16 ## HALLMARK_DNA_REPAIR 6.709712e-24 ## HALLMARK_COMPLEMENT 1.366256e-34 ## HALLMARK_UNFOLDED_PROTEIN_RESPONSE 2.144160e-08 ## HALLMARK_APICAL_SURFACE 1.000000e+00 ## HALLMARK_ANGIOGENESIS 6.462100e-01 ## HALLMARK_HEDGEHOG_SIGNALING 7.856739e-01 ## HALLMARK_KRAS_SIGNALING_UP 3.417063e-01 ## HALLMARK_COAGULATION 4.939474e-05 ## HALLMARK_INFLAMMATORY_RESPONSE 2.460747e-03names(h.sets)[sets.to.use[1:5]] #Compare to the simulate sets we created## [1] "HALLMARK_APICAL_SURFACE" ## [2] "HALLMARK_COMPLEMENT" ## [3] "HALLMARK_DNA_REPAIR" ## [4] "HALLMARK_IL2_STAT5_SIGNALING" ## [5] "HALLMARK_UNFOLDED_PROTEIN_RESPONSE"#Plot the z scores with heatmap3 heatmap3(res.stats,Colv=NA,cexRow=0.5,cexCol=1,scale="row")

恭喜，您已經(jīng)使用singleseqgset執(zhí)行了第一個(gè)基因集富集分析！

你可能還想看

• R語言 - 富集分析泡泡圖

• 拿到基因兩眼一抹黑？沒關(guān)系，先做個(gè)基因富集分析吧！

• GO、GSEA富集分析一網(wǎng)打進(jìn)

• 無需寫代碼的高顏值富集分析神器

• 去東方，最好用的在線GO富集分析工具

往期精品(點(diǎn)擊圖片直達(dá)文字對應(yīng)教程)

后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的一个R包完成单细胞基因集富集分析 （全代码）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。