花花花寫于2020-02-07 閉門不出的第n天，我們將原定于下周一的學(xué)習(xí)小組提前到今天開始了，我也收到公司通知，我們可能五一之前都沒辦法出去做線下培訓(xùn)，要轉(zhuǎn)型線上培訓(xùn)了，如果不是這場災(zāi)難，此刻我應(yīng)該是在廣州的酒店，準(zhǔn)備明天的課程了。

1.準(zhǔn)備輸入數(shù)據(jù)

輸入數(shù)據(jù)是TCGA的表達(dá)矩陣expr、臨床信息meta和group_list。保存為forest.Rdata了，在生信星球公眾號后臺聊天窗口回復(fù)“森林”即可獲得。

load("forest.Rdata")
exprSet = expr[,group_list=="tumor"]

## 必須保證生存資料和表達(dá)矩陣，兩者一致
all(substring(colnames(exprSet),1,12)==meta$ID)
#> [1] TRUE

library(ROCR)
library(genefilter)
library(Hmisc)
library(e1071)

2.構(gòu)建支持向量機(jī)模型

2.1.切割數(shù)據(jù)

用R包c(diǎn)aret切割數(shù)據(jù)，生成的結(jié)果是一組代表列數(shù)的數(shù)字，用這些數(shù)字來給表達(dá)矩陣和meta取子集即可。

#library(caret)
#set.seed(12345679)
#sam<- createDataPartition(meta$event, p = .5,list = FALSE)
load("sam.Rdata")

train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

2.2 train數(shù)據(jù)集建模

x=t(log2(train+1))
y=as.factor(train_meta$event)
model = svm(x,y,kernel = "linear")
summary(model) 
#> 
#> Call:
#> svm.default(x = x, y = y, kernel = "linear")
#> 
#> 
#> Parameters:
#>    SVM-Type:  C-classification 
#>  SVM-Kernel:  linear 
#>        cost:  1 
#> 
#> Number of Support Vectors:  179
#> 
#>  ( 108 71 )
#> 
#> 
#> Number of Classes:  2 
#> 
#> Levels: 
#>  0 1

2.3.模型預(yù)測

用訓(xùn)練集構(gòu)建模型，預(yù)測測試集的生死。不同于其他模型，這個預(yù)測結(jié)果是分類變量，直接預(yù)測生死(0,1)，而不是prob(百分?jǐn)?shù))。

x=t(log2(test+1))
y=as.factor(test_meta$event)
pred = predict(model, x)
table(pred,y)
#>     y
#> pred   0   1
#>    0 151  42
#>    1  36  32

目前留有一個疑問，端詳了一下模型內(nèi)部結(jié)構(gòu)，還不知道如何找到構(gòu)建模型用到的基因。由于公眾號文章一經(jīng)發(fā)出無法修改，此問題后期如需更新，我將在簡書（原文鏈接跳轉(zhuǎn)）更新。

微信公眾號生信星球同步更新我的文章，歡迎大家掃碼關(guān)注！