DESeq2與LASSO對基因表達的預測能力
在這篇文章中，我們將比較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(j in 1:100){#RANDOMLY SPLIT DATA SET INTO TRAIN AND TESTtest_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))]#TRAIN AND EVALUATE MULTIVARIATE LASSO MODELlasso_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)#TRAIN AND EVALUATE UNIVARIATE SPEARMAN CORRELATION MODELrho<-vector(); p<-vector()for(i in 1: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(i in 1: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)#TRAIN AND EVALUATE MANN-WHITNEY U TEST UNIVARIATE MODELwilcox_stat<-vector(); p<-vector(); fc<-vector()for(i in 1: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(i in 1: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)#POPULATE FPR AND TPR VECTORS / 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"))}#AVERAGE MULTIPLE ROC-CURVES OR EACH MODEL AND PLOT MEAN ROC-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"))
在這裡我們可以得出結論，LASSO比這兩種單變量特徵選擇方法具有更大的預測能力。更好看的差異ROC曲線下的面積(AUC ROC)之間的三種方法,以及能夠執行統計測試來解決如何重要ROC曲線之間的差異,讓我們做一個箱線圖的AUC ROC套索,單變量斯皮爾曼相關和Mann-Whitney U測試。
boxplot(roc_obj_multivar_AUC,roc_obj_univar_AUC,roc_obj_wilcox_univar_AUC,names=c("LASSO","SPEARMAN","MANN-WHITNEY U TEST"),col=c("blue","green","red"))mwu<-wilcox.test(roc_obj_wilcox_univar_AUC,roc_obj_multivar_AUC)mtext(paste0("LASSO vs. Mann-Whitney U test: P-value = ",mwu$p.value))
我們可以看到，Spearman correlation和Mann-Whitney U test單變量特徵選擇模型具有相當的AUC ROC指標(儘管Mann-Whitney U test較好)，且兩者的AUC ROC即預測能力均明顯低於multivarite LASSO。
然而，這種比較可能存在偏差。談到單變量模型(Spearman和Mann-Whitney U檢驗)，我們提到只有少數基因在多次檢驗的Bonferroni校正後是顯著的。所以我們在構建單變量模型的預測得分時只使用了少量的基因，而LASSO選擇了更多的基因~30個，請參閱github上的完整代碼。如果LASSO更好的預測能力僅僅是因為它的預測分數使用了更多的特徵呢?為了驗證這一假設，在下一節中，我們將暫時忽略Bonferroni校正，並使用Spearman相關性和Mann-Whitney U檢驗，單獨使用p值排序來確定~30個最具預測性的基因。換句話說，我們將使用與多品種相同數量的基因來構建其預測得分。通過選擇的基因數量來模擬LASSO的相應的單變量模型稱為SPEAR30 (Spearman correlation with ~30 differential expressed genes)和MWU30 (Mann-Whitney U test with ~30 differential expressed genes)。
性別預測：DESeq2與多元方法
在本節中，除了將LASSO與SPEAR30（具有約30個差異表達基因的Spearman相關性）模型和MWU30（具有約30個差異表達基因的Mann-Whitney U檢驗）模型進行比較之外，我們還將添加其他一些流行的單變量和多變量模型。
首先，當進行差異基因表達分析時，DESeq2是業界的一個黃金標準。該工具享有很高的聲譽，在RNAseq社區中非常受歡迎。這是一個單變量工具，即假設基因表達計數為負二項分布，它會執行逐個基因的測試。此外，它採用了方差穩定程序，其中高表達的基因有助於低表達的基因得到正確測試。比較DESeq2預測能力與Mann-Whitney U檢驗和Spearman相關性，這些檢驗本質上利用相同的單變量思想，但是與假定基因表達為負生物學分布的DESeq2相比，都執行非參數類型的檢驗，這將是很有意思的。執行參數測試。並行地，我們將計算DESEQ2_30模型，該模型使用與LASSO選擇的相同數目的基因（按與性別相關的p值排序）建立預測得分。類似於SPEAR30和MWU30模型。
其次，我們將添加另外兩個多元特徵選擇模型，以與LASSO和單變量模型進行比較。這兩個是偏最小二乘判別分析（PLS-DA）和隨機森林，它們都是通用的多元模型。其中一個（PLS-DA）和LASSO是線性的，另一個（Random Forest）是非線性的。在這裡，我們不僅旨在比較單變量或多變量特徵選擇模型，而且還想了解與線性LASSO和PLS-DA相比，非線性隨機森林能否改善預測。
如上一節所述，我們將通過將數據集分為多次訓練和測試來圍繞ROC曲線建立置信區間。它涉及很多編碼，我在這裡不介紹代碼，但是歡迎您在我的github上檢查它們。在下面，我僅介紹經過多次訓練測試後的每個模型的平均ROC曲線，以及AUC的箱線圖。
我們在這裡觀察到一些有趣的事情。首先，與所有多變量模型相比，所有單變量模型的預測能力似乎都更差。即使是最差的多元模型，在這裡是隨機森林(RF)，其AUC ROC也明顯高於最好的單變量模型，在這裡似乎是Mann-Whitney U檢驗(MWU)。
其次，具有與LASSO選擇的基因數量相同的所有單變量模型（DESeq230，SPEAR30和MWU30）無法與所有其他單變量或多變量模型競爭，這暗示單變量模型的預測能力較差的原因不是由於數目不同特徵/基因的選擇，但由於構建預測得分的基因的等級和權重不同而異。
第三，與線性多變量LASSO和PLS-DA模型相比，非線性多變量隨機森林對RNAseq基因表達的預測效果似乎沒有改善。然而，根據我的經驗，許多生命科學問題通常都是這樣，在使用非線性分類器之前，通常值得檢查一下簡單的線性模型。
第四，也是最有趣的是，DESeq2單變量參數預測得分似乎不僅比多變量模型(LASSO, PLS-DA, Random Forest)表現更差，而且比單變量非參數模型，如Spearman相關和Mann-Whitney U檢驗也差。考慮到上述非參數測試的簡單性和DESeq2的卓越聲譽，這是相當出乎意料的。然而，事實證明，至少對於這個特定的數據集，簡單的Spearman和Mann-Whitney非參數測試在預測能力方面優於DESeq2。
LASSO，Ridge和Elastic Net
作為獎勵，在本節中，我們將比較LASSO（L1規範），Ridge（L2規範）和Elastic Net（L1和L2規範的組合）的預測能力的預測得分。 通常，懲罰線性模型族的這三個成員之間的區別並不明顯。 在不贅述的情況下（有很多文獻解釋了這種差異），我們僅強調LASSO是最保守的方法，因此由於生命科學數據的高噪聲水平，在生命科學中通常首選LASSO。 Elastic Net在Ridge的凸優化優勢和LASSO的嚴格性之間提供了很好的平衡，並且在癌症研究中非常受歡迎。 接下來，我們再次對RNAseq基因表達數據集進行100次訓練測試後，得出LASSO，Ridge和Elastic Net的平均ROC曲線，以及AUC ROC指標的箱線圖。
我們可以看到LASSO和Elastic Net給出了幾乎相同的ROC曲線，並且勝過了已知的最寬鬆的Ridge模型，這顯然不利於模型概括，因此對於預測目的無益，因此Ridge可能不是第一個生命科學數據的最佳選擇。
總結
在本文中，我們了解到，與單變量模型相比，多元統計模型似乎具有更好的預測能力。 至少對於本文研究的GTEX骨骼肌RNAseq基因表達數據而言，單變量差異基因表達工具DESeq2不僅比諸如LASSO，PLS-DA和Random Forest等多變量模型，而且與簡單非 參數單變量模型，例如Spearman相關性和Mann-Whitney U檢驗。
本文代碼：github/NikolayOskolkov/UnivariteVsMultivariteModels
作者：Nikolay Oskolkov
deephun翻譯組