概要

本稿では、メタボロームデータ解析において統計学的指標に基づいて代謝物を探索・選定するために開発された手法 PLS-ROG（PLS-rank order of groups） のアプリケーションを紹介する。

原著論文

Yamamoto, H. PLS-ROG: Partial least squares with rank order of groups. J. Chemom. 2017, 31, e2883.

PLS-ROG の概要

PLS ROG（Partial least squares with rank order of groups） は、次元削減手法の 1 つである naïve PLS に、パラメータを 1 つ加えて発展させたものである。

PLS とは

PLS（Partial least squares） は群情報を用いた次元削減手法であり、メタボロームデータ解析によく用いられる。PLS には様々な種類があり、次元削減時に deflation を行うかどうかで得意な分析が異なる。一方で群情報を用いない次元削減手法として PCA（Principal Component Analysis） が知られている。メタボローム解析では、まず PCA を実施して PCA スコアに群間の差が現れなかった際に実施することが多い。（補足1）

Deflation とは

次元削減法では、第 1 成分から第 2 成分、第 3 成分、と順番に成分を決める。ここで n 番目の成分を求めるとき、 n-1 番目までの成分で表現された部分を差し引く操作を deflation という。Deflation を実施する場合、全種類の代謝物の情報をもつ成分は第 1 成分だけとなる。（補足1）（補足2）

PLS-R と PLS-DA

PLS-R（PLS-regression） と PLS-DA （PLS-discrimination analysis） は次元削減時に deflation を行う。これらの手法は第 2 成分以降に代謝物全種類の情報を含まないため、データの可視化や、群の区別に重要な代謝物の抽出にはあまり適さない。代謝物の抽出には VIP（variable importance for prediction） が用いられることがある。VIP は群の区別における重要度を表す指標で、1.0 以上の代謝物は群の区別に重要であると経験則から述べられている。
これらの手法は判別モデルの作成に用いるのがよい。（補足3）

naïve PLS と PLS-ROG

次元削減時に deflation を行わないため、データの可視化や代謝物の抽出に適する。代謝物の抽出では、PCA とよく似た因子負荷量（factor loading） を用いた統計学的仮説検定が可能である。

naïve PLS

メタボロームデータ（説明変数）を X とおき、群情報（目的変数）を Y とおく。naïve PLS では、固有ベクトル $w_x$ 、$w_y$ を用いて潜在説明変数 T と潜在目的変数 S を以下のようにおく。

T = Xw_x \\
S = Yw_y

固有ベクトルには以下のように大きさが 1 になるという制約条件がある。

w'_xw_x = 1 \\
w'_yw_y = 1

ここで T と S の共分散 $cov(S,T)$ が最大になるように次元を削減する。T と S を autoscaling （平均 0、分散 1 への変換）したとき、 $w_x$ は S と $X_p$ （p は代謝物の種類）の相関係数 $corr(S,Xp)$ に比例することが原著論文で示されている。この相関係数は、群情報と関連度の高い代謝物ほど 1 または -1 に近い。このように代謝物の抽出に用いられる固有ベクトルや相関係数は因子負荷量（Factor loading） と呼ばれる（補足1）。さらに相関係数から p 値を求めることができる。naïve PLS の特徴として、p 値という統計学的な指標に基づいて代謝物を探索・選定できることが挙げられる。

PLS-ROG

PLS-ROG は $w_y$ の制約条件が naïve PLS と異なり以下のようになる。詳細は原著論文や開発者の Github を参考にするとよい。
$$
(1-κ)w'_yw_y+κY'P'D'DPYw'_yw_y = 1　　（0≦κ≦1）
$$

少々複雑に見えるが、$Y'P'D'DPY $ は S の平滑化項と考えるとよい。つまり平滑化パラメータ κ を 1 に近づけるほど S が平均化処理される。これに伴って T も変化し、第 1 成分方向に整列するようになる。よって、3 群以上の場合も T の第 1 成分で群を順番付けられるようになるため、PLS with rank order of groups と呼ばれている。

PLS ROG App の紹介

今回は開発者の Github を参考に、平滑化パラメータ κ を変化させながら PLS-ROG を実行できる web アプリケーション PLS ROG App を Shiny from R studio で開発した。ソースコードは Github にも公開している。データをアップロードしたくない場合は、アプリケーションをダウンロードしてローカル環境で実行するとよい。

web アプリに繋いでみよう

こちらから PLS App に接続できる。

Choose TSV File

手元にあるメタボロームデータをアップロードしよう。データ形式は GitHub を参考に、以下のように作成するとよい。

	代謝物A	代謝物B	代謝物C	グループ
サンプルA	30	5	3	グループA
サンプルB	20	7	0	グループA
サンプルC	23	4	0	グループA
サンプルD	24	5	1	グループB
サンプルE	19	8	2	グループB
サンプルF	12	5	0	グループB

Normalization method

PLS-ROG はメタボロームデータ解析のために開発された手法だが、菌叢データの探索にも応用できることがある。菌叢データは組成データに変換して扱われるため、その場合は composition を選択する。

Analytics method

デフォルトでは naïve PLS が選択されている。必要に応じて PLS ROG に変更しよう。

Confidential Interval

確率楕円の信頼区間を設定する。デフォルトでは 0.95 が代入されており、多変量 t 分布における 95% 信頼区間が描画される。

Kappa （smoothing parameter）

PLS-ROG を選択すると出現する。1 に近づけるほど平均化処理される。

Axis1 と Axis2

描画したい成分を選択する。PLS-ROG では(群の数-1)だけ成分が計算されるため、横軸におきたい成分を Axis1 に設定し、縦軸におきたい成分を Axis2 に設定する。デフォルトでは第一成分が横軸に、第二成分が縦軸に設定される。

デモデータで実践

本稿では IBD multi'omics database に公開されている糞便のメタボロームデータをデモデータとして利用する。このデータセットには、クローン病（CD）、潰瘍性大腸炎（UC）、非炎症性腸疾患（nonIBD）の患者から採取されたサンプルがそれぞれ 265 検体、146 検体、135 検体含まれており、こちらの data/ からダウンロードすることもできる。

潜在説明変数 T

まずは naïve PLS を実施し、T の分布をみてみよう。メタボロームの構成と群情報が似ているサンプルほど近くに配置される。

カーソルをプロットに当てるとサンプルの情報が表示される。

では PLS-ROG を実施してみよう。Kappa = 0.90 に設定すると、第1成分 V1 方向に UC、nonIBD、CD が整列する。つまり V1 は UC、nonIBD、CD の違いを表す成分である可能性がある。

潜在目的変数 S

次に潜在目的変数 S の分布をみてみよう。naïve PLS を実施すると次のようになる。

ここで PLS-ROG を Kappa = 0.90 で実施すると平均化されることがわかる。

固有ベクトル Weight X

ここからは PLS-ROG の因子負荷量をみていこう。まずは各代謝物の固有ベクトルを調べると、uridine が第 1 成分 V1 と第 2 成分 V2 の両方で正の重みがつけられていることがわかる。

では生データで uridine の分布を調べてみると、CD、nonIBD、UC のそれぞれで異なる傾向がみられることがわかる。このように固有ベクトルは、代謝物の探索・選定に使える因子負荷量といえる。

相関係数 R

固有ベクトルで顕著な値がみられた uridine の相関係数 R を調べてみよう。第 1 成分の正負こそ逆転しているものの、他の代謝物と比べて顕著な値となっていることがわかる。

Q 値

最後に、相関係数から求められる P 値を Benjamini-Hochberg 法で補正した Q 値を調べてみよう。uridine の Q 値は第 1 成分 V1 と第 2 成分 V2 の両方で 0.05 以下となっており、統計学的にも有意であることがわかる。

菌叢データで実践

IBD multi'omics database には腸生検から採取された菌叢データも登録されている。同様にこちらの data/ からダウンロードできる。Normalization method で composition を選択していることに注意しよう。

まとめ

naïve PLS および PLS-ROG は代謝物の探索や可視化に有効である。
菌叢データや他のオミクスデータに対しても有効である可能性があるが、前処理の仕方などで要検討である。

Source Codes

R コンソールで print(source('app.R')) と実行するとアプリケーションが起動する。

app.R

app.R

rm( list=ls(all=TRUE) ) 
library(shiny)
library(tidytidbits) #isTruthy() に必要
library(ggplot2)
library(plotly)

# Define UI for dataset viewer app ----
ui <- fluidPage(
  navbarPage(
    title = "PLS ROG App",
    tabPanel("Analysis results",
             sidebarLayout(
               sidebarPanel(
                 width = 4,
                 uiOutput("tab"),
                 fileInput("RawData", "Choose TSV File", accept = ".tsv"),
                 selectInput(
                   inputId = "normalization",
                   label = "Normalization method",
                   choices = c("none","composition"),
                   selected = NULL
                 ),
                 selectInput(
                   inputId = "select",
                   label = "Analytics method",
                   choices = c("naïve pls","pls-rog"),
                   selected = NULL
                 ),
                 sliderInput(inputId = "level",
                             label = "Confidence interval",
                             min = 0.01,
                             max = 0.99,
                             value = 0.95),
                 conditionalPanel(
                   condition = "input.select == 'pls-rog'",
                   sliderInput(inputId = "kappa",
                               label = "Kappa (smoothing parameter)",
                               min = 0.00,
                               max = 0.99,
                               value = 0.00)
                 ),
                 selectInput(inputId = "Axis1", 
                             label = "Axes1", 
                             choices = ""),
                 selectInput(inputId = "Axis2", 
                             label = "Axes2", 
                             choices = ""),
               ),
               mainPanel(
                 tabsetPanel(type = "tabs",
                             tabPanel("T",
                                      plotlyOutput(
                                        outputId = "score_X_plots",
                                        width = "700px",
                                        height = "600px",
                                        inline = FALSE
                                      )
                             ),
                             tabPanel("S",
                                      plotlyOutput(
                                        outputId = "score_Y_plots",
                                        width = "700px",
                                        height = "600px",
                                        inline = FALSE
                                      )
                             ),
                             tabPanel("Weight_X",
                                      plotlyOutput(
                                        outputId = "Wx_plots",
                                        width = "600px",
                                        height = "600px",
                                        inline = FALSE
                                      )
                             ),
                             tabPanel("R",
                                      plotlyOutput(
                                        outputId = "R_plots",
                                        width = "600px",
                                        height = "600px",
                                        inline = FALSE
                                      )
                             ),
                             tabPanel("Q",
                                      plotlyOutput(
                                        outputId = "Q_plots",
                                        width = "600px",
                                        height = "600px",
                                        inline = FALSE
                                      )
                             )
                 )
               )
             )
    ),
    tabPanel("Raw Data",
             sidebarLayout(
               sidebarPanel(
                 width = 4,
                 selectInput(inputId = "feature", 
                             label = "select feature", 
                             choices = ""
                 )
               ),
               mainPanel(
                 plotlyOutput(
                   outputId = "boxplot",
                   width = "600px",
                   height = "600px",
                   inline = FALSE
                 ),
               )
             )
    )
  )
)

# Define server logic to summarize and view selected dataset ----
server <- function(input, output, session) {
  url <- a("here", href="https://github.com/Keisuke-H-Ota/PLS_App")
  output$tab <- renderUI({
    tagList("Refer input data format", url)
  })
  observeEvent(input$RawData, {
    # データの読み込み
    mbx <- read.table(file=input$RawData$datapath, sep="\t",header=F, quote='', comment.char="")
    nr <- nrow(mbx) # Number of rows
    nc <- ncol(mbx) # Numbers of columns
    sn <- mbx[2:nr,1] # Sample Name
    fn <- as.vector(t(mbx)[2:(nc-1),1]) # Feature Name
    gr <- mbx[2:nr,nc] # Group
    DF <- mbx[2:nr,2:(nc-1)]
    DF <- data.frame(sapply(DF, function(x) as.numeric(as.character(x))))
    
    # サンプルごとに欠損値を最小値の半分で置換
    for(i in 1:nrow(DF)){
      m <- min(DF[i,],na.rm=TRUE)
      na_index <- which(is.na(DF[i,]))
      DF[i,na_index] <- m/2
    }
    
    # 全サンプルで 0 以下の列を削除する。
    nanidx <- c()
    for(i in 1:ncol(DF)){
      if(sum(DF[,i]) <= 0){
        nanidx <- c(nanidx,i)
      }
    }
    if(length(nanidx)>0){
      DF <- DF[,-nanidx]
      fn <- fn[-nanidx]
    }
    
    # 全 feature で 0 の行を除く
    nanidx <- c()
    for(i in 1:nrow(DF)){
      if(sum(DF[i,])<=0){
        nanidx <- c(nanidx,i)
      }
    }
    if(length(nanidx)>0){
      DF <- DF[-nanidx,]
      sn <- sn[-nanidx]
      gr <- gr[-nanidx]
    }
    
    
    # 前処理
    df <- reactive({
      if(input$normalization=='none'){
        tmp = as.data.frame(DF)
        colnames(tmp) = fn
        return(tmp)
      }else if(input$normalization=='composition'){
        tmp2 = as.data.frame(DF/rowSums(DF))
        colnames(tmp2) = fn
        return(tmp2)}
    })
    
    # 関数の読み込み
    source('plsrog.R')
    
    # pls-rog の実行
    res <- reactive({
      if(input$select=='naïve pls'){
        return(pls(df(), gr))}else{
          return(plsrog(df(), gr, input$kappa))
        }
    })
    
    axis1 <- reactive({
      score_X <- cbind(as.data.frame(res()[[1]]))
      if(!isTruthy(input$Axis1)){
        return(colnames(score_X)[1])
      }else{
        return(input$Axis1)
      }
    })
    
    axis2 <- reactive({
      score_X <- cbind(as.data.frame(res()[[1]]))
      if(!isTruthy(input$Axis2)){
        return(colnames(score_X)[2])
      }else{
        return(input$Axis2)
      }
    })
    
    output$score_X_plots <- renderPlotly({
      score_X <- cbind(as.data.frame(res()[[1]]),sn,gr)
      ng <- ncol(score_X)-2
      observe({
        updateSelectInput(session = session, 
                          inputId = "Axis1", 
                          choices = colnames(score_X)[1:ng],
                          selected = input$Axis1)
      })
      observe({
        updateSelectInput(session = session, 
                          inputId = "Axis2", 
                          choices = colnames(score_X)[1:ng],
                          selected = input$Axis2)
      })
      data = data.frame(x = score_X[,axis1()], y = score_X[,axis2()])
      plot <- ggplot(data, aes(x = x, y = y, color = gr)) + 
        geom_point(aes(color = gr, text = paste("Sample Name:",sn)),alpha = 1.0, size = 1) +
        labs(x = 'Axis1', y = 'Axis2', color = "diagnosis") +
        stat_ellipse(aes(color = gr), type = "t", level = input$level)
      plot
    })
    
    output$score_Y_plots <- renderPlotly({
      score_Y <- cbind(as.data.frame(res()[[2]]),gr)
      data = data.frame(x = score_Y[,axis1()], y = score_Y[,axis2()], color = gr)
      plot <- ggplot(data, aes(x = x, y = y)) + 
        geom_point(aes(color = gr),alpha = 1, size = 2) +
        labs(x = axis1(), y = axis2())
      plot
    })
    
    output$Wx_plots <- renderPlotly({
      Wx <- cbind(as.data.frame(res()[[3]]),fn)
      data = data.frame(x = Wx[,axis1()], y = Wx[,axis2()], fn = fn)
      plot <- ggplot(data, aes(x = x, y = y, fn = fn)) + 
        geom_point(alpha = 1, size = 1) +
        labs(x = axis1(), y = axis2())
      plot
    })
    
    output$R_plots <- renderPlotly({
      R <- res()[[5]]
      colnames(R) <- NULL
      R <- cbind(as.data.frame(R),fn)
      data = data.frame(x = R[,axis1()], y = R[,axis2()], fn = fn)
      plot <- ggplot(data, aes(x = x, y = y, fn = fn)) + 
        geom_point(alpha = 1, size = 1) +
        labs(x = axis1(), y = axis2())
      plot
    })
    
    output$Q_plots <- renderPlotly({
      Q <- res()[[7]]
      colnames(Q) <- NULL
      Q <- cbind(as.data.frame(Q),fn)
      data = data.frame(x = Q[,axis1()], y = Q[,axis2()], fn = fn)
      plot <- ggplot(data, aes(x = x, y = y, fn = fn)) + 
        geom_point(alpha = 1, size = 1) +
        labs(x = axis1(), y = axis2())
      plot
    })
    
    feature <- reactive({
      if(!isTruthy(input$feature)){
        return(fn[1])
      }else{
        return(input$feature)
      }
    })
    
    output$boxplot <- renderPlotly({
      observe({
        updateSelectInput(session = session, 
                          inputId = "feature", 
                          choices = fn,
                          selected = input$feature)
      })
      data = cbind(data.frame(y = df()[,feature()]),gr)
      plot <- ggplot(data, aes(y = y, x = gr)) + 
        geom_boxplot() +
        labs(y = feature(), x = gr)
      plot
    })
  }
  )
}

# Create Shiny app ----
shinyApp(ui, server)

plsrog.R

plsrog.R

plsrog <- function(X,class, kappa){
  
  # data matrix
  X <- as.matrix(X)
  X <- matrix(as.numeric(X),nrow=nrow(X)) # metabolites*samples
  
  # response variable
  Y0 <- factor(class)
  Y <- model.matrix(~ Y0 + 0)
  
  # penalized matrix
  P <- NULL
  p <- colSums(Y)
  for(i in 1:ncol(Y)){
    P <- cbind(P,Y[,i]/p[i])
  }
  P <- t(P)
  
  # differential matrix
  g <- ncol(Y)
  D <- diff(diag(1,g))
    
  # autoscaling
  X <- scale(X)
  Y <- scale(Y,scale=FALSE)
  
  # sample size-1
  N <- nrow(X)-1
  
  # smoothing parameter
  C <- kappa*t(Y)%*%t(P)%*%t(D)%*%D%*%P%*%Y+(1-kappa)*diag(1,g)
  
  # cholesky decomposition
  Rx <- chol(solve(C))
  Ry <- chol(C)
  
  # singular value decomposition
  USVx <- svd(Rx%*%t(Y)%*%X/N)
  USVy <- svd(t(X)%*%Y%*%solve(Ry)/N)
  
  # weght vector (factor loading)
  Wx <- USVx$v
  Wy <- solve(Ry)%*%USVy$v
  
  # PLS score
  T <- X%*%Wx
  S <- Y%*%Wy
    
  # correlation coefficient
  R <- NULL
  for(i in 1:(ncol(Y-1))){
      lambdax <- cov(T[,i],S[,i])
      r <- (sqrt(N)*lambdax*Wx[,i])/as.numeric(sqrt(t(Wy[,i])%*%t(Y)%*%Y%*%Wy[,i]))
      r[is.nan(r)] <- 0
      R <- cbind(R,r)
  }
  
  # statistical test
  P <- NULL
  for(i in 1:(ncol(Y-1))){
    p <- 2*pt(abs(R[,i])*sqrt(nrow(X)-2)/sqrt(1-R[,i]^2), nrow(X)-2, lower.tail=FALSE)
    P <- cbind(P,p)
  }
  
  # Q value
  Q <- NULL
  for(i in 1:(ncol(Y-1))){
    q <- p.adjust(P[,i], method = "BH")
    Q <- cbind(Q,q)
  }
  
  all <- list(T,S,Wx,Wy,R,P,Q)
}

pls <- function(X,class){
  
  # data matrix
  X <- as.matrix(X)
  X <- matrix(as.numeric(X),nrow=nrow(X)) # metabolites*samples

  # response variable
  Y0 <- factor(class)
  Y <- model.matrix(~ Y0 + 0)
  
  # penalized matrix
  P <- NULL
  p <- colSums(Y)
  for(i in 1:ncol(Y)){
    P <- cbind(P,Y[,i]/p[i])
  }
  P <- t(P)
  
  # autoscaling
  X <- scale(X)
  Y <- scale(Y,scale=FALSE)
  
  # ----------------
  #   ordinary PLS
  # ----------------
  
  # (sample size)-1
  N <- nrow(X)-1
  
  # singular value decomposition
  USVx <- svd(t(Y)%*%X/N)
  USVy <- svd(t(X)%*%Y/N)
  
  # weight vector matrix
  Wx <- USVx$v
  Wy <- USVy$v
  
  # score matrix
  T <- X%*%Wx
  S <- Y%*%Wy
  
  # correlation coefficient
  R <- NULL
  for(i in 1:(ncol(Y-1))){
    lambdax <- cov(T[,i],S[,i])
    options(warn=-1)
    r <- (sqrt(N)*lambdax*Wx[,i])/(sqrt(t(Wy[,i])%*%t(Y)%*%Y%*%Wy[,i]))
    r[is.nan(r)] <- 0
    R <- cbind(R,r)
  }
  
  # statistical test
  P <- NULL
  for(i in 1:(ncol(Y-1))){
    p <- 2*pt(abs(R[,i])*sqrt(nrow(X)-2)/sqrt(1-R[,i]^2), nrow(X)-2, lower.tail=FALSE)
    P <- cbind(P,p)
  }
  
  # Q value
  Q <- NULL
  for(i in 1:(ncol(Y-1))){
    q <- p.adjust(P[,i], method = "BH")
    Q <- cbind(Q,q)
  }
  
  all <- list(T,S,Wx,Wy,R,P,Q)
}