前回Structure→Structure Harvester→CLUMPP→distruct の流れに失敗したという記事を書きましたが、今回はその続きです。STRUCTURE の結果をきれいに出すことだけを目標に書いていきます。
今回、K 値はEvanno et al., 2005 のΔK を用いて求め、さらに採用するSTRUCTURE の結果は得られた結果の中で最も尤度が大きかったものを選択しました。念のため、尤度がとりわけ大きかった結果すべてでバープロットを表示したところ、すべて結果が一致したのでおそらくそれ自体は問題ないでしょう。たぶん…。
ということで本題。なんできれいなバープロットにこだわっているかというと、STRUCTURE のソフトで出してくれるバープロットがjpg で、しかもものすごく画質が悪いから。論文には到底使えません。使っているものもあるみたいですが…。
今回使うソフトはＲ のパッケージpophelper です。
なぜかinstall.package が使えなくて序盤から苦労しました。
pophelper/README.md at master · royfrancis/pophelper · GitHub
に載っている詳細のサイトを見たら、全然普段のパッケージと違って道理で…という感じ。こういうことがあるんだなぁ。インストールで困った時はソフトウェアのREADME を読んでみるといいみたい。
まず周辺パッケージのダウンロード。詳細サイトにはMac は、と書いてあったがWindows でもこれをやらないとインストールできませんでした。
# install dependencies
>install.packages(c("ggplot2","gridExtra","label.switching","tidyr","remotes"),repos="https://cloud.r-project.org")
つぎにremote パッケージを使ったGithubからのpophelper のインストール
# install pophelper package from GitHub
>remotes::install_github('royfrancis/pophelper')
ここまできたらいつものようにライブラリの読み込み
# load library for use
>library(pophelper)
これ以降の処理はこのページを参考にしています。
Pophelper 2.3.1 • pophelper (royfrancis.com)
インプットには STRUCTURE run files、TESS 2.3 run files (-TR files in folders)、 BAPS、 BASIC (simple delimited files) 、CLUMPP associated file が使えるみたいです。今回はタイトル通りSTRUCTURE を使います。ほかのインプットのやり方も上記サイトにあります。
【追記】structure の結果ファイルだよ、ということを表現するために「hoge.structure」という名前に勝手に変えていますが、実際に使っているファイルの名前は「ipyrad_NJ_run_run番号_f」みたいなやつです。わかりにくくしてしまってごめんなさい。structureHarvester の結果を見るなどして選んだK 値の結果の中で、一番LnP(K) が大きいファイルを選択しています。
>structure <- "hoge.structure"
インプットの形式の変換
>readQ(structure)
>sfiles <- "hoge.structure"
>readQ(files=sfiles)
>slist <- readQ(files=sfiles)
>readQ(files=sfiles,filetype="structure")
出力クラスチェック
> class(slist)
[1] "list"
最初に読み込まれているかの確認のため、変換されたファイルの先頭を表示する
> head(slist[ [1] ])
Cluster1 Cluster2 Cluster3 Cluster4 Cluster5 Cluster6
1 0.002 0.001 0.002 0 0.001 0.001
2 0.001 0.000 0.000 0 0.001 0.001
3 0.001 0.000 0.001 0 0.001 0.001
4 0.002 0.001 0.001 0 0.001 0.001
5 0.001 0.000 0.001 0 0.001 0.001
6 0.001 0.000 0.001 0 0.001 0.001
Cluster7 Cluster8 Cluster9 Cluster10 Cluster11
1 0.001 0.987 0.001 0.002 0.002
2 0.001 0.995 0.001 0.001 0.001
3 0.001 0.995 0.001 0.001 0.000
4 0.001 0.989 0.002 0.001 0.001
5 0.001 0.994 0.001 0.001 0.001
6 0.001 0.995 0.001 0.001 0.001
こんな感じの結果が出てくる
もしデータが信頼区間を含んでいた場合
>slist <- readQ(files=sfiles, readci=TRUE)
で読み込めるらしい。今回は含んでいないので無理でした。
データが個体番号を含んでいた場合、どの個体がどのクラスターに入るのかも表示できる。
> slist <- readQ(files=sfiles, indlabfromfile=T)
>head(slist[ [1] ])
Cluster1 Cluster2 Cluster3 Cluster4 Cluster5
101_pair 0.002 0.001 0.002 0 0.001
102_pair 0.001 0.000 0.000 0 0.001
4413_pair 0.001 0.000 0.001 0 0.001
104_pair 0.002 0.001 0.001 0 0.001
105_pair 0.001 0.000 0.001 0 0.001
106_pair 0.001 0.000 0.001 0 0.001
Cluster6 Cluster7 Cluster8 Cluster9 Cluster10
101_pair 0.001 0.001 0.987 0.001 0.002
102_pair 0.001 0.001 0.995 0.001 0.001
4413_pair 0.001 0.001 0.995 0.001 0.001
104_pair 0.001 0.001 0.989 0.002 0.001
105_pair 0.001 0.001 0.994 0.001 0.001
106_pair 0.001 0.001 0.995 0.001 0.001
Cluster11
101_pair 0.002
102_pair 0.001
4413_pair 0.000
104_pair 0.001
105_pair 0.001
106_pair 0.001
こんな感じ。
どうも上記のページを見ているとgene helper は同じKで実行した複数のランの結果をまとめたりもできるみたいなんですが、今回は時間の都合上割愛します。
あとで暇ができたら書くかも。
ほかにも色々なファンクションが紹介されているのでSTRUCTURE気になる人は読んでみてね。
いよいよプロッティングです。
> hoge.plot <-plotQ(slist,imgoutput="sep",showlegend=T,showtitle=T,titlelab="Hoge structure", height=20,barbordercolour="white",barbordersize=0,outputfilename="plotq",imgtype="pdf",exportpath=getwd())
Drawing plot ...
書けましたよ、確かに書けましたけれどもさ…
なんだか異様に細長い。
barsizeでバーの太さを変更できるのですが、最大の1 でも全画面表示したら素麺みたいに細い。
ので、普通にバープロットしてウインドウに表示し、それを保存するよくあるやり方で保存することにしました。
> hoge.plot <- plotQ(slist,imgoutput="sep",returnplot=T,exportplot=F,basesize=11)
>hoge.plot
#basesize: base font size, given in pts.
#returnplot: 説明読んでもよくわからないんですが、これを入れないと結果がnull になります。
#exportplot: 検索しても出てこなかったけど、たぶん結果をエクスポートするかどうか、だと思われる
で、あとはウインドウに表示されたものを好きな形式と大きさで保存するだけ！
画質命なのでW860, H476、あとからillustrator で色を編集できるようにsvg で保存。
やった～～～～～～～～～！！！！！！！！！！！！！！！！！！！！！
あと、plot の()の中にshowlegend=T やbarbordercolour="white",barbordersize=0.01 などを追加して、レジェンドを表示したり成分間の境界線を引いたりできます。
plotQ のオプション一覧。
An overview of the components of a plotQ() figure and most of the arguments used to modify them. - Pophelper 2.3.1 • pophelper (royfrancis.com) より引用
おしまい
【追記】うわ、これ縦軸がないじゃないのよ。と思って調べていたら、pophelper のサイトにヒントがありました。
軸はあとからgridExtra パッケージの中のgrid.arrange 関数を使って付け加えるらしい。
これは複数のデータをまとめて表示するのに使える関数みたい。
参考：ggplot2で複数のプロットをまとめて表示する方法あれこれ - Qiita
> library("gridExtra")
> hoge.plot <- plotQ(slist,imgoutput="sep",returnplot=T,exportplot=F,basesize=11,showyaxis=T)
さっきの描画用のスクリプトに太字部分を付け足して実行した後、
> grid.arrange(hoge.plot$plot[ [1] ])
これで
こんな感じの図が描ける。
