典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集

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

現在這套代碼綜合高通量數據中批次效應的鑑定和處理 - 系列總結和更新，已經遷移到高顏值免費在線繪圖工具，新增WGCNA和差異分析 。

線性模型

為了分析發炎和未發炎組織的差異表達，我們需要構建一個線性模型。線性模式是實驗數據分析的常用方法，適用於幾乎任何複雜的實驗設計。後面我們專門出文介紹，推薦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 0

head(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_fitted

ggplot(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.6443

overlap_UC

## [1] 0.6731

total_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] 576

total_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_CD

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

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

library(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的成員。

更多基於基因的分析見（還沒完呢，接著往下看）

基因富集分析

基本原理見（依舊沒完呢，接著往下看）：

如前所述，一般推薦使用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 1

table(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)，按圖示操作。

Session information

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /anaconda2/envs/r-environment/lib/R/lib/libRblas.so
##
## locale:
## [1] zh_CN.UTF-8
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] ggrepel_0.8.0 knitr_1.21
## [3] hexbin_1.27.2 genefilter_1.64.0
## [5] matrixStats_0.54.0 stringr_1.4.0
## [7] dplyr_0.8.0.1 pheatmap_1.0.12
## [9] RColorBrewer_1.1-2 geneplotter_1.60.0
## [11] annotate_1.60.0 XML_3.98-1.16
## [13] lattice_0.20-38 ggplot2_3.1.0
## [15] gplots_3.0.1.1 clusterProfiler_3.10.1
## [17] ReactomePA_1.26.0 topGO_2.34.0
## [19] SparseM_1.77 GO.db_3.7.0
## [21] graph_1.60.0 limma_3.38.3
## [23] arrayQualityMetrics_3.38.0 hugene10sttranscriptcluster.db_8.7.0
## [25] org.Hs.eg.db_3.7.0 AnnotationDbi_1.44.0
## [27] pd.hugene.1.0.st.v1_3.14.1 DBI_1.0.0
## [29] oligo_1.46.0 RSQLite_2.1.1
## [31] Biostrings_2.50.2 XVector_0.22.0
## [33] IRanges_2.16.0 S4Vectors_0.20.1
## [35] ArrayExpress_1.42.0 oligoClasses_1.44.0
## [37] Biobase_2.42.0 BiocGenerics_0.28.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 htmlwidgets_1.3
## [3] beadarray_2.32.0 BiocParallel_1.16.5
## [5] munsell_0.5.0 codetools_0.2-16
## [7] preprocessCore_1.44.0 withr_2.1.2
## [9] colorspace_1.4-0 GOSemSim_2.8.0
## [11] highr_0.7 rstudioapi_0.9.0
## [13] setRNG_2013.9-1 DOSE_3.8.0
## [15] labeling_0.3 urltools_1.7.1
## [17] GenomeInfoDbData_1.2.0 hwriter_1.3.2
## [19] polyclip_1.9-1 bit64_0.9-7
## [21] farver_1.1.0 rprojroot_1.3-2
## [23] xfun_0.5 affxparser_1.54.0
## [25] R6_2.4.0 GenomeInfoDb_1.18.1
## [27] illuminaio_0.24.0 gridSVG_1.6-0
## [29] bitops_1.0-6 fgsea_1.8.0
## [31] gridGraphics_0.3-0 DelayedArray_0.8.0
## [33] assertthat_0.2.0 scales_1.0.0
## [35] ggraph_1.0.2 nnet_7.3-12
## [37] enrichplot_1.2.0 gtable_0.2.0
## [39] Cairo_1.5-9 affy_1.60.0
## [41] rlang_0.3.1 splines_3.5.1
## [43] lazyeval_0.2.1 acepack_1.4.1
## [45] europepmc_0.3 checkmate_1.9.1
## [47] BiocManager_1.30.4 yaml_2.2.0
## [49] reshape2_1.4.3 backports_1.1.3
## [51] qvalue_2.14.1 Hmisc_4.2-0
## [53] tools_3.5.1 ggplotify_0.0.3
## [55] affyio_1.52.0 ff_2.2-14
## [57] ggridges_0.5.1 Rcpp_1.0.0
## [59] plyr_1.8.4 base64enc_0.1-3
## [61] progress_1.2.0 zlibbioc_1.28.0
## [63] purrr_0.3.0 RCurl_1.95-4.11
## [65] prettyunits_1.0.2 rpart_4.1-13
## [67] openssl_1.2.1 viridis_0.5.1
## [69] cowplot_0.9.4 SummarizedExperiment_1.12.0
## [71] cluster_2.0.7-1 magrittr_1.5
## [73] data.table_1.12.0 DO.db_2.9
## [75] openxlsx_4.1.0 triebeard_0.3.0
## [77] reactome.db_1.66.0 hms_0.4.2
## [79] evaluate_0.13 xtable_1.8-3
## [81] gcrma_2.54.0 gridExtra_2.3
## [83] compiler_3.5.1 tibble_2.0.1
## [85] KernSmooth_2.23-15 crayon_1.3.4
## [87] htmltools_0.3.6 Formula_1.2-3
## [89] tidyr_0.8.2 tweenr_1.0.1
## [91] MASS_7.3-51.1 rappdirs_0.3.1
## [93] Matrix_1.2-15 vsn_3.50.0
## [95] gdata_2.18.0 igraph_1.2.4
## [97] GenomicRanges_1.34.0 pkgconfig_2.0.2
## [99] rvcheck_0.1.3 foreign_0.8-71
## [101] xml2_1.2.0 foreach_1.4.4
## [103] BeadDataPackR_1.34.0 affyPLM_1.58.0
## [105] digest_0.6.18 rmarkdown_1.10
## [107] base64_2.0 fastmatch_1.1-0
## [109] htmlTable_1.13.1 gtools_3.8.1
## [111] graphite_1.28.2 jsonlite_1.6
## [113] viridisLite_0.3.0 askpass_1.1
## [115] pillar_1.3.1 httr_1.4.0
## [117] survival_2.43-3 glue_1.3.0
## [119] zip_1.0.0 UpSetR_1.3.3
## [121] iterators_1.0.10 bit_1.1-14
## [123] ggforce_0.2.0 stringi_1.3.1
## [125] blob_1.1.1 latticeExtra_0.6-28
## [127] caTools_1.17.1.1 memoise_1.1.0

References

https://f1000research.com/articles/5-1384/v2

學習津貼

單篇留言點讚數的第一位（點讚數至少為8）可獲得我們贈送的在線基礎課的9折優惠券。

越留言，越幸運。

主編會在每周選擇一位最有深度的留言，評論者可獲得我們贈送的任意一門在線課程的9折優惠券（偷偷告訴你，這個任意是由你選擇哦）。

高顏值免費在線繪圖

往期精品

畫圖三字經 生信視頻 生信系列教程 

心得體會 TCGA資料庫 Linux Python 

高通量分析 免費在線畫圖 測序歷史 超級增強子

生信學習視頻 PPT EXCEL 文章寫作 ggplot2

海哥組學 可視化套路 基因組瀏覽器

色彩搭配 圖形排版 互作網絡

自學生信 2019影響因子 GSEA 單細胞 

後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集