精品伊人久久大香线蕉,开心久久婷婷综合中文字幕,杏田冲梨,人妻无码aⅴ不卡中文字幕

說到富集，富集是將基因根據(jù)一些先驗的知識（也就是常見的注釋）進行分類的過程。我們一般會想到最常見的是GO/KEGG富集，其思路是先篩選差異基因，然后確定這些差異基因的GO/KEGG注釋，然后通過超幾何分布計算出哪些通路富集到了，通常會選擇一個閾值來卡一下，比如p值和FDR等。因此這會涉及到人為的閾值選擇，具有一定的主觀性，而且只能用于差異較大的基因，所以結果可能有一定的局限性。

根據(jù)上述情況，有了GSEA（Gene Set Enrichment Analysis），其思路是發(fā)表于2005年的Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles，主要是要有兩個概念：預先定義的基因集S（基于先驗知識的基因注釋信息）和待測基因集L（一般是表達矩陣）；然后GSEA目的就是為了判斷S基因集中的基因是隨機分布于L（排序后的數(shù)據(jù)集），還是聚集分布在L的頂部或者底部（這也就是富集）。如果待測基因集中的某些基因顯著富集在L的頂部或者底部，這說明這些基因的表達（因為其是根據(jù)表達譜數(shù)據(jù)）對你定義的分組（預先分組）的差異有顯著影響（一致性），從而找到我們關注的基因集；在富集分析的理論中，GSEA可以認為是第二代，即Functional Class Scoring (FCS) Approaches

GSEA的使用

這里不詳細說算法，具體可看GSEA的文章，因為我也是一知半解。。。

下載地址http://software.broadinstitute.org/gsea/downloads.jsp，PS.會驗證下你的郵箱，先注冊下

第一次使用的話，而且數(shù)據(jù)不大的話，建議使用javaGSEA
Desktop Application，Java的圖形化界面，使用較為友好點，根據(jù)你的數(shù)據(jù)大小，選擇不同內存的版本[1-8G不等]，PS.記得看看系統(tǒng)的Java版本，2G內存開始的GSEA版本需要的是64位的Java 1.8版

打開后的界面如下：

數(shù)據(jù)準備

因為GSEA分析一般只作用于人物種的，所以我準備以TCGA的BRCA的mRNA數(shù)據(jù)作為測試數(shù)據(jù)，正好也試下UCSC xena 瀏覽器才是最簡單的TCGA數(shù)據(jù)下載途徑這個方法來下載TCGA數(shù)據(jù)（數(shù)據(jù)還蠻新的，2017年的）

數(shù)據(jù)下載地址：https://xenabrowser.net/datapages/?dataset=TCGA-BRCA%2FXena_Matrices%2FTCGA-BRCA.htseq_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net

其還提供了ID/Gene Mapping的文件（整理好的），正好可以拿來用，因為雖然GSEA有EnsemblID轉化的chip文件，但是感覺有些數(shù)據(jù)有點問題（可能是由于Ensembl的版本一直在更新的緣故），HUGO gene symbol最好

然后用R處理下，將癌組織和對應的癌旁組織的數(shù)據(jù)分別提取出來分別作為兩組的表達矩陣（gct文件）以及或者分組文件（cls文件）

library(dplyr)library(stringr)data <- read.table(file = "TCGA-BRCA.htseq_fpkm.tsv", sep = "\t", header = T, stringsAsFactors = F, check.names = F)data_df <- data[, -1]normal_df <- data_df[, str_sub(names(data_df), 14, 15) %in% 11:19]length(names(normal_df))tumor_df <- data_df[, grepl(paste(str_sub(names(normal_df), 1, 12), collapse = "|"), names(data_df)) & !names(data_df) %in% names(normal_df)]length(names(tumor_df))idmapping <- read.table(file = "gencode.v22.annotation.gene.probeMap", sep = "\t", header = T, stringsAsFactors = F)geneid <- data.frame(id = data$Ensembl_ID, stringsAsFactors = F)geneid2symbol <- left_join(geneid, idmapping, by = "id")all_df <- cbind(Name = geneid2symbol$gene, DESCRIPTION = "na", tumor_df, normal_df)group <- c(rep("Tumor", 118), rep("Normal", 113))group <- paste(group, collapse = " ")group <- c(paste(c(118+113, 2, 1), collapse = " "), "# Tumor Normal", group)write.table(file = "BRCA_fpkm.txt", all_df, sep = "\t", col.names = T, row.names = F, quote = F)write.table(file = "group.cls", group, col.names = F, row.names = F, quote = F)

從上述代碼，我獲得118個癌組織樣本和對應的113個癌旁樣本的表達譜數(shù)據(jù)，并且將Ensembl ID均轉化為了Gene symbol（避免之后用GSEA時，再用chip做ID轉化）；然后可以直接將txt文件作為輸入，也可以將BRCA_fpkm.txt文件做一些處理變成GSEA標準的gct文件，如下圖所示：分別加上前2行內容，第一行照抄就行，第二行則是geneid數(shù)目和樣本數(shù)目

接著是Phenotype labels文件（上述代碼直接出了），即cls文件，格式如下圖所示：第一行231代表樣本數(shù)目，2代表分2組，空格間隔，1照抄；第二行井號注釋說明分組信息；第三行為每個樣本對應的組名，空格分隔

上述文件的詳細格式可參照網站：http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

如果網絡不佳的話，接下來最好將Gene sets file（也就是GSEA軟件上需要輸入的Gene sets database），作者將gene sets都儲存在Signature Database (MSigDb)中，去官網下載即可http://software.broadinstitute.org/gsea/downloads.jsp，比如下載個c2.cp.kegg.v6.1.symbols.gmt文件作為Gene sets database

如果數(shù)據(jù)是芯片數(shù)據(jù)或者需要GSEA的chip文件做ID轉化的話，則也可以先將chip文件下載下來，F(xiàn)TP地址：ftp://ftp.broadinstitute.org/pub/gsea/annotations

軟件使用

因為是windows桌面式軟件，使用就比較簡單了。首先將BRCA_fpkm.txt和group.cls文件導入，如果網速不好的話（比如我自己），那也把上述下載下來的c2.cp.kegg.v6.1.symbols.gmt文件也導入

接下來點擊RUN GSEA，就是幾個指定參數(shù)的選擇了，如下圖所示:

1. Expression dataset：輸入的表達矩陣
2. Gene sets database：gene sets文件，這里是c2.cp.kegg.v6.1.symbols.gmt文件；PS.這個文件也可以自己創(chuàng)建哦
3. Number of permutations：置換檢驗的次數(shù)
4. Phenotype labels：選擇比較組，如果你輸入的文件就只有2個組別的話，這個就很方便選一個就行了；如果你輸入的有三個組別及以上的話，則這里就要跟你的需要選擇兩個組別的比較組，而且GSEA也會根據(jù)你的組別信息去表達矩陣中提取相對應的數(shù)據(jù)
5. Collapse dataset to gene symbols： 如果你已經ID轉化為HUGO gene symbol，那么這里選FALSE，否則選擇TRUE
6. Permutation type：選擇置換的類型，是random phenotype還是random gene sets，一般每組樣本數(shù)目大于7個時，建議選擇phenotype，否則選擇gene sets（這句話一直沒在官網上找到。。。似乎在文章里？）
7. Chip platform：早些時候主要都是芯片數(shù)據(jù)，那時必須要將芯片id轉化為HUGO gene symbol，所以這個參數(shù)一直叫做chip（我猜的。。。），其實就是對ID進行注釋，即ID轉化，選擇你ID對應的chip文件即可，如果已自行轉化了ID的話，則這里空著就行（記得Collapse dataset to gene symbols選擇否）

除了Required field參數(shù)外，下面還有Basic fields和Advanced fields，具體參見官網吧（注：或者鼠標懸浮在對應參數(shù)名稱上，有簡單的參數(shù)介紹哦）

最后點擊RUN，等待左下角的Running變成Success，然后點擊Success即可查看完整的結果，也可以點擊Show results folder，GSEA將所有結果都放在一個文件夾中了！！！

分析結果

來看下文章里最常見的GSEA的結果圖片，如下圖所示：

從圖上，我們一般關注ES值，峰出現(xiàn)在前端還是后端（ES值大于0在前端，小于0在后端）以及Leading-edge subset（即對富集貢獻最大的部分,領頭亞集）；在ES圖中出現(xiàn)領頭亞集的形狀，表明這個功能基因集在某處理條件下具有更顯著的生物學意義；對于分析結果中，我們一般認為|NES|>1，NOM p-val<0.05，F(xiàn)DR q-val<0.25的通路下的基因集合是有意義的

除了上述的結果外，GSEA還提供了Running the Leading Edge Analysis等操作，也可以看看

GSEA的結果解讀我也不是太熟悉，還是得多看看文獻中的解釋說明啦

多于多個樣本的批處理，GSEA也有服務器版本，通過命令行即可操作，適合批處理操作；其還提供了R腳本可供使用（但官網上說似乎并一定可行，需要自己調整？），反正我也正準備都試試看。。。

參考資料：

功能數(shù)據(jù)庫專題-GSEA
GSEA富集分析 – 界面操作
http://software.broadinstitute.org/gsea/doc/desktop_tutorial.jsp
http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html
GSEA使用（初級）

本文出自于http://www.bioinfo-scrounger.com轉載請注明出處