R语言计算两个基因集重叠显著性及可视化实现方法咨询

（图3a、b，扩展数据图3a、b及补充表1）。处理48小时后，超过1/3的转录组出现差异表达（基因数>5000；其中405个基因编码胞外区蛋白，对应Gene Ontology (GO) 编号0005576），该差异基因集与vemurafenib处理5天后体内A375肿瘤的基因表达变化集存在显著重叠（图3a、b及扩展数据图3c）。在vemurafenib处理的Colo800、UACC62黑色素瘤细胞，以及crizotinib处理的H3122肺腺癌细胞中也观测到了类似的大范围基因表达变化（扩展数据图3d）。尽管细胞谱系、致癌驱动因子、靶向治疗药物均不同，我们仍观测到黑色素瘤与肺腺癌细胞的分泌组存在显著重叠（P < 9.11 × 10−5）

实现方案 #

显著性计算逻辑 #

两个基因集的重叠显著性采用领域通用的超几何分布检验，核心输入为4个参数：

• 总背景基因数N：优先使用本次测序实际检测到的所有基因总数，无对应数据时可默认采用人类全蛋白编码基因数（约20000）
• 基因集1大小k：RNA-seq差异基因总数
• 基因集2大小m：ATAC-seq关联的基因总数
• 交集基因数q：两个基因集的共同基因数量

计算得到的p值含义为：随机从总背景基因中抽取k个基因，其中至少有q个与基因集2重叠的概率，p值越小说明重叠的非随机性越强，显著性越高。

代码修改补充 #

在现有代码的基础上，新增超几何检验逻辑，调整画图标注即可实现论文同款展示效果，修改后的完整可用代码如下：

library(reshape2)
library(venneuler)
RNA_seq_cds <- read.csv("~/Downloads/RNA_seq_gene_set.txt", header=TRUE, sep="\t")
ATAC_seq <- read.csv("~/Downloads/ATAC_seq_gene_set.txt", header=TRUE, sep="\t")
RNA_seq <- RNA_seq_cds

cbindPad <- function(...) {
  args <- list(...)
  n <- sapply(args, nrow)
  mx <- max(n)
  pad <- function(x, mx) {
    if (nrow(x) < mx) {
      nms <- colnames(x)
      padTemp <- matrix(NA, mx - nrow(x), ncol(x))
      colnames(padTemp) <- nms
      if (ncol(x) == 0) {
        return(padTemp)
      } else {
        return(rbind(x, padTemp))
      }
    } else {
      return(x)
    }
  }
  rs <- lapply(args, pad, mx)
  return(do.call(cbind, rs))
}

dat <- cbindPad(ATAC_seq, RNA_seq)

vennfun <- function(x) { 
  x$id <- seq(1, nrow(x))
  xm <- melt(x, id.vars="id", na.rm=TRUE)
  xc <- dcast(xm, value~variable, fun.aggregate=length)
  rownames(xc) <- xc$value
  xc$value <- NULL
  xc
}

VennDat <- vennfun(dat)

# ---------- 新增：超几何检验计算重叠显著性 ----------
# 总背景基因数可根据实际检测的基因总数调整
N <- 20000
k <- sum(VennDat[,1] > 0)
m <- sum(VennDat[,2] > 0)
q <- sum(rowSums(VennDat) == 2)

# 计算p值，lower.tail=FALSE对应P(X≥q)
p_value <- phyper(q-1, m, N - m, k, lower.tail = FALSE)
# 转换为科学计数法便于展示
p_sci <- format(p_value, scientific = TRUE, digits = 2)

# ---------- 调整画图逻辑，新增显著性标注 ----------
vd <- venneuler(VennDat)
vd$labels <- paste(c("RNA-seq", "ATAC-seq"), colSums(VennDat))
plot(vd, cex=1.2)

# 标注交集基因数
text(.3, .5, 
     bquote(bold("Common ="~.(q))), 
     col="red", cex=1.2)

# 标注显著性p值，坐标可根据实际显示效果调整
text(.3, .42,
     bquote(bold("P = "~.(p_sci))),
     col="blue", cex=1.2)

可选优化 #

如果需要更高的美观度和自定义空间，可以换用ggVennDiagram包绘制ggplot2体系的韦恩图，内置支持标注交集数和显著性，调整配色、字体也更便捷。

内容的提问来源于stack exchange，提问作者PesKchan