MaxEnt解析手順メモ

前回までの記事で琵琶湖南湖の湖沼図から地形に関するラスターデータを作成しました。今回はこれを利用してMaxEnt解析を試みます。MaxEntはオープンソースのソフトでここからダウンロードできます。なお、ここではバージョン3.4.4を使用しています。サンプルとして、ブルーギルの分布調査¹から抜き出した2018年の分布データを使ってみます。

MaxEnt解析に必要なものはターゲットである種の分布データと地形・植生などその種の分布に関係があると思われる地理空間情報です。

1. 分布データの準備

ダウンロードした分布調査結果の図からブルーギルの分布情報（2018年6月）を抽出し、csv形式で保存しました（ファイルはこちら）。正確なものではありませんが、今回はとりあえずMaxEntを動かすことを目的としているので良しとします。

2. 地理空間情報データの準備

今回は地理空間情報データとしてこれまでの記事で作成した深度・斜度・底質のラスターデータおよびWorldClimからダウンロードした6月（分布調査が6月だったため）の平均気温・風速データを用います（空間解像度: 30 secondsを使用）。

地理空間情報データをMaxEntに読み込ませるためには、それぞれのラスターデータの範囲とセルサイズを一致させた上で、Esri ASCII ラスター形式（.asc）またはDIVA-GISグリッド形式（.grd）にする必要があります。

2.1. データの範囲

ラスターデータから特定の範囲を抜き出すための外枠を用意する必要があります。今回は琵琶湖が舞台なので、湖沼の水域データが活用できそうです。国土交通省の国土数値情報ダウンロードサイトからダウンロードします。

ダウンロードしたファイルを解凍し、QGISで「ベクタレイヤを追加」からシェープファイル（.shp）を読み込みます。QGISに日本の湖沼の水域が表示されます。他の湖沼は邪魔なので、まずは琵琶湖だけを別のレイヤに移します。「シングルクリックによる地物選択」をクリックし、琵琶湖の上でクリックします。すると図1のように琵琶湖全域が選択されます。

図1

続いてメニューの「編集」>「地物をコピー」から「編集」>「新規レイヤへの地物貼り付け」>「一時スクラッチレイヤ」をクリックし、レイヤ名を入力して「OK」を押します（ここでは「LakeBiwa」とした）。これにより琵琶湖だけを別レイヤにコピーできました。最初に読み込んだベクタレイヤを非表示にするか削除し、QGIS上に琵琶湖だけを表示します。

次に琵琶湖のうち、解析に使う南湖のみを切り出します。今はまだ琵琶湖全体が一つの地物となっているので、これを複数に分割します。レイヤを選択した上で、「編集モード切替」をクリックし、「地物を分割」をクリックします（図2）。

図2

この状態で、切りたい部分の開始点をクリックし、終点をクリックすると切断箇所が直線で表示されます。最後に右クリックすることで切断されます（図2の直線は切断後）。切断ができたら「部分の削除」をクリックし、消したい方の地物をクリックすることで地物が削除されます。これにより必要な領域を切り出します（図3）。

図3

これで完成に見えるのですが、このままだとラスターデータを切り抜く際にエラーとなります。これを回避するため、メニューの「ベクタ」>「ジオメトリツール」>「マルチパートをシングルパートに変換」をクリックし、「入力レイヤ」にさきほど切り出した琵琶湖南湖のレイヤを選択して実行します。「出力レイヤ」という名前のレイヤができるので、これを右クリックし、「エクスポート」>「新規ファイルに地物を保存」でシェープファイル形式（.shp）で適当な名前を付けて保存します（ここではSouthLakeBiwa_boundary.shpとした）。またこの際、「座標参照系」はラスターデータと一致させたほうが良いかと思います。

2.2. 地理情報空間データの編集

上で作成した琵琶湖南湖の範囲データを使用してラスターデータを切り抜きます。まずダウンロードした平均気温データ解凍し、「ラスタレイヤを追加」から読み込みます（データは.tif形式）。そして、メニューの「ラスタ」>「抽出」>「マスクレイヤで切り抜く」をクリックします。入力画面が表示されるので、「入力レイヤ」に平均気温データのレイヤ、「マスクレイヤ」に琵琶湖南湖の範囲データ（ここではSouthLakeBiwa_boundary）を選択し、「実行」をクリックします。完了すると「出力ファイル」という名前のレイヤに切り抜いたラスターデータが出力されます（図4）。このレイヤを右クリックし、「エクスポート」>「名前を付けて保存」でファイル形式をGeo TIFF形式（.tif）で保存します（ここではclipped_tavg06.tifとした）。この段階ではまだ（.asc）形式にはしないでください。

図4

繰り返しになりますが、MaxEnt解析を行うには全てのラスターデータの範囲とセルサイズを一致させなければなりません。そこで、先ほど作成したclipped_tavg06.tifを基準として他のラスターデータを変換していきます。QGISには「ラスタを揃える」という機能があり、これが使えるかと思ったのですが、なぜかエラーで動きませんでした。少々面倒ですが、Rを使って変換をしていきます。方法は以下の通りです。

>library(raster)　#使用するパッケージです。なければinstall.packages("raster")でインストールしてください。
>library(rgdal)　#同上。

>s = raster("clipped_tavg06.tif")　#基準とするラスターファイルの読み込み。
>ex = extent(s)　#基準ラスタの範囲を抽出。
>writeRaster(s, "clipped_tavg06.asc", "ascii")　#ここで基準としたclipped_tavg06.tifをEsri ASCII ラスター形式（.asc）で保存。

>r = raster("wc2.1_30s_wind_06.tif")　#変換するラスターファイルの読み込み。
>r = crop(r, ex)　#変換するラスターを基準ラスタの範囲で切り抜く。
>resample(r, s, "ngb", filename="clipped_wind06.asc", format="ascii")　#基準ラスタのセルサイズに調整。"ngb"は調整方法で最近傍法を意味する。他に"bilinear"も選択できる。また、filenameとformatを指定することで、ASCII ラスター形式で結果をファイルに保存している。

上記のうち、下の3行を繰り返すことで、他のラスターファイルも同じようにして範囲とセルサイズを調整します。完了したら作成した.asc形式のファイルをメモ帳などで開き、先頭の5行（NCOLS, NROWS, XLLCONER, YLLCORNER, CELLSIZE）が完全に一致しているか確認してください。できていれば新しいディレクトリを作り、これらのファイルをそこに移動させます。

3. MaxEnt解析の実行

MaxEntを起動し、サンプルに分布データ（.csv）、環境レイヤに.ascファイルを入れたディレクトリを指定します。下の枠にインプットしたデータのリストが表示されます（図5）。

図5

このうち、環境レイヤの底質データ（ここではclipped_sediment）はカテゴリカルデータなので、「Categorical」に設定します。他は「Continuous」のままで大丈夫です。最後に出力ディレクトリを指定し、「RUN」をクリックします。問題なければこのディレクトリに結果のファイルが出力されるはずです。このうちの.htmlファイルを開くと、モデルの評価、環境データごとの寄与度、予測された環境適地の地図などが確認できます。また、出力された.ascファイルをQGISに読み込むことで生息適地の地図を表示することもできます（図6）。

図6

参考文献

1. 田口貴史・石崎大介・根本守仁、令和元年(2019 年)6 月の琵琶湖南湖におけるブルーギルの分布、外来魚駆除対策研究（https://www.pref.shiga.lg.jp/file/attachment/5238608.pdf）

湖沼図から底質ラスターデータの作成

前回の記事では湖沼図の等深線を基に内挿を行うことで水深ラスターデータを作成しました。今回は同じ湖沼図から底質ラスターデータを作成します。

底質データの作成は以下の4つのステップからなります。

1. 湖沼図の入手
2. ジオリファレンス
3. 底質データの書き込み
4. 内挿

このうち1と2については、前々回の記事を参照してください。ここでは3.底質データの書き込みと4.内挿について説明します。

3. 底質データの書き込み

QGISにジオリファリンス済みの地図を表示してください。等深線を書き込んだときは線を用いましたが、底質データは湖沼図上ではある地点の底質をアルファベットで示しているので点で書き込みます。なお、それぞれのアルファベットの意味はここにまとめました。

やることは単純で、各アルファベット上に点を打ち、底質の情報を入力するだけです。

まずは点を打つレイヤを準備します。画面左上にある「新規シェープファイルレイヤ」をクリックします（フロッピーディスクのようなアイコンの真下にあるやつ）。入力画面が開くので、ファイル名を入力し（ここでは、「sediment.shp」とした）、ジオメトリタイプで点（ポイント）を選択します。そして、新規フィールド欄の名前を「Sediment」とし（任意ですが、ここではSedimentとしました）、データ型を「テキストデータ」にしてから「フィールドリストに追加」をクリックします。すると下のフィールドリストに「Sediment」というフィールドが追加されます。下の「OK」をクリックすると「sediment」というレイヤが追加されます。

追加されたレイヤを選択し、「編集モード切替」をクリックし、続いて「点地物を追加する」をクリックします（図1）。

図1

その状態で地図上の底質を示すアルファベットにカーソルを合わせてクリックします（例えば図1の赤い円の中の「M」）。入力画面が現れるので、「Sediment」にその地点のアルファベットを入力してOKをクリックします。すると地図上に点が表示されます。これを繰り返して地図上の全ての底質データをレイヤに入力します。

4. 内挿

前回は水深という数値データを基に内挿を行いましたが、今回は文字データを基に内挿を行います。そのため、今回は「Multilevel b-spline interpolation for categories」を用います。「プロセシングツールボックス」の検索窓にこれを入力して選択します。

パラメータ入力画面が現れるので、以下のように設定して「実行」を押します（図2）。

• Points: sediment（レイヤ名）
• Attribute: Sediment （内挿を行う属性）
• Output extent: 「レイヤから計算」>地図のレイヤを選択（内挿の範囲指定）
• Cell size: 0.001（数字が小さいほど解像度が高い）
• Fit: 0（よくわからないが、試した範囲では「1」にしても結果はほぼ変わらなかった）
• Categories: sediment_cell001（アウトプットファイル名。任意）
• Probability: sediment_cell001_prob（確率のアウトプットファイル名。任意）
  

  

  図2

  
  

完了すると図3のような内装された底質レイヤが得られます。以上で底質ラスターデータの作製は完了です。

図3

＊エラーが出て内挿ができない場合

QGISは内部でSAGAというプログラムを呼び出して「Multilevel b-spline interpolation for categories」を実行しています。しかし、PCにインストールされているSAGAのバージョンがQGIS推奨のバージョンと一致しないとエラーになることがあります。使っているOSがLinuxの場合、以下のコマンドで直接SAGAによる内挿を実行できます（ここでは上記のパラメータ設定と同じ設定をしています）。コマンドを実行してできたファイルのうち拡張子が.sdatとなっているものをラスターとしてQGISに読み込ませれば結果を見ることができます。

saga_cmd grid_spline 7 -POINTS Biwa_Sedimentpoint_shape.shp -FIELD Sediment -TARGET_USER_SIZE 0.001 -TARGET_USER_XMIN 135.853200 -TARGET_USER_XMAX 136.006200 -TARGET_USER_YMIN 34.974000 -TARGET_USER_YMAX 35.173800 -TARGET_USER_FITS 0 -CATEGORIES sediment_cell001 -PROPABILITY sediment_cell001_prob

コマンド先頭の「saga_cmd grid_spline 7」はこの解析が「Multilevel b-spline interpolation for categories」であることを指定しているので変更する必要はありません。「-POINTS」は内挿を行うシェープファイル、「-FIELD」はQGISの入力画面で「Attribute」とされているもの、「-TARGET_USER_SIZE」は「Cell size」、「-TARGET_USER_XMIN」から「-TARGET_USER_YMAX」はQGISの「Output extent」で指定したものを別々に指定、「-TARGET_USER_FITS」は「Fit」、「-CATEGORIES」は「Categories」、「-PROBABILITY」は「Probability」をそれぞれ意味しています。一度QGISで内挿を走らせるとログが出力されるのでそれを参考にすると設定しやすいと思います。

湖沼図からの水深ラスターデータの作成（後編）

 前回の投稿でQGIS上で湖沼図に位置情報を埋め込むことができました。今回はその湖沼図上に等深線を書き込み、それを基にしてラスターデータを作成します。

3. 等深線の書き込み

湖沼図の等深線をなぞります。等深線が明瞭になるように図を拡大してください。マウスのスクロールで拡大・縮小が可能です。また、マウスのスクロールホイールをクリックしながら動かすことで地図上を移動することができます。これは覚えておくと便利です。

画面左上の「新規シェープファイルレイヤ」アイコンをクリックします（図1の赤丸）。すると「新規シェープファイルレイヤ」ウィンドウが表示さます。

図1

ファイル名を入力し（任意ですが、半角英数字にしてください。ここでは「Depth.shp」にしました。）、ジオメトリタイプを「ラインストリング」にします。さらに「新規フィールド」の「名前」欄を「Depth」とし（これも任意ですが、わかりやすいほうが良いです。）、「データ型」を「小数点付き数値」にします。できたら「フィールドリストに追加」をクリックし、下のフィールドリストに追加されたことを確認してください。なお、フィールドリストには最初から「id」というフィールドが入っていますが、必要ないので削除しても構いません。ここまで入力したら「OK」をクリックします。元の画面に戻りますが、左のレイヤ欄に「Depth」というレイヤが追加されたはずです。

追加したレイヤに線を書き込みます。「Depth」レイヤが選択された状態（青く反転している状態）で、「編集モード切替」アイコン（「設定」の下あたりにある黄色い鉛筆のアイコン）をクリックし、次にその2つ右にある「線の地物を追加」をクリックします（図2）。するとカーソルが十字を◯で囲ったような形になったと思います。この状態で地図上で左クリックをするとその場所に頂点が置かれます。さらに別の場所で左クリックするともう一つ頂点が置かれ、最初の頂点との間に線が引かれます。これを続けることで等深線に沿って線を引いていきます。引き終えたら右クリックします。すると地物属性というウィンドウが立ち上がります。ここにはレイヤを作成したときに追加したフィールド（ここでは「Depth」）が表示されるので、今線を引いた場所の水深を入力し、「OK」を押します。この作業を繰り返すことで等深線をシェープファイルレイヤ上に入力していきます。

図2

一度引いた線は「線の地物を追加」の右隣りにある「頂点ツール」で修正することが可能です。「頂点ツール」を選択するとカーソルが十字マークに変わります。この状態で引いた線の上に現れる☓印をクリックするとその場所に新しい頂点を追加することができます。また、すでにある頂点をクリックすることでその頂点の位置を変えることができます。効率的に等深線を引く上での私のおすすめは、とりあえずざっくりと最初から最後まで線を引き、線の形は後から「頂点ツール」で整えることです。端から少しずつ丁寧に線を引いていくと終わる気がしません。

一通り線を引き終えると図3のような感じになります。これでも全部引いたわけではありませんが、面倒なのでとりあえず今回はこれで良しとしました。なお、線の色や太さなどはレイヤを右クリック>「プロパティ」>「シンボロジ」で変更可能です。作成したレイヤは保存しておきましょう。レイヤを右クリック>「エクスポート」>「新規ファイルに地物を保存」でファイル名を入力して「OK」を押します。形式は「ESRI Shapefile」のままで構いません。この形式で保存すると複数の形式のファイルが指定フォルダに現れますが、レイヤを開くときは拡張子がshpになっているものを選んでください。また、このとき「OK」の左にある「保存したファイルを地図に追加する」にチェックが入っていると保存したものと同じものがレイヤに追加されるので注意してください。ファイルとして保存すると、レイヤを削除しても「ベクタレイヤを追加」から保存したshpファイルを選択することで再びレイヤを開くことができます。

図3

4. 内挿

それではいよいよ等深線と等深線の間の深度を補間していきます。「プロセシング」>「ツールボックス」をクリックしてください。右側に「プロセシングツールボックス」というパネルが現れるので、そこの検索窓に「tin interpolation」と入力します。すると「TIN内挿（不規則三角網）」が表示されるので、それをダブルクリックします。ダイアログボックスが表示されるので、「入力ベクタ」に作成した等深線のレイヤを指定し、「内挿対象の属性」で「Depth」を選択してから、すぐ下の緑の「＋」マークをクリックします（図4）。

図4

入力レイヤに関するテーブルが現れるので、「データ型」として「ストラクチャーライン」を選択します。「内挿方法」を「線形」、「領域」は右の「...」マークをクリックして「レイヤから計算」から最初に読み込んだラスターデータを選択します。「ピクセルサイズ」は0.001にします。また、「出力レイヤ」に出力先ファイル名を入力します（ファイル名は任意。ここでは depth_pic0.001.tif ）。これにより、自動的にファイルが保存されます。ここまでセットできたら右下の「実行」をクリックします。完了するとレイヤパネルにラスターレイヤが追加されます（図5）。デフォルトでは水深の深いところほど白く表示されていると思います。

図5

以上で湖沼図からの水深ラスターデータの作成は完了です。

（おまけ） 傾斜量図の作成

作成した水深ラスターデータから傾斜量図を作成できます。「プロセシングツールボックス」の検索窓に「r.slope」と打ち込むと「r.slope.aspect」が出てくるので、これをダブルクリックします。ダイアログボックスが表示されるので、「Elevation」の欄で作成した水深レイヤを指定し（ここでは depth_pic0.001）、「実行」をクリックします。完了すると複数のレイヤが作成されます。このうち、「Slope」が傾斜，「Aspect」が傾斜方位，「Profile curvature」が断面曲率，「Tangential curvature」が接曲率を示しています。

＊「Exception: unknown」と表示され、TIN内挿が動かない場合

自分の環境（Ubuntu 20.04）で上記のようなエラーに遭遇しました。GUI上ではどうにもなりませんでしたが、CUIで実行できました。コマンドは以下の通りです。オプションの設定などはエラーが出たときのログを参照してください。以下のコマンドは上のダイアログボックスの設定と同じです。

qgis_process run qgis:tininterpolation --INTERPOLATION_DATA=Depth.shp::~::0::~::1::~::1 --METHOD=0 --EXTENT=135.833234096,136.000092079,34.974896348,35.166801988 --PIXEL_SIZE=0.001 --OUTPUT=depth_pic0.001.tif

湖沼図からの水深ラスターデータの作成（前編）

 現在では標高についてのデジタルデータはかなり充実してきていますが、湖沼のデータについては利用できるものが限られています。本記事では国土地理院が提供している湖沼図から水深のラスターデータを作成する方法を説明します。ただし、この方法は私が手探りで調べたことをまとめたものですので、もっと良い方法を知っている方がいたら是非教えていただければと思います。

水深データの作成は以下の4つのステップからなります。

1. 湖沼図の入手
2. ジオリファレンス
3. 等深線の書き込み
4. 内挿

1では水深データの基となる湖沼図を国土地理院から入手します。入手した湖沼図はただの画像データなので、2により位置情報を埋め込みます。3は一番手間のかかるステップで、湖沼図に描かれた等深線をなぞって線を描いていきます。現状これを自動化する良い方法が見つかりません。誰か教えてください。本当に。最後に、4では作成した等深線からラスターデータを作成します。以下では琵琶湖の南湖を例として詳しい手順を説明します。

1. 湖沼図の入手

地理院地図にアクセスし、目的の湖沼図を入手します。アクセスしたら最初は標準地図が表示されているはずです。この状態で目的の場所（ここでは琵琶湖の南湖）が画面の中央にくるように移動します。

次に湖沼図へ切り替えます。画面左上にある「地図」と書かれたアイコンをクリックすると、左側にパネルが現れます。ここで「地図の種類」の欄から、「土地の成り立ち・土地利用」>「湖沼図・湖沼データ」>「湖沼図」とクリックしてください。すると「選択中の地図」の欄に湖沼図が追加されたと思います。ここまで来たら標準地図は消したほうが見やすくなると思うので、「☓」をクリックして消しましょう。また、画面右上の「設定」から「中心十字線」を「OFF」にすることで余計なものを省くことができます。

　

重要な点は、同じく「設定」>「グリッド表示」から「緯経度グリッド」を「ON」にしておくことです。これは次のジオリファレンスのステップで利用します。ここまでくると図1のような画面になっていると思います。

図1

最後にこの地図を保存するのですが、注意しなければならないのは地図の解像度が現在表示されている解像度で保存されることです。つまり解像度の高い地図が必要なら、できるだけ拡大してから保存する必要があります。しかし一方で、保存する範囲が広いと解像度を高くしすぎることでデータ量が大きくなり、保存できなくなります。このあたりは必要に応じて調整してください。

保存は「設定」の左にある「共有」から行います。

図2

ここでtwitterアイコンの右にある山のような形のアイコンで画像の保存ができます。保存方法ですが、「範囲を固定」をクリックし、保存したい範囲を指定します。範囲の指定は画面に現れた赤枠でするか、直接数値を入力することで行います。ここでは南湖全体をカバーするため、図2のように指定しました。指定したら「OK」をクリックし、適当な名前を付けて保存します。

2. ジオリファレンス

地理院から入手した地図は現状ではただの画像データです。これに位置情報を埋め込みます。ここではフリーで利用可能なQGISを使ってジオリファレンスを行います（ここではバージョン3.22を使っています）。QGISを開いたら「ラスタ」>「ジオリファレンサ」を選択します。するとジオリファレンサのウィンドウが開くので、「ファイル」>「ラスタを開く」を選択し、先ほど保存した湖沼図を開きます。そして、緯度経度のグリッド線の交点が明瞭になる程度に画像を拡大します。

図3

画面上部の黄色い歯車のマークのすぐ隣にある「点を追加」が選択されていることを確認したら、画像の角あたりにある適当なグリッド線の交点をクリックします（ここでは図3の◯で囲んだ位置）。すると図3のように「地図座標の入力」というウィンドウが開くので、クリックした交点の座標を10進数で入力します。湖沼図に表示されているのは60進数なので、ここのサイトなどで変換してください。入力欄の下は投影法の選択欄です。ここは「EPSG4326」にしておきます。

同じ手順で画像の四隅それぞれに位置情報を入力します。入力し終えたら、左上の緑の三角形のマーク（あるいは「ファイル」>「ジオリファレンスを開始」）をクリックします。「変換の設定」というウィンドウが開きますが、「出力ラスタ」の欄で出力ファイルの名前を指定し、下の「完了時にQGISにロードする」にチェックが入っていることを確認したら「OK」を押します。うまくいけば画面の上部にジオリファレンスが成功したとの表示が出るはずです。出なければ再度三角形マークを押してみましょう（なぜか自分の環境ではもう一度押さないと出ない…）。できたらジオリファレンサのウィンドウを閉じます。元のQGISのウィンドウに湖沼図がロードされているはずです。ここで湖沼図にマウスカーソルを持ってくると、その位置の緯度経度がウィンドウ下部の座標欄に表示されているでしょう。

念の為、座標が合っているかを確認します。国土数値情報から湖沼図のシェープファイルをダウンロードし、解凍したら、「レイヤ」>「レイヤを追加」>「ベクタレイヤを追加」で解凍したシェープファイル（.shp）をロードします。その結果、レイヤにシェープファイルが追加され、図4のように湖沼図の上に琵琶湖の領域が色付きで示されると思います。この色付き部分と湖沼図上の琵琶湖の位置が一致していたらジオリファレンスは成功です。

図4

長くなるので残りは後半に続きます。

Rによる主成分分析（作図編）

少し間が空きましたが、前回の投稿（Rによる主成分分析）ではggbiplotを用いて主成分分析結果を作図しました。今回はやり方を変えて、グループごとに枠で囲む方法を説明します。使用するデータは前回と同様です。

データの準備

プロットする前にデータを少し編集する必要があります。前回に引き続き、xという変数にアヤメのデータ、pcaという変数にその主成分分析結果が入っているとします。以下にそれぞれの要約を示してみます（図1）。

図1

次にプロット用データを作成し、変数df.pcaへ代入します。今回は第一主成分と第二主成分でプロットするとします。pcaから第一主成分と第二主成分の主成分得点、xからそれぞれの行に対応する種名を抜き出してデータフレームにまとめています。ここで、左辺のPC1、PC2、Speciesは任意です。[]内の数字を3や4にすることで、第三主成分、第四主成分を抜き出すこともできます。

df.pca <- data.frame(PC1=pca$x[,1], PC2=pca$x[,2], Species=x$Species)

グループ分け

グループ分けにはdplyrライブラリを使います。なければインストールしてください。

library(dplyr)
hull <- df.pca %>%
+ dplyr::group_by(Species) %>%
+ dplyr::do(.[chull(.[c('PC1','PC2')]),])

少々ややこしいですが、打ち間違えないよう注意してください。

プロット

プロットにはggplot2を用います。まずデフォルトでプロットしてみましょう。一度プロットを変数pへ代入してから表示します。

library(ggplot2)
p <- ggplot(df.pca, aes(x = PC1, y = PC2, shape = Species)) + geom_point()
p

これで以下のような図が表示されると思います（図2）。

図2

ggplot2の良いところはこの変数pに書き加えていくことで図を編集することができることです。もっと表示をシンプルにしてみましょう（図3）。

p + theme_test()

図3

次にプロットのシンボルを変更します（図4）。

p + theme_test() + scale_shape_manual(values=c(1,2,3))

図4

このような感じで図を編集することが可能です。一度ここまでの編集内容を変数p2へ保存しておきます。

p2 <- p + theme_test() + scale_shape_manual(values=c(1,2,3))

その他、編集方法についての詳細はggplot2のドキュメントを参照してください。

各グループを枠線で囲う

それではこれまで作った図に枠線を追加します。グループ分けのところで作成した変数hullを利用します（図5）。なお、alphaは枠内の塗りつぶしの透過度で0が透明です。colorは枠線の色を示します。1は黒です。

p2 + geom_polygon(data=hull, alpha=0, color=1)

図5

最後に、全ての枠線が実線というのもアレなので、グループごとに線の種類を変えましょう（図6）。

p2 + geom_polygon(data=hull, alpha=0, color=1, aes(linetype=Species))

図6

以上です。

PopARTによるネットワーク解析

 系統地理学の分野では、ミトコンドリアの各ハプロタイプの関係をネットワークで表現することがあります。ネットワークの作成によく使われるソフトとしてはTCS(http://bioresearch.byu.edu/tcs/)がありますが、きれいな図をつくるために色をつけるといったことはできません。

ここではPopART(http://popart.otago.ac.nz/documentation.shtml)を用いたネットワーク解析の手順を紹介します。

まずインプットファイルとしてNexus形式のアラメントデータを準備します。すでにFasta形式のデータがあればMEGA(https://www.megasoftware.net)などで簡単に変換できます。

続いてPopARTのアイコンをダブルクリックして起動します。場所は研究室のWindowsパソコンの「デスクトップ」 > 「Phylogenetic_Analyses_TOOLs」 >  「popart-1.7」フォルダの中です。

起動後、画面左上の「nex」と書かれたアイコンをクリックしてインプットファイルを選択してください。うまくいけば画面左下に「Loaded file (ファイル名)」と表示されるはずです。この状態で左のパネルの「Data Viewer」の下にある「Alignment」タブを開くとインプットした配列を見ることができます。

ファイルのインプットが終わったらネットワーク解析に移ります。画面上の「Network」から解析に用いる手法を選ぶことができます。ここでは左下の「TCS Network」を選択します。作成されたネットワークは以下のようになります。

各円がハプロタイプで、それぞれがネットワーク状に繋がっているのがわかるかと思います。ハプロタイプを繋ぐ枝の上にある棒の数はハプロタイプ間で異なる塩基の数です。ちなみに、「Statistics」から塩基多様度などの指標を出力することも可能です。

解析自体は以上で終了ですが、PopARTではネットワークを好みに合わせて編集することができます。円の色を変えるには、「Edit」 > 「Set default vertex color」をクリックします。「Edit」 > 「Set default vertex size」では円のサイズを変更することができます。また、ネットワークの形そのものは直接円をドラッグ&ドロップすることで変えることができます。好みの図ができたら「File」 > 「Export graphics」から保存しましょう。PopARTの一通りの使い方は以上です。

グループ分けの全ての可能な組み合わせ

複数の要素をグループ分けする際の可能な全ての組み合わせが知りたいことがありました。例えば要素A,B,Cがあった場合、グループ分けの全ての組み合わせは以下の5通りになります。

(ABC), (AB)(C), (AC)(B), (BC)(A), (A)(B)(C)

要素の数が4つまでであれば一つ一つ自力で書いていってもいいのですが、5つ以上になると組み合わせが50を超え、面倒な上に間違う可能性が高くなってしまいます。何かいい方法はないかとググってみたものの、探し方が悪いせいか見つかりませんでした。いまいちスマートではないのですが、何とか方法を思いついたので以下にそのPythonスクリプトを載せます。この方法ではitertoolsとmore_intertoolsというライブラリを使います。itertoolsはpythonに標準で入っていますが、more_itertoolsの方はpipかcondaでインストールする必要があります。

import itertools
import more_itertools

d = ['A','B','C']
n = len(d)
l = list(more_itertools.powerset(d))

m = []
for h in l[1:]:
    m.append(''.join(list(h)))
    
r = []
for i in range(1,n+1):
    r += list(itertools.combinations(m, i))
      
print([j for j in r if sorted("".join(list(j))) == d])

手順としては、最初にmore_itertools.powerset(d)で、各要素の全部分集合を生成します。次にグループ数ごとに、itertools.combinations(m, i)を使って部分集合の全ての組み合わせのリストを作成します。当然このリストの中には要素の数が異なっていたり重複が含まれていたりするので、最後にフィルタリングをかけます。このスクリプトを実行した結果は以下のようになります。

#要素がA,B,Cのとき
[('ABC',), ('A', 'BC'), ('B', 'AC'), ('C', 'AB'), ('A', 'B', 'C')]

#要素がA,B,C,Dのとき
[('ABCD',), ('A', 'BCD'), ('B', 'ACD'), ('C', 'ABD'), ('D', 'ABC'), ('AB', 'CD'), ('AC', 'BD'), ('AD', 'BC'), 
('A', 'B', 'CD'), ('A', 'C', 'BD'), ('A', 'D', 'BC'), ('B', 'C', 'AD'), ('B', 'D', 'AC'), ('C', 'D', 'AB'), 
('A', 'B', 'C', 'D')]