典型醫學設計實驗GEO數據分析 (step-by-step) - 數據獲取到標準化介紹了實驗的設計、數據獲取、數據標準化和注釋，下面是如何利用Limma和線性模型鑒定差異基因，并進行GO富集分析。

線性模型

為了分析發炎和未發炎組織的差異表達，我們需要構建一個線性模型。線性模式是實驗數據分析的常用方法，適用于幾乎任何復雜的實驗設計。后面我們專門出文介紹，推薦Mike Love和Michael Irizzary的genomics class和可汗學院的線性代數免費公開課。

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

在構建線性模型時，采用縮寫UC和CD代表疾病類型，non_infl.和infl.代表發炎與否。

# 獲得所有的個體信息 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")

因為要找發炎和未發炎組織的差異表達基因，所以理論上只需要比較這兩個變量。但因為每個獨立的個體有兩套芯片檢測結果 (發炎和未發炎組織)，這是一個配對樣品實驗設計。下游分析時需要考慮個體差異的影響。如果一個實驗特征對結果可能有系統性影響，需要對引入這個特征作為阻斷因子 (bolcking factors)。

為了與文章一致，我們為每個疾病分別構建了一個設計矩陣。(也可以針對完整數據集設計一個聯合模型，但兩種疾病可能特征差別很大，放在一起可能不太合適，從典型醫學設計實驗GEO數據分析 (step-by-step) - 數據獲取到標準化文中的PCA結果也可以看出來)

# 獲得CD疾病的個體 i_CD <- individual[disease == "CD"] # 獲得兩種組織類型和個體的矩陣 # 0 + 表示不設置截距項，所有樣品都有自己的回歸系數 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

檢查下設計矩陣:

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

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

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

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

單個基因差異表達分析測試

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

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

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")

我們擬合線性模型計算回歸系數。

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

設計矩陣 (design_palmieri_CD)與回歸系數向量(crat_coef)相乘獲得擬合后的表達值。

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

設計矩陣中每一行只有值為1的變量用于計算擬合的表達值，crat_fitted的每一項代表每個樣品各個變量回歸系數的和。

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

可視化發炎和未發炎組織擬合后的表達值:

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")

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

這是因為一個病人的兩個樣本的相同blocking variable在設計矩陣中都是激活的，只能在病人內比較。這些blocking variable校正擬合后的組織特異性表達值趨向每個病人的表達值，因此最終的估計是個體內差異的平均值，也就是對應基因的差異倍數 (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基因差異表達檢測

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

$$t=\frac{\bar{d}}{s/\sqrt{n}}$$

個體間表達值的差異的平均值。配對t-test的方差來源于配對樣品。這與標準t-test不同，因此只要配對樣品的表達是相關的 ，配對t檢驗就有更高的統計檢出力 (https://en.wikipedia.org/wiki/Paired_difference_test)。

我們通過blocking variable改善了常規t-test的統計檢出力。下面就對擬合值進行t-test分析，檢測發炎和未發炎的組織中該基因的差異是否顯著不同于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

統計檢測獲得p-value值為 0，因此可以說CRAT基因在炎癥和非炎癥組織中表達差異顯著。

這里算出的p-value跟limma輸出的有些差異，這是因為limma還會做一個方差校正 (variance moderation)。

設定比較分組進行統計檢驗

需要比較發炎和未發炎組織的基因表達差異，使用limma包的makeContrasts函數構建只含有一對比較I-nI的比較矩陣。

對所有數據擬合線性模型，并應用?contrasts.fit()鑒定差異表達基因。

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))

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

提取差異表達基因

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

作為質控，我們繪制了p-value分布的直方圖 (這是結果檢測的很重要標準，后續會專門介紹)。p-value在零假設 (null hypotheses)處應該符合均勻分布，而在0值附近有一個峰，富集差異基因對應的低p-value。

如果某個數據集的p-value分布與這里展示的不同，后續的分析可能會造成誤導 (quality loss)。如果p-value的分布比較發散，可能是有批次效應或有其它blocking factor沒有在設計矩陣中考慮進去。需要嘗試在設計矩陣中增加可能的批次因子或blocking factor再重復執行上述分析。如果仍然沒有解決，應用于多重假設的經驗貝葉斯/ null estimation methods可能會有幫助。(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")

多重假設檢驗校正

在原文中，p-value<0.001是顯著性篩選標準，使用這個標準，UC疾病中獲得947差異表達基因。

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

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

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

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

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

原始p-value 0.001校正后是0.0222, 與我們上面根據p-value估算出的FDR一致。(原文說有一個量級的差異，沒看出來：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.)

這里為了與原文一致，還是繼續使用p-value<0.001作為篩選標準。自己分析時，需要根據adj.P.Val進行篩選。

原文的結果可以從http://links.lww.com/IBD/A795獲得并存儲在當前目錄的palmieri_DE_res.xlsx文件中。原始文章雖然使用p-value進行的篩選，但結果中也包含了FDR值。(生信寶典：這個地方實際是不推薦用Excel查看或存儲結果的，具體見Excel改變了你的基因名，30% 相關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中的數據到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。我們鑒定出更多的差異基因可能是因為考慮到個體作為blocking factor和使用limma做了方差校正。

火山圖可視化差異基因

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

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)

輸出差異分析數據，自己繪制火山圖

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"))

導出結果，用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)

根據log2FolcChange排序基因

#head(res_output) # 整理數據，按差異倍數排序，增加橫軸，選擇log2差異倍數大于5的標記基因名字 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")# 標記基因名字 p + geom_text_repel(aes(label=ID))

我們可以對這些標記的基因做些功能然所，如上調基因S100A8，在genecards.org中搜索后發現是一個pro-inflammatory complex的成員。

更多基于基因的分析見

• 科研小萌新，掌握這些技巧，輕松玩轉各個基因！

• 生信寶典之傻瓜式 (二) 如何快速查找指定基因的調控網絡

• 生信寶典之傻瓜式 (三) 我的基因在哪里發光 - 如何查找基因在發表研究中的表達

• 生信寶典之傻瓜式 (四) 蛋白蛋白互作網絡在線搜索

• 生信寶典之傻瓜式 (五) 文獻挖掘查找指定基因調控網絡

• 生信寶典之傻瓜式 (六) 查找轉錄因子的靶基因

基因富集分析

基本原理見：

• GO、GSEA富集分析一網打進

• GSEA富集分析 - 界面操作

• 無需寫代碼的高顏值富集分析神器

• 去東方，最好用的在線GO富集分析工具

• 沒錢買KEGG怎么辦？REACTOME開源通路更強大

• 超簡便的國產lncRNA預測工具LGC

• 我想做信號通路分析，但我就是不想學編程

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

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

選擇合適的背景基因集

這里我們使用genefilter::genefinder篩選與差異基因表達模式分布相近的一批背景基因集以免因為表達值的篩選對后續的富集分析進行誤導 (差異基因是表達的基因的一部分，嚴格來講用全部注釋基因集做背景注釋集不太妥當)。這個方法于GOrilla的算法類似。

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

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)# 驗證扣除結果，應該為空 intersect(back_genes, DE_genes_CD)## character(0)length(back_genes)## [1] 9756

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

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包進行富集分析，其優勢是會考慮GO的層級結構，只輸出最特異的基因集。

運行topGO

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

gene_IDs <- rownames(table_CD) # 獲取差異基因和與其表達模式相近的背景基因 in_universe <- gene_IDs %in% c(DE_genes_CD, back_genes) in_selection <- gene_IDs %in% DE_genes_CD # 定義一個有名字的列表，同時包含差異基因和背景基因 all_genes <- factor(as.integer(in_selection[in_universe])) names(all_genes) <- gene_IDs[in_universe] # 差異基因為1，背景基因為0 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數據集，注釋數據來源于所用的芯片。?nodeSize參數設置GO term中允許的最小基因數，少于這個數目的注釋不用于后續分析。

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

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

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中的數目 # Expected: 期望的差異基因落在term中的數目 # 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)

富集分析表格結果

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富集分析結果

上一步輸出的文件top10_enriched_go.txt，導入高顏值免費在線繪圖網站 (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>

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。