https://www.sohu.com/a/281843011_120055884
http://www.reibang.com/p/cff8fe8b548a
之前介紹的一種是用R包做的圖，利用的是log2foldchange 和ensemble ID
這次介紹的是GESA軟件操作
原始分析數(shù)據(jù)
原始的FPKM值為
基因名（需要是ENTREZID格式）
構(gòu)建其他物種
GSEA軟件下載：https://www.gsea-msigdb.org/gsea/index.jsp 直接選擇你操作系統(tǒng)對(duì)應(yīng)的版本截型，下載安裝即可。
準(zhǔn)備上傳文件
1.基因表達(dá)量表（.gct）
通常用.gct為后綴纽窟。文件第一行以“#1.2”開(kāi)頭玷室；文件第二行的第一列為基因個(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分割奠货。
GSEA 分析
導(dǎo)入數(shù)據(jù)
打開(kāi)GSEA軟件，點(diǎn)擊左側(cè)菜單欄的Load data選項(xiàng)座掘，將數(shù)據(jù)按如下圖所示方法導(dǎo)入
設(shè)置參數(shù)
按照下圖標(biāo)記的批注設(shè)置合適的參數(shù)递惋。
Expression dataset（表達(dá)文件）： 選擇上一步上傳的表達(dá)gct文件
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
跳到相應(yīng)的分析結(jié)果
如果是其他物種,http://www.reibang.com/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ù)文件
大家可以看到缤言，這個(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()
另外一種方法订晌，構(gòu)建斑馬魚(yú)物種的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
