yokaのblog

琵琶湖で微生物の研究してます

RによるNMDSを用いた微生物群集構造解析

今日も解析で苦労したので備忘メモ。
※今後、自分の理解に合わせて勝手に加筆・修正するかもしれません。

以下のものは「とりあえず動かすところまで」を目標に書いたものです。もし間違いがあった場合はご指摘くださると有難いです。

なお、今回勉強するにあたっては、以下の記事を大いに参考にしました。英語がOKであれば、こっちを読んだ方が正確で手っ取り早いかと思います。
NMDS Tutorial in R – sample(ECOLOGY)
NON-METRIC MULTIDIMENSIONAL SCALING (MDS)
Ordination analysis [Analysis of community ecology data in R]

目的

 Rを使って、16Sアンプリコンシーケンスの分析から得られたOTUテーブル(site = 50 / OTU = 800 程度)から、NMDS(non metric multidimensional scaling)を用いてサイト間の群集の類似性(β多様性)を比較し、さらに各サイトの環境要因のメタデータ(温度・栄養塩・酸素濃度・・・)の情報を載せて、どの要因がサイト間の群集の差を強く説明しているのかを一覧できる散布図を作る。

こんな図を作りたい(論文でよく見るやつ)↓

f:id:yokazaki:20160629164818j:plain

なぜNMDS?

 複雑なものを2次元のプロットに落として視覚化するタイプの多変量解析にはPCAやCAなどの手法もある(参考※pdfです)。NMDSは微生物生態学の論文では良く見るけど、他の分析と比べて一体何がいいのだろうか?
 一番の利点は、rankに基づいて距離を計算しその距離関係に基づいてプロットを描くというノンパラメトリックな手法であるため、より幅広い事例に応用できてロバストな結果が出やすい、ということだ。さらに目的に応じて距離の計算方法を自由に変えられる、ということも利点とされている。(後者に関しては、以下で紹介する手法のデフォルト設定値(Bray-Curtis指数による計算)がすでに生態学ではほとんどの場合最適とされる計算方法になっているため、今回は触れない(とくに設定を触る必要は無い)。)
 逆に欠点としては、NMDSが反復的に計算して最適解を見つけてくる手法であるため、他の手法のように1点の最適解に落ち着かず、計算の度に結果が微妙に変わることや、マシンパワーを必要とすることが挙げられる。これに関しては、以下で示すように、ランダム試行の出発点となるseed値と、計算の反復回数を固定することで、再現性を担保した結果が得られるようにすることはできる。また各軸がどれだけデータのばらつきに寄与しているかが示されるPCAなどと違って、NMDSでは軸自体には何の意味も無いため、点が近いか遠いかの位置関係だけによって結果が解釈されるのも注意点だ。

準備するもの

Rにパッケージ"vegan"を入れておく。
以下2つのデータを行名・列名を付けた状態でデータフレームにしてRに読み込ませる。

  • ①OTU_table (各行がsite,各列がOTU)
  • ②Environmental_parameters (各行がsite,各列が環境要因)

Rarefactionを行う場合

 通常アンプリコンシーケンスを行うときは、各siteのDNA量が均等になるようにサンプルを準備するため、本来はどのsiteも同じ数だけDNAリードが取れてくるべきなのだけど、このゆらぎが結構大きいようで、往々にしてサンプル間のリード数が最大3倍くらい違うシーケンス結果が出てくる。リード数が違うまま異なるsiteのサンプル同士を比較すると、当然バイアスがかかる。なので、そのばらつきを平準化するために、各siteからリードをランダムにリサンプリングするのがRarefactionだ。この手法、微生物生態学では広く行われてきた作業だけど、得られたリードの一部しか使わない(往々にして半分以上捨てることになる)、言い換えると(ランダムに選ぶとはいえ)「とったデータを無かったことにする行為」になるため、不正確な分析につながっているという問題点も指摘されていたりするので、(特にレアなものを対象にしたい場合など)データの扱いには慎重になる必要がある。
 なお、古くから一般的に行われているRarefactionは、「最もDNAリード数が少ないサンプルに合わせて他のサンプルのリード数を削る」というAbundanceベースのRarefactionなのだけど、最近になって「母集団の多様性推計値をどれくらいカバーできているかで各サンプルのリード数を合わせる」というCoverageベースのRarefactionも提唱されていて、こちらのほうがサンプル間の多様性の違いを正確に比較できるという報告もある(文献)。まずは両方やってみる、ということで、以下、それぞれの方法を説明。

AbundanceベースのRarefaction

 veganに入っているrrarefyという関数でランダムリサンプリングができるので、リード数が最も少ないサンプルに合わせてサンプリング数を指定してやれば良い。

set.seed(123) #再現性をとるためにランダム変数を固定(数字は何でもいい)
OTU_aburared<-rrarefy(OTU_table,min(rowSums(OTU_table))) #最小サンプルに合わせてリサンプリング
write.csv(OTU_aburared, "abundance_rared.csv") #CSVに出力したい場合

CoverageベースのRarefaction

※以下に関してはメタバーコーディング・環境DNAバーコーディング解析の技法の「4.メタゲノムデータを用いた群集統計解析法」の掲載資料を大いに参考にさせていただきました。

 Coverageベースでやる場合、最後は同じくrrarefyでリサンプリングをしてくるのだけど、その前にRarefactionカーブの傾き(=カバー率)が一定の値(=一番カバー率が低いサンプルの傾きと同じ傾き)以下になるDNAリード数をサンプルごとに求めなくてはならない。veganのrareslopeという関数を用いて、Rarefactionカーブの傾きを求められるので、各サンプル、1リード増やすごとの傾きを、forを使ってとってくるところがポイント。

#各サンプル、1リード増えるごとの傾きをリスト形式で一括出力(要計算時間)
rareslopelist<-list()
for(i in 1:nrow(OTU_table)){
  rareslopelist[[i]]<-rareslope(OTU_table[i,],1:(sum(OTU_table[i,])-1))
}#最後の1リードを入れずに「-1」としているのは、最後の1リードについては傾きが計算できない(0になる)ため

#最後のリードの傾きが最も大きい(=カバー率が最も低い)サンプルを探す
getmincov<-c()
  for(i in 1:nrow(OTU_table)){
    getmincov[i]<-rareslopelist[[i]][length(rareslopelist[[i]])]
}#最後の1リードの傾きだけをベクトルでとってくる

(1-max(getmincov))*100 #カバレッジの%を確認

#指定したカバレッジに到達した(=傾きが指定値を下回る)瞬間のリード数をサンプルごとに採ってくる
cvrfun<-function(x){min(which(x<=max(getmincov)))+1} #関数を設定。上記で1を引いた分を足し戻す
cvrrare<-unlist(lapply(rareslopelist,cvrfun)) #lapply+unlistでベクトル形式にして一括で値を取得

set.seed(123) #再現性をとるためにランダム変数を固定(数字は何でもいい)
OTU_covrared<-rrarefy(OTU_table,cvrrare) #得られたリード数に沿って各サンプルからリサンプリング

write.csv(OTU_covrared,"covarage_rared.csv")  #CSVに出力したい場合

(補足)α多様性の確認

diversity関数で各siteの群集そのものの多様性(α多様性)の計算も可能。変数を指定することで、複数の多様性指数(ShannonとかSimpsonとか)を計算してくれる。

diversity(OTU_aburared,index="shannon")

NMDSをかける

分析にかけるデータフレームができたら、"vegan"のmetaMDS関数を使って、NMDSプロットに必要な「距離の計算」を行う。

set.seed(123) #NMDSの再現性を担保するため、ランダム変数を固定(数字はなんでもいい)
nmds<-metaMDS(abundance_rared,k=2,trymax=20) 
#"abundance_rared"の部分はphyseq-classのデータでもOK
#kは次元数。基本は2次元にしたいので2を使う
#trymaxは計算の最大試行回数。増やすとより良い解が見つかりやすくなるが、
#計算時間もかかるため、検討段階では少なめに、本番では大きな値を使ってやる
#この回数に達する前にプログラムが最適と判断して計算を止めることもある

ここで、NMDSの出来を評価するのが"stress"という値だ。これが小さければ小さいほど、データを上手く2次元に落とし込めたということを意味しているらしい(原理は完全には理解していないので不正確な表現かもしれません)。この値は、

nmds$stress

で見ることができる。で、どの程度であればOKか?という話なのだけど、Web上の情報で良く出てくるのは「>0.3だとダメダメ、>0.2だとまあまあ、>0.1だと良い、>0.05だと最高!」という感じの指標なのだけど、サンプルサイズや、site間の環境の違いの大きさなどにも影響するため、一概には言うべきではない、という主張もあったりして、統一見解はなさそう。ので、できるだけ低いstress値が出た結果を、その値と一緒に載せるのが外に出すときの一般的なやり方っぽい。
 なお、set.seedの値をいじれば分かるけど、最初の出発点によって結構結果が変わってくる。出発点によっては局所最適みたいなところに落ち着いてしまって、計算回数をいくら増やしてもstressが減らないということも起こる。ので、本番の分析をやるときは、

stress.value<-NA
for(i in 1:100){
 set.seed(i)
 nmds<-metaMDS(abundance_rared,k=2,trymax=50)
 stress.value[i]<-nmds$stress
}
grep(min(stress.value),stress.value)

とかして、低いStress値を叩きだすseed値を探してくるのが良さそう。
 ここまでくれば、あとは作図関数に投げるだけだ。veganに入っているordiplotorditorpを使えば、metaMDSで生成したデータ(上記でいう"nmds")を、縮尺など調整して綺麗にプロットしてくれる。

win.graph(30,30) #正方形の作画領域を生成
ordiplot(nmds,type="n") #枠を書く
orditorp(nmds,display="sites",air=0.1,cex=1) #点を書く
#displayをsitesとすればsiteの、speciesとすればOTUのプロットを描いてくれる
#airは重なって表示されるラベルをずらすための変数
#変数labelsで表示ラベルを、selectで表示/非表示プロットを指定可能
#普通のplotの変数がそのまま使え、好きな色や形の点を打てる

環境変数をNMDS上に載せる

 ここまでで各site間の組成がどの程度違うのか(β多様性)を、視覚化することができた。近いものほど似ており、離れているものほど違う。この段階で誰の目にもはっきり分かるくらい傾向が出ていれば、これ以上の説明はいらないだろう。だけど、一目で結果が分かりにくい時、各siteの環境条件のデータを付け加えることで、「温度が強く効いている」とか「酸素濃度はそれほど影響がなさそう」とかが見えてくることがある。そこで冒頭に示したもう一つのデータセット「②Environmental_parameters」をveganのenvfit関数を使ってNMDS上に反映させる。

envNMDS<-envfit(nmds,Environmental_parameters)
#データにNAがある場合は変数"na.rm=T"を加える

これで先ほど作った"nmds"上に"Environmental_parameters"の値を載せた場合の各環境変数の説明力が計算された。このデータを見れば、各環境変数が、どの程度NMDS上の各siteの距離関係のばらつきを上手く説明できているかを示す表が得られる。

envNMDS
#             NMDS1    NMDS2     r2 Pr(>r)    
#env1       0.77127 -0.63650 0.4564  0.001 ***
#env2      -0.67797  0.73509 0.6592  0.001 ***
#env3      -0.07703 -0.99703 0.1305  0.175    
#env4      -0.15870  0.98733 0.5373  0.001 ***
#env5       0.56661  0.82399 0.3727  0.004 ** 

当てはまりの強さは、P値で示されており、上の例ではenv1,env2,env4がうまくsite間のばらつきを説明できていて、env3が説明できていない環境要因だということが分かる。
これを先ほどのNMDS上にプロットする。普通のplot関数で"envNMDS"を読んで上記のNMDS上に上書きするだけで、矢印まで描いてくれるようになっている。

win.graph(30,30) 
ordiplot(nmds,type="n") 
orditorp(nmds,display="sites",air=0.1,cex=1) 
#ここまでは上記と同じ
plot(envNMDS,p.max=0.01) #p.maxでその値以下のP値を叩きだした環境要因のみプロットする

これで冒頭の図と同じものができる。なお、各環境要因の矢印だけど、

向きと相対的な長さだけが意味を持ち、矢印の始点と終点、絶対的な長さには意味が無い

ということに注意。まず矢印の向きだけど、「この向きがsite間のその環境要因の傾きが最も急な方向である」ということを示している。そして矢印の長さは「その傾きの急さ」を示しており、上記のr2の値の平方根に比例してスケーリングされるようになっている。そしてこの「長さ」は各矢印間の相対値で見るべきで、絶対的な長さは、プロット上で最も見栄えが良いように自動的に決定されているに過ぎない。つまり、矢印はNMDS上のプロットの縮尺とは独立した縮尺で描かれていて、その位置関係を比較することはできないということだ。
 この説明では少しわかりにくいかもしれないけど、同じくveganについているordisurfという関数を用いることで、この傾きの「傾斜」と「向き」を等高線的に表示することもできて、これを見れば矢印がどういう意味なのかももう少し分かり易いのではないかと思う。

ordisurf(nmds,Environmental_parameters$env1,add=T)
#"add=T"によって、先述ののNMDSプロットに上書きできる。
#インプットデータは元のデータフレームから1列のみ抜き出したものにする

この関数の原理については、僕はきちんと理解できていないのだけど、以下のページが大変詳しい。
What is ordisurf() doing...?
 なお、これらの情報を調べていると「ノンパラメトリックな手法であるNMDS上にenvfitで出したリニアな回帰の結果を載せていいのか?」という議論を所々で見かけるのだけど、veganの開発元でも本件はFAQとして取り上げられていて、「距離の計算方法はノンパラメトリックであるが、NMDSのプロット自体は距離スコアに応じた線形のデータなので、envfitordisurfで作ったデータを載せることは問題ない」という見解のよう。一方で、この方法で論文のレフェリーに突っ込まれたという話も出てきたりして、自分の直感を超えた統計手法を使うときは、目的をはっきりさせて原理としっかり理解してからでないと怖いなぁということで、まだまだ勉強途上ですが、ここまでの理解をまとめてみた次第です。