VSE—连接遗传位点和功能基因组的桥梁

后GWAS时代人们最关心的问题是如何解释显著遗传位点的生物学功能。一般来讲，由GWAS找到的显著位点大部分落在不编码蛋白的基因组区域，然而我们对这些区域的了解目前还十分有限。随着ENCODE, FANTOM等数据库的逐渐完善，基因组中的非编码区段正在被各种新的组学技术揭开神秘面纱。具有某些特定功能的基因组非编码区段是否与GWAS易感位点在位置上具有显著的交互关系就显得尤为重要。

为了帮助研究以上科学问题，我给各位读者介绍一款不错的工具——Variant Set Enrichment (VSE). 首先，这是一个非常有名的R包，曾在《NatureGenetics》被多次引用。其核心算法是比较真实的遗传位点集合是否比随机匹配的模拟集合更加富集与某一特定的功能基因组区段上。

下面我就带领大家利用示例数据学习使用VSE。

（1）      上传遗传位点数据

首先需要说明的是GWAS报道的显著位点以及与其有连锁不平衡的位点都有可能发挥生物学功能，因此我们需要将tag SNP以及其LD SNP一起纳入分析。VSE支持直接从rAggr网站（http://raggr.usc.edu/）上获取的具有连锁关系的SNP集合。一旦完成以上工作，数据导入就很简单了。示例代码如下：

library("VSE")

bca.ld <-loadLd(file.path(system.file("extdata","ld_BCa_raggr.csv", package = "VSE")), type = "raggr")

bca.avs <- makeAVS(bca.ld)

avs.size <- avsSize(bca.avs)

（2）      构建匹配的随机模拟位点

这一步中，VSE要根据第一步得到的真实的位点信息去基因组中随机产生匹配的模拟位点。只用一行代码就可以轻松搞定。

bca.mrvs.200 <- makeMRVS(bca.avs, bgSize=200, mc.cores = 8)

###该步骤有两个参数，其中bgSzie表示随机模拟的次数，mc.cores表示用到的计算机节点数，Windows下只能是1

（3）      载入功能基因组区段信息

还是经典的bed文件格式即可。VSE自带了5个来自MCF7细胞系的组蛋白Chip-seq的数据，下载地址为（ www.hansenhelab.org/VSE/sample_regions/）。读者也可以通过自带函数“loadSampleRegions”下载。

# Downloading sample regions

sampleSheet_path <- loadSampleRegions()

# Loading sample sheet

samples <- read.csv(sampleSheet_path, header = TRUE)

（4）      用热图直观描绘遗传位点与功能区段的交互情况

在进行统计分析之前，VSE可以给出一个近似热图的交互矩阵（如图1）。不禁让人窃喜：又能在文章中多放一个图了。

bca.intersect <- intersectMatrix(bca.avs, regions = samples, col = c("white", 

    "grey10"), scale = "none", margins = c(5, 5), cexRow = 1, cexCol = 0.5, 

    Rowv = NA, Colv = NA)

图1：遗传位点与基因组区段的交互热图

 

（5）      富集分析

VSE需要SNP集合，模拟集合以及基因组坐标三个输入来进行富集分析。一行代码即可实现：

bca.vse <- variantSetEnrichment(bca.avs, bca.mrvs.200, samples)

在展示P值之前，VSE还提供QQ plot对null distribution的正态性进行检验。只有符合正态分布的结果才可行。

par.original <- par(no.readonly = TRUE)

par(mfrow = c(ceiling(length(samples$Peaks)/3), 3), mai = c(1, 1, 0.5, 0.1))

VSEqq(bca.vse)

par(par.original)

如图2，富集分析的正态性进行检验结果

最后，我们可以放心大胆地展示结果了。VSE提供图和表两种形式的展示结果的方式。

（i）表

bca.vse.res <- VSESummary(bca.vse)

bca.vse.res

结果如下：

 

（ii）图

VSEplot(bca.vse, las = 2, pch = 20, cex = 1, cex.main = 0.6, padj = 0.05, main = "BCa AVS in MCF7 genomic features")

如图3，红色的点对应的区域既是遗传位点显著富集的区域（Bonferroni adjusted P-value < 0.01）

图3：最终富集结果图

 

最后，给各位读者一些使用注意事项

（1）      VSE对tag SNP的数量非常敏感，低于15会导致不准确的结果。

（2）      确保使用0.8作为LD的阈值。

（3）      用户需尽量选择可靠来源的功能基因组数据。

（4）      随机次数决定结果的显著程度。随机次数越大，分布的正态性越好，结果越可信。