以下是我们想要绘制的图片

图片.png

目前发现了两种protocol，BWERF和GENIE3，都是基于R进行运算的。
先开始说BWERF吧！
1、由于BWERF需要大量的计算，建议在多核服务器环境下运行BWERF。
2、打开文件 "use_bwerf.R" 修改参数.
3、运行命令

>  nohup Rscript use_bwerf.R & 

附件(复制r脚本即可)

1、两个示例输入文件

pathway_data.csv & TF_data.csv

pathway_data.csv

通路基因表达矩阵，该文件的行列分别表示通路基因名称和样品。

图片.png

TF_data.csv

TFs表达矩阵，该文件的行列分别表示转录因子名称和样品。

图片.png

2、use_bwerf.R

###################parameters need to change################
#number of cores
ncore = 14

#pathway gene file name, the file must be .csv file
pathway_file = "./example_data/pathway_data.csv"

#TF file name, the file must be .csv file
tf_file = "./example_data/TF_data.csv"

#elimination ratio
ratio = 1 / 10

#the format of gene expression data, either "rowissample" or "colissample"
format = "colissample"

#output file name
out_file = "./output/bwerf_output.csv"

#the manner to cut-off TFs, either "assign" or "auto" 
mode = "assign"

#the number of layers of ML-GRN
num.layers = 3

#the numbers of TFs of eachlayer
nums.tf.eachlayer = c(20,10,5)
##############################################################


source("bwerf.R")


if (format == "rowissample") {
  pathway.data = read.csv(pathway_file,header = TRUE)
  tf.data = read.csv(tf_file,header = TRUE)
}else{
  pathway.data = read.csv(
    pathway_file,stringsAsFactors = FALSE,header = FALSE,row.names = 1
  )
  tf.data = read.csv(
    tf_file,stringsAsFactors = FALSE,header = FALSE,row.names = 1
  )
  pathway.data = t(as.matrix(pathway.data))
  tf.data = t(as.matrix(tf.data))
}
ans = multiLayer(
  tf.data, pathway.data, num.layers = num.layers, ratio = ratio, mode = mode, nums.tf.eachlayer =
    nums.tf.eachlayer, verbose = TRUE
)
write.csv(ans,file = out_file,row.names = FALSE)

3、bwerf.R

# add comments on version BWERF_2017_01_10_version2

#install.packages("randomForest","foreach","doMC","mixtools")
#
library(randomForest)
library(foreach)
library(doMC)
library(mixtools)
library(reshape2)

registerDoMC(cores=ncore)

#function:
#           Use backward elimination method to compute the importance of predictors.
#Parameters:
#           x: the predictor matrix (TF genes), each column contains the expression value of a gene
#           y: response varible vector (a pathway gene)
#           ratio: the ratio of TF genes being elimated in each elimation
#     repeats: As there are some randomness in the result of randomForest, we run randomForest "repeats" times.
#           verbose: if verbose = TRUE, the progress of elimation will be shown
#Output:
#           a named vector that contains the importance of the TFs
rf.be=function(x,y,ratio=1/10,repeats=30,verbose=FALSE){
  n.TF = ncol(x)
  IMP = numeric(n.TF)
  names(IMP) = colnames(x)
  while (ncol(x) > 10) {
    if (verbose) {
      print(paste(ncol(x)," TFs left for training random forest model"))
    }
    imp = rep(0,ncol(x))
    names(imp) = colnames(x)
    for (i in 1:repeats) {
      fit = randomForest(
        x = x,y = y,ntree = 1000,mtry = floor(sqrt(ncol(x))),importance = TRUE
      )
      imp = imp + importance(fit)[,1]
    }
    imp = imp / repeats
    imp = sort(imp)
    num.to.el = max(floor(ncol(x) * ratio),1)
    IMP[names(imp)] = imp
    x = x[,colnames(x) %in% names(imp)[(num.to.el + 1):length(imp)]]
  }
  return(IMP)
}
#function:
#           For each pathway gene, use rf.be function to compute importance of TFs to the pathway gene.
#         use multiple cores parallelization to accelerate.
#Parameters:
#           tf.data: each columne contains the expression values of a TF gene
#           pathway.data: each columne contains the expression values of a pathway gene
#           ratio: the ratio parameter used by rf.be function
#           verbose: if verbose=TRUE, the progress will be shown
#Output:    A matrix each row corresponds to a TF, each columne corresponds to a pathway gene,
#           the value in the matrix is the important of the TF to the pathway gene.  
im.mat=function(tf.data,pathway.data, ratio=1/10,verbose=FALSE){
  num.samples <- nrow(pathway.data)
  num.pathway <- ncol(pathway.data)
  num.tf <- ncol(tf.data)
  im.list =
    foreach(i = 1:num.pathway) %dopar% {
      if (verbose) {
        print(paste("pathway ",i,"/",num.pathway,"is processed",sep = ""))
      }
      im = rf.be(tf.data,pathway.data[,i], ratio = ratio,repeats=30, verbose = verbose)
      list(im = im,index.pathway = i)
    }
  im.res = matrix(0,num.tf,num.pathway)
  rownames(im.res) = colnames(tf.data)
  colnames(im.res) = colnames(pathway.data)
  for (i in 1:num.pathway) {
    j = im.list[[i]][[2]]
    a = im.list[[i]][[1]]
    im.res[names(a),j] = a
  }
  return(im.res)
}

#function:
#           Build one layer of gene regulatory network.
#Parameters:
#           tf.data: each columne contains the expression values of a TF gene
#           pathway.data: each columne contains the expression values of a pathway gene
#           ratio: the ratio parameter used by rf.be function
#     mode: either "assign" , "auto" or "noCut"
#           num.alayer: the number of TFs user want to keep.
#Output: A list that contains one layer of GRN   
oneLayer=function(tf.data,pathway.data,num.alayer=20,ratio=1/10,mode="assign",verbose=TRUE){
  im.matrix = im.mat(tf.data,pathway.data,ratio = ratio,verbose = verbose)
  num.tf = dim(tf.data)[2]
  num.pathway = dim(pathway.data)[2]
  im.tf = apply(im.matrix,1,sum)
  im.tf = im.tf[order(im.tf,decreasing = TRUE)]
  im.tf = im.tf[im.tf > 0]
  if (mode == "assign") {
    num.cut = num.alayer
  }else{
    mixmdl = normalmixEM(im.tf,k = 3)
    k = which.max(mixmdl$mu)
    n = 1
    while (which.max(mixmdl$posterior[n,]) == k) {
      n = n + 1
    }
    n = n - 1
    if (n > num.alayer * 1.5) {
      mixmdl = normalmixEM(im.tf[1:n],k = 3)
      k = which.min(mixmdl$mu)
      n = 1
      while (which.max(mixmdl$posterior[n,]) != k) {
        n = n + 1
      }
      n = n - 1
    }
    num.cut = n
  }
  tf.selected = names(im.tf)[1:num.cut]
  sorted.indices <- order(im.matrix, decreasing = TRUE)
  edge.list = data.frame(from = character(),to = character(),im = numeric())
  indx = 0
  num.tf.selected = 0
  set.tf.selected = character()
  while (num.tf.selected < num.cut) {
    indx = indx + 1
    i = sorted.indices[indx]
    im = im.matrix[i]
    row.col <- lin.to.square(i, num.tf)
    row <- row.col[1]
    col <- row.col[2]
    if (colnames(tf.data)[row] %in% tf.selected) {
      edge.list = rbind(edge.list,data.frame(
        from = colnames(tf.data)[row],to = colnames(pathway.data)[col],im = im
      ))
      if (!(colnames(tf.data)[row] %in% set.tf.selected)) {
        set.tf.selected = c(set.tf.selected, colnames(tf.data)[row])
        num.tf.selected = num.tf.selected + 1
      }
    }
  }
  edge.list[,1] = as.character(edge.list[,1])
  edge.list[,2] = as.character(edge.list[,2])
  return(edge.list)
}

#function:
#     for a matrix, change the linear order to square order.
#parameters:
#     i:  the linear order
#     nrow: the number of rows of the matrix.
#Output:
#     a vector contains two elements, containing the row order and the column order. 
lin.to.square <- function(i, nrow) {
  col <- ((i - 1) %/% nrow) + 1
  row <- ((i - 1) %% nrow) + 1
  return(c(row, col))
}


#function:
#     Build a multilayered gene regulatory network.
#parameters:
#           tf.data: each columne contains the expression values of a TF gene
#           pathway.data: each columne contains the expression values of a pathway gene
#           ratio: the ratio parameter used by rf.be function
#     mode: either "assign" or "auto"
#           nums.tf.eachlayer: the numbers of TFs in each layer.
#Output:
#     A list that contains ML-GRN    
multiLayer=function(tf.data,pathway.data,num.layers=3,ratio=1/10,mode="assign",nums.tf.eachlayer=c(20,10,5),verbose=FALSE){
  result = data.frame(
    layer = numeric(),from = character(),to = character(),im = numeric()
  )
  for (i in 1:num.layers) {
    if (verbose) {
      print(paste("current layer: ",i))
      print("---------------------------------")
    }
    alayer = oneLayer(
      tf.data,pathway.data,num.alayer = nums.tf.eachlayer[i],mode = mode,ratio =
        ratio, verbose = verbose
    )
    select.tf = unique(alayer[,1])
    n = nrow(alayer)
    result = rbind(result,data.frame(
      layer = rep(i,n),from = alayer[,1],to = alayer[,2],im = alayer[,3]
    ))
    pathway.data = tf.data[,colnames(tf.data) %in% select.tf]
    tf.data = tf.data[,!(colnames(tf.data) %in% select.tf)]
  }
  result
}

toys_data.R

###################parameters need to change################
library(randomForest)
library(reshape2)
#sample size
n=100

#elimination ratio
ratio = 1 / 10

#number of noise variables
sp=1000

set.seed(100)
###########################################################

##########################################################

set.seed(100)
x1=rnorm(n,1,1)
x2=rnorm(n,1,1)
x3=rnorm(n,1,1)

x4=rnorm(n,3,1)
x5=rnorm(n,3,1)
x6=rnorm(n,3,1)

y=1*x1+1*x2+1*x3+1*x4+1*x5+1*x6+rnorm(n,0,0.1)
  
  
x.noise=matrix(rnorm(sp*n,0,1),n,sp)

x=cbind(x1,x2,x3,x4,x5,x6,x.noise)
colnames(x)=paste("x",1:ncol(x),sep="")


#############Don't use backward elimination as GENIE3#####
par(mfrow=c(1,2))
a=proc.time()
res2=numeric()
for(i in 1:30){
  fit=randomForest(x=x,y=y)
  im=importance(fit,scale=TRUE)
  res2=cbind(res2,im[,1])
}
res2=t(res2)
m2=apply(res2, 2, mean)
m2=sort(m2,decreasing = TRUE)

names2=names(m2)[1:15]
ans2=melt(res2[,colnames(res2) %in% names2])
ans2$Var2=factor(ans2$Var2,levels = names2)

col=rep("blue",15)
indx=which(names2 %in% c("x1","x2","x3","x4","x5","x6"))
col[indx]="red"
boxplot(value~Var2,data=ans2,xlab="",ylab="importance score",
        col=col,xaxt="n",main="GENIE3")
axis(1, at=seq(1,15,by=1), labels=FALSE)
text(x=seq(1,15,by=1), y=par()$usr[3]-0.05*(par()$usr[4]-par()$usr[3]),
     labels=names2, srt=45, adj=1, xpd=TRUE,cex=0.6)
b=proc.time()
b-a
###########################################################
################Backward elimination#####################
a=proc.time()
while(ncol(x)>15){
  res=numeric()
  for(i in 1:30){
    fit=randomForest(x=x,y=y)
    im=importance(fit,scale=TRUE)
    res=cbind(res,im[,1])
  }
  res=t(res)
  m=apply(res, 2, mean)
  m=sort(m)
  n.to.el=floor(ncol(x)/10)
  x=x[,colnames(x) %in% names(m)[(n.to.el+1):ncol(x)]]
}
m=sort(m,decreasing = TRUE)
names1=names(m)[1:15]
ans1=melt(res[,colnames(res) %in% names1])
ans1$Var2=factor(ans1$Var2,levels = names1)

col=rep("blue",15)
indx=which(names1 %in% c("x1","x2","x3","x4","x5","x6"))
col[indx]="red"
boxplot(value~Var2,data=ans1,xlab="",ylab="importance score",
        col=col, xaxt="n",main="BWERF")
axis(1, at=seq(1,15,by=1), labels=FALSE)
text(x=seq(1,15,by=1), y=par()$usr[3]-0.05*(par()$usr[4]-par()$usr[3]),
     labels=names1, srt=45, adj=1, xpd=TRUE,cex=0.6)
b=proc.time()
b-a

再说GENIE3~
我直接搬运一个作者的博客，我怕他会删掉。
正文开始~
在网上随意浏览的时候,无意间发现另一种可以构建调控网络的软件,叫做GENIE3。这是一个R包,直接在R里运行。这款R包基于你的表达矩阵来推测基因之间的调控关系。虽然共表达网络分析(WGCNA)可以帮助我们鉴定在疾病或生物功能中起重要作用的基因,但是从共表达网络中推论因果关系仍然很难。GENIE3通过确定最能解释每个靶基因表达的转录因子表达模式,整合转录因子信息来构建调控网络。但这个软件也有缺点, A limitation ofGENIE3 is that TF information is required for it to perform better than random chance.(摘自:https://www.jianshu.com/p/a324ab116bcd差异共表达网络-Co-expression networks)
GENIE3官网说明书:https://bioconductor.riken.jp/packages/release/bioc/vignettes/GENIE3/inst/doc/GENIE3.html
NOTE1:GENIE3的input表达矩阵:矩阵的每一列代表一个细胞的不同基因,每一行代表一个基因在不同细胞里的表达量。矩阵的元素可以是count,也可以是其它指标,例如TPM (Transcripts PerMillion)、FPKM/RPKM等等。输入矩阵应避免进行标准化或归一化处理,这样会人为的引入多余的协方差。(https://links.jianshu.com/go?to=%255Bhttp%3A%2F%2Fchenyin.top%2Fbioinfo%2F20190315-94a8.html%255D生信分析工具-SCENIC)
NOTE2:表达矩阵不能是一个data.frame。必须是一个matrix。如果你用class()检查你的表达矩阵,显示是data.frame,你必须用as.matrix()把你的矩阵转换一下,否则会报错。(https://links.jianshu.com/go?to=https%3A%2F%2Fwww.biostars.org%2Fp%2F362357%2FQuestion: genie3package 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值应该是下面这３种情况中的一种：

[图片上传失败...(image-231493-1669645500383)]

参数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没有任何统计学上的意义，只提供了可能存在的调控关系。所有没有标准的阈值来筛选，完全根据你的数据情况来进行判断。

写在最后

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

需要示例数据和脚本的可以私信我~