花花花写于2020-02-07 闭门不出的第n天,我们将原定于下周一的学习小组提前到今天开始了,我也收到公司通知,我们可能五一之前都没办法出去做线下培训,要转型线上培训了,如果不是这场灾难,此刻我应该是在广州的酒店,准备明天的课程了。
1.准备输入数据
输入数据是TCGA的表达矩阵expr、临床信息meta和group_list。保存为forest.Rdata了,在生信星球公众号后台聊天窗口回复“森林”即可获得。
load("forest.Rdata")
exprSet = expr[,group_list=="tumor"]
## 必须保证生存资料和表达矩阵,两者一致
all(substring(colnames(exprSet),1,12)==meta$ID)
#> [1] TRUE
library(ROCR)
library(genefilter)
library(Hmisc)
library(e1071)
2.构建支持向量机模型
2.1.切割数据
用R包caret切割数据,生成的结果是一组代表列数的数字,用这些数字来给表达矩阵和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数据集建模
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.模型预测
用训练集构建模型,预测测试集的生死。不同于其他模型,这个预测结果是分类变量,直接预测生死(0,1),而不是prob(百分数)。
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
目前留有一个疑问,端详了一下模型内部结构,还不知道如何找到构建模型用到的基因。由于公众号文章一经发出无法修改,此问题后期如需更新,我将在简书(原文链接跳转)更新。
微信公众号生信星球同步更新我的文章,欢迎大家扫码关注!
我们有为生信初学者准备的学习小组,点击查看◀️
想要参加我的线上线下课程,也可加好友咨询🔼
如果需要提问,请先看生信星球答疑公告