在这篇文章中，我们将比较LASSO、PLS、Random Forest等多变量模型与单变量模型的预测能力，如著名的差异基因表达工具DESeq2以及传统的Mann-Whitney U检验和Spearman相关。使用骨骼肌RNAseq基因表达数据集，我们将展示使用多变量模型构建的预测得分，以优于单变量特征选择模型。
骨骼肌RNAseq基因表达数据
在这里，我们将量化几种特征选择方法的预测能力:a)单变量(逐个)特征选择，b)多变量(一起)特征选择。出于演示目的，我们将使用来自GTEX人体组织基因表达联盟的骨骼肌RNAseq基因表达数据集(为简单起见，随机抽取1000个基因)。
版本6的GTEX骨骼肌数据集包括157个样本。我们装载基因表达矩阵X，去除低表达基因。为了与组学整合的后选择特征相一致，我们将使用性别作为我们将要预测的响应变量Y。让我们通过降维(如PCA)来快速可视化这些样本。
library("mixOmics") X<-read.table("GTEX_SkeletalMuscles_157Samples_1000Genes.txt", header=TRUE,row.names=1,check.names=FALSE,sep="\t") X<-X[,colMeans(X)>=1] Y<-read.table("GTEX_SkeletalMuscles_157Samples_Gender.txt", header=TRUE,sep="\t")$GENDERnames(Y)<-rownames(X) pca.gtex<-pca(X, ncomp=10) plotIndiv(pca.gtex, group=Y, ind.names=FALSE, legend=TRUE)
正如我们所看到的，PCA图并没有根据骨骼肌基因表达数据显示出男性和女性之间的明显区别。我们只能假设基因性染色体(X, Y)可以提供一个完美的男性和女性之间的分离,但是我们没有观察到这种情况,因为:1)大多数基因来自常染色体,2)随机抽样 大约有1000个基因，因此我们可能没有从X和Y染色体上取样许多基因。在下一节中，我们将把数据集分割成训练和测试子集，然后在训练集上实现单变量和多变量特征选择(训练)模型，并使用平衡假阳性(FPR)和真阳性(TPR)率的roc曲线技术对测试集上的模型进行评估。
性别预测：LASSO与单变量方法
为了评估单变量和多变量模型的预测能力，我们需要对独立的数据集进行训练和评估。为此，让我们将数据集划分为训练（80％的样本）和测试（20％的样本）子集。为了确保我们的评估不基于数据的特定“幸运”拆分，我们将多次对训练和测试子集执行随机拆分，多次计算每次ROC曲线并平均多个ROC曲线。这样，我们将一举两得：为每个模型的ROC曲线建立置信区间，并使ROC曲线平滑且美观，否则，除非您在测试中有大量样本，否则它们将会出现各种问题。
在本篇文章中,我将添加另一种非ivarite特征选择方法，即Mann-Whitney U检验，该检验应在很大程度上与Spearman相关性相比较，因为两者都是非参数和基于秩的单变量方法，因此不假定数据的特定分布。为了将通过单变量方法单独选择的基因组合到预测得分中，我们将使用它们的表达与性别之间的个体关联性的p值对它们进行排名，并通过Bonferroni程序校正多次测试。经过多次测试调整后，这通常会在训练数据集上产生一些差异表达的基因。此外，可以通过计算其在训练数据集上的作用/权重的乘积（Spearman rho和男性与女性之间基因表达的倍数变化的对数），将差异显着表达（男性与女性之间）的基因折叠成预测得分。以及来自测试数据集的样本的基因表达值。让我们对其进行编码，并比较模型之间的ROC曲线。
roc_obj_univar_FPR_list<-list(); roc_obj_univar_TPR_list<-list() roc_obj_wilcox_univar_FPR_list<-list();roc_obj_wilcox_univar_TPR_list<-list() roc_obj_multivar_FPR_list<-list(); roc_obj_multivar_TPR_list<-list() roc_obj_univar_AUC<-vector(); roc_obj_wilcox_univar_AUC<-vector() roc_obj_multivar_AUC<-vector() for(jin1:100) { #RANDOMLYSPLITDATASETINTOTRAINANDTESTtest_samples<-rownames(X)[sample(1:length(rownames(X)), round(length(rownames(X))*0.2))] train_samples<-rownames(X)[!rownames(X)%in%test_samples] X.train<-X[match(train_samples,rownames(X)),] X.test<-X[match(test_samples,rownames(X)),] Y.train<-Y[match(train_samples,names(Y))];Y.test<-Y[match(test_samples,names(Y))] #TRAINANDEVALUATEMULTIVARIATELASSOMODELlasso_fit<-cv.glmnet(as.matrix(X.train), Y.train,family="binomial",alpha=1) coef<-predict(lasso_fit, newx=as.matrix(X.test), s="lambda.min") roc_obj_multivar<-rocit(as.numeric(coef),Y.test) #TRAINANDEVALUATEUNIVARIATESPEARMANCORRELATIONMODELrho<-vector(); p<-vector() for(iin1:dim(X.train)[2]) { corr_output<-cor.test(X.train[,i],as.numeric(Y.train),method="spearman") rho<-append(rho,as.numeric(corr_output$estimate)) p<-append(p,as.numeric(corr_output$p.value)) } output_univar<-data.frame(GENE=colnames(X.train),SPEARMAN_RHO=rho,PVALUE=p) output_univar$FDR<-p.adjust(output_univar$PVALUE,method="BH") output_univar<-output_univar[order(output_univar$FDR,output_univar$PVALUE, -output_univar$SPEARMAN_RHO),] output_univar<-output_univar[output_univar$FDR<0.05,] my_X<-subset(X.test,select=as.character(output_univar$GENE)) score<-list() for(iin1:dim(output_univar)[1]) { score[[i]]<-output_univar$SPEARMAN_RHO[i]*my_X[,i] } roc_obj_univar<-rocit(colSums(matrix(unlist(score), ncol=length(Y.test), byrow=TRUE)),Y.test) #TRAINANDEVALUATEMANN-WHITNEYUTESTUNIVARIATEMODELwilcox_stat<-vector(); p<-vector(); fc<-vector() for(iin1:dim(X.train)[2]) { wilcox_output<-wilcox.test(X.train[,i][Y.train=="Male"], X.train[,i][Y.train=="Female"]) wilcox_stat<-append(wilcox_stat,as.numeric(wilcox_output$statistic)) fc<-append(fc,mean(X.train[,i][Y.train=="Male"])/mean( X.train[,i][Y.train=="Female"])) p<-append(p,as.numeric(wilcox_output$p.value)) } output_wilcox_univar<-data.frame(GENE=colnames(X.train),MWU_STAT=wilcox_stat, FC=fc,PVALUE=p) output_wilcox_univar$LOGFC<-log(output_wilcox_univar$FC) output_wilcox_univar$FDR<-p.adjust(output_wilcox_univar$PVALUE,method="BH") output_wilcox_univar<-output_wilcox_univar[order(output_wilcox_univar$FDR, output_wilcox_univar$PVALUE),] output_wilcox_univar<-output_wilcox_univar[output_wilcox_univar$FDR<0.05,] my_X<-subset(X.test,select=as.character(output_wilcox_univar$GENE)) score<-list() for(iin1:dim(output_wilcox_univar)[1]) { score[[i]]<-output_wilcox_univar$LOGFC[i]*my_X[,i] } roc_obj_wilcox_univar<-rocit(colSums(matrix(unlist(score),ncol=length(Y.test), byrow=TRUE)),Y.test) #POPULATEFPRANDTPRVECTORS/MATRICESroc_obj_univar_FPR_list[[j]]<-roc_obj_univar$FPRroc_obj_univar_TPR_list[[j]]<-roc_obj_univar$TPRroc_obj_wilcox_univar_FPR_list[[j]]<-roc_obj_wilcox_univar$FPRroc_obj_wilcox_univar_TPR_list[[j]]<-roc_obj_wilcox_univar$TPRroc_obj_multivar_FPR_list[[j]]<-roc_obj_multivar$FPRroc_obj_multivar_TPR_list[[j]]<-roc_obj_multivar$TPRroc_obj_univar_AUC<-append(roc_obj_univar_AUC,roc_obj_univar$AUC) roc_obj_wilcox_univar_AUC<-append(roc_obj_wilcox_univar_AUC, roc_obj_wilcox_univar$AUC) roc_obj_multivar_AUC<-append(roc_obj_multivar_AUC,roc_obj_multivar$AUC) print(paste0("FINISHED ",j," ITERATIONS")) } #AVERAGEMULTIPLEROC-CURVESOREACHMODELANDPLOTMEANROC-CURVESroc_obj_univar_FPR_mean<-colMeans(matrix(unlist(roc_obj_univar_FPR_list), ncol=length(Y.test)+1,byrow=TRUE)) roc_obj_univar_TPR_mean<-colMeans(matrix(unlist(roc_obj_univar_TPR_list), ncol=length(Y.test)+1,byrow=TRUE)) roc_obj_wilcox_univar_FPR_mean<-colMeans(matrix( unlist(roc_obj_wilcox_univar_FPR_list), ncol=length(Y.test)+1,byrow=TRUE)) roc_obj_wilcox_univar_TPR_mean<-colMeans(matrix( unlist(roc_obj_wilcox_univar_TPR_list), ncol=length(Y.test)+1,byrow=TRUE)) roc_obj_multivar_FPR_mean<-colMeans(matrix(unlist(roc_obj_multivar_FPR_list), ncol=length(Y.test)+1,byrow=TRUE)) roc_obj_multivar_TPR_mean<-colMeans(matrix(unlist(roc_obj_multivar_TPR_list), ncol=length(Y.test)+1,byrow=TRUE)) plot(roc_obj_univar_FPR_mean,roc_obj_univar_TPR_mean,col="green", ylab="SENSITIVITY (TPR)",xlab="1 - SPECIFISITY (FPR)",cex=0.4,type="o") lines(roc_obj_multivar_FPR_mean,roc_obj_multivar_TPR_mean,cex=0.4,col="blue", type="o",lwd=2) lines(roc_obj_wilcox_univar_FPR_mean,roc_obj_wilcox_univar_TPR_mean,cex=0.4, col="red",type="o",lwd=2) lines(c(0,1),c(0,1),col="black") legend("bottomright", inset=.02, c("Multivarite: LASSO","Univarite: SPEARMAN", "Univarite: MANN-WHITNEY U TEST"), fill=c("blue","green","red"))
