Go to Qiita Advent Calendar 2024 Top
LoginSignup
10
19

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@keisuke-ota(太田 圭祐)

MetaboanalystR を用いたメタボロームデータ解析

Last updated at Posted at 2020-10-17

概要

本稿では、Metaboanalyst の R パッケージ を用いたメタボロームデータ解析の流れを紹介する。

目的

本稿では、IBD multiomics database に登録されているメタボロームデータの一部をデモデータとして用い、CD (クローン病患者) 群と nonIBD (非炎症性腸疾患患者) 群の糞便中に含まれる菌叢代謝物を比較する。今回の代謝物は LC-MS スペクトルから定量されたが、NMR や MS で測定した場合も本稿の解析が可能である。

R コマンド

ライブラリとデータの読み込み

まず MetaboAnalystR を読み込む。

Rコンソール
library("MetaboAnalystR")

次に mSetオブジェクトを作る。MetaboanalystR では、mSet に次々と属性を加えていくことで解析が進行する。

Rコンソール
mSet<-InitDataObjects("pktable", "stat", FALSE)

メタボロームデータを csv ファイルで読み込む。

Rコンソール
mSet<-Read.TextData(mSet, "MetabolomeData.csv", "colu", "disc");

csv ファイルのフォーマットは以下のとおり。
一行目に被験者の ID、二行目に被験者の診断名、三行目以降に代謝物の LC-MS スペクトル強度が記載されている。
詳しくは 公式HP を参考にすると良い。

MetabolomeData.csv
Metabolite,MSM9VZHX,MSM9VZJ3,MSMA26AL,PSMA265T,MSM9VZHF,HSMA33JH
diagnosis,nonIBD,CD,CD,CD,nonIBD,CD
1-3-7-trimethylurate,264995.0,7268.0,165010.0,183495.0,4448779.0,171060.0
1-methylguanine,81422.0,156158.0,103739.0,53031.0,126482.0,50352.0
1-methylguanosine,8369.0,17360.0,9266.0,2041.0,6577.0,1725.0
1-methylhistamine,118680.0,178849.0,1271.0,0.0,12265.0,569462.0
1-methylnicotinamide,505048.0,212683.0,1786.0,0.0,4796.0,82993.0

欠損値のフィルタリングと推定

RemoveMissingPercent では、全サンプルの何 % で欠損値だった代謝物をフィルタリングするのかを決める。以下の例では 50% 以上のサンプルで取れてこなかった代謝物が除去される。
ImputeVarでは、欠損値を推定する。method ="min" の場合、全サンプルの最低値に置換される。

Rコンソール
mSet<-SanityCheckData(mSet) # データの状態をチェック
mSet<-RemoveMissingPercent(mSet, percent=0.5)
mSet<-ImputeVar(mSet, method="min")
mSet<-PreparePrenormData(mSet)

サンプルの標準化

SumNorm（組成比変換）、LogNorm（対数変換）、AutoNorm（平均値との差を標準偏差で割る） という3種類のオプションがある。標準化の方法によって解析結果は異なるため、注意深く検討する必要がある。

Rコンソール
mSet<-Normalization(mSet, "SumNorm", "LogNorm", "AutoNorm", ratio=FALSE, ratioNum=20)
mSet<-PlotNormSummary(mSet, "norm_0_", "png", 72, width=NA)
mSet<-PlotSampleNormSummary(mSet, "snorm_0_", "png", 72, width=NA)

標準化結果は、norm_0_dpi72.png と snorm_0_dpi72.png にされる。標準化後の分布が正規分布に近いほど、t 検定の妥当性が上がる。

norm_0_dpi72.png

PCA（主成分分析）

以下のコマンドで PCA を実施する。

Rコンソール
mSet<-PCA.Anal(mSet)
mSet<-PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet<-PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet<-PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA, 1,2,0.95,0,0)
mSet<-PlotPCALoading(mSet, "pca_loading_0_", "png", 72, width=NA, 1,2);
mSet<-PlotPCABiplot(mSet, "pca_biplot_0_", "png", 72, width=NA, 1,2)
mSet<-PlotPCA3DLoading(mSet, "pca_loading3d_0_", "json", 1,2,3)

結果は、pca_pair_0_dpi72.png に表示される。
pca_pair_0_dpi72.png

各軸の寄与度は、pca_scree_0_dpi72.png の青線に表示される。緑色は累積寄与度を表す。
pca_scree_0_dpi72.png

PLS および PLS-DA

PLS は、PCA で代謝物の違いが見られなかった場合に実施されることがある。

Rコンソール
# PLS
mSet<-PLSR.Anal(mSet, reg=TRUE)
mSet<-PlotPLSPairSummary(mSet, "pls_pair_0_", "png", 72, width=NA, 5)
mSet<-PlotPLS2DScore(mSet, "pls_score2d_0_", "png", 72, width=NA, 1,2,0.95,0,0)
mSet<-PlotPLS3DScoreImg(mSet, "pls_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
mSet<-PlotPLSLoading(mSet, "pls_loading_0_", "png", 72, width=NA, 1, 2);
mSet<-PlotPLS3DLoading(mSet, "pls_loading3d_0_", "json", 1,2,3)
# PLS-DA
mSet<-PLSDA.CV(mSet, "L",5, "Q2")
mSet<-PlotPLS.Classification(mSet, "pls_cv_0_", "png", 72, width=NA)
mSet<-PlotPLS.Imp(mSet, "pls_imp_0_", "png", 72, width=NA, "vip", "Comp. 1", 15,FALSE)

PLS の結果は、pls_score2d_0_dpi72.png や pls_score3d_0_dpi72.png に表示される。各軸は、CD 群か nonIBD 群かという群情報からの寄与も受けるため、PCA よりも分類されやすい。

pls_score2d_0_dpi72.png

pls_score3d_0_dpi72.png

PLS-DA では、CD 群か nonIBD 群かを区別する際の重要度を表す VIP score が計算される。pls_imp_0_dpi72.png から、どちらの群でどのような代謝物が増加あるいは減少しているのかが分かる。
pls_imp_0_dpi72.png

Fold change

Fold change (FC) は、2 群間で代謝物を比較した際の、差の大きさを表す。

Rコンソール
mSet<-FC.Anal.unpaired(mSet, 2.0, 0)
mSet<-PlotFC(mSet, "fc_0_", "png", 72, width=NA)

結果は fc_0_dpi72に表示される。 $Log_2(FC)$ が 1 以上の代謝物はピンク色でプロットされる。

fc_0_dpi72.png

t 検定

標準化した結果、正規分布に近づけば t 検定を実施することができる。

Rコンソール
mSet<-Ttests.Anal(mSet, F, 0.05, FALSE, TRUE, FALSE)
mSet<-PlotTT(mSet, "tt_0_", "png", 72, width=NA)

結果は tt_0_dpi72.png に表示される。

tt_0_dpi72.png

Random Forest 分析

PLS-DA と同様に、CD 群と nonIBD 群を区別する際に重要な代謝物を選定することができる。

Rコンソール
mSet<-RF.Anal(mSet, 500,7,1)
mSet<-PlotRF.Classify(mSet, "rf_cls_0_", "png", 72, width=NA)
mSet<-PlotRF.VIP(mSet, "rf_imp_0_", "png", 72, width=NA)
mSet<-PlotRF.Outlier(mSet, "rf_outlier_0_", "png", 72, width=NA)

結果はrf_imp_0_dpi72.pngに表示される。MeanDecreaseAccuracyは、その代謝物を除いたときに、CD 群と nonIBD 群の分類精度がどれくらい落ちるのかを表す指標である。VIP score と同様に、この値が大きいほど 2 群に分類する際に重要な代謝物といえる。
rf_imp_0_dpi72.png

Streamlit を用いた欠損値処理と標準化条件の検討

メタボローム解析の結果は、欠損値の処理方法および標準化条件によって異なることがある。
ここでは、欠損値処理と標準化条件の検討方法を提案する。

網羅的解析の実施

以下の script を実行し、欠損値処理と標準化を変えながら網羅的に解析する。
具体的には、RemoveMissingPercent に 0.1, 0.5, 0.9 のいずれかを代入し、Normalization に SumNorm, LogNorm, AutoNorm または NULL を代入する。
今回の場合、RemoveMissingPercent で 3 通り、Normalization で $2^3$ 通り、合計 24 通りのパラメータで解析する。
実行すると、各代謝物に対して Fold change (FC) 、 Q value (t検定)、 VIP score (PLS-DA)、MDA (Random Forest 解析) が計算される。
成果物は results.csv に出力される。

Metaboanalyst.R
library("MetaboAnalystR")
rm( list=ls(all=TRUE) )

df <- data.frame()

file <- paste0('MetabolomeData.csv')

for (p in c(0.1, 0.5, 0.9)){
	for (n1 in c('NULL','SumNorm')){
		for (n2 in c('NULL','LogNorm')){	
			for (n3 in c('NULL','AutoNorm')){
				mSet<-InitDataObjects("pktable", "stat", FALSE)
				mSet<-Read.TextData(mSet, file, "colu", "disc");
				mSet<-SanityCheckData(mSet)

				mSet<-RemoveMissingPercent(mSet, percent=p)
				mSet<-ImputeVar(mSet, method="min")
				mSet<-PreparePrenormData(mSet)

				mSet<-Normalization(mSet, n1, n2, n3, ratio=FALSE, ratioNum=20)

				mSet<-Volcano.Anal(mSet, FALSE, 2.0, 0, 0.75,F, 0.1, TRUE, "raw")
				mSet<-PLSR.Anal(mSet, reg=TRUE)
				mSet<-PLSDA.CV(mSet, "L",5, "Q2")
				mSet<-RF.Anal(mSet, 1000,7,1)

				results = cbind (as.data.frame(mSet$analSet$volcano$fc.log)[mSet$dataSet$orig.var.nms,],
					as.data.frame(mSet$analSet$volcano$p.log)[mSet$dataSet$orig.var.nms,],
					as.data.frame(mSet$analSet$plsda$vip.mat[,'Comp. 1'])[mSet$dataSet$orig.var.nms,],
					as.data.frame(mSet$analSet$rf.sigmat)[mSet$dataSet$orig.var.nms,],
					as.data.frame(mSet$dataSet$orig.var.nms))
				colnames(results) <- c('FC', 'Q_value','VIP','RF_MDA', 'metabolites')

				results <- transform(results,Normalization_method = rep(paste(n1,n2,n3), length(mSet$dataSet$orig.var.nms)))
				results <- transform(results,RemoveMissingPercent = rep(p, length(mSet$dataSet$orig.var.nms)))
				rownames(results) <- NULL

				df <- rbind(df, results)

				rm(mSet)
			}	
		}	
	}
}

write.csv(df, paste0('/results.csv'), fileEncoding = "CP932")

Streamlit を用いた条件検討

Streamlit を用いて、24 通りの解析結果をまとめて表示する。
streamlit run ファイル名.py で、ブラウザにアプリケーションが表示される。Streamlit については、こちらを参考にされるとよい。

st_condition_search.py

import os
import streamlit as st
import altair as alt
import pandas as pd

cwd = os.getcwd()
path = ['results.csv']
file = os.path.join(cwd, *path)

df = pd.read_table(file, sep=',' ,header=0)

axis_list = ['FC', 'Q_value', 'VIP', 'RF_MDA']
metabo_list = sorted(set(list(df['metabolites'])))
RemoveMissingPercent_list = sorted(set(list(df['RemoveMissingPercent'])))
Normalization_method_list = sorted(set(list(df['Normalization_method'])))

x_axis = st.selectbox('Select x axis',  axis_list)
y_axis = st.selectbox('Select y axis',  axis_list)
metabolites = st.multiselect('Select metabolites',metabo_list)
RemoveMissingPercent = st.multiselect('Select RemoveMissingPercent', RemoveMissingPercent_list)
Normalization_method = st.multiselect('Select normalization method', Normalization_method_list)

def select_all_elements(l,ALL):
    if not l:
        l = ALL
    else:
        pass
    return(l)

metabolites = select_all_elements(metabolites,metabo_list)
RemoveMissingPercent = select_all_elements(RemoveMissingPercent,RemoveMissingPercent_list)
Normalization_method = select_all_elements(Normalization_method,Normalization_method_list)

df = df[(df['metabolites'].isin(metabolites))&
        (df['RemoveMissingPercent'].isin(RemoveMissingPercent))&
        (df['Normalization_method'].isin(Normalization_method))] 

plot = alt.Chart(df).mark_circle(size=40).encode(
    x=x_axis,
    y=y_axis,
    tooltip=[x_axis, y_axis, 'metabolites', 'Normalization_method', 'RemoveMissingPercent'],
    color='Normalization_method',
    opacity='RemoveMissingPercent'
).properties(
    width=700,
    height=500
).interactive()

st.write(plot)

以下の散布図で、各代謝物の値が表示される。各 plot の色は標準化条件を表し、不透明度 (opacity) は、RemoveMissingPercent を表す。

スクリーンショット 2020-10-18 11.42.19.png

軸の種類は、ドロップダウンで選択できる。
スクリーンショット 2020-10-18 11.51.46.png

カーソルをプロットに当てると、標準化条件と解析結果が表示される。

スクリーンショット 2020-10-18 11.53.01.png

2-hydroxyibuprofen の FC が著しく高いことに興味がある場合、他の代謝物のプロットを消すこともできる。同様に、ドロップダウンから 2-hydroxyibuprofen を選択すると良い。

スクリーンショット 2020-10-18 11.59.36.png

特定の条件の解析結果のみを表示したい場合も、ドロップダウンで同様に選択することができる。

スクリーンショット 2020-10-18 12.01.20.png

スクリーンショット 2020-10-18 12.01.39.png

参考

10
19
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
10
19

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?