機器學(xué)習(xí)，在生信文章中出現(xiàn)的頻率越來越高，成為不可或缺的一部分分析內(nèi)容，今天小云想帶著大家學(xué)習(xí)一下生存相關(guān)機器學(xué)習(xí)的兩種算法-superpc和Ridge,接下來小云為大家對這兩種算法做一下簡單介紹，superpc這個算法的中文名稱叫做監(jiān)督主成分分析，Ridge回歸也是對參數(shù)進行約束的一種求解方程，Ridge回歸是在均方損失函數(shù)后面添加L2正則化項，能夠參數(shù)無限接近零；通過訓(xùn)練集數(shù)據(jù)進行模型構(gòu)建，然后利用測試數(shù)據(jù)進行模型預(yù)測，并計算模型的預(yù)測效能的C指數(shù)，通過該指數(shù)來判斷模型的好壞，選擇最適合的模型，如果覺得推文不錯，點贊加關(guān)注奧，話不多說，

小云開始今天的分享啦，代碼如下：

1.?安裝需要的R包

install.packages(“tidyverse”)

BiocManager::install("survcomp")

install.packages("superpc")

install.packages(“survival”)

install.packages("glmnet")

2.?載入需要的R包

library(superpc)

library(survival)

library(survcomp)

library(tidyverse)

library(glmnet)

3.?讀取數(shù)據(jù)

#TCGA.txt-訓(xùn)練集數(shù)據(jù)，第一列為樣本信息，第二列表示生存狀態(tài)，第三列表示生存時間，其他列為基因名

tcga <- read.table("TCGA.txt", header = T,sep = "\t", quote = "", check.names = F

#GSE57303-測試集數(shù)據(jù)，第一列為樣本信息，第二列表示生存狀態(tài)，第三列表示生存時間，其他列為基因名.

GSE57303 <- read.table("GSE57303.txt", header = T, sep = "\t", quote = "", check.names = F)

4.?所需數(shù)據(jù)格式整理

# 生成包含2個數(shù)據(jù)集的列表

mm <- list(TCGA = tcga,

???????????GSE57303 = GSE57303)

# 數(shù)據(jù)標準化

mm <- lapply(mm,function(x){

??x[,-c(1:3)] <- scale(x[,-c(1:3)])

??return(x)})

result <- data.frame()

# TCGA作為訓(xùn)練集

est_data <- mm$TCGA

#GSE57303作為測試集

w<-mm$GSE57303

#獲得基因名

pre_var <- colnames(est_data)[-c(1:3)]

#提取訓(xùn)練集需要的列數(shù)據(jù)

est_dd <- est_data[, c('OS.time', 'OS', pre_var)]

#提取測試集需要的列數(shù)據(jù)

w <- w[, c('OS.time', 'OS', pre_var)]

5.superpc算法進行模型構(gòu)建

#訓(xùn)練集數(shù)據(jù)列表

data <- list(x = t(est_dd[, -c(1,2)]), y = est_dd$OS.time, censoring.status = est_dd$OS, featurenames = colnames(est_dd)[-c(1, 2)])

set.seed(123)

#利用訓(xùn)練集進行模型構(gòu)建

fit <- superpc.train(data = data,type = 'survival', s0.perc = 0.5) #default

cv.fit <- superpc.cv(fit, data, n.threshold = 20, #default

?????????????????????n.fold = 10,

?????????????????????n.components = 3,

?????????????????????min.features = 5,

?????????????????????max.features = nrow(data$x),

?????????????????????compute.fullcv = TRUE,

?????????????????????compute.preval = TRUE)

#測試集數(shù)據(jù)列表

test <- list(x = t(w[,-c(1,2)]), y = w$OS.time, censoring.status = w$OS, featurenames = colnames(w)[-c(1,2)])

#利用訓(xùn)練集進行模型預(yù)測

ff <- superpc.predict(fit, data, test, threshold = cv.fit$thresholds[which.max(cv.fit[["scor"]][1,])], n.components = 1)

#將數(shù)據(jù)類型改為數(shù)字

rr <- as.numeric(ff$v.pred)

rr2 <- cbind(w[,1:2], RS = rr)

#計算模型的預(yù)測效能的C指數(shù)

cc<-data.frame(Cindex=as.numeric(summary(coxph(Surv(OS.time,OS)~ RS,rr2))$concordance[1])) %>%

??rownames_to_column('ID')

cc$Model <- paste0('SuperPC')

result <- rbind(result, cc)

6.Ridge算法進行模型構(gòu)建

#模型構(gòu)建

x1 <- as.matrix(est_dd[, pre_var])

x2 <- as.matrix(Surv(est_dd$OS.time, est_dd$OS))

set.seed(seed)

fit = glmnet(x1, x2, family = "binomial", alpha = 0, lambda = NULL)

cvfit = cv.glmnet(x1, x2,

??????????????????nfold = 10, #例文描述：10-fold cross-validation

??????????????????family = "binomial",

??????????????????type.measure = "class"

)

#模型預(yù)測

rs = as.numeric(predict(fit, type = 'response', newx = as.matrix(w[, -c(1,2)]), s = cvfit$lambda.min))

rr2 <- cbind(w[,1:2], RS = rs)

##計算模型的預(yù)測效能的C指數(shù)

cc <- data.frame(Cindex = as.numeric(summary(coxph(Surv(OS.time,OS) ~ RS, rr2))$concordance[1])) %>%

??rownames_to_column('ID')

cc$Model <- paste0('Ridge')

Result <- rbind(result, cc)

write.table(Result,file="superpc_ridge.txt",row.names=F,col.names=T,sep="\t",quote=F)

小云最終順利的利用supersc和Ridge兩種算法進行了模型構(gòu)建和預(yù)測，并計算了預(yù)測效能C指數(shù)，來判斷模型的好壞，機器學(xué)習(xí)相關(guān)的其他分析內(nèi)容歡迎嘗試本公司新開發(fā)的云平臺生物信息分析小工具，零代碼完成分析，云平臺網(wǎng)址：http://www.biocloudservice.com/home.html，主要包括lasso回歸模型篩選特征基因(http://www.biocloudservice.com/116/116.php),隨機森林的十折交叉驗證(http://www.biocloudservice.com/646/646.php),svm-ref支持向量機的機器學(xué)習(xí)(http://www.biocloudservice.com/372/372.php)等機器學(xué)習(xí)相關(guān)的小工具，

今天小云的分享就到這里，歡迎大家和小云一起討論學(xué)習(xí)，下期再見奧。