Jimmy老師去年解析了CIBERSORT使用SVM算法實現去卷積,我決定亦步亦趨的跟著Jimmy老師的代碼學習
背景知識首先有一些背景知識需要了解（特別是一些算法），但是我的理解方法特別粗暴，不知道Jimmy老師會不會打我。當然了，如果是原始的CIBERSORT R腳本 https://rdrr.io/github/singha53/amritr/src/R/supportFunc_cibersort.R 其實懂得使用即可。
最開始接觸非負矩陣分解其實是在鹼基突變的signature知識點：
鹼基突變類型先說6種鹼基組合與96種組合：
以下解釋來源自wiki百科https://en.wikipedia.org/wiki/Mutational_signatures
首先要明白，鹼基突變共有六類鹼基取代：C> A，C> G，C> T，T> A，T> C，T> G。為什麼只有6種呢？因為G> T取代被認為等同於C> A取代，因為不可能區分最初發生在哪條DNA鏈（正向或反向）上。因此，C> A和G> T替換都計為「 C> A」類的一部分。出於相同的原因，G> C，G> A，A> T，A> G和A> C突變被計為「 C> G」，「 C> T」，「 T> A」，「 T> C」和「 T> G」類。
從5'和3'相鄰鹼基（也稱為側翼鹼基對或三核苷酸上下文）中獲取信息會導致96種可能的突變類型（例如A [C> A] A，A [C> A] T等）。腫瘤的突變目錄是通過將96種突變類型之一中的每個單核苷酸變體（SNV）分類（同義詞：鹼基對取代或置換點突變）並計算這96種突變類型中每種突變的總數來創建的（見圖）。
以下解釋來源自知乎專題https://zhuanlan.zhihu.com/p/27460660
NMF的基本思想可以簡單描述為：對於任意給定的一個非負矩陣V，NMF算法能夠尋找到一個非負矩陣W和一個非負矩陣H，使得滿足 ，從而將一個非負的矩陣分解為左右兩個非負矩陣的乘積。如下圖所示，其中要求分解後的矩陣H和W都必須是非負矩陣。
img矩陣V分解為左矩陣W和右矩陣H，可理解為原始矩陣V的列向量是H中的所有列向量的加權和，對應的權重係數則是W的列向量的元素，所有H稱為基矩陣，W稱為係數矩陣。
NMF在人臉識別的應用中和PCA還有VQ分解不同。VQ分解是用一張完整的圖像直接代表源臉部圖像；PCA是將幾個完整人臉加減壓成一張臉；而NMF是取甲的眼睛，乙的鼻子，丙的嘴巴直接拼成一張臉，也就是說NMF分解後的基矩陣H是每張人臉的一個特徵部分，例如眼睛，鼻子，嘴巴等，然後權重矩陣對這些特徵進行加權處理從而得到一張完整的人臉。如下圖所示3種矩陣分解方式的區別。
imgimg**個人理解：**通過這2副圖粗略的講解，30個signatures就是第一幅圖的左邊的五官，第二幅圖右邊的數字，也就是基矩陣H，通過權重矩陣W對這些特徵進行加權處理從而得到一張完整的人臉，一個非負矩陣。
ssGSEAGSEA分析，jimmy老師在《生信技能樹》公眾號多次講解：
但實際上，絕大部分讀者並沒有去細看這個統計學原理，也不需要知道gsea分析的nes值如何計算 、
去卷積算法這個就更麻煩了，見這位仁兄的「打屁股理論」https://www.zhihu.com/question/48279880
image-20200418235612174通俗易懂，而CIBERSORT的反卷積算法就是，你的mrna就是信號系統的輸出（一個基因在樣本中的表達量），而lm22就是衝擊響應（細胞分數權重），可以反推信號系統的輸入（該基因在不同細胞亞群表達水平），也就是Jimmy老師文中所說的一個細胞在樣本中的表達量是該基因在不同細胞亞群表達水平與細胞分數權重的線性組合，
SVM釋義來源自https://www.jiqizhixin.com/articles/2018-10-17-20與https://tangshusen.me/2018/10/27/SVM/，SVM就是一種二類分類模型，他的基本模型是的定義在特徵空間上的間隔最大的線性分類器，SVM的學習策略就是間隔最大化。
我們憑直觀感受應該覺得答案是H3。首先H1不能把類別分開，這個分類器肯定是不行的；H2可以，但分割線與最近的數據點只有很小的間隔，如果測試數據有一些噪聲的話可能就會被H2錯誤分類(即對噪聲敏感、泛化能力弱)。H3以較大間隔將它們分開，這樣就能容忍測試數據的一些噪聲而正確分類，是一個泛化能力不錯的分類器。
但是很多情況是線性不可分，這時候就需要引入罰分機制C與鬆弛變量ξ；
如果是非線性，則需要通過kerenel將數據變化成高維空間，在新空間裡用線性方法從訓練數據中學習得到模型
還是不懂就看後面的函數輸出結果進行結合理解吧，反正絕大部分情況下，不懂原理似乎是並不會影響使用
R包e1071這裡用e1071的SVM函數進行運算，需要對其參數進行一些了解,這裡用詳細的中文解釋這個函數算法https://rpubs.com/skydome20/R-Note14-SVM-SVR http://blog.fukuball.com/lin-xuan-tian-jiao-shou-ji-qi-xue-xi-ji-fa-machine-learning-techniques-di-3-jiang-xue-xi-bi-ji/
svm(...
type = 決定svm是要用來classification(類別)、還是regression(連續)。
scale = 將資料正規化成(平均值, 標準差) = (0,1) 的分佈。
kernel = 將資料映射到特徵空間的kernel-fun，用來處理「非線性可分」的問題。
*cost = 在Lagrange formulation中的大C，決定給被誤差/分錯的資料「多少」懲罰值。
*epsilon = margin of tolerance。越大，表示在容忍範圍內的誤差/分錯的資料，不會被懲罰；反之，越接近0，每一個誤差/分錯的資料都會被懲罰。
*gamma = 在kernel-fun裡面的參數(linear-fun除外)。
...
)
# 更多詳情請參考R的官方Help文件以及原始的CIBERSORT R腳本 https://rdrr.io/github/singha53/amritr/src/R/supportFunc_cibersort.R
函數分解step1.載入數據這裡就直接使用Jimmy老師的第一步了，我覺得有困難的再加注釋
rm(list = ls())
options(stringsAsFactors = F)
# 首先讀取兩個文件
sig_matrix <-"LM22-ref.txt" # CIBERSORT 內置資料庫挖掘
mixture_file <- "mRNA2.txt" # 約80M，TCGA資料庫
# 兩個表達矩陣需要取交集
#read in data
X <- read.table(sig_matrix,header=T,sep="\t",row.names=1,check.names=F)
Y <- read.table(mixture_file, header=T, sep="\t", check.names=F)
Y <- Y[!duplicated(Y[,1]),]###去重複基因名
rownames(Y)<-Y[,1]
Y<-Y[,-1]
X <- data.matrix(X)###convert data as matrix
Y <- data.matrix(Y)
Y[1:4,1:4]##check data
X[1:4,1:4]
dim(X)
dim(Y)
X <- X[order(rownames(X)),]###行名字母排序
Y <- Y[order(rownames(Y)),]###行名字母排序
#anti-log if max < 50 in mixture file
if(max(Y) < 50) {Y <- 2^Y} ###如果Y矩陣中最大值<50，則變為2的y次方，也就是原始Y是被log2的
QN = F #QN = Quantile normalization of input mixture (default = TRUE)
#quantile normalization of mixture file
if(QN == TRUE){
tmpc <- colnames(Y)
tmpr <- rownames(Y)
Y <- normalize.quantiles(Y)#preprocessCore的函數，正態化數據
colnames(Y) <- tmpc
rownames(Y) <- tmpr
}
#intersect genes
Xgns <- row.names(X)
Ygns <- row.names(Y)
YintX <- Ygns %in% Xgns ###y中取x
Y <- Y[YintX,] ###取共有子集
XintY <- Xgns %in% row.names(Y)###x中取y
X <- X[XintY,]
dim(X)
dim(Y)
#standardize sig matrix
X <- (X - mean(X)) / sd(as.vector(X)) ###標準化數據
Y[1:4,1:4]
X[1:4,1:4]
boxplot(X[,1:4])
save(X,Y,file = 'input.Rdata')第一步就是取交集並將mrna矩陣zscore化
step2. 核心算法rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
Y[1:4,1:4]
X[1:4,1:4]
dim(X)
dim(Y)
# 下面的演示是為了搞清楚 CoreAlg 函數
# 並不需要保存任何信息
# 從表達矩陣Y裡面，隨機挑選LM22矩陣基因數量的表達量值
Ylist <- as.list(data.matrix(Y)) ###將Y矩陣中每一個數值作為list的一個元素
yr <- as.numeric(Ylist[ sample(length(Ylist),dim(X)[1]) ])###從Ylist隨機挑選nrow（x）個元素
# yr 這個時候是一個假設的樣本
#standardize mixture,就是scale 函數
yr <- (yr - mean(yr)) / sd(yr) ##標準化樣本
boxplot(yr)
# 每次隨機挑選的yr，都是需要走後面的流程
# 一切都是默認值的支持向量機
# 這裡的X是LM22矩陣，不同的免疫細胞比例組合成為不同的yr
# 這裡的yr是隨機的，反推免疫細胞比例
out=svm(X,yr)##e1071包的函數
out
out$SV ##SV=支持向量
# 需要修改的參數包括：type="nu-regression",kernel="linear",nu=nus,scale=F
###參數含義，採用nu-regression,nu回歸，採用線性方法從訓練數據中學習得到模型，nu參數為nus值，不scale
###SVM參數https://stats.stackexchange.com/questions/94118/difference-between-ep-svr-and-nu-svr-and-least-squares-svr
###https://blog.csdn.net/csqazwsxedc/article/details/52230092
svn_itor <- 3
y=yr
#try different values of nu
res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
model
}
#Execute In a parallel way the SVM
####Windows沒有辦法用mclapply開多核的，可以用parlapply
library(parallel)
if(Sys.info()['sysname'] == 'Windows') {
out <- mclapply(1:svn_itor, res, mc.cores=1)
}else {
out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
}
# 運行了Support Vector Machines，函數是 svm {e1071}
###windows開多核
library(parallel)
clnum<-detectCores()
cl <- makeCluster(getOption("cl.cores", clnum))
clusterExport(cl, c("X","y"),envir=environment())
clusterEvalQ(cl,library(e1071))
out <- parLapply(cl,1:svn_itor,res)
####
out
#Initiate two variables with 0
nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)
t <- 1
while(t <= svn_itor) {
# 得到兩個向量之間矩陣乘法的權重，此時應該只得到一個數字。
# 這樣做是乘以係數
# 支持向量是數據集的點，它們靠近分隔類別的平面
# 現在的問題是，我沒有任何類別（離散變量，例如「運動」、「電影」），但我有一個連續變量
mySupportVectors <- out[[t]]$SV ###不同nu參數的支持向量
# 係數定義
myCoefficients <- out[[t]]$coefs ###不同nu參數的係數
weights = t(myCoefficients) %*% mySupportVectors ####2個矩陣的乘積
# 設置權重和相關性
weights[which(weights<0)]<-0 ##小於0的乘積為0
w<-weights/sum(weights) ##相關性
# 根據對應的權重與參考集相乘
u <- sweep(X,MARGIN=2,w,'*') ###sweep類似於apply，多了一個STATS，代表是運算的參數
# 統計每行總和
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k - y)^2))) ###標準差
corrv[t] <- cor(k, y) ##相關性
t <- t + 1 ###t從1開始循環，直到t=3
}
#pick best model
rmses <- nusvm
corrv
mn <- which.min(rmses) ###去標準差最小的nu值為best model
mn
#[1] 1
model <- out[[mn]]
# 從nus為0.25,0.5,0.75的3個模型裡面挑選一個即可
#get and normalize coefficients
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
# w 就是計算後的22種免疫細胞的比例
w <- (q/sum(q))
mix_rmse <- rmses[mn]
mix_r <- corrv[mn]
# 會返回這個隨機的y的免疫細胞組成情況，就是權重w
newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
newList
# 根據對應的權重與參考集相乘
u <- sweep(X,MARGIN=2,w,'*')
k <- apply(u, 1, sum)
plot(y,k)
sqrt((mean((k - y)^2)))
cor(k, y)
# 通常這個預測結果都慘不忍睹
# 每次把表達矩陣通過去卷積拆分成為LM22的免疫細胞比例結果
# 並且包裝成為函數，如下：
#' CIBERSORT R script v1.03 (last updated 07-10-2015)
#' Note: Signature matrix construction is not currently available; use java version for full functionality.
#' Author: Aaron M. Newman, Stanford University (amnewman@stanford.edu)
#Core algorithm
CoreAlg <- function(X, y){
#try different values of nu
svn_itor <- 3
res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-e1071::svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
model
}
if(Sys.info()['sysname'] == 'Windows') out <- parallel::mclapply(1:svn_itor, res, mc.cores=1) else
out <- parallel::mclapply(1:svn_itor, res, mc.cores=svn_itor)
nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)
#do cibersort
t <- 1
while(t <= svn_itor) {
weights = t(out[[t]]$coefs) %*% out[[t]]$SV
weights[which(weights<0)]<-0
w<-weights/sum(weights)
u <- sweep(X,MARGIN=2,w,'*')
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k - y)^2)))
corrv[t] <- cor(k, y)
t <- t + 1
}
#pick best model
rmses <- nusvm
mn <- which.min(rmses)
model <- out[[mn]]
#get and normalize coefficients
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
w <- (q/sum(q))
mix_rmse <- rmses[mn]
mix_r <- corrv[mn]
newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}上面就是CIBERSORT函數的運算過程，不想了解算法過程的就當成1個函數來使用就好
step3. 運行CIBERSORTrm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
Y[1:4,1:4]
X[1:4,1:4]
dim(X)
dim(Y)
library(preprocessCore)
library(parallel)
library(e1071)
source("cibersort.R")
itor <- 1
Ylist <- as.list(data.matrix(Y))
dist <- matrix()
# 就是把 CoreAlg 函數運行1000次
perm=1000
while(itor <= perm){
print(itor) # 列印進度
#random mixture
yr <- as.numeric(Ylist[ sample(length(Ylist),dim(X)[1]) ])
#standardize mixture
yr <- (yr - mean(yr)) / sd(yr)
#run CIBERSORT core algorithm
result <- CoreAlg(X, yr)
mix_r <- result$mix_r
#store correlation
if(itor == 1) {dist <- mix_r}
else {dist <- rbind(dist, mix_r)}
itor <- itor + 1
}
####
newList <- list("dist" = dist)###獲取1000次計算相關係數
nulldist=sort(newList$dist) ###w值排序
# 這個nulldist 主要是用來計算P值
if(F){
P=perm
#empirical null distribution of correlation coefficients
if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)} ###取最小的非負p值
print(nulldist)
}
save(nulldist,file = 'nulldist_perm_1000.Rdata')
step.4 計算混合樣本對每個單樣本計算1000次，得到單個樣本的1000次計算中最小的非零p值，以及對應的權重，相關性與最小標準差
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
Y[1:4,1:4]
X[1:4,1:4]
dim(X)
dim(Y)
library(preprocessCore)
library(parallel)
library(e1071)
source("cibersort.R")
header <- c('Mixture',colnames(X),"P-value","Correlation","RMSE")
print(header)
# [1] "Mixture" "B cells naive"
# [3] "B cells memory" "Plasma cells"
# [5] "T cells CD8" "T cells CD4 naive"
# [7] "T cells CD4 memory resting" "T cells CD4 memory activated"
# [9] "T cells follicular helper" "T cells regulatory (Tregs)"
# [11] "T cells gamma delta" "NK cells resting"
# [13] "NK cells activated" "Monocytes"
# [15] "Macrophages M0" "Macrophages M1"
# [17] "Macrophages M2" "Dendritic cells resting"
# [19] "Dendritic cells activated" "Mast cells resting"
# [21] "Mast cells activated" "Eosinophils"
# [23] "Neutrophils" "P-value"
# [25] "Correlation" "RMSE"
load(file = 'nulldist_perm_1000.Rdata')
print(nulldist)
fivenum(print(nulldist))
#[1] -0.079400079 -0.009724867 0.019402797 0.056885990 0.604810742
#minimum, lower-hinge, median, upper-hinge, maximum
output <- matrix()
itor <- 1
mix <- dim(Y)[2] ###取mrna中的樣本
pval <- 9999
# 表達矩陣的每個樣本，都需要計算一下LM22的比例
#iterate through mix
while(itor <= mix){
##################################
## Analyze the first mixed sample
##################################
y <- Y[,itor]
#標準化樣本數據集
y <- (y - mean(y)) / sd(y)
#執行SVR核心算法
result <- CoreAlg(X, y)
#獲得結果
w <- result$w
mix_r <- result$mix_r
mix_rmse <- result$mix_rmse
#計算p-value
if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}
#輸出output
out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
if(itor == 1) {output <- out}
else {output <- rbind(output, out)}
itor <- itor + 1
}
head(output)
#save results
write.table(rbind(header,output), file="CIBERSORT-Results.txt", sep="\t", row.names=F, col.names=F, quote=F)
#return matrix object containing all results
obj <- rbind(header,output)
obj <- obj[,-1]
obj <- obj[-1,]
obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
rownames(obj) <- colnames(Y)
colnames(obj) <- c(colnames(X),"P-value","Correlation","RMSE")
obj
save(obj,file = 'output_obj.Rdata')
#return matrix object containing all results
obj <- rbind(header,output)
obj <- obj[,-1]
obj <- obj[-1,]
obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
rownames(obj) <- colnames(Y)
colnames(obj) <- c(colnames(X),"P-value","Correlation","RMSE")
obj[1:4,1:4]
# B cells naive B cells memory Plasma cells
# TCGA-DK-AA74-01A-11R-A39I-07 0.0125034 0 0.00000000
# TCGA-DK-A3IM-01A-11R-A20F-07 0.0000000 0 0.00000000
# TCGA-GU-A42P-01A-11R-A23W-07 0.1414380 0 0.03075805
# TCGA-4Z-AA7W-01A-11R-A39I-07 0.0000000 0 0.03125399
# T cells CD8
# TCGA-DK-AA74-01A-11R-A39I-07 0.11688532
# TCGA-DK-A3IM-01A-11R-A20F-07 0.04647556
# TCGA-GU-A42P-01A-11R-A23W-07 0.02934341
# TCGA-4Z-AA7W-01A-11R-A39I-07 0.37869645
save(obj,file = 'output_obj.Rdata')
step5. 結果可視化最後一步就是結果可視化了，這裡Jimmy老師用了3種可視化方法：熱圖，柱狀圖以及箱圖
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
Y[1:4,1:4]
X[1:4,1:4]
dim(X)
dim(Y)
library(preprocessCore)
library(parallel)
library(e1071)
load(file = 'output_obj.Rdata')
# Step3:將CIBERSORT_Result挑選整理，去除沒有差異表達的細胞。
library(dplyr)
library(tidyr)
library(tidyverse)
cibersort_raw <- read.table("CIBERSORT-Results.txt",header = T,sep = '\t') %>%
rename("Patients" = "Mixture") %>%
select(-c("P.value","Correlation","RMSE"))
# 通過管道符一步步先將CIBERSORT_Results讀入R語言中，並將其第一列列名「Mixture」修改為「Patiens」並去除了後三列。
#並賦值給cibersort_raw。
cibersort_tidy <- cibersort_raw %>%
remove_rownames() %>%
column_to_rownames("Patients")
# 將cibersort_raw第一列變為列名後賦值給cibersort_tidy。
flag <- apply(cibersort_tidy,2,function(x) sum(x == 0) <
dim(cibersort_tidy)[1]/2)
# 篩選出0值超過樣本的一半的一些細胞
cibersort_tidy <- cibersort_tidy[,which(flag)] %>%
as.matrix() %>%
t()
# 留下在大部分樣本中有所表達的細胞。
bk <- c(seq(0,0.2,by = 0.01),seq(0.21,0.85,by=0.01))
# breaks用來定義數值和顏色的對應關係。
# Step4:將CIBERSORT_Result進行可視化
#1）熱圖
library(pheatmap)
library(RColorBrewer)
pheatmap(
cibersort_tidy,
breaks = bk,
cluster_cols = T,
scale = "row",
cluster_row = T,
border_color = NA,
show_colnames = F,
show_rownames = T,
color = c(colorRampPalette(colors = c("blue","white"))(length(bk)/2),
colorRampPalette(colors = c("white","red"))(length(bk)/2)
))
#調整參數讓熱圖更加美觀。
#柱狀圖可視化細胞佔比預測
library(RColorBrewer)
mypalette <- colorRampPalette(brewer.pal(8,"Set1"))###8種顏色
cibersort_barplot <- cibersort_raw %>%
gather(key = Cell_type,value = Proportion,2:23)
#使用RColorBrewer包配置需要的色彩方案，使用gather函數中的key-value對應關係重建細胞名稱和比例的對應關係並賦值給cibersort_barplot
#cibersort_barplot$Patient1 <- factor(cibersort_barplot$Patient,
# levels = str_sort(unique(cibersort_barplot$Patient),
# numeric = T))
ggplot(cibersort_barplot,aes(Patients,Proportion,fill = Cell_type)) +
geom_bar(position = "stack",stat = "identity") +
labs(fill = "Cell Type",x = "",y = "Estiamted Proportion") + theme_bw() +
theme(axis.text.x = element_blank()) + theme(axis.ticks.x = element_blank()) +
scale_y_continuous(expand = c(0.01,0)) +
scale_fill_manual(values = mypalette(23))
#調整參數讓柱狀圖更加美觀。
#直觀箱線圖
ggplot(cibersort_barplot,aes(Cell_type,Proportion,fill = Cell_type)) +
geom_boxplot(outlier.shape = 21,coulour = "black") + theme_bw() +
labs(x = "", y = "Estimated Proportion") +
theme(axis.text.x = element_blank()) + theme(axis.ticks.x = element_blank()) +
scale_fill_manual(values = mypalette(23))
#調整參數讓柱狀圖更加美觀。這裡有個比較新奇的函數是library(RColorBrewer)，是樣本配色更加好看
mypalette <- colorRampPalette(brewer.pal(8,"Set1"))
#這裡表示按照Set1的配色的前8種
mypalette_1<-brewer.pal(8,"Set1")
image(1:8,1,as.matrix(1:8),col=mypalette_1,xlab="Greens (sequential)",
ylab="",xaxt="n",yaxt="n",bty="n")
set1裡面配色最多9種，如果設置為10的話mypalette_1<-brewer.pal(10,"Set1")
image(1:10,1,as.matrix(1:10),col=mypalette_1,xlab="Greens (sequential)",
ylab="",xaxt="n",yaxt="n",bty="n")會發現還是只有9種
image-20200419210132730最後按照Jimmy老師的可視化函數的結果如下
熱圖
柱狀圖
小結CIBERSORT就是將已知的LM22作為參考組，與混合樣本表達矩陣取交集後隨機從合樣本表達矩陣中抽樣對抽樣結果運用SVM算法的nu regression+line kernel，採取不同nu值獲得最佳的model。通過model運算1000次得到隨機樣本的權重W，相關係數r與標準差rmse，並將結果進行運算得到p值將1000次的p值進行排序，採用最小非零p值的結果作為預測結果
箱型圖