前言

今天來討論下利用Tree-Based的方法來建立基因之間的調(diào)控關(guān)系酒唉，最近看到一篇17年發(fā)表在Nature Methods的文章《SCENIC: single-cell regulatory network inference and clustering》臀玄，其中有一個步驟是利用樹模型來預(yù)測轉(zhuǎn)錄因子作用的靶基因流椒，而這一步調(diào)用的是GENIE3來做預(yù)測的，所以這一次我們來看看GENIE3是如何進(jìn)行工作的
附上GitHub地址：https://github.com/vahuynh/GENIE3/tree/master/GENIE3_R

步驟

首先我們先來看一下輸入的數(shù)據(jù)長什么樣子：

行名為不同的sample婿牍，而列名為基因名哥牍，表格可以理解為標(biāo)準(zhǔn)化后的表達(dá)量，接下來就可以上代碼了：

# transpose expression matrix to (samples x genes)
#將普通的基因表達(dá)譜進(jìn)行轉(zhuǎn)換环鲤，sample為行名纯趋，基因為列名，
expr.matrix <- t(expr.matrix)
# setup weight matrix
num.samples <- dim(expr.matrix)[1]
num.genes <- dim(expr.matrix)[2]
gene.names <- colnames(expr.matrix)
weight.matrix <- matrix(0.0, nrow=num.genes, ncol=num.genes)
rownames(weight.matrix) <- gene.names
colnames(weight.matrix) <- gene.names
# get number of input genes, names of input genes
# 這里的input.gene.names指的是所有的基因
input.gene.names <- gene.names
# compute importances for every target gene
# 遍歷每一個基因冷离，target.gene.name是其中的一個基因
for (target.gene.idx in seq(from=1, to=num.genes)) {
  target.gene.name <- gene.names[target.gene.idx]
  # remove target gene from input genes
  these.input.gene.names <- setdiff(input.gene.names, target.gene.name)
  # 對遍歷到的基因做差集吵冒，setdiff是差集的意思，these.input.gene.names則是除去target.gene.name以外的所有基因
  num.input.genes <- length(these.input.gene.names)
  # 以除去target.gene.name基因以外酒朵，剩下基因的表達(dá)量作為決策變量
  x <- expr.matrix[,these.input.gene.names]
  # target.gene.name基因的表達(dá)量作為響應(yīng)變量
  y <- expr.matrix[,target.gene.name]
  # set mtry, mtry為特征值變量
  mtry <- round(sqrt(num.input.genes))
  # Normalize output
  y <- y / sd(y)
  # 隨機(jī)森林回歸建模
  rf <- randomForest(x, y, mtry=mtry, ntree=1000, importance=FALSE)  
  im <- importance(rf)[,'IncNodePurity']
  im.names <- names(im)
  weight.matrix[im.names, target.gene.name] <- im
}
result = weight.matrix / num.samples

注意我們進(jìn)行建模的變量的輸入和輸出：

# 以除去target.gene.name基因以外桦锄，剩下基因的表達(dá)量作為決策變量，有若干列表達(dá)數(shù)據(jù)
x <- expr.matrix[,these.input.gene.names]
# target.gene.name基因的表達(dá)量作為響應(yīng)變量蔫耽，為一列數(shù)據(jù)
y <- expr.matrix[,target.gene.name]

#建模
rf <- randomForest(x, y, mtry=mtry, ntree=1000, importance=FALSE)

其中结耀，隨機(jī)森林回歸的決策變量是除去目標(biāo)基因以外的所有基因表達(dá)矩陣，響應(yīng)變量是該目標(biāo)基因的表達(dá)向量匙铡。那么決策變量的列名（一串基因名）可以理解為回歸的變量图甜；而響應(yīng)變量可以理解為在各個sample中，該基因的表達(dá)量

這里在randomForest中有一個比較重要的參數(shù)叫IncNodePurity鳖眼，維基上給的解釋是：

Mean Decrease Gini (IncNodePurity) - This is a measure of variable importance based on the Gini impurity index used for the calculating the splits in trees. The higher the value of mean decrease accuracy or mean decrease gini score, the higher the importance of the variable to our model.

是衡量在建樹過程中黑毅，該變量作為判別回歸變量對模型是否重要的指標(biāo)（IncNodePurity）

事實上我們在建樹的時候，將y當(dāng)作響應(yīng)變量钦讳，將x當(dāng)作回歸的決策變量矿瘦，比方說第一個target.gene.name為TBX3，第一個指標(biāo)為GATA5=0.13572435愿卒，那么GATA5表達(dá)值的改變會對TBX3的表達(dá)值的改變產(chǎn)生比較大的影響缚去，在生物學(xué)上可以理解為這兩個基因有互作，即改變GATA5的表達(dá)量琼开，那么會對TBX3的表達(dá)量產(chǎn)生重要影響易结，可解釋為這兩個基因有互作，且是GATA5調(diào)控TBX3

最終的result為：

它反應(yīng)了基因之間的權(quán)重關(guān)系柜候，值得注意的是該矩陣并不是實對稱矩陣搞动，所以行基因為from gene，而列基因為to gene渣刷，比方說TBX3 -> GATA5 和 GATA5 -> TBX3的權(quán)重是不一樣的鹦肿，最終整理為：

那么接下來，根據(jù)權(quán)重閾值可以進(jìn)行篩選辅柴，選出有互作的基因Candidates

小tip

(1). 隨機(jī)森林回歸的概念

其實真正區(qū)分回歸還是分類問題的關(guān)鍵是響應(yīng)變量的類型狮惜，若響應(yīng)變量為數(shù)量型變量（數(shù)值型變量高诺，連續(xù)型變量）則屬于回歸問題；若響應(yīng)變量為是分類型變量（因子型變量碾篡，定性變量）虱而，則屬于分類問題
我們簡單介紹一些隨機(jī)森林回歸的大致原理：
由上面的代碼為例子，

  #以其中第一個基因為例子
  target.gene.name <- gene.names[1]
  # remove target gene from input genes
  these.input.gene.names <- setdiff(input.gene.names, target.gene.name)
  # 對遍歷到的基因做差集开泽，setdiff是差集的意思牡拇，these.input.gene.names則是除去target.gene.name以外的所有基因
  num.input.genes <- length(these.input.gene.names)
  # 以除去target.gene.name基因以外，剩下基因的表達(dá)量作為決策變量
  x <- expr.matrix[,these.input.gene.names]
  # target.gene.name基因的表達(dá)量作為響應(yīng)變量
  y <- expr.matrix[,target.gene.name]
  # set mtry, mtry為特征值變量
  mtry <- round(sqrt(num.input.genes))
  # Normalize output
  y <- y / sd(y)
  # 隨機(jī)森林回歸建模
  rf <- randomForest(x, y, mtry=mtry, ntree=1000, importance=FALSE)  
  im <- importance(rf)[,'IncNodePurity']
  im.names <- names(im)

建模的數(shù)據(jù)類型如下圖所示穆律，隨機(jī)森林回歸的本質(zhì)是若干個樹回歸

假設(shè)我們擁有D = { (x1,y1), (x2,y2), ..., (xn,yn) }惠呼，i 代表上面數(shù)據(jù)框的每一行； xi 是一個 m 維的向量峦耘，即 xi 含有 m 個 features
那么樹回歸是將數(shù)據(jù) D 里面的元素按照決策變量的條件由父節(jié)點分散到子節(jié)點剔蹋，假設(shè)一共有 N 個葉子，每個葉子分得一定量的數(shù)據(jù)辅髓，而最終的結(jié)果顯示在終子節(jié)點（沒有往下分叉的節(jié)點為終子節(jié)點）上

注意：這里的決策變量的條件指的是按照每一個決策變量的取值范圍（以GATA5為例泣崩，小于 0.05 的為一類，介于 0.05 到 0.1的為一類洛口，大于 0.1 的又為一類矫付，其他決策變量以此類推），那么按照決策變量的條件篩選出的 yi 將會被分配到不同的節(jié)點里面第焰，理論上按照不同的決策變量篩選的 yi 值應(yīng)該會具有相同的特征买优，即它們的值（ yi ）相差不會太大

那么接下來根據(jù)決策變量的條件將 D 里面的數(shù)據(jù)平均分配給這些節(jié)點，每個節(jié)點擬合一個值（比方說平均值）挺举，最終要求：

即每一個節(jié)點分得若干個數(shù)據(jù)杀赢，用每個節(jié)點的平均值代替每個節(jié)點的特征值，那么每個終子節(jié)點里面的數(shù)據(jù)與該終子節(jié)點的均值MSE最小湘纵，則達(dá)到我們的要求
每個模型里面的均方差達(dá)到最小脂崔，那么這個回歸樹才算合格，而隨機(jī)森林回歸則是由若干個回歸樹組成瞻佛，選擇最優(yōu)的那一個

最后我們不妨用模型預(yù)測一下：

predict = predict(rf,x)

貌似還是有些差別的脱篙。娇钱。伤柄。

(2). IncNodePurity的概念

根據(jù)前面所敘述的那樣，IncNodePurity是基于基尼系數(shù)計算的值文搂，而基尼系數(shù)越大适刀，代表分出的類不確定性較大，分類效果不好煤蹭。（IncNodePurity代表基尼系數(shù)的下降）
“IncNodePurity”即increase in node purity笔喉，值越大表示該變量的重要性越大取视。
IncNodePurity的大小是衡量某一個回歸變量對隨機(jī)森林回歸模型貢獻(xiàn)的一個指標(biāo)，數(shù)值越大常挚，貢獻(xiàn)越大作谭，反之越小

varImpPlot(rf,main = "The importance")

上圖描述的就是各個回歸變量對模型的貢獻(xiàn)
當(dāng)然，lncMSE是衡量回歸變量重要性的一種指標(biāo)奄毡，lncMSE越高折欠，該變量越重要，反之越不重要
lncMSE即Increase in MSE(%)值吼过，更高的MSE％值意味著更重要的變量

除此之外還有幾個衡量分類指標(biāo)重要性的因素：Mean Decrease Accuracy 和 Mean Decrease Gini 為隨機(jī)森林模型中的兩個重要指標(biāo)锐秦。其中，mean decrease accuracy 表示隨機(jī)森林預(yù)測準(zhǔn)確性的降低程度盗忱，該值越大表示該變量的重要性越大酱床；mean decrease gini 計算每個變量對分類樹每個節(jié)點上觀測值的異質(zhì)性的影響，從而比較變量的重要性趟佃。該值越大表示該變量的重要性越大

GENIE3預(yù)測基因調(diào)控網(wǎng)絡(luò)

GENIE3預(yù)測基因調(diào)控網(wǎng)絡(luò)

前言

步驟

小tip

(1). 隨機(jī)森林回歸的概念

(2). IncNodePurity的概念