前言
前面蹂季，我們介紹過了差異基因的功能富集分析，今天疏日，我們對這部分的內(nèi)容作一些補充
主要介紹一下 GEVA偿洁、ssGSEA 和單基因的富集分析
GSVA
我們知道，GSEA 富集分析方法是針對兩組樣本來進(jìn)行評估的沟优，也就是說對基因列表的排列方式是根據(jù)基因與表型的相關(guān)度（例如涕滋，FC 值）來計算的，無法對單個樣本使用
其富集分?jǐn)?shù)（Enrichment Score挠阁，ES）的計算方式為
依次判斷基因列表中的基因是否在基因集合中宾肺，如果在基因集合中，則 ES 加上該基因與表型的相關(guān)度侵俗，如果不是集合中的基因锨用，則減去對應(yīng)的值，最后可以計算出一個最大 ES隘谣。
Gene Set Variation Analysis（GSVA）與 GSEA 的原理類似增拥，只是計算每個基因集合在每個樣本中的 enrichment statistic（ES，或 GSVA score）寻歧，其算法流程如下
不同于 GSEA 之處在于掌栅，對于不同的數(shù)據(jù)類型（只支持 log 表達(dá)值或原始的 read counts 值），假設(shè)了不同的累積密度函數(shù)（cumulative density function码泛，CDF）
- 芯片數(shù)據(jù)：正態(tài)分布密度函數(shù)
-
RNA-seq數(shù)據(jù)：泊松分布密度函數(shù)
而且猾封，GSVA 是為每個樣本的每個基因計算對應(yīng)的 CDF 值，然后根據(jù)該值對基因進(jìn)行排序噪珊，這樣晌缘，每個樣本都有一個從大到小排序的基因列表
對于某一基因集合，計算其在每個樣本中的 ES 值痢站，也就是評估基因集合在基因列表中的富集情況磷箕。
例如，我們有一個排序后的樣本
基因集合包含：B瑟押、E搀捷、H，我們可以繪制這樣一張 K-S 分布圖
x 軸為排序后的基因順序，依照這一順序嫩舟，如果基因在集合內(nèi)氢烘，則累積和會加上該基因的值（與基因的順序有關(guān)，排名越靠前值越大）家厌，否則累積和不變播玖。
將基因列表分為基因集內(nèi)核基因集外兩個集合，就可以繪制兩個分布（紅色和綠色曲線）饭于，分別計算兩個分布之間的最大間距蜀踏，以基因集內(nèi)的分布值更大視為正間距（即紅色曲線更高），兩個最大間距之和即為該基因集的 ES 值
這樣掰吕，就把基因水平的表達(dá)矩陣轉(zhuǎn)換成了基因集水平的評分矩陣果覆，可以使用差異表達(dá)基因識別算法，尋找顯著差異的基因集殖熟，從而達(dá)到類似功能富集的作用
1. GSVA 分析
先獲取基因表達(dá)矩陣局待，我們使用 TCGA 肺腺癌和肺鱗癌各 10 個樣本的 read counts 數(shù)據(jù)
library(TCGAbiolinks)
# 獲取表達(dá)矩陣
get_count <- function(cancer, n = 10) {
query <- GDCquery(
project = cancer,
data.category = "Gene expression",
data.type = "Gene expression quantification",
platform = "Illumina HiSeq",
file.type = "results",
sample.type = c("Primary Tumor"),
legacy = TRUE
)
# 選擇 n 個樣本
query$results[[1]] <- query$results[[1]][1:n,]
GDCdownload(query)
# 獲取 read count
exp.count <- GDCprepare(
query,
summarizedExperiment = TRUE,
)
return(exp.count)
}
luad.count <- get_count("TCGA-LUAD")
lusc.count <- get_count("TCGA-LUSC")
dataPrep_luad <- TCGAanalyze_Preprocessing(
object = luad.count,
cor.cut = 0.6,
datatype = "raw_count"
)
dataPrep_lusc <- TCGAanalyze_Preprocessing(
object = lusc.count,
cor.cut = 0.6,
datatype = "raw_count"
)
# 合并數(shù)據(jù)并使用 gcContent 方法進(jìn)行標(biāo)準(zhǔn)化
dataNorm <- TCGAanalyze_Normalization(
tabDF = cbind(dataPrep_luad, dataPrep_lusc),
geneInfo = TCGAbiolinks::geneInfo,
method = "gcContent"
)
# 分位數(shù)過濾
dataFilt <- TCGAanalyze_Filtering(
tabDF = dataNorm,
method = "quantile",
qnt.cut = 0.25
)
# 將數(shù)據(jù)拆分
luad.exp <- subset(dataFilt, select = luad.count$barcode)
lusc.exp <- subset(dataFilt, select = lusc.count$barcode)
我們使用 GSVA 包提供的 gsva 函數(shù)來將基因表達(dá)矩陣轉(zhuǎn)換為基因集分?jǐn)?shù)矩陣
library(GSVA)
library(GSEABase)
# 讀取從 GSEA 官網(wǎng)下載的通路數(shù)據(jù)
c2gmt <- getGmt("~/Downloads/data/pathway/c2.cp.v7.2.symbols.gmt")
# 刪選出常用的這三個數(shù)據(jù)庫中的通路
gene.set <- c2gmt[grep("^KEGG|REACTOME|BIOCARTA", names(c2gmt)),]
# gsva 分析，read counts 使用泊松分布菱属，通路至少包含 10 個基因
gs.exp <- gsva(dataFilt, gene.set, kcdf = "Poisson", min.sz = 10)
雖然 GSVAdata 包提供了通路數(shù)據(jù) c2BroadSets 是基因 ID钳榨，但我們的基因表達(dá)數(shù)據(jù)的行是基因 Symbol，所以通路信息也必須是 Symbol 格式纽门，要進(jìn)行格式轉(zhuǎn)換薛耻，比較麻煩
所以我們使用 GSEABase 包提供的 getGmt 函數(shù)來讀取從 GSEA 官網(wǎng)下載的 C2 通路信息
得到結(jié)果如下，共包含 1511 條通路
然后赏陵，使用差異基因識別方法
2. 差異分析
我們使用 limma 分析差異通路
DEA.gs <- TCGAanalyze_DEA(
mat1 = gs.exp[, colnames(luad.exp)],
mat2 = gs.exp[, colnames(lusc.exp)],
metadata = FALSE,
pipeline = "limma",
Cond1type = "LUAD",
Cond2type = "LUSC",
fdr.cut = 0.05,
logFC.cut = 0.5,
)
通過設(shè)置 FDR = 0.05饼齿，logFC = 0.5 共篩選出 40 條差異通路
查看通路的火山圖
我們可以一起看下基因的火山圖
DEA.gene <- TCGAanalyze_DEA(
mat1 = luad.exp,
mat2 = lusc.exp,
metadata = FALSE,
pipeline = "limma",
Cond1type = "LUAD",
Cond2type = "LUSC",
fdr.cut = 0.01,
logFC.cut = 1
)
總共識別出 804 個差異表達(dá)基因
ssGSEA
single sample Gene Set Enrichment Analysis (ssGSEA) 是針對單個樣本進(jìn)行 GSEA 分析，其基因列表的排序方式和 ES 的計算方式都是依賴于樣本中基因的表達(dá)值瘟滨，而不再是依賴基因與表型的相關(guān)度
使用方式也很簡單候醒，只要在 gsva 函數(shù)中指定 method = "ssgsea"，例如
res.ssgsea <- gsva(dataFilt, gene.set, method = "ssgsea", kcdf = "Poisson", min.sz = 10)
也可以進(jìn)行差異分析
DEA.ssgsea <- TCGAanalyze_DEA(
mat1 = res.ssgsea[, colnames(luad.exp)],
mat2 = res.ssgsea[, colnames(lusc.exp)],
metadata = FALSE,
pipeline = "limma",
Cond1type = "LUAD",
Cond2type = "LUSC",
fdr.cut = 0.05,
logFC.cut = 0.1,
)
或者繪制熱圖
annotation_col <- data.frame(sample = rep(1:2, each = 10))
rownames(annotation_col) <- colnames(res.ssgsea)
pheatmap(
res.ssgsea[rownames(DEA.ssgsea),],
show_colnames = F,
# 不展示行名
cluster_rows = F,
# 不對行聚類
cluster_cols = F,
# 不對列聚類
annotation_col = annotation_col,
# 加注釋
cellwidth = 5,
cellheight = 5,
# 設(shè)置單元格的寬度和高度
fontsize = 5
)
單基因富集分析
單基因富集分析并不是說拿單個基因來進(jìn)行富集分析杂瘸，單個基因怎么能進(jìn)行富集分析呢？一個基因根本沒法進(jìn)行統(tǒng)計檢驗伙菊。
其實败玉，這里說的單基因并不是拿單個基因來富集，而是基于單個基因來進(jìn)行富集分析镜硕，這個“基于”运翼，就是以單個基因為基礎(chǔ)，向外擴展兴枯，抓取與其相關(guān)的基因血淌，然后用這些相關(guān)的基因來進(jìn)行功能富集
所以，要理解這個單基因富集分析的意思，這樣一說就已經(jīng)很明了了悠夯。針對單個基因我們可以做什么癌淮？
主要有兩種做法：
定性分組：我們可以根據(jù)給定基因的表達(dá)值對樣本進(jìn)行分組，然后識別在兩組樣本之間差異表達(dá)的基因沦补，最后用這些差異表達(dá)基因來進(jìn)行功能富集
定量相關(guān)：通過計算其他基因與目標(biāo)基因表達(dá)之間的相關(guān)性乳蓄，將具有顯著相關(guān)的基因作為一個集合，也可以進(jìn)行富集分析
1. 定性分組
我們以 CCDC134 基因為例夕膀，以該基因表達(dá)值的中位值來對樣本進(jìn)行分組
gene <- "CCDC134"
gene.exp <- dataFilt[gene,]
label <- if_else(gene.exp < median(gene.exp), 0, 1)
group.low <- dataFilt[,label == 0]
group.high <- dataFilt[,label == 1]
識別兩組樣本之間的差異表達(dá)基因
DEGs <- TCGAanalyze_DEA(
mat1 = group.low,
mat2 = group.high,
metadata = FALSE,
pipeline = "limma",
Cond1type = "CCDC134_Low",
Cond2type = "CCDC134_High",
fdr.cut = 0.01,
logFC.cut = 1,
)
共識別出 873 個差異表達(dá)基因
2. 定量相關(guān)
我們對其他基因與 CCDC134 基因進(jìn)行相關(guān)性檢驗虚倒，由于基因較多，我們使用并行的方式來計算
library(future.apply)
batch_cor <- function(exp, gene){
y = as.numeric(exp[gene,])
gene_list = rownames(exp)
gene_list = gene_list[rownames(exp) != gene]
do.call(rbind, future_lapply(gene_list, function(x){
ct <- cor.test(as.numeric(exp[x,]), y, type='spearman')
data.frame(key = gene, gene = x, cor = ct$estimate,p.value = ct$p.value )
}))
}
plan(multiprocess)
system.time(res.cor <- batch_cor(dataFilt, gene))
對結(jié)果進(jìn)行過濾产舞，篩選出顯著相關(guān)且相關(guān)系數(shù)的絕對值大于 0.6 的基因魂奥，共篩選出 232 個基因
cor.genes <- filter(res.cor, p.value < 0.05 & abs(cor) > 0.6)
3. 富集分析
格式化識別出的差異基因
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
gene.id <- bitr(
rownames(DEGs), fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db
)
go <- enrichGO(
gene = gene.id,
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = T
)
dotplot(go)
gene_info <- DEGs %>%
rownames_to_column(var = "SYMBOL") %>%
inner_join(., gene.id[,1:2], by = "SYMBOL") %>%
# 必須降序
arrange(desc(logFC))
# 構(gòu)造輸入數(shù)據(jù)格式
geneList <- gene_info$logFC
names(geneList) <- as.character(gene_info$ENTREZID)
go2 <- gseGO(
geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "ALL",
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.1,
verbose = FALSE
)
兩種富集方法都沒有富集到 go 通路，對于相關(guān)基因也是沒有富集到通路的易猫。選的這個基因不行
