初稿という既成事実

　月例の調査やサンプル処理は進めつつも、最近はデスクワーク中心の毎日。次の論文に載せるデータを作りながら、そのデータを使って来月の国際学会での発表ポスターを作成中。

　前回の論文はかなりスピード重視で書いたけど、今回のはあえて焦らずに分析を深める方向で進めている。とにかく早くアウトプットが欲しかった前回と違って、今回はスピードにそこまで貪欲ではないということや、前回よりもデータ量もストーリーもでかいので分析に時間がかかるということもある。だけどそれ以上に、前回の論文執筆・査読対応を通じて

丁寧な分析というのは思ったよりも評価されるものであり、初稿の段階で分析を深めておくことがむしろ、労力に対して質の高いアウトプットを出すことにつながるのではないか

という思いを強くしたことが理由としては大きい。

　前回の論文は、会社員時代の仕事の仕方を引きずっていたこともあって、

判断に迷うところはとりあえず仮決めで進めて、早めに審判を仰ぎ、軌道修正が必要かどうかを早く知るほうが効率がいい

という考えで、あまり深い方向に突っ走らず、できるだけ最小公倍数的なモノを書いて査読に回すことを目指してすすめていた。その結果、根本の考え方のエラーを指摘され、分析まるごとやり直しになって、二度手間を回避することができたパートもあった。だけど同時に面白かったのが、自分は「最小公倍数で仮決めだ」と思っていた分析でも、評価されて結局最後まで生き残ったパートも結構たくさんあった、ということだ。必要であればもっと深い分析に突っ走ることもできたけど、そこまで行く前に出しておいて意外に大丈夫なんだな、ということを思った。もちろん、掲載先のジャーナルのレベルにも大きく依るのだろうけど、結果として「とにかく早く出したい」という僕の欲望に応える結末になった。

　これは決して「手抜きで不完全な論文を出してしまった」ということではない。言いたいのは「もっと深い方向に突っ走ることもできたけど、敢えて深入りせずに置いておいたら、それで意外とうまくいった」ということだ。仮に深い分析をして出していたところで、「不要なこだわり」と判定されて、ボツになっていた可能性だって大いにあるだろう。だけど一通り論文が出てみて「最初からできる限り深い分析をしておいても良かったな」というのもちょっと思った。

会社にいた時も思っていたけど、

「仮決めで仕事を進めること」で二度手間を回避するチャンスも増えるけど、誰かが修正してくれるだろうと思ったら誰も指摘してくれずにそれがそのまま既成事実になってしまうリスクも増える

ということが、研究論文にも当てはまると感じた。いくら業界のプロが査読するとはいえ、「査読⇔修正」のやり取りはせいぜい数回までだし、時間が限られている中で、初稿がベースになって議論が進んでいくのはまあ当たり前のことだ。

　だから今回は、最初から分析できるところまで分析しつくした、最高クオリティの初稿を書き上げたいと考えている。判断に迷うところも、仮決めで進めるのではなく、全ての可能性を盛り込んだ最大公約数的な分析にしていきたい。その結果、出したものの多くが不要になるかもしれないけれど、深く分析したこと自体は査読者はちゃんと評価してくれて、無駄になることはなさそうだという自信が出てきたことと、初稿で最高クオリティのものを送っておかないと、後からその既成事実を上回るものに書き換えることは難しいということが分かったからだ。

　なので、当初の予定からかなり後ろ倒しになっているけど、緻密な議論・分析をすべく、じっくりと時間をかけて次の作品の準備を進めている。「通ればいい」という欲望ではなく「驚きとワクワクを与えたい」という欲望をモチベーションにして、恥ずかしくないものを世の中に出したいと思う。

7月びわこ

今日は琵琶湖調査の日でした。ちょっと波風があるけど、良い天気。

気温28度でまだマシな暑さ。水温は26度、透明度は6mに増えて、水中も真夏モードに入りつつある。

いつもは深層ばかり採水しているのだけど、今回は別の目的で久しぶりに表層から採水。やっぱり深層と全然違って粒子が多いので、処理の大変さが全然違う。いつもの半分の濾過量でフィルターが水漏れを起こしてしまった。

さっそく細菌数を数えてみたのだけど、7.1×10^6 cells/mlで、深層の7倍近い現存量。これは琵琶湖にしてはかなり多い数だけど、盛夏に向かっていく今の時期には過去にも報告されているレンジでもある。梅雨が終わって、表層に残された最後の栄養の取り合いで細菌達も数を増やすのに必死なのかもしれない。ウイルスは4.8×10^7 particle/mlで、これも深層の倍以上の値。

変化が激しい表層の調査は色々緊張する。深層の安定性や再現性があってロバストなデータが出る感じは、やっぱりいいよなぁと思った。

琵琶湖に生息する微生物の図鑑

自分の研究の論文が出版されたので、簡単に紹介。

Vertical partitioning of freshwater bacterioplankton community in a deep mesotrophic lake with a fully oxygenated hypolimnion (Lake Biwa, Japan) - Okazaki - 2016 - Environmental Microbiology Reports - Wiley Online Library

一言でいえば

琵琶湖の水中で季節・深さごとに出現する細菌を網羅した図鑑を作成

3行で言えば

・琵琶湖沖で15か月・3水深にわたって、水中の細菌群集を網羅的に調査

・特に、研究例が少ない、深層に生息する細菌に注目

・深層にのみ生息する、10以上の「変わった系統」の細菌を特定

という内容。

研究する意味

「琵琶湖の細菌の研究なんてもうやりつくされているんじゃないの？」と思われそうだけど、実はまだほとんど何もわかっていない。単に研究者がいなかった、というのも大きな理由だけど、やはり一番大きな障壁は、 

ほとんどの細菌は培養ができないので、研究もできない

ということだ。

琵琶湖の水中には1ml（ティースプーン1杯）あたり100万を超える細菌が生息している。顕微鏡でみれば、それこそ星の数ほどの細菌が見られる。

　だけど、環境中にはこんなにたくさんいるのに、なぜか研究室に持ち帰って人工的な環境で飼おうとすると、ほとんどの細菌は増えてくれない（正確には、少数の増えやすい奴ばかり増えてしまう）。その理由はまだよく分かっていない。人工的な環境では再現できないような、環境中に存在する絶妙な何かが欠けているのだろう。で、培養ができなければ、個々の細菌系統の性質を詳しく調べることもできない。だから、難培養の細菌系統を培養する、というのは大仕事で、日夜多くの研究者たちがリソースを割いているし、重要な系統を単離することに成功すれば、それだけでNatureに掲載されるレベルのインパクトを世界に与えることができる。

　琵琶湖の細菌も例外ではなく、細菌全体の存在数やその季節変動を報告した研究は昔から例があったものの、培養できないことが障壁となって、「どのような系統が生息しているか」ということは解明されずにいた。

　そこに訪れた大きな転機が「次世代シーケンサー（NGS）」の登場だ。簡単に言えば、環境中に存在するDNAを集めてきて、片っ端から塩基配列を読むことで、どの遺伝子を持った細菌がどれくらいいるかを、培養することなしに推定することができる技術だ。細菌の系統は特定の遺伝子（16S rRNA遺伝子）の塩基配列を元に分類でき、データベースに情報が蓄積されている。だから、この遺伝子の塩基配列の情報さえ手に入れば、それが培養できない系統であっても、どのような分類に属し、過去にどのような環境で見つかってきたものなのかは分かる。

 このNGSの登場は本当に革命的で、水中に限らず、あらゆる環境に生息する細菌の生態が急速に明らかになりつつある。微生物学者にとっては「ボーナスステージ」と言っていいような状況で、様々な動物の腸内、家のドアや電車のつり革、歴史的絵画の表面に至るまで、調べれば調べるだけ新しい発見が出てくる事態で、論文もバンバン出版されている状況だ。

 その一方で、「NGSを使ってまだ調べられていない環境の細菌の組成を調べました、おしまい」というだけの、単なる記載的な研究が増えているという批判があるのも事実だ。

新しいところを調べれば、新しいものが出てくるのは当たり前。それが何を意味しているのか、どう面白いのかが肝要だ

という流れになりつつあるというのが、NGS技術が一通り普及しきった今の業界の状況だ。

何をしたのか？

　このNGSを琵琶湖で初めて使い、どのような細菌が生息するのかを網羅したのがこの研究だ。だけど、単に表層の水を採って細菌群集を調べるだけの研究であれば、すでに世界中の湖沼で行われている研究の二番煎じになってしまう。

そんな中で、本研究は、

①これまでに研究が行われてこなかった、深層に焦点を当てたこと

②15か月にわたる継続的な調査を行い、季節変動を明らかにしたこと

という二点を大きなウリにしている。

特に①の「大型淡水湖の深層にはどのような細菌が生息しているのか」ということに関しては、世界的にもほとんど報告が無い。手掛かりになるような断片的な情報は存在するものの、NGSで季節変化を含めて網羅した研究は、本研究が初めてになる。そのような意味で、本研究は

湖の深層に生息する細菌の組成を、世界で初めて網羅的にみる 

のが目的と言える。

何が分かったのか？

　過去に琵琶湖沖に毎月通って採り溜めていた3水深×15か月のDNAサンプルを使って、NGSで16S rRNA遺伝子のシーケンスを行い、とれてきた塩基配列をもとに各サンプルの細菌群集を分析した。すると、表層と深層では生息する細菌系統が大きく異なることが分かった。しかも、これまで表層で良く知られている細菌とは、「門」のレベル（一番大きな生物の分類の区分）で異なる系統が深層にたくさん生息していることが分かった。

　これらの系統の塩基配列をDNAのデータベースで調べてみると、実はすでに、世界中の深い湖からの報告が散発的に存在することが分かった。このことから、これらの系統は、琵琶湖にだけ生息するのではなく、世界中の湖の深水層に普遍的に生息する系統であると考えられる。

　本研究では、これらの過去に報告されてきた知見をまとめながら、琵琶湖で季節・深度ごとに採られたデータを加えて、湖沼の深水層に生息する細菌系統を初めて網羅的に整理することができた。

何がすごいのか？今後の発展は？

　湖の深層に生息する細菌は、湖全体の物質循環や生態系の中で重要な役割を果たしている。本研究が「深層には表層とは全く違う系統が生息している」ことを明らかにしたことで、今後これらの深層独自の細菌系統の生態の解明を進めていかなければならないことが分かった。

　ただ、「じゃあ彼らは生態系内でどのような役割を果たしているのか？」とことについては、現時点では分からない。今回見つかった「門」の多くは、多数の難培養性系統で構成される、ほとんど正体が分かっていない分類群だからだ。残念ながら本研究は「変なのがたくさんいるからもっと研究を進めましょう」というところで終わってしまっている。

　なので、今後は今回明らかにした個々の細菌系統に焦点を当て、どのような遺伝子や代謝活性を持っているか調べて、彼らが生態系で何をしているのか、なぜそこにたくさんいるのか、ということを明らかにしていく必要がある。

　また、今回は琵琶湖のみでの調査だったけど、データベースの情報などから、他の湖にも似たような系統が存在している可能性が高い。今回琵琶湖で見つかった発見が、どの程度一般化できる内容なのか、今後調査地の数を増やして検討していく必要がある。

　前者については、データをとり始めており、後者についてはすでに結果が出揃って、論文の執筆に着手しているところだ。

　また、さらに長い目で見たとき、僕はこの湖の深層の細菌の研究を、単なる湖の生態学の研究で終わらせずに、「生物学の新発見の宝庫」に育て上げたいと思っている。断言できる根拠はないけれど、環境の特殊性や、調査の容易さ、出現する細菌系統の異質さなどを考えると、この環境を深堀りしまくることで、生態学や微生物学、もしかすると生命科学全般に影響を及ぼすような、面白い現象を発見できるのではないか、という気がしている。本研究は、人生を賭けるその長い道のりの、最初の出発点となる予定の研究だ。

ここに至るまで

　しかし、長かった。会社を辞めて研究に戻ってきてから、1年以上アウトプットが無かったことがとてもストレスだった。ライバルが似たような研究をしているのも知っていたし、何よりも「色々やっているのにまだ何も成果になっていない」ということが悔しくて、この1年少し、早く論文を出したい一心だった。

　本研究は、去年の4月、研究に戻ってきた直後に思いついて、修士課程の時に採り溜めていたサンプルを使ってすぐに着手した研究だ。「とりあえず論文1本出そう」くらいの気持ちで、そこからほぼ滞りなく、最速で進めてきたつもりだけど、

・4-5月　予備実験

・6-7月　シーケンス外注

・8-9月　結果分析、一部サンプルの失敗判明　　

・10-11月シーケンス再外注

・12月　　分析

・1月　　論文執筆⇒英文校閲

・2月　　論文投稿⇒瞬殺リジェクト⇒すぐに再投稿

・3月　　レビューコメント付きリジェクト（姉妹誌への投稿はOK）

・4月　　リバイス⇒姉妹誌に投稿

・6月　　マイナーリビジョン⇒すぐに再投稿⇒アクセプト

・7月　　オンライン公開

　という流れで、結局1年以上もかかってしまった。特に、初稿を書き終えた1月から、アクセプトまでに半年も待たされたことは、とてつもなくストレスだった（それでもかなり早いほうだと言われるけど）。また、2回目のリジェクトのあと、査読コメントに従って大幅に書き直しをすることになった4月は本当に苦しかった。一度書いたものを全て分析しなおして、書き直す大変さは並大抵ではない。初稿ができた段階で、8割くらい終わっているかと思ったけど、時間的にも労力的にも、初稿ができた段階でようやく半分くらいだなぁ、と思った。

　会社員時代と同じ仕事スピードで大学院生をやれば、サクサク早く効率的にアウトプットが出せるだろうと思っていた頃もあった。だけどやっぱ研究は、厳密さが全然違うし、一人でやらなければならないから、必要な労力も時間も全然違う。種を撒いてから収穫まで、1年以上は余裕でかかる。だから、去年の秋に種を撒いて、今論文を書いているネタが公開されるのは来年になるだろうし、今種を撒いているやつが出荷されるのは再来年の春とかになるかもしれない。本当に気が遠くなるし、のろまだと思う。でも、しょうがない。ライバルもみんな同じだ。効率的に収穫できるように、できるだけ計画的に種まきをすすめるのみだ。

　しかし、やっぱり、自分の仕事が自分の名前で論文になるのはとても嬉しい瞬間だ。一本目の論文が出た時も「これでいつ死んでも世の中に名前が残せる」みたいなことを思ったけど、またこの感動を味わえて嬉しい。研究に戻ってきてよかったです。

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)を用いてサイト間の群集の類似性（β多様性）を比較し、さらに各サイトの環境要因のメタデータ（温度・栄養塩・酸素濃度・・・）の情報を載せて、どの要因がサイト間の群集の差を強く説明しているのかを一覧できる散布図を作る。

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

なぜ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バーコーディング解析の技法の「４．メタゲノムデータを用いた群集統計解析法」の掲載資料を大いに参考にさせていただきました。

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

#各サンプル、１リード増えるごとの傾きをリスト形式で一括出力（要計算時間）
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に入っているordiplot、orditorpを使えば、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プロットに上書きできる。
#インプットデータは元のデータフレームから１列のみ抜き出したものにする

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

湖沼の水質の調べ方

　今書いている論文で、日本中の湖沼の水質データが必要になって色々調べたのだけど、結構苦労したので方法を備忘にメモ。

　まず、データベースなのだけど、おそらく日本国内の水環境を最も網羅的に調査しているのが、環境省の公共用水域水質測定調査だと思う。大型湖沼はもちろん、河川や小さめの池まで、あらゆる水環境について、年1回以上の調査が行われていて、水温、ｐH、溶存酸素、全リン、全窒素等の基礎的データに加え、サンプルによっては微量金属や有害有機物の数値まで測定されている。湖沼によっては、都道府県や個々の研究者によって、もっと細かい測定項目のデータが測定・公開されている場合もあるけど、湖沼間のデータを比較したい場合などは、網羅的に統一された項目・方法で測定されている本調査の結果を使うのがやっぱり便利だ。

　で、このデータへのアクセスなんだけど、普通に「公共用水域　水質」とかでググっても全然出てこない。出てこないどころか、リンク切れだらけの前のバージョンのウェブサイトや、環境省の元データを転載して作られた中途半端なデータベースなど、ノイズ情報がヒットしまくる。しょうがないので、僕は頑張って探した本家のサイトをブックマークに登録することにした↓

水環境総合情報サイト

　で、ここにたどり着いてからも、色々分かりにくい。まず、非常に似た作りで複数の異なる情報が並列で公開されているため、自分が間違ったところにいることに気づかず、グルグルしてしまうことが起こる。

　ここでは大きく3つの情報が公開されている。その入り口が、左側にある「水質調査データの公開」の項目にある「①公共用水域水質測定データ」「②広域総合水質測定データ」「③水浴場水質測定データ」だ。今回見たいデータベースは①だ。なのでまずこれをクリックする。

　すると日本地図の上に色々表示された画面に出る（いちいち書かないけど、この画面の操作性も最悪だ）。今回欲しいのはデータそのもの（数値が載ったエクセルファイル）なので、地図には目もくれずダウンロードページに飛ぶ。

　が、ここでまた罠がある。右上に「ダウンロード」という項目があるので、そこをクリックすると、色々ダウンロードできそうなページに飛ぶ。だけど、なぜか僕のネット環境からのアクセスだと、どのボタンを押しても一切反応せず、何もダウンロードすることができない！なんなんだこれは。

　ということで正解は「測定結果検索」だ。ここをクリックすると、データを検索できるページに飛ぶ。このページの使い方もどこにも説明がなく、未だに半数以上の項目の使い方が分かっていないのだけど、とりあえず目的に到達するための道順だけ書く。

　まず、「データ種別」だけど、月ごとの一番細かいデータが欲しいなら「公共用水域（検体値）」だ。年間平均のデータだけで良い場合は「公共用水域（年間値）」を選ぶ。

　次に「測定地域」だけど、知りたい湖や川の名前がもう決まっているのであれば、その名称で「水域名」で検索（「総括表を表示」をクリック）をかければ一発だ。だけどここにも罠があって、完全一致でないと検索結果を拾ってこない仕様になっているので、入力する文字を間違えるとデータが出てこない。例えば「琵琶湖」で調べても出てこないのだけど、これは調査対象に琵琶湖が含まれていないのではなく、琵琶湖は北湖と南湖で別の扱いになっているので、例えば北湖を調べたければ「琵琶湖（１）（琵琶湖大橋北）」と入れなければヒットしないのだ。なんという使いにくさ。なので、対象の水域が検索で引っかからなくても、都道府県別の検索などもやってみながら、本当にその水域のデータが存在しないのかどうかは確認したほうがいい。

　で、ついに対象湖沼のデータが表示され、左上にダウンロードボタンが現れる。が、ここも罠だ。このままダウンロードを押してしまうと、デフォルトの表示項目（今ブラウザ上に表示されている項目）しかエクセルになって出てこない。本当はもっといろんな項目が測定されているはずなのにだ。全ての項目のデータをゲットするためには「表示物質の設定」から全項目を選択し、画面を更新してからダウンロードボタンを押す必要がある。これでようやく、全リン・全窒素等を含んだ、その湖の対象期間の全ての水質調査結果が手に入る。

　データはCSV形式だ。湖の場合、データは測定地点、測定月、測定深度の順にソートされたものが出てくる。このデータも非常に見づらい。まず各測定地点が「ST-1」みたいな書き方しかされてなくて、一体どの地点のデータなのかがどこにも書かれていない。これはググりながら環境省の過去の報告書を漁ってきて自分で調べるしかなく、すごく時間が掛かった。また、各項目の単位もどこにも書かれていない。これは、検索画面の「測定値」のプルダウンに載っていることを発見して、事なきを得た。これでようやく日本中の湖沼の水質データの比較ができる。

　・・・以上、こんなことに丸一日近く費やしてしまったので、今後のためにメモっておく。労力もコストもかかった、継続的・網羅的な素晴らしい調査なのに、データベースが使いにくすぎて、本当に残念だ。

　ちなみに日本中の湖沼のデータを比較していると、都道府県ごとに調査の性格（？）みたいなのが見えてきて面白い。北海道は地点が多いせいか、広く薄く採られている印象。富士五湖は測定項目も頻度も充実している。猪苗代湖・中禅寺湖・池田湖は深いところまできちんとデータが採られているのが有難かった。琵琶湖は地点数・項目数・測定頻度とも他を圧倒していて、さすがのこだわりを感じた。

早期取得断念

　今所属している、京大理学研究科（生物科学専攻）の博士の学位の最低取得用件は「査読付き国際誌に一報以上掲載」だ。僕はすでにその基準を満たしていて、できれば早く学位をとって、研究費申請や共同研究の自由度を上げたい、と思い、3年未満での早期取得を目指していたのだけど・・・どうやらこれは難しそうなことが分かった。

　すでにD論のストーリーは固まっていて、指導教官からも前向きな意見を頂いていたので、今執筆中のD論のシメとなる論文が良いところに通れば、すでに持っている中堅どころ（＞IF３）のジャーナルの2報と合わせて、十分に学位の申請ができると思っていた。だが、本部に確認したところ、前例と照らし合わせると、早期取得は「とてつもない科学的貢献」が必要とのことで、超一流誌（Natureクラス）に載るようなレベルの研究業績がなければ、認めることはできない、との回答だった。

　「そう言われたのでしょうがないんでNatureに載せました」ってなれば超かっこいいんだけど、今やっている研究は、将来Natureになるかもしれないネタを見つけてくるための下地固めなので、残念ながらそこまでのインパクトは無くて、うまくいって分野内のトップジャーナル（ISME J）レベルかなぁと思っている。

　にしても、3年以上いれば「査読付国際誌1報（どこでもいい）」で学位がもらえるのに、3年未満だったら「超一流でなければだめ」というのはなんか腑に落ちない。それにその「超一流」ってのも客観的な指標ではない。調べてみると、よその大学では「これこれ以上のジャーナルに何報以上」みたいな客観的な基準が書かれていたりして、多くの大学では僕はその基準を満たせるように思われた。けどでも確かに、僕ぐらいの業績の人は京大にはゴロゴロしているし、その中でも前例が超一流だけ、という事実を突きつけられると、僕は黙るしかない。

　まぁ、博士の学位は免許みたいなもので「めんどいからさっさととってしまいたい」くらいの価値にしか思っていないので、そんなものにこだわるのもカッコ悪いし、逆に「3年間在籍してればくれるんでしょ」くらいの気持ちで、あまり気にすることなく研究に集中できるようになった、ということでポジティブにとらえて頑張ろうと思います。

　・・・しかし、もう1年、自分の名前で研究費の応募が出せないのはつらいな・・・DC１の外部資金に応募できない縛り、撤廃してくれないかな。新しい研究アイデアを人に評価してもらうチャンスがない、というのは若手研究者を育成する制度としては間違っているような気がするのですが。

6月びわこ

今日は月例の琵琶湖調査の日。

湖の北端が見渡せるほど済んだ空気。最近ずっと天気に恵まれている気がする・・・

透明度は5.6m、水温は21.2℃で、深層の温度は約8.2℃。早速顕微鏡で細菌を観察したけど、今年も例年と同じ変化が起こり始めていて、すごくワクワクしている。これだけ毎年予測可能な微生物の季節的挙動ってあまりないと思う。この生態系を調べ尽すことで、必ず面白い発見が出てくると信じている。