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

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

底質データの作成は以下の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')]

GenBankファイルから遺伝子を切り出すPythonスクリプト

これまでGenBankファイルから遺伝子を切り出すため、とあるPerlスクリプトを使っていたのですが、Perlのアップデートに伴って使えなくなっていました。そこで練習がてらGenBank形式のファイルからタンパク質遺伝子の塩基配列を取り出すPythonスクリプトを作成しました。

以下よりダウンロードできます。

GBstrip.py

使い方:

python GBstrip.py inputfile

*Python version 3.8での動作を確認しています。なお、inputfileの名前は任意です。

追記（2021/12/13）
GBstrip.pyを改良し、タンパク質遺伝子以外にもtRNAとrRNAも切り出せるようにしました。また、複数のaccessionがインプットのGenBankファイルに含まれる場合、accessionごとに別々のファイルに出力するようにしました。