當(dāng)前位置：首頁 >

典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集

發(fā)布時(shí)間：2025/3/15 78 豆豆

生活随笔收集整理的這篇文章主要介紹了典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - 數(shù)據(jù)獲取到標(biāo)準(zhǔn)化介紹了實(shí)驗(yàn)的設(shè)計(jì)、數(shù)據(jù)獲取、數(shù)據(jù)標(biāo)準(zhǔn)化和注釋，下面是如何利用Limma和線性模型鑒定差異基因，并進(jìn)行GO富集分析。

線性模型

為了分析發(fā)炎和未發(fā)炎組織的差異表達(dá)，我們需要構(gòu)建一個(gè)線性模型。線性模式是實(shí)驗(yàn)數(shù)據(jù)分析的常用方法，適用于幾乎任何復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)。后面我們專門出文介紹，推薦Mike Love和Michael Irizzary的genomics class和可汗學(xué)院的線性代數(shù)免費(fèi)公開課。

我們這里主要用limma包構(gòu)建線性模型進(jìn)行差異表達(dá)分析。這個(gè)包可以同時(shí)比較很多實(shí)驗(yàn)組并且盡量維持其易用性。首先對每個(gè)基因的表達(dá)擬合一個(gè)線性模型，然后用經(jīng)驗(yàn)貝葉斯 (Empirical Bayes)或其他方法進(jìn)行殘差分析獲得合適的t統(tǒng)計(jì)量，并針對小樣本實(shí)驗(yàn)的方差估計(jì)進(jìn)行優(yōu)化，使得分析結(jié)果更加可靠。

在構(gòu)建線性模型時(shí)，采用縮寫UC和CD代表疾病類型，non_infl.和infl.代表發(fā)炎與否。

# 獲得所有的個(gè)體信息 individual <- as.character(Biobase::pData(palmieri_final)$Characteristics.individual.) # 組織信息的空格替換為下劃線 tissue <- str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.phenotype., " ", "_") # 簡化組織的描述信息tissue <- ifelse(tissue == "non-inflamed_colonic_mucosa", "nI", "I") # 疾病信息替換下劃線，并簡化描述 disease <- str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.disease., " ", "_") disease <- ifelse(str_detect(Biobase::pData(palmieri_final)$Factor.Value.disease., "Crohn"), "CD", "UC")

因?yàn)橐野l(fā)炎和未發(fā)炎組織的差異表達(dá)基因，所以理論上只需要比較這兩個(gè)變量。但因?yàn)槊總€(gè)獨(dú)立的個(gè)體有兩套芯片檢測結(jié)果 (發(fā)炎和未發(fā)炎組織)，這是一個(gè)配對樣品實(shí)驗(yàn)設(shè)計(jì)。下游分析時(shí)需要考慮個(gè)體差異的影響。如果一個(gè)實(shí)驗(yàn)特征對結(jié)果可能有系統(tǒng)性影響，需要對引入這個(gè)特征作為阻斷因子 (bolcking factors)。

為了與文章一致，我們?yōu)槊總€(gè)疾病分別構(gòu)建了一個(gè)設(shè)計(jì)矩陣。(也可以針對完整數(shù)據(jù)集設(shè)計(jì)一個(gè)聯(lián)合模型，但兩種疾病可能特征差別很大，放在一起可能不太合適，從典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - 數(shù)據(jù)獲取到標(biāo)準(zhǔn)化文中的PCA結(jié)果也可以看出來)

# 獲得CD疾病的個(gè)體 i_CD <- individual[disease == "CD"] # 獲得兩種組織類型和個(gè)體的矩陣 # 0 + 表示不設(shè)置截距項(xiàng)，所有樣品都有自己的回歸系數(shù) design_palmieri_CD <- model.matrix(~ 0 + tissue[disease == "CD"] + i_CD) colnames(design_palmieri_CD)[1:2] <- c("I", "nI") rownames(design_palmieri_CD) <- i_CD i_UC <- individual[disease == "UC"] design_palmieri_UC <- model.matrix(~ 0 + tissue[disease == "UC"] + i_UC ) colnames(design_palmieri_UC)[1:2] <- c("I", "nI") rownames(design_palmieri_UC) <- i_UC

檢查下設(shè)計(jì)矩陣:

head(design_palmieri_CD[, 1:6])## I nI i_CD183 i_CD2114 i_CD2209 i_CD2255 ## 164 0 1 0 0 0 0 ## 164 1 0 0 0 0 0 ## 183 0 1 1 0 0 0 ## 183 1 0 1 0 0 0 ## 2114 0 1 0 1 0 0 ## 2114 1 0 0 1 0 0head(design_palmieri_UC[, 1:6])## I nI i_UC2424 i_UC2987 i_UC2992 i_UC2995 ## 2400 0 1 0 0 0 0 ## 2400 1 0 0 0 0 0 ## 2424 0 1 1 0 0 0 ## 2424 1 0 1 0 0 0 ## 2987 0 1 0 1 0 0 ## 2987 1 0 0 1 0 0

設(shè)計(jì)矩陣 (design matrix)中行代表每個(gè)芯片，列代表囊括入線性模型的變量，包含是否發(fā)炎，個(gè)體信息。i_UC2424是病人2424的變量，UC代表病人所患疾病。矩陣中的0和1代表所屬關(guān)系 (也稱激活狀態(tài))。

在線性模型中，第一個(gè)個(gè)體 (設(shè)計(jì)矩陣第一行)會(huì)作為其他個(gè)體的基準(zhǔn)，不會(huì)包含到樣品變量中。這里~0是去除截距項(xiàng)，每個(gè)樣品都計(jì)算回歸系數(shù)。

除了像上面把個(gè)體作為blocking factor的方式，也可以構(gòu)建混合模型 (mixed model)，增加random effect，以后再詳細(xì)講述。

單個(gè)基因差異表達(dá)分析測試

先選擇一個(gè)基因查看其分布和擬合出的線性模型，這里選擇的是PROBEID為8164535，gene symbol為CRAT的基因。

首先看下其表達(dá)，整體在未發(fā)炎樣品中高。而且不同病人間基因的絕對表達(dá)豐度差別挺大。如果我們沒有在設(shè)計(jì)矩陣中考慮到這個(gè)問題，線性模型就會(huì)把這些數(shù)據(jù)視為一個(gè)整體。考慮到每個(gè)個(gè)體的基準(zhǔn)表達(dá)水平不同，最終獲得的差異倍數(shù)會(huì)有較高的方差。

tissue_CD <- tissue[disease == "CD"] crat_expr <- Biobase::exprs(palmieri_final)["8164535", disease == "CD"] crat_data <- as.data.frame(crat_expr) colnames(crat_data)[1] <- "org_value" crat_data <- mutate(crat_data, individual = i_CD, tissue_CD)crat_data$tissue_CD <- factor(crat_data$tissue_CD, levels = c("nI", "I"))ggplot(data = crat_data, aes(x = tissue_CD, y = org_value)) +geom_boxplot(aes(fill=tissue_CD)) +geom_line(aes(group = individual, color = individual)) +ggtitle("Expression changes for the CRAT gene") +theme(legend.position = "none")

我們擬合線性模型計(jì)算回歸系數(shù)。

crat_coef <- lmFit(palmieri_final[,disease == "CD"],design = design_palmieri_CD)$coefficients["8164535",]crat_coef## I nI i_CD183 i_CD2114 i_CD2209 i_CD2255 i_CD255 i_CD2826 ## 6.76669 7.19173 0.12382 -0.22145 0.55759 -0.39905 0.29204 -1.07285 ## i_CD2853 i_CD2978 i_CD321 i_CD3262 i_CD3266 i_CD3271 i_CD3302 i_CD3332 ## -0.78285 -0.11633 0.01692 -0.62480 -0.46209 -0.61732 -0.30257 -0.09709

設(shè)計(jì)矩陣 (design_palmieri_CD)與回歸系數(shù)向量(crat_coef)相乘獲得擬合后的表達(dá)值。

crat_fitted <- design_palmieri_CD %*% crat_coef rownames(crat_fitted) <- names(crat_expr) colnames(crat_fitted) <- "fitted_value"crat_fitted## fitted_value ## 164_I_.CEL 7.192 ## 164_II.CEL 6.767 ## 183_I.CEL 7.316 ## 183_II.CEL 6.891 ## 2114_I.CEL 6.970 ## 2114_II.CEL 6.545 ## 2209_A.CEL 7.749 ## 2209_B.CEL 7.324 ## 2255_I.CEL 6.793 ## 2255_II.CEL 6.368 ## 255_I.CEL 7.484 ## 255_II.CEL 7.059 ## 2826_I.CEL 6.119 ## 2826_II.CEL 5.694 ## 2853_I.CEL 6.409 ## 2853_II.CEL 5.984 ## 2978_I.CEL 7.075 ## 2978_II.CEL 6.650 ## 321_I.CEL 7.209 ## 321_II.CEL 6.784 ## 3262_I.CEL 6.567 ## 3262_II.CEL 6.142 ## 3266_I.CEL 6.730 ## 3266_II.CEL 6.305 ## 3271_I.CEL 6.574 ## 3271_II.CEL 6.149 ## 3302_I.CEL 6.889 ## 3302_II.CEL 6.464 ## 3332_I.CEL 7.095 ## 3332_II.CEL 6.670

設(shè)計(jì)矩陣中每一行只有值為1的變量用于計(jì)算擬合的表達(dá)值，crat_fitted的每一項(xiàng)代表每個(gè)樣品各個(gè)變量回歸系數(shù)的和。

例如對病人樣品?2114_I.CEL?6.9703: 對應(yīng)的激活變量的回歸系數(shù)的和 “nI” (7.1917) and “i_CD2114” (-0.2215)。

可視化發(fā)炎和未發(fā)炎組織擬合后的表達(dá)值:

crat_data$fitted_value <- crat_fittedggplot(data = crat_data, aes(x = tissue_CD, y = fitted_value)) +geom_boxplot(aes(fill=tissue_CD)) +geom_line(aes(group = individual, color = individual)) +ggtitle("Expression changes for the CRAT gene") +theme(legend.position = "none")

可以看到，所有病人中發(fā)炎組織和未發(fā)炎組織的CRAT基因表達(dá)差別一致，并且是由變量I(6.7667) 和?nI?(7.1917)的回歸系數(shù)的差值決定。

這是因?yàn)橐粋€(gè)病人的兩個(gè)樣本的相同blocking variable在設(shè)計(jì)矩陣中都是激活的，只能在病人內(nèi)比較。這些blocking variable校正擬合后的組織特異性表達(dá)值趨向每個(gè)病人的表達(dá)值，因此最終的估計(jì)是個(gè)體內(nèi)差異的平均值，也就是對應(yīng)基因的差異倍數(shù) (log2 fold change)。(These blocking variables correct the fitted tissue specific expression values towards the expression levels of the individual patients. Therefore the final estimate is like an average of all the within-individual differences.)

CRAT基因差異表達(dá)檢測

為了檢測基因是否是差異表達(dá)，需要執(zhí)行零假設(shè) (null hypothesis)為基因在發(fā)炎和未發(fā)炎的樣品中表達(dá)無差異的t-test。我們的這種設(shè)計(jì)概念上類似于配對t-test，統(tǒng)計(jì)量t的計(jì)算如下：

$\frac{\barozvdkddzhkzd }{s/\sqrt{n}}$

個(gè)體間表達(dá)值的差異的平均值。配對t-test的方差來源于配對樣品。這與標(biāo)準(zhǔn)t-test不同，因此只要配對樣品的表達(dá)是相關(guān)的，配對t檢驗(yàn)就有更高的統(tǒng)計(jì)檢出力 (https://en.wikipedia.org/wiki/Paired_difference_test)。

我們通過blocking variable改善了常規(guī)t-test的統(tǒng)計(jì)檢出力。下面就對擬合值進(jìn)行t-test分析，檢測發(fā)炎和未發(fā)炎的組織中該基因的差異是否顯著不同于0。

crat_noninflamed <- na.exclude(crat_data$org_value[tissue == "nI"]) crat_inflamed <- na.exclude(crat_data$org_value[tissue == "I"]) res_t <- t.test(crat_noninflamed, crat_inflamed, paired = TRUE) res_t## ## Paired t-test ## ## data: crat_noninflamed and crat_inflamed ## t = 6.8, df = 14, p-value = 8e-06 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.2919 0.5581 ## sample estimates: ## mean of the differences ## 0.425

統(tǒng)計(jì)檢測獲得p-value值為 0，因此可以說CRAT基因在炎癥和非炎癥組織中表達(dá)差異顯著。

這里算出的p-value跟limma輸出的有些差異，這是因?yàn)閘imma還會(huì)做一個(gè)方差校正 (variance moderation)。

設(shè)定比較分組進(jìn)行統(tǒng)計(jì)檢驗(yàn)

需要比較發(fā)炎和未發(fā)炎組織的基因表達(dá)差異，使用limma包的makeContrasts函數(shù)構(gòu)建只含有一對比較I-nI的比較矩陣。

對所有數(shù)據(jù)擬合線性模型，并應(yīng)用?contrasts.fit()鑒定差異表達(dá)基因。

contrast_matrix_CD <- makeContrasts(I-nI, levels = design_palmieri_CD)palmieri_fit_CD <- eBayes(contrasts.fit(lmFit(palmieri_final[,disease == "CD"],design = design_palmieri_CD),contrast_matrix_CD))contrast_matrix_UC <- makeContrasts(I-nI, levels = design_palmieri_UC)palmieri_fit_UC <- eBayes(contrasts.fit(lmFit(palmieri_final[,disease == "UC"],design = design_palmieri_UC),contrast_matrix_UC))

這里應(yīng)用eBayes()函數(shù)對模型進(jìn)行經(jīng)驗(yàn)貝葉斯方差校正 (empirical Bayes variance moderation)獲得校正后的t-統(tǒng)計(jì)量。這是因?yàn)樾酒瑪?shù)據(jù)中樣本量較少，方差估計(jì)困難。通過組合一個(gè)先驗(yàn)方差 (prior variance)和每個(gè)基因的方差可以改善方差估計(jì)，獲得校正后的方差 (也稱?moderation)。經(jīng)驗(yàn)貝葉斯 (Empirical Bayes)就是從數(shù)據(jù)中估算先驗(yàn)方差。eBayes()步驟獲得的方差是個(gè)體方差向先驗(yàn)值趨近后的結(jié)果 (individual variances are shrunken towards the prior value)。

提取差異表達(dá)基因

topTable函數(shù)可用于提取差異基因。我們獲得克羅恩病和潰瘍性結(jié)腸炎的差異分析結(jié)果，并把基因按照絕對t-統(tǒng)計(jì)量排序。

作為質(zhì)控，我們繪制了p-value分布的直方圖 (這是結(jié)果檢測的很重要標(biāo)準(zhǔn)，后續(xù)會(huì)專門介紹)。p-value在零假設(shè) (null hypotheses)處應(yīng)該符合均勻分布，而在0值附近有一個(gè)峰，富集差異基因?qū)?yīng)的低p-value。

如果某個(gè)數(shù)據(jù)集的p-value分布與這里展示的不同，后續(xù)的分析可能會(huì)造成誤導(dǎo) (quality loss)。如果p-value的分布比較發(fā)散，可能是有批次效應(yīng)或有其它blocking factor沒有在設(shè)計(jì)矩陣中考慮進(jìn)去。需要嘗試在設(shè)計(jì)矩陣中增加可能的批次因子或blocking factor再重復(fù)執(zhí)行上述分析。如果仍然沒有解決，應(yīng)用于多重假設(shè)的經(jīng)驗(yàn)貝葉斯/ null estimation methods可能會(huì)有幫助。(If this does not help, empirical Bayes / null estimation methods for multiple testing are useful.)

table_CD <- topTable(palmieri_fit_CD, number = Inf) head(table_CD)hist(table_CD$P.Value, col = brewer.pal(3, name = "Set2")[1],main = "inflamed vs non-inflamed - Crohn's disease", xlab = "p-values")

table_UC <- topTable(palmieri_fit_UC, number = Inf) head(table_UC)hist(table_UC$P.Value, col = brewer.pal(3, name = "Set2")[2],main = "inflamed vs non-inflamed - Ulcerative colitis", xlab = "p-values")

多重假設(shè)檢驗(yàn)校正

在原文中，p-value<0.001是顯著性篩選標(biāo)準(zhǔn)，使用這個(gè)標(biāo)準(zhǔn)，UC疾病中獲得947差異表達(dá)基因。

nrow(subset(table_UC, P.Value < 0.001))## [1] 947

然而，這個(gè)列表中很難界定哪些是假陽性結(jié)果。采用p-value，我們可以說至多有 21041 (total number of tests) * 0.001 = 21.041 假陽性基因。那么在我們鑒定的差異基因列表中，有最多 2.22% 的基因可能是假陽性。

因此，如果同時(shí)做了特別多比較，采用原始的p-values作為篩選標(biāo)準(zhǔn)就有些寬松了，需要做多重假設(shè)檢驗(yàn)校正。分子生物學(xué)中最流行的校正方式是假陽性率 (false discovery rate, FDR)。采用簡單的p-value閾值獲得的差異基因列表的FDR可能會(huì)比較高。

另一方面，在p-value分布直方圖中有一個(gè)差異基因構(gòu)成的明顯的峰。因此我們期待我們的差異基因列表中FDR比較低。

tail(subset(table_UC, P.Value < 0.001))

原始p-value 0.001校正后是0.0222, 與我們上面根據(jù)p-value估算出的FDR一致。(原文說有一個(gè)量級的差異，沒看出來：The adjusted p-value for a raw p-value of 0.001 in the table is 0.0222, which is an order of magnitude lower than the FDR we can infer from p-values alone.)

這里為了與原文一致，還是繼續(xù)使用p-value<0.001作為篩選標(biāo)準(zhǔn)。自己分析時(shí)，需要根據(jù)adj.P.Val進(jìn)行篩選。

原文的結(jié)果可以從http://links.lww.com/IBD/A795獲得并存儲在當(dāng)前目錄的palmieri_DE_res.xlsx文件中。原始文章雖然使用p-value進(jìn)行的篩選，但結(jié)果中也包含了FDR值。(生信寶典：這個(gè)地方實(shí)際是不推薦用Excel查看或存儲結(jié)果的，具體見Excel改變了你的基因名，30% 相關(guān)Nature文章受影響，NCBI也受波及。)

#fpath <- system.file("extdata", "palmieri_DE_res.xlsx", package = "maEndToEnd") fpath <- "palmieri_DE_res.xlsx" palmieri_DE_res <- sapply(1:4, function(i) openxlsx::read.xlsx(fpath, cols = 1, sheet = i, startRow = 4))names(palmieri_DE_res) <- c("CD_UP", "CD_DOWN", "UC_UP", "UC_DOWN") palmieri_DE_res <- lapply(palmieri_DE_res, as.character) paper_DE_genes_CD <- Reduce("c", palmieri_DE_res[1:2]) paper_DE_genes_UC <- Reduce("c", palmieri_DE_res[3:4])overlap_CD <- length(intersect(subset(table_CD, P.Value < 0.001)$SYMBOL, paper_DE_genes_CD)) / length(paper_DE_genes_CD)overlap_UC <- length(intersect(subset(table_UC, P.Value < 0.001)$SYMBOL,paper_DE_genes_UC)) / length(paper_DE_genes_UC) overlap_CD## [1] 0.6443overlap_UC ## [1] 0.6731total_genenumber_CD <- length(subset(table_CD, P.Value < 0.001)$SYMBOL) total_genenumber_UC <- length(subset(table_UC, P.Value < 0.001)$SYMBOL)total_genenumber_CD## [1] 576total_genenumber_UC## [1] 947

繪制Venn圖，拷貝all_de_for_venn.txt中的數(shù)據(jù)到www.ehbio.com/ImageGP，按圖示選擇，即可獲得Venn圖。

allList <- rbind(cbind(list=paper_DE_genes_CD, type="paper_DE_genes_CD"),cbind(list=paper_DE_genes_UC, type="paper_DE_genes_UC"),cbind(list=subset(table_CD, P.Value < 0.001)$SYMBOL, type="our_DE_genes_CD"),cbind(list=subset(table_UC, P.Value < 0.001)$SYMBOL, type="our_DE_genes_UC")) head(allList)## list type ## [1,] "SEC61B" "paper_DE_genes_CD" ## [2,] "PRKD1" "paper_DE_genes_CD" ## [3,] "AVPR1A" "paper_DE_genes_CD" ## [4,] "VTRNA1-3" "paper_DE_genes_CD" ## [5,] "LGALS1" "paper_DE_genes_CD" ## [6,] "FKBP14" "paper_DE_genes_CD"write.table(allList,"all_de_for_venn.txt",sep="\t", row.names=F, col.names=F, quote=F)

我們在CD中找到576差異基因 (原文是298)，在UC中找到947差異基因 (原文是520)。二者之間的吻合度還是比較好的，CD中共有差異基因比例0.6443，UC中共有差異基因比例0.6731。我們鑒定出更多的差異基因可能是因?yàn)榭紤]到個(gè)體作為blocking factor和使用limma做了方差校正。

火山圖可視化差異基因

為了更好的可視化效果，只在火山圖標(biāo)注差異倍數(shù)大于1的p-value值最小的100個(gè)基因的名字。

volcano_names <- ifelse(abs(palmieri_fit_CD$coefficients)>=1, palmieri_fit_CD$genes$SYMBOL, NA)limma::volcanoplot(palmieri_fit_CD, coef = 1L, style = "p-value", highlight = 100, names = volcano_names,xlab = "Log2 Fold Change", ylab = NULL, pch=16, cex=0.35)

輸出差異分析數(shù)據(jù)，自己繪制火山圖

library(ggrepel) res_output <- table_CDres_output$level <- ifelse(res_output$adj.P.Val<=0.05, ifelse(res_output$logFC>=1, "UP", ifelse(res_output$logFC<=1*(-1), "DW", "NoDiff")) , "NoDiff")res_output <- res_output[order(res_output$P.Value),]top100_p <- res_output$P.Value[100]res_output$ID <- ifelse((res_output$SYMBOL %in% volcano_names) & (res_output$P.Value<top100_p), res_output$SYMBOL,NA) res_output$padj <- (-1) * log10(res_output$adj.P.Val) res_output$padj <- replace(res_output$pad, res_output$pad>5, 5.005)boundary = ceiling(max(abs(res_output$logFC))) p = ggplot(res_output, aes(x=logFC,y=padj,color=level)) +#geom_point(aes(size=baseMean), alpha=0.5) + theme_classic() +geom_point(size=1, alpha=0.5) + theme_classic() +xlab("Log2 transformed fold change") + ylab("Negative Log10 transformed FDR") +xlim(-1 * boundary, boundary) + theme(legend.position="top", legend.title=element_blank()) +geom_text_repel(aes(label=ID)) #ggsave(p, filename=paste0(file_base1,".volcano.pdf"),width=13.5,height=15,units=c("cm"))

導(dǎo)出結(jié)果，用ImageGP (www.ehbio.com/ImageGP)繪圖

colnames(res_output) <- str_replace_all(colnames(res_output), '[ .][ .]*', '_') write.table(res_output, "for_volcano.txt", row.names=F, sep="\t", quote=F)

根據(jù)log2FolcChange排序基因

#head(res_output) # 整理數(shù)據(jù)，按差異倍數(shù)排序，增加橫軸，選擇log2差異倍數(shù)大于5的標(biāo)記基因名字 res_output_line <- res_output[order(res_output$logFC),] res_output_line$x <- 1:nrow(res_output_line) res_output_line$ID <- ifelse((res_output_line$SYMBOL %in% volcano_names) & (res_output_line$P.Value<top100_p), res_output_line$SYMBOL,NA) head(res_output_line)## PROBEID SYMBOL ## 7919055 7919055 HMGCS2 ## 7994252 7994252 AQP8 ## 8063590 8063590 PCK1 ## 8109194 8109194 SLC26A2 ## 8101675 8101675 ABCG2 ## 7962559 7962559 SLC38A4 ## GENENAME logFC ## 7919055 3-hydroxy-3-methylglutaryl-CoA synthase 2 -2.231 ## 7994252 aquaporin 8 -2.184 ## 8063590 phosphoenolpyruvate carboxykinase 1 -1.816 ## 8109194 solute carrier family 26 member 2 -1.736 ## 8101675 ATP binding cassette subfamily G member 2 (Junior blood group) -1.727 ## 7962559 solute carrier family 38 member 4 -1.688 ## AveExpr t P.Value adj.P.Val B level padj x ID ## 7919055 9.103 -3.953 0.0009322 0.03573 -0.5639 DW 1.447 1 NA ## 7994252 9.309 -3.025 0.0072765 0.07721 -2.4367 NoDiff 1.112 2 NA ## 8063590 7.886 -3.618 0.0019661 0.04395 -1.2468 DW 1.357 3 NA ## 8109194 10.058 -3.480 0.0026722 0.04936 -1.5269 DW 1.307 4 NA ## 8101675 7.657 -4.046 0.0007566 0.03430 -0.3727 DW 1.465 5 NA ## 7962559 4.910 -4.570 0.0002370 0.03004 0.6902 DW 1.522 6 NAlibrary(ggrepel) p <- ggplot(res_output_line, aes(x=x, y=logFC)) + geom_point(aes(color=logFC)) + scale_color_gradient2(low="dark green", mid="yellow", high= "dark red", midpoint = 0) + theme_classic() + geom_hline(yintercept = 0, linetype="dotted")# 標(biāo)記基因名字 p + geom_text_repel(aes(label=ID))

我們可以對這些標(biāo)記的基因做些功能然所，如上調(diào)基因S100A8，在genecards.org中搜索后發(fā)現(xiàn)是一個(gè)pro-inflammatory complex的成員。

更多基于基因的分析見

科研小萌新，掌握這些技巧，輕松玩轉(zhuǎn)各個(gè)基因！
生信寶典之傻瓜式 (二) 如何快速查找指定基因的調(diào)控網(wǎng)絡(luò)
生信寶典之傻瓜式 (三) 我的基因在哪里發(fā)光 - 如何查找基因在發(fā)表研究中的表達(dá)
生信寶典之傻瓜式 (四) 蛋白蛋白互作網(wǎng)絡(luò)在線搜索
生信寶典之傻瓜式 (五) 文獻(xiàn)挖掘查找指定基因調(diào)控網(wǎng)絡(luò)
生信寶典之傻瓜式 (六) 查找轉(zhuǎn)錄因子的靶基因

基因富集分析

基本原理見：

GO、GSEA富集分析一網(wǎng)打進(jìn)
GSEA富集分析 - 界面操作
無需寫代碼的高顏值富集分析神器
去東方，最好用的在線GO富集分析工具
沒錢買KEGG怎么辦？REACTOME開源通路更強(qiáng)大
超簡便的國產(chǎn)lncRNA預(yù)測工具LGC
我想做信號通路分析，但我就是不想學(xué)編程

如前所述，一般推薦使用FDR(也稱adj.P.Val)作為篩選差異基因的閾值，可以更好的估計(jì)假陽性率和比較解釋結(jié)果。在后續(xù)富集分析中，我們使用FDR<0.1作為篩選標(biāo)準(zhǔn)。

DE_genes_CD <- subset(table_CD, adj.P.Val < 0.1)$PROBEID

選擇合適的背景基因集

這里我們使用genefilter::genefinder篩選與差異基因表達(dá)模式分布相近的一批背景基因集以免因?yàn)楸磉_(dá)值的篩選對后續(xù)的富集分析進(jìn)行誤導(dǎo) (差異基因是表達(dá)的基因的一部分，嚴(yán)格來講用全部注釋基因集做背景注釋集不太妥當(dāng))。這個(gè)方法于GOrilla的算法類似。

對每個(gè)差異表達(dá)的基因，使用genefinder函數(shù)鑒定與其有相似表達(dá)的基因。genefinder函數(shù)返回一個(gè)列表有兩個(gè)元素：背景基因的索引和背景基因與差異基因表達(dá)分布的距離度量。

back_genes_idx <- genefilter::genefinder(palmieri_final, as.character(DE_genes_CD), method = "manhattan", scale = "none")# 提取背景基因的PROBIDs back_genes_idx <- sapply(back_genes_idx, function(x)x$indices) head(back_genes_idx)## 7928695 8123695 8164535 8009746 7952249 8105348 8018558 8007008 8072876 ## 8162586 7935180 8084589 7982377 8091411 8081890 8154245 8040362 7993126 ## 7982878 8120927 7956009 7907859 7901549 8008263 8138834 8169504 7901140 # 提取背景基因 back_genes <- featureNames(palmieri_final)[back_genes_idx]# 從中扣除差異基因 back_genes <- setdiff(back_genes, DE_genes_CD)# 驗(yàn)證扣除結(jié)果，應(yīng)該為空 intersect(back_genes, DE_genes_CD)## character(0)length(back_genes)## [1] 9756

繪制所有基因、差異基因和背景基因的表達(dá)分布密度曲線，以看其表達(dá)分布模式是否相近，對后續(xù)的富集分析是否有影響。整體分布模式相仿，只是差異基因 (fore,紫色的線)有些向右偏移，說明鑒定出的差異基因整體表達(dá)較高，在背景基因集中難找到類似的分布。

multidensity(list(all = table_CD[,"AveExpr"] ,fore = table_CD[DE_genes_CD , "AveExpr"],back = table_CD[rownames(table_CD) %in% back_genes, "AveExpr"]),col = c("#e46981", "#ae7ee2", "#a7ad4a"),xlab = "mean expression",main = "DE genes for CD-background-matching")

我們采用topGO包進(jìn)行富集分析，其優(yōu)勢是會(huì)考慮GO的層級結(jié)構(gòu)，只輸出最特異的基因集。

運(yùn)行topGO

獲取差異基因和與其表達(dá)模式相近的背景基因, 定義一個(gè)有名字的列表，同時(shí)包含差異基因 (level 1)和背景基因(level 0)作為topGO的一個(gè)基因列表輸入。

gene_IDs <- rownames(table_CD) # 獲取差異基因和與其表達(dá)模式相近的背景基因 in_universe <- gene_IDs %in% c(DE_genes_CD, back_genes) in_selection <- gene_IDs %in% DE_genes_CD # 定義一個(gè)有名字的列表，同時(shí)包含差異基因和背景基因 all_genes <- factor(as.integer(in_selection[in_universe])) names(all_genes) <- gene_IDs[in_universe] # 差異基因?yàn)?，背景基因?yàn)? head(all_genes) # 7928695 8123695 8164535 8009746 7952249 8105348 # 1 1 1 1 1 1 # Levels: 0 1table(all_genes) # all_genes # 0 1 # 9756 2490

初始化topGO數(shù)據(jù)集，注釋數(shù)據(jù)來源于所用的芯片。?nodeSize參數(shù)設(shè)置GO term中允許的最小基因數(shù)，少于這個(gè)數(shù)目的注釋不用于后續(xù)分析。

top_GO_data <- new("topGOdata", ontology = "BP", allGenes = all_genes, nodeSize = 10, annot = annFUN.db, affyLib = "hugene10sttranscriptcluster.db")

這里應(yīng)用?topGO的兩個(gè)算法進(jìn)行富集測試，常規(guī)的Fisher檢驗(yàn)和elim算法。elim算法考慮GO的層級結(jié)構(gòu)致力于富集最特異的條目。這是一個(gè)自底向上的富集方式，先對最特異的通路進(jìn)行富集分析，如果顯著，分析其上一層注釋時(shí)去除這部分基因，以此類推。

result_top_GO_elim <- runTest(top_GO_data, algorithm = "elim", statistic = "Fisher") result_top_GO_classic <- runTest(top_GO_data, algorithm = "classic", statistic = "Fisher")

篩選top 100富集的通路，顯著富集的通路?raw p-value == 2，不顯著富集的通路?raw p-value == 1。

res_top_GO <- GenTable(top_GO_data, Fisher.elim = result_top_GO_elim,Fisher.classic = result_top_GO_classic,orderBy = "Fisher.elim" , topNodes = 100)genes_top_GO <- printGenes(top_GO_data, whichTerms = res_top_GO$GO.ID,chip = "hugene10sttranscriptcluster.db", geneCutOff = 1000)res_top_GO$sig_genes <- sapply(genes_top_GO, function(x){str_c(paste0(x[x$'raw p-value' == 2, "Symbol.id"],";"), collapse = "")})# Annotated: 背景基因集中term總注釋基因 # Significant：差異基因落在term中的數(shù)目 # Expected: 期望的差異基因落在term中的數(shù)目 # Rank in Fisher.classic：按Fisher.classic pvalue排序 # Fisher.elim Fisher.classic：富集p-value (未做multiple test adjust) head(res_top_GO)# 替換下列名字，不含空格等特殊字符colnames(res_top_GO) <- str_replace_all(colnames(res_top_GO), '[ .][ .]*', '_')write.table(res_top_GO[1:10,], "top10_enriched_go.txt", row.names=F, sep="\t", quote=F)

富集分析表格結(jié)果

GO_ID Term Annotated Significant Expected Rank_in_Fisher_classic Fisher_elim Fisher_classic sig_genes GO:0030593 neutrophil chemotaxis 61 39 13.4 68 1.9e-10 2.0e-12 PIK3CD;PDE4B;S100A9;FCER1G;TGFB2;CSF3R;S100A8;NCKAP1L;C3AR1;LGALS3;DAPK2;CCL22;CX3CL1;CCL2;CCL18;CCL3;CCR7;VAV1;C5AR1;DYSF;IL1RN;CXCR2;IL1B;CXCR1;CXADR;ITGB2;RAC2;CXCL8;CXCL6;CXCL1;CXCL13;CXCL5;CXCL3;CXCL2;CXCL9;CXCL10;CXCL11;RIPOR2;PIK3CG; GO:0051897 positive regulation of protein kinase B ... 107 53 23.51 104 2.6e-10 2.6e-10 PIK3CD;F3;CHI3L1;TCF7L2;PIK3AP1;FGFR2;ERBB3;P2RX4;KITLG;FGF7;CD19;CX3CL1;ERBB2;CSF3;MIR21;PIK3R5;CCL3;CCR7;VAV1;MYDGF;INSR;ITGB1BP1;IGFBP5;IRS1;ITSN1;OSM;CD86;CX3CR1;CD80;PIK3CB;PIK3R1;SEMA5A;PDGFRB;GCNT2;TNF;RAMP3;EGFR;PIK3CG;HGF;NRG1;BAG4;FGFR1;ARFGEF1;ANGPT1;PTK2;TEK;TGFBR1;MYORG;ENG;MIR221;TNF;TNF;GPX1;

可視化GO富集分析結(jié)果

上一步輸出的文件top10_enriched_go.txt，導(dǎo)入高顏值免費(fèi)在線繪圖網(wǎng)站 (www.ehbio.com/ImageGP<www.ehbio.com imagegp="" style=“margin: 0px; padding: 0px; max-width: 100%; box-sizing: border-box !important; overflow-wrap: break-word !important;”>)，按圖示操作。</www.ehbio.com>

總結(jié)

以上是生活随笔為你收集整理的典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

上一篇：哈佛大学教授刘小乐：我与生物信息学的不解
下一篇：别再吼孩子骂孩子了，他们的脑子真的会受伤