今天魄幕，和大家分享一篇生信文章的解讀與復(fù)現(xiàn)——隨機(jī)森林與人工神經(jīng)網(wǎng)絡(luò)聯(lián)合診斷心衰模型的構(gòu)建與分析。

Part1文獻(xiàn)解讀

1摘要

心力衰竭是一個(gè)全球性的健康問題仅偎，影響到全世界約 2600 萬人孕似。由于傳統(tǒng)的心力衰竭診斷技術(shù)在實(shí)踐中存在諸多局限性，有必要開發(fā)新的診斷模型來補(bǔ)充現(xiàn)有的診斷方法娘荡。

近年來，隨著基因測(cè)序技術(shù)的進(jìn)步和進(jìn)步驶沼，人們發(fā)現(xiàn)了更多與心力衰竭相關(guān)的基因炮沐。利用基因表達(dá)綜合(GEO)數(shù)據(jù)庫中現(xiàn)有的基因表達(dá)數(shù)據(jù)，篩選出心力衰竭患者差異表達(dá)基因(DEGS)回怜，并通過隨機(jī)森林分類器鑒定出 6 個(gè)關(guān)鍵基因（HMOX 2, SERPINA 3, LCN6, CSDC 2, FREM 1大年，和 ZMAT1)。

這些基因中玉雾，CSDC 2, FREM 1翔试，和 ZMAT1 從未與心力衰竭有關(guān)聯(lián)。利用人工神經(jīng)網(wǎng)絡(luò)成功地建立了心力衰竭的診斷模型复旬，并在公共數(shù)據(jù)集上驗(yàn)證了該模型的有效性垦缅。

2前言

心力衰竭(HF)是所有類型心臟病常見的一種慢性疾病。心衰實(shí)質(zhì)上是由心臟功能異常引起的病理生理狀態(tài)赢底，心臟不能滿足正常心臟壓力下正常代謝所需的抽吸速度失都。HF 可分為兩類：一種是伴有射血分?jǐn)?shù)降低的 HF柏蘑，另一種是保留射血分?jǐn)?shù)(HFpEF)的 HF幸冻，這兩類 HF 的發(fā)生和發(fā)展機(jī)制明顯不同。

HFpEF 常發(fā)生于壓力超負(fù)荷肥大疾病咳焚。與 HFrEF 相比洽损，HFpEF 更有可能降低心臟儲(chǔ)備。在 HFpEF 的發(fā)病機(jī)制中革半，心肌細(xì)胞本身的凋亡程度較小碑定，而其特征性變化是異常成纖維細(xì)胞的增殖和細(xì)胞基質(zhì)蛋白的積累。這是 HFpEF 與 HFrEF 最顯著的區(qū)別又官。

對(duì)于臨床常用的 HF 診斷技術(shù)有幾個(gè)限制延刘。腦鈉尿肽/N 端脯氨酸型鈉尿肽在各種非心衰疾病如肺動(dòng)脈高壓、肝硬化腹水六敬、急性或慢性腎功能衰竭碘赖、感染和炎癥中也可能升高，但 HFpEF 患者正常。超聲心動(dòng)圖是另一種常用的心功能評(píng)估技術(shù)普泡，它更多地依賴于專家的個(gè)人操作能力和診斷經(jīng)驗(yàn)播掷，使檢查的可重復(fù)性差。此外撼班，單純測(cè)量 EF 值歧匈，難以識(shí)別 HFpEF 患者。因此砰嘁，有必要開發(fā)新的診斷模型來補(bǔ)充現(xiàn)有的診斷方法件炉。

近年來，第二代測(cè)序技術(shù)的迅速發(fā)展為識(shí)別與多種疾病相關(guān)的標(biāo)記基因提供了基礎(chǔ)般码，為建立新的與基因相關(guān)的 HF 診斷模型奠定了堅(jiān)實(shí)的基礎(chǔ)妻率。在本研究中，我們?cè)?GEO 數(shù)據(jù)庫中篩選了 HF 與正常心肌標(biāo)本之間的差異表達(dá)基因(DEGS)板祝。在這些 DEG 數(shù)據(jù)的基礎(chǔ)上宫静，我們采用隨機(jī)森林算法對(duì) HF 中的關(guān)鍵基因進(jìn)行了識(shí)別。然后券时，我們將這些關(guān)鍵基因輸入到人工神經(jīng)網(wǎng)絡(luò)中孤里，以構(gòu)建 HF 的遺傳診斷模型.

3方法&結(jié)果

GEOQuery 軟件包用于下載數(shù)據(jù)以獲得芯片數(shù)據(jù)集的表達(dá)譜和臨床表型數(shù)據(jù):GSE57345, GSE42955，和 GSE84796 和 rna-seq 數(shù)據(jù)集:GSE141910 和 GSE116250橘洞。從 GEO 數(shù)據(jù)庫中獲取相應(yīng)平臺(tái)芯片探針的相應(yīng)標(biāo)注信息捌袜。在芯片探針 ID 和基因符號(hào)轉(zhuǎn)換過程中，發(fā)現(xiàn)多個(gè)探針對(duì)應(yīng)于 1 個(gè)基因符號(hào)炸枣，在這種情況下虏等，平均探針表達(dá)作為基因表達(dá)水平。使用 org.Hs.eg.db 包 (3.7.0 版)對(duì) rna-seq 表達(dá)式配置文件進(jìn)行基因 ID 轉(zhuǎn)換适肠。

R 軟件包 limma 對(duì) GSE57345 數(shù)據(jù)集 136 個(gè)正常和 177 個(gè)心衰樣本進(jìn)行了差異分析霍衫。limma 軟件包使用經(jīng)典的貝葉斯數(shù)據(jù)分析來篩選 DEGS。DEGS 的顯著性標(biāo)準(zhǔn)設(shè)置在 P 值小于 0.05侯养， logFoldChang (LogFC) 大于 1.5敦跌。用 phatmap 軟件包繪制 DEGS 的熱圖。

用 R 軟件包繪制聚類剖面圖對(duì)相關(guān)基因進(jìn)行 GO 功能富集分析和 KEGG 富集分析逛揩，使用了 Benjamini–Hochberg 校正方法柠傍，閾值設(shè)置為 P 值<0.01，q 值<0.01辩稽。為了避免 GO 富集結(jié)果中的冗余惧笛，我們對(duì) GO 富集項(xiàng)進(jìn)行了去重復(fù)，消除了基因重疊>0.75（詳見附錄 2）逞泄。確定三種顯著富集 GO 術(shù)語(P < 0.05) 和顯著富集通路 (P < 0.05)患整。

在這些結(jié)果中静檬，與 HF 相關(guān)的生物學(xué)過程包括細(xì)胞外基質(zhì)的組織、心臟收縮并级、巨噬細(xì)胞活化和細(xì)胞-基質(zhì)粘附拂檩。所涉及的細(xì)胞成分包括含有膠原的細(xì)胞外基質(zhì)。分子功能包括整體結(jié)合和其他重要功能嘲碧。圖 3B 顯示了部分 GO 富集的術(shù)語和顯著表達(dá)的差異基因稻励。我們還對(duì) DEGS 進(jìn)行了 kegg 通路富集分析，展示了所涉及的重要的富集生物途徑的結(jié)果和相應(yīng)的 DEG愈涩。

將 281 個(gè) DEG 輸入到隨機(jī)森林分類器中望抽。為了找到最優(yōu)參數(shù) mtry(即指定節(jié)點(diǎn)中二叉樹的最優(yōu)變量數(shù))，我們對(duì) 1-281 個(gè)變量中的所有可能數(shù)進(jìn)行了遞歸隨機(jī)森林分類履婉，并計(jì)算了模型的平均錯(cuò)誤率煤篙。圖 4A 顯示所有變量被選中時(shí)的平均錯(cuò)誤率。最后毁腿，選取 6 作為變量數(shù)的參數(shù)辑奈。變量數(shù)目越少越好，帶外誤差(out-of-band error)越小越好已烤。引用模型誤差與決策樹數(shù)之間的關(guān)系圖(圖 4B)選取 2000 棵樹作為最終模型的參數(shù)鸠窗，模型誤差穩(wěn)定。

在建立隨機(jī)森林模型的過程中胯究，從降低精度和減小均方誤差的角度出發(fā)稍计，對(duì)輸出結(jié)果的變量重要性(基尼系數(shù)法)進(jìn)行了測(cè)量。然后裕循，我們確定了 6 個(gè)重要性大于 2 的 DEGS 作為后續(xù)分析的候選基因臣嚣。圖 4C 表明在這六個(gè)變量中，HMOX 2 和 CSDC 2 是最重要的剥哑，其次是 ZMAT1, SERPINA3, FREM1硅则，和 LCN6⌒浅郑基于這六個(gè)重要變量抢埋，我們對(duì) GSE57345 數(shù)據(jù)集進(jìn)行了 k-均值無監(jiān)督聚類弹灭。圖 4D 結(jié)果表明督暂，在 GSE57345 數(shù)據(jù)集 313 例樣本中，這 6 個(gè)基因可用于患者和正常人的鑒別穷吮。其中逻翁，ZMAT1 和 FREM 1 基因在正常組織中呈低表達(dá)，在疾病組織中呈高表達(dá)捡鱼。另一方面八回，SERPINA3, LCN6, HMOX2，和 CSDC2 屬于另一組， 在正常標(biāo)本中高表達(dá)缠诅，在疾病標(biāo)本中低表達(dá)溶浴。

我們使用了另一個(gè)數(shù)據(jù)集 GSE141910 建立了基于 neuralnet 包的人工神經(jīng)網(wǎng)絡(luò)模型。第一步 是對(duì)數(shù)據(jù)進(jìn)行預(yù)處理管引，以實(shí)現(xiàn)數(shù)據(jù)的規(guī)范化士败。其次，選擇最小-最大值方法[0褥伴，1]谅将，在訓(xùn)練神經(jīng)網(wǎng)絡(luò)之前按下分離縮放數(shù)據(jù)。在開始計(jì)算之前重慢，對(duì)最大和最小數(shù)據(jù)值進(jìn)行標(biāo)準(zhǔn)化處理饥臂，并將隱藏層數(shù)設(shè)置為 5 層。在參數(shù)的選擇上似踱，不存在使用多少層和多個(gè)神經(jīng)元的固定規(guī)則隅熙。神經(jīng)元的數(shù)量應(yīng)該介于輸入層大小和輸出層大小之間，通常是輸入層大小的三分之二核芽。為了更有效地評(píng)價(jià)神經(jīng)網(wǎng)絡(luò)模型的結(jié)果猛们，我們選擇了一種 5 折交叉驗(yàn)證方法。數(shù)據(jù)集隨機(jī)分為訓(xùn)練集和驗(yàn)證集狞洋。訓(xùn)練集的目的是確定候選指標(biāo)的權(quán)重弯淘。利用驗(yàn)證集驗(yàn)證了基于基因表達(dá)和基因權(quán)重構(gòu)建的模型評(píng)分的分類效率。得到疾病神經(jīng)網(wǎng)絡(luò)模型分類評(píng)分吉懊。

5 折交叉驗(yàn)證結(jié)果使用受試者工作曲線(ROC)顯示模型分類性能(圖 5A)庐橙。此外，還使用了一個(gè)混淆矩陣來評(píng)估預(yù)測(cè)的性能(表 1)借嗽。

五折交叉驗(yàn)證結(jié)果的 ROC 曲線下面積接近 1(平均 AUC>0.99)态鳖，表明模型具有較強(qiáng)的魯棒性(Robust)。因此恶导，接下來浆竭，我們利用所有的數(shù)據(jù)來構(gòu)建神經(jīng)網(wǎng)絡(luò)模型。

根據(jù)神經(jīng)網(wǎng)絡(luò)模型的輸出結(jié)果(補(bǔ)充檔案 4 和圖 5B)惨寿，可以看出邦泄，整個(gè)培訓(xùn)是按照 1423 步驟進(jìn)行的。終止條件為誤差函數(shù)的絕對(duì)偏導(dǎo)數(shù)<0.01(達(dá)到閾值)裂垦。輸出結(jié)果表明顺囊，該模型的權(quán)重范圍為 ?4.67~4.53。重量預(yù)測(cè)為 4.527373(HMOX 2)蕉拢，?4.7670777(CSDC 2), 1.478590 (ZMAT1), 2.332519 (SERPINA 3)特碳，?4.522891(FREM 1)诚亚，以及 1.940819(LCN6).

使用三個(gè)獨(dú)立的驗(yàn)證數(shù)據(jù)集 GSE116250, GSE42955，和 GSE84796 在進(jìn)行最大和最小標(biāo)準(zhǔn)化數(shù)據(jù)處理后午乓，計(jì)算出三個(gè)評(píng)分站宗，并對(duì)其分類效率進(jìn)行評(píng)價(jià)，并對(duì) AUC 進(jìn)行比較益愈。這三個(gè)評(píng)分結(jié)果 如下：1)neuraHF份乒，將本研究中所識(shí)別的 DGS 與神經(jīng)網(wǎng)絡(luò)中的權(quán)重相乘得到的分?jǐn)?shù)；2)CD8K腕唧；3)TP 53或辖，這是文獻(xiàn)中報(bào)道的與 HF 疾病相關(guān)的特征基因。

圖 6 顯示了三個(gè)獨(dú)立驗(yàn)證數(shù)據(jù)集的三個(gè)分?jǐn)?shù)的比較枣接。在 GSE116250 數(shù)據(jù)集(圖 6A neuraHF颂暇、 CD8K 和 TP 53 的 AUC 分別為 0.991、0.683 和 0.597但惶。neuraHF 的敏感性為 100%耳鸯，特異性為 96%。在 GSE 42955 數(shù)據(jù)集(圖 6B 神經(jīng)高頻膀曾、CD8K 和 TP 53 的 AUC 分別為 0.858县爬、0.517 和 0.65。neuraHF 的敏感性為 80%添谊，特異性為 95.8%财喳。在核查結(jié)果中 GSE 84796 (圖 6C neuraHF、 CD8K 和 TP 53 的 AUC 分別為 0.871斩狱、0.586 和 0.486耳高。neuraHF 的敏感性和特異性分別為 85.7%和 80%。

4討論

在本研究中所踊，我們首次計(jì)算了與 HF 相關(guān)的 DEGS泌枪，并通過隨機(jī)森林分類器獲得了 6 個(gè)重要候選 DEGS。利用神經(jīng)網(wǎng)絡(luò)模型確定了相關(guān)基因的預(yù)測(cè)權(quán)重秕岛，構(gòu)建了與心衰相關(guān)的 neuraHF 分類模型碌燕，并在三個(gè)獨(dú)立樣本數(shù)據(jù)集上評(píng)價(jià)了模型評(píng)分的分類效率。與其他兩種 HF 相關(guān)的生物標(biāo)志物 相比继薛，AUC 的分類效率較高修壕，且具有較好的分類效率。然而惋增，由于 rna-seq 在構(gòu)建神經(jīng)網(wǎng)絡(luò)模型時(shí)所預(yù)測(cè)的權(quán)重在理論上更適用于 rna-seq 數(shù)據(jù)的疾病分類叠殷，GSE116250 數(shù)據(jù)集在驗(yàn)證結(jié)果中顯示出最佳性能改鲫。同時(shí)诈皿，由于 GEO 數(shù)據(jù)庫中缺乏 HFpEF 的基因數(shù)據(jù)，HFpEF 的遺傳特征沒有被納入診斷模型的構(gòu)建中，從而影響了 HFpEF 模型的診斷效果嫉父。

更有趣的是富俄，在構(gòu)建 HF 診斷模型的分析過程中，我們首次發(fā)現(xiàn)了三個(gè)關(guān)鍵基因(CSDC 2, FREM 1截歉，和 ZMAT1)可能在 HF 的發(fā)病機(jī)制中起一定作用胖腾。

心臟標(biāo)本的獲取困難可能會(huì)減少 HF 的潛在應(yīng)用。然而瘪松，我們目前的研究并不打算完全取代現(xiàn)有的診斷和治療方法咸作，而是為了補(bǔ)充這些方法。通常宵睦，目前的診斷標(biāo)準(zhǔn)和程序是基于來自 HFrEF 患者的數(shù)據(jù)记罚。然而，仍不清楚這些是否完全適用于 HFpEF 患者壳嚎。例如桐智，使用這些非侵入性方法很難診斷 HFpEF 的輕微癥狀。然而烟馅，從我們的研究中得出的診斷模型可以通過及時(shí)的心臟活檢來確定心力衰竭的可能性说庭。因此，我們的方法具有一定的臨床應(yīng)用價(jià)值郑趁。顯然刊驴，根據(jù)我們目前的研究結(jié)果，該模型的精度還有待進(jìn)一步研究寡润。

5總結(jié)

亮點(diǎn)：

1.使用了機(jī)器學(xué)習(xí)來篩選 hub gene

2.使用了深度學(xué)習(xí)來構(gòu)建診斷模型

3.使用了多個(gè)數(shù)據(jù)集進(jìn)行數(shù)據(jù)分析

4.文章思路清晰缺脉，邏輯明了

不足：

1.文章里包含明顯錯(cuò)誤如文中實(shí)際 logFC 是>0.5 而不是 1.5,GSE141910 在 GEO 中僅包含 366 個(gè)樣本而非 399 個(gè)樣本

2.流程圖未描述 kegg,帶外誤差的英文為 out-of-bag error

3.本文臨床意義不大

Part2文章復(fù)現(xiàn)

step1 start GEO

#數(shù)據(jù)下載

rm(list?=?ls())

library(GEOquery)

gse?="GSE57345"

Sys.setenv("VROOM_CONNECTION_SIZE"=?2621440*2)

eSet?<-?getGEO(gse,

destdir?='.',

getGPL?=?F)

#(1)提取表達(dá)矩陣exp

exp?<-?exprs(eSet[[1]])

exp[1:4,1:4]

#boxplot(exp)

#exp?=?log2(exp+1)

#(2)提取臨床信息

pd?<-?pData(eSet[[1]])

#(3)調(diào)整pd的行名順序與exp列名完全一致

p?=?identical(rownames(pd),colnames(exp));p

if(!p)?exp?=?exp[,match(rownames(pd),colnames(exp))]

#(4)提取芯片平臺(tái)編號(hào)

gpl?<-?eSet[[1]]@annotation

save(eSet,gse,pd,exp,gpl,file?="step1output.Rdata")

因?yàn)?GSE57345 是常規(guī)的微陣列芯片，所以我們按照標(biāo)準(zhǔn)流程代碼來就好悦穿，得到了 GSE57345 的表達(dá)矩陣攻礼。

step2 group id

#1.group_list(實(shí)驗(yàn)分組)和ids(芯片注釋)

rm(list?=?ls())

load("step1output.Rdata")

library(stringr)

library(dplyr)

gr?=?pd$characteristics_ch1.1

k1?=?str_detect(gr,"yes");table(k1)

k2?=?str_detect(gr,"no");table(k2)

group?=?ifelse(k1,"disease",'normal');table(group)

group?=?factor(group,levels?=?c("normal","disease"))

#2.探針注釋----

#library(devtools)

#install_github("jmzeng1314/idmap1")

library(idmap1)

temp=getIDs("GPL11532")

rowname?<-?data.frame(probe_id=rownames(exp))

gene?<-?left_join(rowname,temp,by='probe_id')

gene?<-?na.omit(gene)

exp?<-?exp[match(gene$probe_id,rownames(exp)),]

#?多個(gè)探針對(duì)應(yīng)一個(gè)基因的話，取平均表達(dá)量最大的那個(gè)探針

result?<-?data.frame(ex=rowMeans(exp),

symbol=gene$symbol)

result?<-?arrange(result,symbol,desc(ex))

de?<-?result[!duplicated(result$symbol),]

exp?<-?exp[match(rownames(de),rownames(exp)),]

exp_sy?<-?data.frame(symbol=de$symbol,

exp)

save(eSet,pd,exp_sy,group,gpl,file?="step2output.Rdata")

將樣本分組栗柒，用曾老師的 idmap1 包給探針注釋礁扮，這里開始和文章不一樣了。我們經(jīng)常會(huì)遇到多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況瞬沦，但是處理方法有很多種太伊，我們選擇取平均表達(dá)量最大的那個(gè)探針。得到 18837 個(gè)基因逛钻。

step3 pca

rm(list?=?ls())

load(file?="step1output.Rdata")

load(file?="step2output.Rdata")

exp?=?exp_sy[,-1]

group_list?=?group

dat=as.data.frame(t(exp))

library(FactoMineR)#畫主成分分析圖需要加載這兩個(gè)包

library(factoextra)

#?pca的統(tǒng)一操作走起

dat.pca?<-?PCA(dat,?graph?=?FALSE)

pca_plot?<-?fviz_pca_ind(dat.pca,

geom.ind?="point",#?show?points?only?(nbut?not?"text")

col.ind?=?group_list,#?color?by?groups

#palette?=?c("#00AFBB",?"#E7B800"),

addEllipses?=?TRUE,#?Concentration?ellipses

legend.title?="Groups"

)

pca_plot

save(pca_plot,file?="pca_plot.Rdata")

做主成分分析圖僚焦，看看該芯片中對(duì)照組和實(shí)驗(yàn)組的差異性。

step4 DEG

rm(list?=?ls())

load(file?="step2output.Rdata")

#差異分析曙痘，用limma包來做

#需要表達(dá)矩陣和group_list芳悲，不需要改

exp?=exp_sy[,-1]

group_list?=?group

library(limma)

design=model.matrix(~group_list)

fit=lmFit(exp,design)

fit=eBayes(fit)

deg=topTable(fit,coef=2,number?=?Inf)

#為deg數(shù)據(jù)框添加幾列

#1.加probe_id列立肘，把行名變成一列

library(dplyr)

deg?<-?mutate(deg,probe_id=rownames(deg))

head(deg)

#2.加symbol列，火山圖要用

deg?<-?cbind(deg,symbol=exp_sy$symbol)

head(deg)

#3.加change列,標(biāo)記上下調(diào)基因

logFC_t=0.5

P.Value_t?=?0.05

k1?=?(deg$P.Value?<?P.Value_t)&(deg$logFC<?-logFC_t)

k2?=?(deg$P.Value?<?P.Value_t)&(deg$logFC>?logFC_t)

change?=?ifelse(k1,

"down",

ifelse(k2,

"up",

"stable"))

table(change)

deg?<-?mutate(deg,change)

ex_deg?<-?deg[k1|k2,]

save(group_list,deg,logFC_t,P.Value_t,ex_deg,

file?="step4output.Rdata")

做差異分析名扛，logfc 取 0.5谅年，結(jié)果如下，得到 426 個(gè)差異基因肮韧。

step5 random forest

rm(list?=?ls())

library(tidyverse)

load(file?="step2output.Rdata")

load(file?="step4_output.Rdata")

df?<-?exp_sy[match(ex_deg$symbol,exp_sy$symbol),]

rownames(df)?<-?df$symbol

df?<-?data.frame(t(df[,-1]))?%>%

cbind(group,.)

library(randomForest)

set.seed(1423)

fit.forest?<-?randomForest(group~.,?data=df,

ntree=3000)

fit.forest

plot(fit.forest)

#找出使模型準(zhǔn)確率達(dá)到最優(yōu)所需要的樹的數(shù)量

which.min(fit.forest$err.rate[,1])

set.seed(1423)

fit.forest_2?<-?randomForest(group~.,?data=df,

ntree=49)

fit.forest_2

(127+173)/313

#?模型正確率為95.85%

varImpPlot(fit.forest_2)

result_df?<-?data.frame(importance(fit.forest_2,type=2))?%>%

arrange(desc(.))

impo_df?<-?filter(result_df,MeanDecreaseGini>4)

impo_df

#?選取了比較重要的幾個(gè)基因放入后續(xù)分析融蹂，

#?但是我覺得更好的方法應(yīng)該是使用Boruta包進(jìn)行特征選擇。

#?Kursa?M.,?Rudnicki?W.?(2010),?Feature?Selection?with?the?Boruta?Package,

#?Journal?of?Statistical?Software,?36(11),?1?-?13

impo_exp?<-?exp_sy[match(rownames(impo_df),exp_sy$symbol),]

rownames(impo_exp)?<-?impo_exp[,1]

n=impo_exp[,-1]

dim(n)

library(pheatmap)

annotation_col=data.frame(group=group_list)

rownames(annotation_col)=colnames(n)

library(ggplotify)

heatmap_plot?<-?as.ggplot(pheatmap(n,show_colnames?=F,

show_rownames?=?T,

cluster_cols?=?T,

scale?='row',

annotation_col=annotation_col))

heatmap_plot

rownames(impo_exp)

#?"ALG3"?????"AGFG1"????"ADAMTS15"?"ACE"??????"ADGRD1"

save(list?=?ls(),file?="step5output.Rdata")

我們?cè)O(shè)隨機(jī)種子為1423弄企，然后做隨機(jī)森林超燃，確定最佳隨機(jī)樹數(shù)目為49

看看變量重要性統(tǒng)計(jì)圖

選取了基尼指數(shù)平均減少量大于4的變量作為后續(xù)分析的重要變量。

并畫了個(gè)熱圖拘领。

step6 neuralnet

rm(list?=?ls())

filelist?<-?dir('gse141910')

file_path?<-?paste('gse141910',filelist,sep?='/')

data_ne?<-??read.csv(file_path[1])

for(nin2:length(file_path))?{

a?<-?read.csv(file_path[n])

data_ne?<-?merge(data_ne,a)

}

#?轉(zhuǎn)換基因id

library(clusterProfiler)

library(org.Hs.eg.db)

keytypes(org.Hs.eg.db)##?查看org.Hs.eg.db包含的數(shù)據(jù)

#?browseVignettes("clusterProfiler")?#查看教程

eg?=?bitr(data_ne$X,?fromType="ENSEMBL",

toType="SYMBOL",?OrgDb="org.Hs.eg.db")

impo_sy?<-?c("ALG3","AGFG1","ADAMTS15","ACE","ADGRD1")

nt?<-?eg[match(impo_sy,eg$SYMBOL),]

nt_exp?<-?data_ne[match(nt$ENSEMBL,data_ne$X),]?%>%

.[,-1]?%>%

scale()

rownames(nt_exp)?<-?impo_sy

library(tidyverse)

k?<-?str_detect(colnames(nt_exp),'P')

table(k)

#?240名病人淋纲，126名正常人

nt_gr?<-?rbind(nt_exp,group=k)?%>%

t(.)?%>%

data.frame()

nt_gr$group<-?factor(nt_gr$group)

library(neuralnet)

fit?<-?neuralnet(group~.,nt_gr,

hidden?=?5)

fit$result.matrix

plot(fit)

#?我也不知道做出來的結(jié)果對(duì)不對(duì)，反正看得不是很懂院究，

#?不知道那篇文章中的權(quán)重是怎么得出來的洽瞬。

#?但是如果真要應(yīng)用深入學(xué)習(xí)的方法來研究問題時(shí)，

#?應(yīng)該多嘗試一些不同的算法或r包來比較一下結(jié)果业汰。

#?比如TensorFlow,MXNet伙窃，H2O包等。

save(list?=?ls(),file?="step6output.Rdata")

GSE141910是rna-seq芯片样漆，處理原始數(shù)據(jù)可能有點(diǎn)麻煩为障，所以我就直接用了芯片上傳者處理好的表達(dá)矩陣來做后續(xù)分析。

用Y叔的clusterProfiler包放祟，和org.Hs.eg.db數(shù)據(jù)庫來對(duì)探針進(jìn)行注釋鳍怨，并提取出之前的5個(gè)重要變量。

然后構(gòu)建神經(jīng)網(wǎng)絡(luò)模型跪妥，隱藏神經(jīng)元數(shù)量設(shè)為5鞋喇，和文章里的圖比起來就丑好多了

。

step7 auc

rm(list?=?ls())

load('step6output.Rdata')

library(ggplot2)

library(pROC)

source("ROC_Plot.R")

source("Theme_Publication.R")

HUB?<-?nt_gr

colnames(HUB)

#?"ALG3","AGFG1","ADAMTS15","ACE","ADGRD1","group"

x1?<-?HUB$ALG3

x2?<-?HUB$AGFG1

x3?<-?HUB$ADAMTS15

x4?<-?HUB$ACE

x5?<-?HUB$ADGRD1

y?<-?HUB$group

R1?<-?roc(y,x1,ci=T);?R2?<-?roc(y,x2,ci=T);?R3?<-?roc(y,x3,ci=T)

R4?<-?roc(y,x4,ci=T);R5?<-?roc(y,x5,ci=T)

dat1?<-?data.frame(Sensitivities?=?R1$sensitivities,

FalsePositiveRate?=?(1-R1$specificities),

Model?=?rep("ARG1",length(R1$sensitivities)))

dat2?<-?data.frame(Sensitivities?=?R2$sensitivities,

FalsePositiveRate?=?(1-R2$specificities),

Model?=?rep("CXCL1",length(R2$sensitivities)))

dat3?<-?data.frame(Sensitivities?=?R3$sensitivities,

FalsePositiveRate?=?(1-R3$specificities),

Model?=?rep("LRG1",length(R3$sensitivities)))

dat4?<-?data.frame(Sensitivities?=?R4$sensitivities,

FalsePositiveRate?=?(1-R4$specificities),

Model?=?rep("MMP9",length(R4$sensitivities)))

dat5?<-?data.frame(Sensitivities?=?R5$sensitivities,

FalsePositiveRate?=?(1-R5$specificities),

Model?=?rep("OSCAR",length(R5$sensitivities)))

dat?<-?rbind(dat1,dat2,dat3,dat4,dat5)

colnames(dat)[2]?<-?("1-Specificities")

ggplot(dat,aes(`1-Specificities`,Sensitivities,colour?=?Model))?+?geom_roc_plot()+

scale_colour_Publication()?+?theme_Publication()+

theme(panel.border?=?element_rect(colour?='black'),legend.title?=?element_blank(),

legend.position?=?c(0.9,0.15),legend.direction?='vertical')

#做表

sy?<-?(c("ALG3","AGFG1","ADAMTS15","ACE","ADGRD1"))

auc?<-?data.frame(rn=?c('low','auc','up'))

roc?<-?list(R1,R2,R3,R4,R5)

for(nin1:length(sy))?{

auc1?<-?data.frame(roc[[n]]$ci)

auc?<-?cbind(auc,auc1)

}

rownames(auc)?<-?auc[,1]

auc?<-?auc[,-1]

colnames(auc)?<-?sy

write.csv(auc,'auc.csv')

save(list?=?ls(),file?="step7output.Rdata")

#?這里只是舉了個(gè)畫auc的例子眉撵，并沒有像文章一樣去獨(dú)立驗(yàn)證了侦香，

#?因?yàn)槲也恢郎窠?jīng)網(wǎng)絡(luò)模型怎么在其他芯片里驗(yàn)證并做roc曲線。

我并沒有像原文一樣去做auc纽疟，原文的思想是用GSE57345和GSE141910兩個(gè)大樣本芯片作為訓(xùn)練集罐韩，再用三個(gè)小樣本GSE42995，GSE84796污朽，GSE116250作為驗(yàn)證集散吵。最近有點(diǎn)其他事耽誤了，就沒有做這一步，但是還是找了一下之前做auc的代碼放上來了矾睦。

還有富集分析忘了做了晦款，大家可以參考一下其他推文：批量GO-KEGG富集分析注釋何必自己寫腳本，看Y叔的神器