https://www.sohu.com/a/281843011_120055884
http://www.itdecent.cn/p/cff8fe8b548a
之前介紹的一種是用R包做的圖，利用的是log2foldchange 和ensemble ID
這次介紹的是GESA軟件操作
原始分析數(shù)據(jù)
原始的FPKM值為
基因名（需要是ENTREZID格式）

image.png

GSEA軟件下載：https://www.gsea-msigdb.org/gsea/index.jsp 直接選擇你操作系統(tǒng)對(duì)應(yīng)的版本，下載安裝即可。
準(zhǔn)備上傳文件

1.基因表達(dá)量表（.gct）

通常用.gct為后綴。文件第一行以“#1.2”開頭；文件第二行的第一列為基因個(gè)數(shù)、第二列為樣品個(gè)數(shù)；文件的第三行為表達(dá)譜的矩陣的title信息，第一列為基因symbol/探針號(hào)，第二列為基因／探針的描述信息，第三列以后為樣品id。接下來(lái)的行對(duì)應(yīng)每個(gè)基因／探針在每個(gè)樣品中的表達(dá)信息。文件以tab作為分隔符。

image.png

第一列不能有重復(fù)

image.png

2. 樣品表型分類文件(cls)——必需文件

樣品表型分類文件需以.cls為后綴。文件第一行為三個(gè)數(shù)字，第一個(gè)是樣品的總數(shù)，第二個(gè)是樣品分為幾類，第三個(gè)數(shù)字通常為1。第二行也通常三個(gè)字符串，第一個(gè)為#，第二個(gè)為分類1的名稱，第三個(gè)位分類2的名稱。第三行為每個(gè)樣品的分類信息，0代表分類1，1則代表分類2。文件以空格或者tab分割。

image.png

GSEA 分析

導(dǎo)入數(shù)據(jù)
打開GSEA軟件，點(diǎn)擊左側(cè)菜單欄的Load data選項(xiàng)，將數(shù)據(jù)按如下圖所示方法導(dǎo)入

image.png

設(shè)置參數(shù)
按照下圖標(biāo)記的批注設(shè)置合適的參數(shù)。

image.png

image.png

image.png

Gene sets database （功能基因集數(shù)據(jù)庫(kù)）：GSEA包含了MSigDB數(shù)據(jù)庫(kù)中的功能基因集，可以從中選擇感興趣的通路、癌癥標(biāo)記、轉(zhuǎn)錄因子數(shù)據(jù)庫(kù)等。
Number of permutations（擾動(dòng)／隨機(jī)次數(shù)）：通常設(shè)置1000，此參數(shù)不可過(guò)小。
Phenotypes labels(樣品表型分類文件)：選擇上一步上傳的表型cls文件
Collapse dataset to gene symbols：通常true

如果已經(jīng)變?yōu)榱薵ene symobols 則不需要進(jìn)行轉(zhuǎn)換
如果不是,則需要進(jìn)行轉(zhuǎn)換
Permutation type(擾動(dòng)類型)： 通常選擇phenotype，如果樣品數(shù)目較少可以選擇gene_set。
Chip platform(芯片類型)：如果表達(dá)gct文件的第一列為芯片探針id則此處需要選擇對(duì)應(yīng)的芯片平臺(tái)，如果是基因symbol則無(wú)需選擇。

image.png

image.png

image.png

gsea分析結(jié)束后

點(diǎn)擊相應(yīng)的success

image.png

跳到相應(yīng)的分析結(jié)果

image.png

如果是其他物種,http://www.itdecent.cn/p/cff8fe8b548a

可以構(gòu)建心得基因集
<meta charset="utf-8">

基因集文件的準(zhǔn)備相對(duì)比較復(fù)雜，我用到的是KEGG的通路基因集，也可以用GO的的基因集，甚至可以自己創(chuàng)造一個(gè)，只需要滿足上圖所示的文件格式。下面介紹一下如何從KEGG網(wǎng)站上下載某物種（小鼠）的基因集信息。

首先，在KEGG數(shù)據(jù)庫(kù)中下載小鼠的KEGG通路數(shù)據(jù)庫(kù)文件

image.png

image.png

image.png

image.png

image.png

大家可以看到，這個(gè).keg文件和我們需要的.gmt文件還有不小的差異，那么，我們需要按如下方法，先在Linux下解析keg文件，然后用R制作.gmt文件：

#解析信息
perl -alne '{if(/^C/){/PATH:mmu(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' mmu00001.keg > mmu00001_kegg2gene.txt
grep '^C' mmu00001.keg | grep "mmu" | sed 's/^C    //g' | sed 's/\([0-9]*\)\s/\1\t/' > mmu00001_kegg2description.txt

#轉(zhuǎn)換ID同時(shí)制作.gmt文件
#安裝轉(zhuǎn)換ID需要用到的各個(gè)R包
#安裝data.table
install.packages("data.table")
#安裝org.Xx.eg.db系列包中小鼠資源包
BiocManager::install("org.Mm.eg.db")
#安裝clusterProfiler包，安裝時(shí)需要“科學(xué)上網(wǎng)”
BiocManager::install("clusterProfiler")

#調(diào)用各個(gè)R包
library(data.table)
library(org.Mm.eg.db)
library(clusterProfiler)

#讀入數(shù)據(jù)
path2gene_file<-"mmu00001_kegg2gene.txt"
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
dt <- tmp
dt$geneid <- rownames(dt)

# transform id  
map_dt <- bitr(dt$V2, fromType = "ENTREZID",toType = c( "SYMBOL"),OrgDb = org.Mm.eg.db)
dt_merge <- merge(map_dt,dt, by.y = "V2", by.x = "ENTREZID")

# 可以把轉(zhuǎn)換ID后的表輸出
#write.table(dt_merge, file = "example42.txt", sep = "\t", col.names = T, row.names = T, quote = F)

# first column is kegg ID, second column is entrez ID
GeneID2kegg_list<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
kegg2symbol_list<- tapply(dt_merge$SYMBOL,as.factor(dt_merge$V1),function(x) x)

#輸出數(shù)據(jù)，但描述列為NA
write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
  sink( gmt_file )
  for (i in 1:length(geneSet)){
    cat(names(geneSet)[i])
    cat('\tNA\t')
    cat(paste(geneSet[[i]],collapse = '\t'))
    cat('\n')
  }
  sink()
}

write.gmt()

image.png

另外一種方法，構(gòu)建斑馬魚物種的bp基因集

setwd("I:/ZHS_Seq_data/RUN")
library(msigdbr)
msigdbr_species() #可以看到，這個(gè)包涵蓋了20個(gè)物種

zebrafish <- msigdbr(species = "Danio rerio")
zebrafish[1:5,1:5]
table(zebrafish$gs_cat) #查看目錄，與MSigDB一樣，包含9個(gè)數(shù)據(jù)集
#例如，本例中，我們要分析GO，因?yàn)閙ouse文件包含了所有的基因集，所以要查看GO在哪里，然后將需要的文件提出來(lái)。
table(zebrafish$gs_subcat)
zebrafish_GO_bp = msigdbr(species = "Danio rerio",
                      category = "C5", #GO在C5
                      subcategory = "GO:BP") %>% 
  dplyr::select(gs_name,gene_symbol)#這里可以選擇gene symbol，也可以選擇ID，根據(jù)自己數(shù)據(jù)需求來(lái)，主要為了方便
head(zebrafish_GO_bp)

zebrafish_GO_bp_Set = zebrafish_GO_bp %>% split(x = .$gene_symbol, f = .$gs_name)#后續(xù)gsva要求是list，所以將他轉(zhuǎn)化為list
write.csv(zebrafish_GO_bp_Set,"zebrafish_GO_bp_Set.csv")
###names(zebrafish_GO_bp_Set)
##zebrafish_GO_bp_Set

write.gmt <- function(geneSet=zebrafish_GO_bp_Set,gmt_file='kegg2symbol.gmt'){
  sink( gmt_file )
  for (i in 1:length(geneSet)){
    cat(names(geneSet)[i])
    cat('\tNA\t')
    cat(paste(geneSet[[i]],collapse = '\t'))
    cat('\n')
  }
  sink()
}


write.gmt()

構(gòu)好文件如下，第二列為NA

image.png