GENIE3——一款预测基因调控网络的R包

在网上随意浏览的时候，无意间发现另一种可以构建调控网络的软件，叫做GENIE3。这是一个R包，直接在R里运行。这款R包基于你的表达矩阵来推测基因之间的调控关系。

虽然共表达网络分析（WGCNA）可以帮助我们鉴定在疾病或生物功能中起重要作用的基因，但是从共表达网络中推论因果关系仍然很难。GENIE3通过确定最能解释每个靶基因表达的转录因子表达模式，整合转录因子信息来构建调控网络。但这个软件也有缺点，A limitation of GENIE3 is that TF information is required for it to perform better than random chance。（摘自：差异共表达网络-Co-expression networks）

GENIE3官网说明书：
https://bioconductor.riken.jp/packages/release/bioc/vignettes/GENIE3/inst/doc/GENIE3.html

NOTE1:
GENIE3的input表达矩阵：矩阵的每一列代表一个细胞的不同基因，每一行代表一个基因在不同细胞里的表达量。矩阵的元素可以是count，也可以是其它指标，例如TPM (Transcripts Per Million)、FPKM/RPKM等等。输入矩阵应避免进行标准化或归一化处理，这样会人为的引入多余的协方差。（生信分析工具-SCENIC）
NOTE2:
表达矩阵不能是一个data.frame。必须是一个matrix。如果你用class()检查你的表达矩阵，显示是data.frame，你必须用as.matrix()把你的矩阵转换一下，否则会报错。（Question: genie3 package Error in (function (classes, fdef, mtable)）。

以下内容仅对官网说明书进行翻译

准备表达矩阵

如果你有自己的表达矩阵，可以直接使用。作者这里只是单纯的随机生成了一个表达矩阵用来演示而已。

> exprMatr <- matrix(sample(1:10, 100, replace=TRUE), nrow=20)
> rownames(exprMatr) <- paste("Gene", 1:20, sep="")
> colnames(exprMatr) <- paste("Sample", 1:5, sep="")
> head(exprMatr)

##       Sample1 Sample2 Sample3 Sample4 Sample5
## Gene1       4       1       8      10       5
## Gene2       2       4       7       6      10
## Gene3       1       4       8       3      10
## Gene4       2       8       1      10       4
## Gene5      10       3       8       9       4
## Gene6      10       8       1       4       7

#run GENIE3
##Run GENIE3 with the default parameters
> library(GENIE3)
> set.seed(123)# For reproducibility of results
> weightMat <- GENIE3(as.matrix(exprMatr_f2),nCores=4,nTrees=50) #nTrees: the default value is 1000
> dim(weightMat)
[1] 12279 12279

运行GENIE3

（１）用默认参数运行GENIE3

> library(GENIE3)
> set.seed(123) # 设置seed，以便重复你的结果
> weightMat <- GENIE3(exprMatr)
> dim(weightMat)
## [1] 20 20
> weightMat[1:5,1:5]
##              Gene1     Gene10      Gene11     Gene12      Gene13
## Gene1  0.000000000 0.12287968 0.066797270 0.11860905 0.009444409
## Gene10 0.123264175 0.00000000 0.035959286 0.08078616 0.026713401
## Gene11 0.118580165 0.10273149 0.000000000 0.01047697 0.070549368
## Gene12 0.070181914 0.03900421 0.008827685 0.00000000 0.018925424
## Gene13 0.007990827 0.01520731 0.015534427 0.01199503 0.000000000

这一步输出一个包含假设的权重links，weight值越高，调节关系的可能性就越大。weightMat[i,j]是指从第i个基因到第j个基因的link的weight。

（２）用一个小基因集来进行GENIE3分析

上面用的是所有的基因来进行GENIE3，但是当你真正分析数据的时候，会有几万个基因，即使你进行了过滤，仍然会有上万个基因。这就使得你的GENIE3运行的非常的慢，生成的link数也会几十万或者上百万。你很难从中提取有效的信息。这时，你可以设置一组基因，可以是你感兴趣的，也可以是差异基因，或者其他方法得到的一个基因list。来作为GENIE3的分析目标。

# Genes that are used as candidate regulators
> regulators <- c(2, 4, 7) #这里假设作者用了３个基因来进行后续分析
# Or alternatively:
> regulators <- c("Gene2", "Gene4", "Gene7")
> weightMat <- GENIE3(exprMatr, regulators=regulators)

在何种情况下，在你的表达矩阵里，除了这３个基因，其他基因的weight值都是０．

（３）更改基于Tree的方法，以及其他参数

GENIE3是基于回归树方法进行计算的。这些“Tree”可以利用随机森林法（Random Forest method）或额外树法（Extra-Trees method）来学习。你可以使用参数设置来指定用哪一种方法。(tree.method="RF" 是Random Forests, 也是默认参数； tree.method="ET" 是 Extra-Trees)。

这里我查了一下这两种算法哪一种更好，查到很多文章，但是没一篇能看懂。。。唯一看懂的一句话就是：在高维数据集中时，ET会更差一些。（摘自：随机森林和极端随机树之间的区别）这里还需要再了解一下。

这两种算法都各自有两个参数：K 和ntrees。K是在每个树节点上(node)随机选择的用于最佳分割确定的候选regulators的数量。如果说假设p是上面你的候选基因的数量，那么k值应该是下面这３种情况中的一种：

参数ntrees指定每个集成生成的树的数量。它可以设置为任何正整数(默认值是1000)。

# Use Extra-Trees (ET) method
# 7 randomly chosen candidate regulators at each node of a tree
# 5 trees per ensemble
> weightMat <- GENIE3(exprMatr, treeMethod="ET", K=7, nTrees=50)

（４）平行运行GENIE3

为了减少运行时间，GENIE3可以在多个核上运行。参数ncores指定你想要使用的核的数目。

> set.seed(123) # For reproducibility of results
> weightMat <- GENIE3(exprMatr, nCores=4, verbose=TRUE)

请注意！设置set.seed可以让你在不同次运行中得到相同的结果，但你要保证你的nCore数也是一样的才行。例如，使用set.seed(123)和nCores=1的运行和使用set.seed(123)但nCores=4可能会有不同的结果。

获取调控网络

（１）获取全部调控网络

> linkList <- getLinkList(weightMat)
> dim(linkList)
## [1] 57  3
> head(linkList)
##   regulatoryGene targetGene    weight
## 1          Gene7      Gene4 0.6824708
## 2          Gene2      Gene3 0.6422339
## 3          Gene4     Gene14 0.6307626
## 4          Gene4     Gene16 0.6239141
## 5          Gene7      Gene2 0.6222118
## 6          Gene4     Gene11 0.5670418

这一步会生成一个linkList表，包含每一个link的weight值。每一行对应一个调控关系。第一列是调控因子，第二列是目标基因，最后一列是该link的weight值。

实际上这个表很像WGCNA分析输出的node/edge表，你可以把它放在cytoscape或者VisANT里进行调控网络的可视化。

注意，这一步每次运行获得的link会略有不同。这是由于随机森林和额外树方法的内在随机性导致的。你可以通过增加nTree的值来减少这种variance。

（２）获取top link

上面得到的是所有候选基因的可能调控网络，是一个非常大的表。如果你的基因list有几百个基因，你会得到十几万或者几十万的link，很难从中获得有效信息。这时你可以选择只获取top link，即最有可能的调控关系网络。

> linkList <- getLinkList(weightMat, reportMax=5)

（３）根据weight值筛选link

> linkList <- getLinkList(weightMat, threshold=0.1)

对于weight值的重要说明！这里的weight没有任何统计学上的意义，只提供了可能存在的调控关系。所有没有标准的阈值来筛选，完全根据你的数据情况来进行判断。

写在最后

看了一遍下来，感觉这个软件的随机性很大。所以除了用软件进行预测以外，更重要的还是要用实验证明你的预测和推断。