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

Original 易生信生信宝典 2022-03-28

典型医学设计实验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        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)，按图示操作。

GEO是当今最大、最全的公共基因数据资源库，包括基因的表达、突变、修饰等信息，涵盖几乎所有的疾病，且单个实验检测样品数目较多。TCGA数据库包含11,000个病人的33种肿瘤的7个不同层面的基因数据 (包括基因表达、CNV，SNP，DNA甲基化，miRNA，外显子组等)和临床数据，意在解析癌症发生的分子机制、肿瘤的亚型和治疗靶点等。

这两个来源的数据都是对外公开的，是学习、研究和应用的一个资源宝库。从2006年TCGA计划启动以来，基于TCGA数据发表的文章呈指数增加，一大部分来源于对TCGA数据的再次挖掘。因此学习利用生物信息技术挖掘GEO/TCGA公共数据中疾病的分子特征、合适的检测指标具有重要的临床和科研价值。本课程将从GEO/TCGA的表达、突变数据入手，探索公共数据挖掘的基本套路，分享数据分析和可视化的思路和代码，以便应用于自己的研究。

课程涉及主要内容如下：

每节课1小时一个主题，理论结合实战，学懂原理，实战实操，全是老司机多年经验和代码的无私分享。利用自己电脑，轻松实现整套分析。如果有基础，可以多理解代码内容，做更多定制。如果基础弱一些，只需修改几个备注的变量，即可完成全部分析。下面是课程安排，如11代表第一天第一节课，26代表第二天第六节课，(实际上课会有调整，理论和实战穿插调和)，41为两周后的线上集中视频答疑。

编号	主题	简介
11	GEO挖掘案例介绍和结果解读	学习套路和宏观把控
12	R语言基础知识介绍	基础变量和数据操作
13	ggplot2绘图基础	热图、火山图等常见图绘制
14	GEO数据库使用介绍	数据查找、下载、清洗、批次校正、可视化
15	芯片全套分析	差异基因、GO GSEA(时间序列)富集分析、可视化
16	WGCNA共表达网络	KEGG/Reactome通路可视化
21	实战重现2018 纯公共数据Science文章	文章整体解读和亮点分析
22	实战重现2018 纯公共数据Science文章	表达模式评估，样品差异分析
23	实战重现2018 纯公共数据Science文章	不同来源数据校正和比较的原理和操作
24	实战重现2018 纯公共数据Science文章	WGCNA共表达模块分析、网络可视化、模块性状关联
25	实战重现2018 纯公共数据Science文章	基因表达和突变数据关联
26	常见图快速绘制、解读和排版	见图
31	TCGA数据	案例介绍和结果解读
32	实战重现2018 JAMA文章	TCGA数据下载和整理
33	实战重现2018 JAMA文章	TCGA数据表达分析全套
34	实战重现2018 JAMA文章	突变模式, 突变负荷和生存分析
35	实战重现2018 JAMA文章	突变与临床因素相关性分析以森林图展示
36	考试、圆桌论坛	自评学习效果、知识点回顾
41	答疑-线上	答疑、考试内容串讲

R基础知识和图形绘制

GEO数据分析

GEO数据库使用介绍, 数据查找、下载、清洗、批次校正、可视化
芯片全套分析, 差异基因、GO GSEA(时间序列)富集分析、可视化
WGCNA共表达网络、网络可视化、模块性状关联, KEGG/Reactome通路可视化
实战重现纯公共数据Science文章, 整体解读, 亮点分析, 结果重现

TCGA数据分析

TCGA数据，案例介绍和结果解读
TCGA数据下载和整理，临床数据获取和整理
TCGA数据表达分析全套，差异基因富集分析等
实战重现2018 JAMA文章, 突变模式, 突变负荷和生存分析, 突变与临床因素相关性分析以森林图展示, 新队列数据验证结果
TCGA数据库在线使用
UCSC XENA - 集大成者(TCGA, ICGC)
ICGC数据库使用

希望大家报名时，给出自己想重复的文章或结果，我们综合评估，优先照顾，定制专属课程。如果当时不能满足的，也会在后续讨论群提供解决方案，毕竟我们为所有线下学员提供免费绘图。

(如果基础薄弱，报名付款成功后，可免费领取基础程序课，做好准备工作, 让程序成为我们的得力工具而不是学习新知识的绊脚石。)

往期精彩回顾

主讲教师

陈同，博士，2015毕业于中科院遗传与发育生物学研究所，生物信息专业博士，在Cell Stem Cell(IF=23.2，第一作者兼封面文章)，Nucleic Acids Research，Stem Cells and Development等高水平杂志以第一作者或主要作者发表文章，运营有数万人关注的《生信宝典》微信公众号，给你不一样的学习生信体验。

刘永鑫，博士。2008年毕业于东北农大微生物学专业。2014年中科院遗传发育所获生物信息学博士学位，2016年博士后出站留所工作，任宏基因组学实验室工程师，目前主要研究方向为宏基因组学数据分析与可重复计算。发表论文10余篇，SCI收录7篇。2017年7月创办“宏基因组”公众号，目前分享宏基因组、扩增子原创文章185篇，关注人数2.3万人，累计阅读近300万次。

授课模式

本课程以讲解流程和实际操作为主，采用独创四段式教学，封装好的代码全部分享，随处可用：

第一阶段 3天集中授课；
第二阶段自行练习2周；
第三阶段在线直播答疑；
第四阶段培训视频继续学习；
实现教-练-答-用四个环节的统一协调。

培训时间

2019-04-19 到 2019-04-21 (线下讲解实战)
每天早9点到晚6点，半封闭式教学 (最后1小时为集中讨论时间，最后一天会稍微提前一些，多留出时间讨论，也方便老师乘车返回)
报到时间：培训当天

授课地点

北京市西城区鼓楼附近 (具体开课前一周通知)

课程价格

截止 2019-04-12前 4500 元/人
名额有限，每次课程报名满40人后自动关闭报名通道
提供易汉博基因科技实习机会或工作机会
课程连报有优惠，具体见报名网站

课程福利

座位按报名并缴费或预付款成功顺序从前到后龙摆尾式排序
赠送程序基础课一份 (http://bioinfo.ke.qq.com)
多人 (N，10>N>1) 组团报名并同时缴费，每人还可减免N-1百元 (最高500)
赠送金士顿U盘一个（32G含培训数据和脚本）
附推荐语分享对应的招生信息到朋友圈，截图发到train@ehbio.com 可获得200元生信宝典腾讯课堂课程优惠券（可拆分供多个课程使用）

注意事项 *

需自备笔记本电脑，推荐使用win10系统，4G以上内存 (推荐8G)。课程实践根据需要会提供云计算平台
培训班所有数据，文档为内部资料，仅供参阅，未经允许不得翻印外传登刊
上课期间禁止录音，录像
成功付款的学员，若临时有紧急事情不能到来的，可申请延期，更换后续培训班；也可申请退款
若开课2周 (含) 前申请退款可退还85%费用；开课3个工作日 (含) 前申请退款退还70%的费用 (若已开发票需承担相应手续费)
不可先延期再退款

更多课程的详细介绍，请扫描下方二维码。

复制以下链接http://www.ehbio.com/Training/ 或点击阅读原文跳转报名页，成为实验中不可或缺的人，赶快报名吧！

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

R统计和作图

References

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

“家属和记者取得联系”：记者的退场意味深长

劲爆！为了姜萍两位女CEO互揭老底！

中石化一副总被曝出轨人妻，本人嚣张回应：旧情复燃尔

（待会删）大家低调浏览

又一女明星涉毒被判刑！自称为了“刺激大脑”创作，央视网发文痛批

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

线性模型

单个基因差异表达分析测试

设定比较分组进行统计检验

提取差异表达基因

多重假设检验校正

火山图可视化差异基因

根据log2FolcChange排序基因

基因富集分析

选择合适的背景基因集

运行topGO

可视化GO富集分析结果

R基础知识和图形绘制

GEO数据分析

TCGA数据分析

往期精彩回顾

授课模式

培训时间

授课地点

课程价格

课程福利

Session information

R统计和作图

References

易生信系列培训课程，扫码获取免费资料

更多阅读

您可能也对以下帖子感兴趣

“家属和记者取得联系”：记者的退场意味深长

劲爆！为了姜萍两位女CEO互揭老底！

中石化一副总被曝出轨人妻，本人嚣张回应：旧情复燃尔

（待会删）大家低调浏览

又一女明星涉毒被判刑！自称为了“刺激大脑”创作，央视网发文痛批

生成图片，分享到微信朋友圈

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

线性模型

单个基因差异表达分析测试

设定比较分组进行统计检验

提取差异表达基因

多重假设检验校正

火山图可视化差异基因

根据log2FolcChange排序基因

基因富集分析

选择合适的背景基因集

运行topGO

可视化GO富集分析结果

R基础知识和图形绘制

GEO数据分析

TCGA数据分析

往期精彩回顾

授课模式

培训时间

授课地点

课程价格

课程福利

Session information

R统计和作图

References

易生信系列培训课程，扫码获取免费资料

更多阅读

您可能也对以下帖子感兴趣