这是2014年发表在nature communication上的一张图，大意如下：
a：作者想用181个假基因预测癌症亚型，将306个样本分成训练集（223个），验证集（83个）
b: 用训练集的数据（n=223）进行分析，采用五折交叉验证避免过拟合，选用随机森林，支持向量机，逻辑回归三种机器学习方法进行分类。
c: 三种机器学习方法的AUC曲线，发现SVM的效果最好。（AUC=0.962）
d: 用整个训练集（n=223）对验证集（n=83）进行分析，用SVM进行分类，AUC为0.922。
（原文题目：The Pan-Cancer analysis of pseudogene expression reveals biologically and clinically relevant tumour subtypes）
我把这个思路搬运到我的课题里，复现图如下：
相似度：75%。a和b是用AI画的，所以今天分享下c和d的代码
一 差异分析
假设我们有了这样的一个表达矩阵，叫Expr_all。以及表示表型数据的temp (character)。差异分析你可以用DESeq，但本次我就用T检验。
temp
##转置
Expr_all<-t(Expr_all)
##将每列变成numeric
Expr_all[,1]<-as.numeric(Expr_all[,1])
for(i in 2:535){Expr_all[,i]<-as.numeric(Expr_all[,i])}
##T检验
pvalue<-sapply(1:266,function(i){t.test(Expr_all[which(temp=="cancer"),i],Expr_all[which(temp=="normal"),i])[[3]]})
##将T检验得到的P值进行矫正，降低假阳性
pvalue1<-p.adjust(pvalue1)
##挑出P<0.05的基因
t_expr<-Expr_all[,which(pvalue1<0.05)]
t_expr<-as.data.frame(t_expr)
t_expr<-t(t_expr)
t_expr<-cbind(temp,t_expr)
得到的t_expr矩阵如下
t_expr
二 构建三种机器学习模型
1. 划分训练集和验证集
library(randomForest)
library(caret)
library(pROC)
library(caret)
inTrain<-createDataPartition(y=t_expr[,1],p=0.25,list=F)
test<-t_expr[inTrain,]
train<-t_expr[-inTrain,]
2. 构建五折交叉验证RF模型
folds<-createFolds(y=train[,1],k=5)
fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
fold_test<-train[folds[[i]],]
fold_train<-train[-folds[[i]],]
model<-randomForest(temp~.,data=fold_train,proximity=T,importance=T)
model_pre<-predict(model,newdata = fold_test,type="prob")
fc<-append(fc,as.numeric(fold_test$temp))
mod_pre<-append(mod_pre,model_pre[,1])
}
df<-cbind(fc,as.numeric(mod_pre))
3. 构建五折交叉验证SVM模型
fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
fold_test<-t_expr[folds[[i]],]
fold_train<-t_expr[-folds[[i]],]
model<-svm(temp~.,data=fold_train,probability=T)
model_pre<-predict(model,newdata = fold_test,decision.values = TRUE, probability = TRUE)
fc<-append(fc,as.numeric(fold_test$temp))
mod_pre<-append(mod_pre,as.numeric(attr(model_pre, "probabilities")[,2]))
}
df<-cbind(df,cbind(fc,mod_pre))
4. 构建五折交叉验证LR模型
max=0
num=0
fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
fold_test<-train[folds[[i]],]
fold_train<-train[-folds[[i]],]
model<-glm(temp~.,family=binomial(link=logit),data=fold_train)
model_pre<-predict(model,type='response',newdata=fold_test)
fc<-append(fc,fold_test$temp)
mod_pre<-append(mod_pre,as.numeric(model_pre))
}
df<-cbind(df,cbind(fc,mod_pre))
5. 画三者的AUC曲线
pdf("/Users/baiyunfan/desktop/ROC.pdf",height=6,width=6)
mycol <- c("slateblue","seagreen3","dodgerblue","firebrick1","lightgoldenrod","magenta","orange2")
x<-plot.roc(df[,1],df[,2],
smooth=F,
lwd=2,
ylim=c(0,1),
xlim=c(1,0),
legacy.axes=T,
main="",
col=mycol[2])
x<-plot.roc(df[,3],df[,4],
smooth=F,
add=T,
lwd=2,
ylim=c(0,1),
xlim=c(1,0),
legacy.axes=T,
main="",
col=mycol[3])
x<-plot.roc(df[,5],df[,6],
smooth=F,
add=T,
lwd=2,
ylim=c(0,1),
xlim=c(1,0),
legacy.axes=T,
main="",
col=mycol[4])
legend.name <- c(paste("RF","AUC",0.9836,sep=" "),paste("SVM","AUC",0.9903,sep=" "),paste("LR","AUC",0.9986,sep=" "))
legend("bottomright",
legend=legend.name,
col = mycol[2:4],
lwd = 2,
bty="n")
dev.off()
发现LR的AUC最高，于是单拿出来用整个train对test画
6. 用整个训练集对验证集建立LR模型以及AUC曲线
model<-glm(temp~.,data=train,family=binomial(link=logit))
model_pre<-predict(model,type='response',newdata=test)
pdf("/Users/baiyunfan/desktop/RF.pdf",height=6,width=6)
x<-plot.roc(test[,1],model_pre,
smooth=F,
lwd=2,
ylim=c(0,1),
xlim=c(1,0),
legacy.axes=T,
main="",
col=mycol[4])
legend.name <- paste("LR","AUC",0.9557,sep=" ")
legend("bottomright",
legend=legend.name,
col = mycol[4],
lwd = 2,
bty="n")
dev.off()
AUC掉下来了，0.9557，大家来探讨下为什么吧～
网友评论