測(cè)序上游分析系列：
mRNA-seq轉(zhuǎn)錄組二代測(cè)序從raw reads到表達(dá)矩陣：上中游分析pipeline
miRNA-seq小RNA高通量測(cè)序pipeline：從raw reads羹铅，鑒定已知miRNA-預(yù)測(cè)新miRNA舞肆，到表達(dá)矩陣【一】
miRNA-seq小RNA高通量測(cè)序pipeline：從raw reads，鑒定已知miRNA-預(yù)測(cè)新miRNA迁匠，到表達(dá)矩陣【二】
其他文章系列：ggplot2作圖篇：
（1）基于ggplot2的RNA-seq轉(zhuǎn)錄組可視化：總述和分文目錄
（2）測(cè)序結(jié)果概覽：基因表達(dá)量rank瀑布圖调榄，高密度表達(dá)相關(guān)性散點(diǎn)圖卿城，注釋特定基因及errorbar的表達(dá)相關(guān)性散點(diǎn)圖繪制
（3）單/多個(gè)基因在組間同時(shí)展示的多種選擇：分組小提琴圖部服、分組/分面柱狀圖产场、單基因蜂群點(diǎn)圖拼圖的繪制
（4）配對(duì)樣本基因表達(dá)趨勢(shì)：優(yōu)化前后的散點(diǎn)連線圖+拼圖繪制
TCGA轉(zhuǎn)錄組數(shù)據(jù)在完成差異分析后鹅髓，我們通常希望系統(tǒng)地獲取這些成百上千的差異基因的功能信息，幫助我們分析下游實(shí)驗(yàn)的思路京景。面對(duì)大量的差異基因窿冯，逐個(gè)查詢基因功能是不切實(shí)際的。所以我們需要借助基因功能富集分析工具醋粟，獲取有意義的功能基因集靡菇。
基因功能富集分析的基本思路是一致的：開發(fā)者將具有相似功能的基因劃分到同一個(gè)基因集中重归，然后使用統(tǒng)計(jì)檢驗(yàn)分析測(cè)序獲得的差異基因是否顯著聚于這些已知功能的基因集中。若檢驗(yàn)結(jié)果顯著厦凤，則可說(shuō)明差異基因可能參與了這一基因集所注釋的生物功能鼻吮，為我們的下游研究提供了系統(tǒng)精準(zhǔn)的思路。目前常用的基因富集分析方法有超幾何檢驗(yàn)富集分析较鼓、Gene Set Enrichment Analysis (GSEA)和Gene Set Variation Analysis (GSVA) （或功能類似的single sample Gene Set Enrichment Analysis (ssGSEA)）椎木。
如何從TCGA數(shù)據(jù)庫(kù)中獲得差異基因？我們從這里開始：
- 從TCGA數(shù)據(jù)庫(kù)下載并整合清洗高通量腫瘤表達(dá)譜-臨床性狀數(shù)據(jù)
-
TCGA數(shù)據(jù)整合后進(jìn)行DESeq2差異表達(dá)分析和基于R的多種可視化
下面使用的對(duì)象名稱和內(nèi)容與前文保持一致博烂，仍使用之前教程選擇的的TCGA-LUSC白色人種肺鱗癌表達(dá)譜數(shù)據(jù)（整合后cancer=344, normal=42）香椎。
進(jìn)行基因富集分析，我們需要的R對(duì)象是：
（1）dds_DE: this object created by DESeq2 contains results of differential expression analysis. We need to extract log2FoldChange, gene symbol and other available information from it.
（2）condition_table: this table contains the grouping info of all samples. It will be used to create a new grouping table compatible with Limma package.
（3）mRNA_exprSet: the expression matrix containing columns 'gene_name' and 'gene_id' representing GENE SYMBOL and ENSEMBL ID. They will be used for ID transfer.
注：這里的mRNA_exprSet相比第一篇教程整合得到的矩陣多了一行'gene_id'禽篱。因?yàn)榭紤]到基因的各種ID之間并沒(méi)有非常完美的對(duì)應(yīng)畜伐，尤其是SYMBOL（基因的各種名字版本更新或許不齊），所以既然在整合的時(shí)候已經(jīng)進(jìn)行過(guò)一次SYMBOL-ENSEMBL ID轉(zhuǎn)換躺率，盡可能還是用相同批次的數(shù)據(jù)玛界。其實(shí)差別不大：
#the last filtering step I used in the tutorial before:
mRNA_exprSet <- mRNA_exprSet[,c('gene_name', condition_table$TCGA_IDs)]
#We just don't need to remove 'gene_id' column:
mRNA_exprSet <- mRNA_exprSet[,c('gene_name', 'gene_id', condition_table$TCGA_IDs)]
1. 超幾何檢驗(yàn)GO、KEGG基因富集分析
這是相對(duì)簡(jiǎn)單粗暴一些的基因富集分析方法悼吱。不需要輸入基因的表達(dá)值慎框，只需要通過(guò)DESeq2和自己設(shè)定的閾值（如|log2FoldChange| > 2 & FDR < 0.05）篩選得到的差異基因名。進(jìn)行統(tǒng)計(jì)檢驗(yàn)后返回顯著富集的功能基因集后添。
GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）是兩種常用的數(shù)據(jù)庫(kù)笨枯。GO中的功能數(shù)據(jù)集分為3個(gè)sub-classes: Biological Process (BP), Molecular Function (MF) 和 Cellular Component (CC)，分別含有生物學(xué)功能活動(dòng)遇西、分子功能和細(xì)胞組分定位的基因集合（GO terms）馅精。KEGG包含的基因集合均為在某特定生物學(xué)功能過(guò)程中發(fā)揮作用的分子信號(hào)通路（KEGG pathways）。
（1）需要R包努溃。
DESeq2 (獲得差異基因信息), clusterProfiler（ID轉(zhuǎn)換+富集分析+作圖一站式神包Ａ蛩弧）。
install.packages('DESeq2')
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler") #this package can be downloaded from Bioconductor.
library(DESeq2)
library(clusterProfiler)
（2）提取差異基因列表梧税。
根據(jù)我們指定的硬閾值提取出差異表達(dá)基因的基因SYMBOL（做差異分析時(shí)輸入的矩陣中rownames為SYMBOL）沦疾。
res_DE <- results(dds_DE, alpha=0.05, contrast=c("sample_type","cancer","normal"))
resSig <- subset(res_DE,padj < 0.05)
resSigAll <- subset(resSig, (log2FoldChange < (-2)|log2FoldChange > 2))
gene_list_for_GO_KEGG <- resSigAll@rownames
我們一般會(huì)遇到的基因ID類型有：ENSEMBL GENE ID （如ENSG00000186092），GENE SYMBOL（如OR4F5）和ENTREZID（如79501）第队。絕大部分的基因ID都是能夠一一對(duì)應(yīng)的哮塞，但自然是不排除有例外...
clusterProfiler包的KEGG分析函數(shù)只接受基因的ENTREZID而不是SYMBOL，所以我們進(jìn)行基因ID的轉(zhuǎn)換凳谦。而且一般由于SYMBOL版本繁雜忆畅，也建議用ENTREZID進(jìn)行富集分析比較穩(wěn)妥，防止過(guò)多SYMBOL不能和數(shù)據(jù)庫(kù)匹配而導(dǎo)致分析結(jié)果的不可信尸执。ENSEMBL和SYMBOL之間的轉(zhuǎn)換可以直接借助第一篇數(shù)據(jù)整合教程中的gtf文件來(lái)完成家凯。但其中不含對(duì)于富集分析較為關(guān)鍵的ENTREZID缓醋。這里我們使用bitr函數(shù)可以進(jìn)行SYMBOL, ENSEMBL ID和ENTREZID的在線轉(zhuǎn)換。
可想而知绊诲，編號(hào)之間的對(duì)應(yīng)關(guān)系更好送粱。所以我習(xí)慣用ENSEMBL ID和ENTREZID互換。這里比較懶掂之，普通的GO/KEGG分析就直接用SYMBOL轉(zhuǎn)換了抗俄，只是可能轉(zhuǎn)換丟失率會(huì)稍高一些。但丟失非常少部分的基因是合理而且問(wèn)題不大的世舰。
#from/toType = 'SYMBOL'/'ENTREZID'/'ENSEMBL'
gene_list_for_GO_KEGG <- clusterProfiler::bitr(gene_list_for_GO_KEGG, fromType = 'SYMBOL',
toType = 'ENTREZID', OrgDb = 'org.Hs.eg.db')
于是便獲得了具有對(duì)應(yīng)信息的table和提示有基因轉(zhuǎn)換丟失Warning message：
（3）GO單種subclass的富集分析动雹。
接下來(lái)就可以用enrichGO函數(shù)分析和可視化了。先只做Biological Process亞類的分析跟压。（還是懶胰蝠，就直接用SYMBOL做GO分析了）
注：函數(shù)keyType指輸入的ID類型，OrgDb為數(shù)據(jù)庫(kù)種屬裆馒，ont為subclass類型（BP/CC/MF/All）
GO_BP <- enrichGO(gene = gene_list_for_GO_KEGG$SYMBOL, keyType = 'SYMBOL', OrgDb = 'org.Hs.eg.db', ont = 'BP', pvalueCutoff = 0.05)
barplot(GO_BP, showCategory = 20, title = 'Top20 GO-BP terms', x = 'Count', color = 'p.adjust') #ont='BP'/'MF'/'CC', x='Count'/'GeneRatio'
dotplot(GO_BP, title = 'Top20 GO-BP terms', showCategory = 20, color = 'p.adjust')
結(jié)果解釋：一般可視化top10/20的富集基因集合姊氓。barplot x軸設(shè)定為在該基因集合中的enriched gene counts，顏色代表FDR喷好，排序根據(jù)FDR。dotplot橫軸設(shè)定GeneRatio = m/n（m為在某顯著子基因集中的輸入差異基因個(gè)數(shù)读跷，n為能在總BP基因集中檢測(cè)到的輸入差異基因數(shù)）梗搅，氣泡大小為絕對(duì)enriched gene counts數(shù)，顏色代表FDR效览，排序是根據(jù)gene ratio无切。個(gè)人比較喜歡氣泡圖，多一個(gè)維度的數(shù)據(jù)丐枉，而且排序根據(jù)富集程度而不是FDR顯著程度哆键。（類比一下普通RT-PCR實(shí)驗(yàn)：我們通常只需要p值顯著<0.05，相比之下更關(guān)心的絕對(duì)是fold change）
不過(guò)本人認(rèn)為clusterProfiler包dotplot默認(rèn)的x軸排序方法GeneRatio仍舊不是最好的辦法瘦锹，以enrichment score（存在于某子基因集中的輸入差異基因/該子基因集的總基因數(shù)）或許比較客觀籍嘹，之后可以使用ggplot2自己實(shí)現(xiàn)一下。
對(duì)GO三個(gè)subclasses的基因集合同時(shí)做富集分析：
GO_all <- enrichGO(gene = gene_list_for_GO_KEGG$SYMBOL,ont = 'All' , keyType = 'SYMBOL', OrgDb = 'org.Hs.eg.db', pvalueCutoff = 0.05)
dotplot(GO_all, title = 'Top10 GO terms of each sub-class', showCategory = 10, color = 'p.adjust', split='ONTOLOGY')+ facet_grid(ONTOLOGY~.,scale="free")
可以獲得分面的富集圖弯院，每個(gè)subclass選取Top10可視化：
（4）KEGG富集分析辱士。
KEGG <- enrichKEGG(gene_list_for_GO_KEGG$ENTREZID, organism = 'hsa', pvalueCutoff = 0.05)
dotplot(KEGG, title = 'Top20 KEGG pathways', showCategory = 20, color = 'p.adjust')
品味一下KEGG的富集通路，是不是內(nèi)味和GO略微不同呢~~
2. GSEA富集分析
相比于第一種簡(jiǎn)單粗暴的用硬閾值截取+往籃子里塞雞蛋的方法听绳，GSEA同時(shí)考慮了基因在整個(gè)表達(dá)譜中所處的FoldChange rank以及同一基因集中的基因在表達(dá)譜rank中的距離颂碘。通俗來(lái)講，GSEA基于如下假設(shè)：一個(gè)基因集中的基因如果在表達(dá)譜中所處的rank越極端（高/低FoldChange）椅挣，而且基因之間的距離越短（rank相近）头岔，則這個(gè)基因集越可能顯著塔拳。所以GSEA需要輸入的是【所有納入差異分析的完整基因list和它們的FC值，并已經(jīng)預(yù)先排序（pre-ranked）】峡竣。
（1）需要R包：
DESeq2 (獲得差異基因信息)蝙斜，clusterProfiler（ID轉(zhuǎn)換+富集分析+作圖一站式神包！）澎胡，dplyr（表格操作）孕荠，ggplot2（修飾clusterProfiler的作圖），pheatmap（熱圖繪制）攻谁。
install.packages('dplyr)
install.packages('DESeq2')
install.packages('ggplot2')
install.packages('pheatmap')
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler") #this package can be downloaded from Bioconductor.
library(DESeq2)
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(pheatmap)
（2）獲戎晌椤（幾乎）所有基因的log2FC數(shù)據(jù)并進(jìn)行ID轉(zhuǎn)換。
首先使用results函數(shù)獲得所有納入差異分析的基因信息戚宦。（這函數(shù)一定得設(shè)置0-1間的陽(yáng)性錯(cuò)誤率个曙，那就多打幾個(gè)9吧，盡可能把完整基因譜都提取出來(lái)）
GSEA還是用ENTREZID好好做受楼，所以祭出了之前用到的ENSEMBL ID和SYMBOL互換的mRNA_exprSet[,c(1,2)]垦搬，至少可以保證SYMBOL到ENSEMBL ID的百分之百轉(zhuǎn)換率。隨后使用bitr將ENSEMBL對(duì)應(yīng)到ENTREZID艳汽，去除重復(fù)后最終與log2FC數(shù)據(jù)合并猴贰。
dds_for_GSEA <- results(dds_DE, contrast=c("sample_type","cancer","normal"), alpha = (0.9999999999999999))
lFC_for_GSEA <- data.frame('gene_name'=dds_for_GSEA@rownames, 'lFC'=dds_for_GSEA$log2FoldChange, stringsAsFactors = F)
ENS_for_GSEA <- dplyr::filter(mRNA_exprSet, gene_name %in% dds_for_GSEA@rownames)[,c('gene_name','gene_id')] %>%
dplyr::inner_join(lFC_for_GSEA, by='gene_name')
ENTREZ <- clusterProfiler::bitr(ENS_for_GSEA$gene_id, fromType = 'ENSEMBL', toType = c('ENSEMBL','ENTREZID'), OrgDb = 'org.Hs.eg.db')
names(ENTREZ) <- c('gene_id', 'ENTREZID')
ENTREZ <- ENTREZ[!duplicated(ENTREZ$gene_id),]
ENTREZ <- ENTREZ[!duplicated(ENTREZ$ENTREZID),] %>% dplyr::inner_join(ENS_for_GSEA, 'gene_id')
（3）建立輸入對(duì)象并分析表達(dá)譜在GO BP和KEGG中的富集程度。
提取對(duì)應(yīng)上ENTREZID的log2FC vector并以ID作為vector各個(gè)元素的名字河狐。按lFC降序排序后輸入到專門用于GO數(shù)據(jù)庫(kù)GSEA分析的gseGO函數(shù)米绕，和專用于KEGG數(shù)據(jù)庫(kù)GSEA分析的gseKEGG函數(shù)中進(jìn)行分析。
注：nPerm指permutation重復(fù)檢驗(yàn)富集是否隨機(jī)的次數(shù)馋艺，建議500栅干，1000或更多。keyType中’ncbi-geneid‘就是指ENTREZID捐祠。
input_GSEA <- ENTREZ$lFC
names(input_GSEA) <- ENTREZ$ENTREZID
input_GSEA <- sort(input_GSEA, decreasing = T)
GSEGO_BP <- gseGO(input_GSEA, ont = 'BP', OrgDb = 'org.Hs.eg.db', nPerm = 1000, pvalueCutoff = 0.05)
GSEA_KEGG <- gseKEGG(input_GSEA, organism = 'hsa', keyType = 'ncbi-geneid', nPerm = 1000, pvalueCutoff = 0.05)
（4）GSEA結(jié)果的可視化和解讀碱鳞。
將分析獲得的對(duì)象壓縮到扁平的dataframe，判定富集基因集踱蛀，并根據(jù)NES排序窿给。
GSEA中判斷基因集是否富集一般取決于如下參數(shù)：
（1）NES（normalized enrichment score）: 絕對(duì)值>1為富集，嚴(yán)苛可>1.5星岗。
（2）p值：<0.05填大，嚴(yán)苛可<0.01。
（3）q值：<0.25俏橘，嚴(yán)苛可<0.1/0.05/0.01允华。
找出有興趣進(jìn)行可視化的顯著基因集，并使用gseaplot和gseaplot2分別可視化。（我這里用ggplot2稍微加了個(gè)字幕~）
#further filter out the significant gene sets and order them by NES scores.
#GO
GSEA_BP_df <- as.data.frame(GSEGO_BP) %>% dplyr::filter(abs(NES)>1 & pvalue<0.05 & qvalues<0.25)
GSEA_BP_df <- GSEA_BP_df[order(GSEA_BP_df$NES, decreasing = T),]
p <- gseaplot(GSEGO_BP, GSEA_BP_df$ID[1], by='all', title = GSEA_BP_df$Description[1], color.vline = 'gray50', color.line='red', color='black') #by='runningScore/preranked/all'
p <- p+annotate(geom = 'text', x = 0.87, y =0.85, color='red', fontface = 'bold', size= 4,
label=paste0('NES= ', round(GSEA_BP_df$NES[1], 1), '\n', 'p.adj= ', round(GSEA_BP_df$p.adjust[1], 2)))+
theme(panel.grid = element_line(colour = 'white'))
p
#KEGG
GSEA_KEGG_df <- as.data.frame(GSEA_KEGG) %>% dplyr::filter(abs(NES)>1 & pvalue<0.05 & qvalues<0.25)
GSEA_KEGG_df <- GSEA_KEGG_df[order(GSEA_KEGG_df$NES, decreasing = T),]
p <- enrichplot::gseaplot2(GSEA_KEGG, geneSetID = GSEA_KEGG_df$ID[1], pvalue_table = F, ES_geom = 'line')+
annotate('text', x = 0.87, y =0.85, color='red', fontface = 'bold', size= 4,
label=paste0('NES= ', round(GSEA_KEGG_df$NES[1], 1), '\n', 'p.adj= ', round(GSEA_KEGG_df$p.adjust[1], 2)))+
labs(title = GSEA_KEGG_df$Description[1])+theme(plot.title = element_text(hjust = 0.5))
p
兩種風(fēng)格任選其一即可靴寂。左圖的富集結(jié)果是相當(dāng)相當(dāng)好了磷蜀。基因集中大部分基因都聚集在表達(dá)譜的正頂端百炬，且基因rank距離近褐隆，曲線呈明顯的偏分布正峰。
（5）查看GSEA顯著基因集中基因在樣本中的具體表達(dá)情況（熱圖可視化）剖踊。
得到了顯著的基因集庶弃，但最終我們的落腳點(diǎn)還是在具體基因上。所以或許需要對(duì)我們感興趣的基因集中參與富集的基因的表達(dá)情況進(jìn)行可視化德澈。
這里仍用GO BP富集NES最高的keratinization演示歇攻。熱圖沒(méi)有進(jìn)行聚類，col按樣本類別分類梆造，row按log2FC降序排列缴守，比較直觀清晰。
#to visualize the expression trends of genes within a significant GSEA geneset.
ENTREZID_within_a_set <- unlist(strsplit(GSEA_BP_df$core_enrichment[1], split = '/', fixed = T))
symbols_within_a_set <- dplyr::inner_join(ENTREZ, data.frame('ENTREZID'=ENTREZID_within_a_set, stringsAsFactors = F), by='ENTREZID')[,c('gene_name','lFC')]
symbols_within_a_set <- symbols_within_a_set[order(symbols_within_a_set$lFC, decreasing = T),]
sub_expr_norm <- expr_norm[symbols_within_a_set$gene_name,]
annotation_col <- data.frame('Sample_type'=condition_table$sample_type)
rownames(annotation_col) <- condition_table$TCGA_IDs
pheatmap::pheatmap(t(scale(t(sub_expr_norm))), show_colnames = F, show_rownames = T,
annotation_col = annotation_col, cluster_cols = F, cluster_rows = F,
color = colorRampPalette(c("navy", "white", "firebrick3"))(100), fontsize = 4,
annotation_colors = list(Sample_type=c(cancer = 'orange', normal = 'blue')),
legend = T, legend_breaks = c(-2, 0, 2), breaks = unique(seq(-2, 2, length=100)),
lagend_labels = c('≤2','0','≥2'), main = paste0('Heatmap of enriched genes within ', GSEA_BP_df$Description[1]))
（6）自建基因集合進(jìn)行GSEA分析镇辉。
前面使用了函數(shù)內(nèi)置好的基因集合進(jìn)行GSEA屡穗。如果我們需要自定義基因集合怎么辦？GSEA函數(shù)可以勝任忽肛，它接受表達(dá)情況和自定義基因集作為輸入村砂，進(jìn)行GSEA分析。
作為測(cè)試麻裁，我去GSEA官網(wǎng)Molecular Signatures Database下載了hallmark gene sets的entrez gene id基因集文件（.gmt）箍镜。
用文本編輯器打開.gmt文件，發(fā)現(xiàn)就是我們理解的基因集格式：一個(gè)功能描述+網(wǎng)站廣告后煎源，緊接所有包含在本功能集合中的基因ID。
而GSEA函數(shù)要求的database輸入也是一定的：包含兩列的表格：第一列為功能基因集的名稱香缺；第二列為基因ID手销。
所以接下來(lái)是geneset格式的建立。
original_gmt <- readLines('h.all.v7.0.entrez.gmt')
strsplit_name <- function(gmt.list_layer){
GSgenes <- as.character(unlist(strsplit(gmt.list_layer, split = '\t', fixed = T)))
data.frame('Genesets'=rep(GSgenes[1], (length(GSgenes)-2)), 'Genes'=GSgenes[c(-1,-2)], stringsAsFactors = F)
}
database_list <- lapply(original_gmt, strsplit_name)
db_for_GSEA <- do.call(rbind, database_list)
db_for_GSEA$Genesets <- as.character(db_for_GSEA$Genesets)
db_for_GSEA$Genes <- as.character(db_for_GSEA$Genes)
同前图张，分析并可視化锋拖。gseaplot2可以同時(shí)可視化多個(gè)genesets敌买。不過(guò)感覺(jué)用處不是太大逻澳。
GSEA_test <- clusterProfiler::GSEA(input_GSEA, nPerm = 1000, pvalueCutoff = 0.05, TERM2GENE = db_for_GSEA)
GSEA_test_df <- as.data.frame(GSEA_test) %>% dplyr::filter(abs(NES)>1 & pvalue<0.05 & qvalues<0.25)
GSEA_test_df <- GSEA_test_df[order(GSEA_test_df$NES, decreasing = T),]
enrichplot::gseaplot2(GSEA_test, GSEA_test_df$Description[c(1:2, (length(GSEA_test_df$Description)-1):length(GSEA_test_df$Description))], color = c('red','orange','blue','navy'))
3. GSVA/ssGSEA分析
ssGSEA是為無(wú)重復(fù)的樣本進(jìn)行g(shù)eneset enrichment analysis準(zhǔn)備的，所以不同于上方以組別為單位（cancer vs normal）的GSEA分析涡相，通過(guò)ssGSEA适袜，每個(gè)樣本都可以得到相應(yīng)基因集的評(píng)分柄错。GSVA的原理和作用類似，所以GSVA和ssGSEA被寫入了同一個(gè)R包中，性能等同售貌。
所以我們甚至可以對(duì)每個(gè)樣本的pathway scores進(jìn)行差異分析给猾，從而獲得組間差異表達(dá)的pathways。這又是另外一種船新體驗(yàn)了有沒(méi)有~現(xiàn)在興起的新玩法還包括使用ssGSEA對(duì)腫瘤測(cè)序結(jié)果進(jìn)行免疫細(xì)胞特征基因的評(píng)分颂跨，從而鑒定腫瘤免疫浸潤(rùn)的程度和特征敢伸。
（1）需要R包：
GSVA包（進(jìn)行GSVA/ssGSEA分析），limma包（差異pathway的篩選）恒削，pheatmap包（熱圖繪制）池颈。
source("https://bioconductor.org/biocLite.R")
biocLite(c("GSVA",'limma'))
install.packages('pheatmap')
library(GSVA)
library(limma)
library(pheatmap)
（2）自建用于GSVA/ssGSEA分析的基因集合。
GSVA包不提供預(yù)設(shè)基因集钓丰，所以我們得自己建房子躯砰。本例中使用的原始基因集.gmt文件依舊來(lái)自于GSEA官網(wǎng)Molecular Signatures Database。使用的基因集為C2: curated gene sets-Canonical pathways-KEGG subset of CP斑粱。
然后進(jìn)行一番整理弃揽，最后成為GSVA需要的list類型的輸入基因集：list的每一層為一個(gè)基因集，層名為基因集合的名稱则北。（唉還是動(dòng)用了for循環(huán)它老人家矿微，俺仍需努力鴨。
#prepare the database imported from .gmt file.
original_gmt_GSVA <- readLines('c2.cp.kegg.v7.0.entrez.gmt')
strsplit_no_name <- function(gmt.list_layer){
as.character(unlist(strsplit(gmt.list_layer, split = '\t', fixed = T)))[-2]
}
database_list_GSVA <- lapply(original_gmt_GSVA, strsplit_no_name)
for (layers in 1:length(database_list_GSVA)) {
names(database_list_GSVA)[layers] <- database_list_GSVA[layers][[1]][1]
database_list_GSVA[layers][[1]] <- database_list_GSVA[layers][[1]][-1]
}
注：GSVA的基因集要求和GSEA好不一樣尚揣！看圖說(shuō)話涌矢！
（3）基因表達(dá)譜的選擇和ID轉(zhuǎn)換。
GSVA分析需要輸入【expression matrix】而不僅僅是差異基因名/Log2FC數(shù)據(jù)進(jìn)行分析快骗。而且其對(duì)表達(dá)譜的分布有一定要求娜庇。
官方文檔說(shuō)明為：【對(duì)于RNA-seq，如使用原始整數(shù)counts矩陣作為輸入方篮，其分布近似泊松名秀，GSVA函數(shù)的kcdf參數(shù)選擇'Poisson'；若為log2轉(zhuǎn)換后的FPKM或已經(jīng)用DESeq2標(biāo)準(zhǔn)化藕溅，近似符合高斯分布（正態(tài)分布）的數(shù)據(jù)匕得，則kcdf參數(shù)選擇'Gaussian'〗肀恚】
對(duì)于芯片數(shù)據(jù)直接選擇'Gaussian'（過(guò)時(shí)不看）汁掠。
關(guān)于數(shù)據(jù)的分布，這里用表達(dá)譜中隨便選取的一個(gè)基因做一下樣本間基因表達(dá)水平的密度圖展示：
圖1：原始counts分布（近似Poisson）集币；圖2：DESeq2標(biāo)準(zhǔn)化+vst轉(zhuǎn)換（近似Gaussian）考阱；圖3：僅FPKM標(biāo)準(zhǔn)化（近似Poisson）；圖4：log2(FPKM+1)（近似Gaussian）鞠苟。
其實(shí)目前認(rèn)為原始數(shù)據(jù)的分布更符合負(fù)二項(xiàng)而不是泊松分布乞榨。Anyway秽之，R包說(shuō)明文檔上說(shuō)行就行吧，總之是一種近似姜凄。
再品味一下上面的參數(shù)選擇政溃，有感覺(jué)了嗎~
我們首先選擇DESeq2+vst轉(zhuǎn)換表達(dá)矩陣進(jìn)行GSVA：ID轉(zhuǎn)換。用到了GSEA時(shí)建立的ENTREZ轉(zhuǎn)換表态秧。
#prepare the normalized counts matrix (rownames are ENTREZIDs of genes)
expr_norm_for_GSVA <- assay(vst(dds_DE, blind = T))
expr_norm_for_GSVA <- data.frame('gene_name'=rownames(expr_norm_for_GSVA), expr_norm_for_GSVA, stringsAsFactors = F)
expr_norm_for_GSVA <- dplyr::inner_join(ENTREZ, expr_norm_for_GSVA, by='gene_name')[,c(-1,-3,-4)]
rownames(expr_norm_for_GSVA) <- expr_norm_for_GSVA$ENTREZID
colnames(expr_norm_for_GSVA) <- gsub(colnames(expr_norm_for_GSVA), pattern = '\\.',replacement='-')
expr_norm_for_GSVA <- expr_norm_for_GSVA[,-1]
（4）GSVA分析董虱。
我們使用GSVA方法來(lái)分析（如果想用ssGSEA， method輸入'ssGSEA'）申鱼。GSVA輸入的表達(dá)矩陣必須是matrix格式愤诱，別忘了轉(zhuǎn)換。kcdf參數(shù)這里顯然就是'Caussian'了捐友。最后一個(gè)參數(shù)如果電腦好可以設(shè)大一點(diǎn)淫半。
#the input data must be in matrix format.
#We should set the distribution as 'Poisson' when using RNA-seq integer counts data rather than the default 'Gaussian'.
#If the matrix is derived from microarray or the RNA-seq integer counts have been transformed to log2 form or..., we can directly set kcdf as 'Gaussian'.
es <- GSVA::gsva(as.matrix(expr_norm_for_GSVA), database_list_GSVA,
mx.diff=FALSE, verbose=FALSE, method='gsva', kcdf='Gaussian',
parallel.sz=1)
（5）limma包篩選差異通路。
limma包匣砖，和edgeR科吭，DESeq2并駕齊驅(qū)為三大差異分析R包。
limma篩選所需要的東西不多：
（1）es: the pathway score we got from GSVA.
（2）condition_table_for_limma: the matrix containing the grouping info. It contains 2 columns (1 for normal, another for cancer) in this tutorial. Samples (recorded by rownames) affilifated to this group should be 1, otherwise they are tagged with 0.
（3）contrast_matrix: tell limma which the reference group is.
接下來(lái)就是線性擬合以及用經(jīng)驗(yàn)貝葉斯法計(jì)算p值猴鲫。最后我們提取|Log2FC|>0.5对人，p值<0.05的顯著pathways。
#cancer=1, normal=0.
condition_table_for_limma <- model.matrix(~0+ condition_table$sample_type)
colnames(condition_table_for_limma) <- levels(condition_table$sample_type)
rownames(condition_table_for_limma) <- condition_table$TCGA_IDs
#generate the contrast matrix. cancer=1, normal=-1. It means that genes in normal group will be at denominator level.
contrast_matrix <- makeContrasts('cancer-normal',levels = condition_table$sample_type)
lmfit <- lmFit(es, condition_table_for_limma)
lmfit <- contrasts.fit(lmfit, contrast_matrix)
lmfit <- eBayes(lmfit)
#get significant DE pathways.
sigPathways <- topTable(lmfit, coef=T, p.value=0.05, adjust="BH", lfc = 0.5)
（6）熱圖繪制拂共。
annotation_col <- data.frame('Sample_type'=condition_table$sample_type)
rownames(annotation_col) <- condition_table$TCGA_IDs
pheatmap::pheatmap(t(scale(t(es[rownames(sigPathways),]))), show_colnames = F, show_rownames = T,
annotation_col = annotation_col, cluster_cols = F,
color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
annotation_colors = list(Sample_type=c(cancer = 'orange', normal = 'blue')))
結(jié)果說(shuō)明這組肺鱗癌患者的腫瘤細(xì)胞可能存在細(xì)胞周期失控牺弄、細(xì)胞分裂活躍、DNA損傷和修復(fù)高度活動(dòng)的問(wèn)題宜狐。這都是老生常談的腫瘤特征了势告，也間接說(shuō)明這套方法學(xué)是非常有效的。
*Optional：不同數(shù)據(jù)類型進(jìn)行GSVA分析的比較
根據(jù)之前GSVA對(duì)數(shù)據(jù)的要求抚恒，我們知道多種數(shù)據(jù)類型都可以做這套分析咱台。接下來(lái)我使用原始counts矩陣，以及z score of log2(FPKM+1)矩陣做GSVA俭驮。
#I just show the difference in GSVA analysis. The codes for data processing (like ID conversion, and log2/z_score calculation) have been omitted.
#Integer counts
es_counts <- GSVA::gsva(as.matrix(mRNA_expr_for_DESeq), database_list_GSVA,
mx.diff=FALSE, verbose=FALSE, method='gsva', kcdf='Poisson',
parallel.sz=4)
#z_score(log2(FPKM+1))
es_FPKM_zscore <- GSVA::gsva(mRNA_exprSet_FPKM_log2zscore, database_list_GSVA,
mx.diff=FALSE, verbose=FALSE, method='gsva', kcdf='Gaussian',
parallel.sz=1)
最后得到的差異通路比較如下吵护。雖然略有差異，但是大部分差異通路以及聚類情況都是高度統(tǒng)一的表鳍，說(shuō)明不同數(shù)據(jù)類型之間的可重復(fù)性。大家隨意選用祥诽。
寫了一天...已die...如果覺(jué)得好用就點(diǎn)個(gè)贊8譬圣！
