研究メモブログの目次

研究メモ

プログラミング初心者が研究で使ったソフトウェア(主に遺伝子解析)や情報などを綴るブログです。

Structure解析とその周辺ソフトウェアの使い方(注意:CLUMPP実行までで終わっています)

Structure 解析は、 このサイト(How to STRUCTURE)を見れば万事オッケー!と思っていたわたしですが、いざやってみたらこのサイトで紹介されているインプットファイルがうまく使えませんでした。たぶんわたしがポンコツなせいなのですが。

したがって、あきらかに対数尤度LnP(D)が特定のK値で大きくなっていてばらつきも小さい場合や、上記サイトで紹介されているインプットファイルの使い方がわかる場合は以下の説明を読む必要はありません。上記サイトがとても親切ですのでそちらへどうぞ!

 ちと時間がないので、Structure 解析の結果を出すところまでは上記サイトで説明があるのでそちらで見てください。

今回自分の扱ったデータはK値ごとの対数尤度がばらついていてどのK が尤もらしいのかわかりませんでした。

Structure 解析の結果をまとめてくれるソフトウェアがあるらしいという話は研究室の後輩のK 君にちらっと聞いていたのでそのソフトウェアを調べてみた。

CLUMPP: CLUster Matching and Permutation Program

あまり英語は得意でないので日本語版チュートリアルがないかな~、と最初に検索してみた結果、紹介してくれているサイトがありました。

Structureとその仲間たち - サイコロにもてあそばれる日々

更新が2012年なのでやや不安でしたが、これを参考に大体できたので、自分がやったことを以下で紹介していきます。

※自分の場合は、同じK値での標準偏差が大きいのでStructure をもう一度かけ直す必要があるのですが、結局そのあともう一回同じ工程を繰り返すと思うので、その前にここに方法をメモしています。結果を見て、そもそもパラメータが収束してなくない!?と思った人がいたら、そういうことなので大目に見てくださいまし。

 

一連の流れとしては

Structure Software: V. 2.3.4

↓  最適K 値決定 & CLUMPP 用インプット作成

Structure Harvester

↓ 

CLUMPP: CLUster Matching and Permutation Program

↓ 図示

distruct: graphical display of population structure

 

という感じです。

 

まず、上記サイトで紹介されているStructure Harvesterに行きます。

このサイトのInstructions を見てもらえればわかるのですが、input はStructure のResults をZip. 化したものです。

Structure を実行すると、プロジェクトが入っているフォルダの中にもう一つフォルダが作られます。

f:id:nemunemu_nyanko:20210102193512p:plain

この画像でいうとNJというフォルダです。

f:id:nemunemu_nyanko:20210102193713p:plain

この中に入るとResults というフォルダがあります。このフォルダを右クリックでZip 化します。Structure 使い慣れている人にはわかりきった話だと思うのですが、わたしはこの段階で「Resultsどこ!?」となって数分失ったので一応書いておきます。

Zip ファイルができたらサイトの

f:id:nemunemu_nyanko:20210102193924p:plain

参照のところからアップロードします。

するとExampleと同じようなHTML が出てきます。

下に表示される表を見てどのK 値が最適そうか考えましょう。

f:id:nemunemu_nyanko:20210102211648p:plain

複数K 値でL(K) が大きくて困ったのですが、An overview of STRUCTURE: applications, parameter settings and supporting software. (p15) を見ると

複数のK値が似たようなLn Pr(X|K)の推定値を持つ場合には、最小のK値がデータに最も適していると合理的な確信を持って推測することができます(DeepL 訳)」

とあるので、今回はK=5のようです。このグラフの下に各K についてCLUMPP 用ファイルが保存できるリンクがあるのでここで保存したいK値のファイルを保存します。

【追記】 これなんですが、Evanno et al., 2005 ΔK を用いた方がいい可能性があります。ΔK はStructure Harvester の結果の中にあるので、L(K) にばらつきがあるときはそれも見てみましょう。それから、このΔK も、STRUCTURE での計算を増やさないとあまり参考にならない値になってしまうので、最初にも書きましたが、L(K) がばらついていた場合、MCMC を増やすか、K ごとの計算回数を増やすかして再解析したほうがいいです。

 

ちなみにですが、全結果のダウンロードはHTML ページの一番上からできます。

f:id:nemunemu_nyanko:20210102212146p:plain

 

あと、これはおまけですが、このHTMLの一番下に、Structure Results をExcel 化できる生データ一覧が出てくるので、コピーしてテキストファイルとして保存しておくとかなり便利だと思いました。

Structure が出してくるサマリー.txtはExcel で取り込める形式ではないので(ほかにいい保存形式があるかもしれませんが…)。

f:id:nemunemu_nyanko:20210102212257p:plain

こんな感じです。

 

CLUMPP 用ファイルをダウンロードしたらCLUMPP ソフトウェアをインストールしないといけないのですが、インストール方法は公式サイトの一番初めに書いてあるので省略します。

tar.gz を解凍&実行するだけです。

 

CLUMPP とは

CLUMPP は、R 個の複製から得られた Q 個の行列間の類似度尺度 G を最大化するプログラム (An overview of STRUCTURE より)

K :グループの数を表す整数

C :調査される集団または個体(個体行列 - DATATYPE=0; 母集団行列 - DATATYPE=1)

R:整列されたQ行列または実行の数を表す整数

Q 行列:個人または集団の数に対応するC行と、個別の仮定された集団グループに対応するK列を持つK値(STRUCTUREにおける仮定された集団の数)の分析から得られたC x K係数行列

H:G  の値を使用して求められるペア比較によって類似度の中央値

H’:類似度の 2 番目の関数 G'  から計算される類似度の中央値

 

以下の使い方は解凍したCLUMPP フォルダの中にあるCLUMPP_Manual.pdf を拾い読みした内容なので、疑わしいと思ったら自分で読んでみてください。

 

まずはじめにインプットファイルであるparamfile を作成します。これがとても面倒くさいです。必要最低限のパラメータのことしか書かないのでadditional option については各自調べてください。よければ教えてください(他力本願)。

で、特にこのブログを参考にしている初心者の方々には一からファイルの中身を書くのは難易度高めだと思うので、解凍フォルダの中のexample の中のparamfile を書き換えるといいです。

 

確実に必要なパラメータは以下

Parameter settings
----------------------- Main parameters -------------------------
DATATYPE = 0
INDFILE = Hoge.indfile
POPFILE =
OUTFILE = Hoge.outfile
MISCFILE = Hoge.miscfile
K = 5
C = 67
R = 10
M = 1
W = 0
S = 1
- Additional options for the Greedy and LargeKGreedy algorithms -

ここはどのアルゴリズムを選んだかで変わります。M = 1ならいらない。
GREEDY_OPTION =
REPEATS =
PERMUTATIONFILE =
----------------------- Optional outputs ------------------------
PRINT_PERMUTED_DATA = 1
PERMUTED_DATAFILE = 1
PRINT_EVERY_PERM = 0
EVERY_PERMFILE =
PRINT_RANDOM_INPUTORDERFILE =
RANDOM_INPUTORDERFILE =
----------------------- Advanced options ------------------------
OVERRIDE_WARNINGS = 1
ORDER_BY_RUN = 1

 

説明していきます(ガバガバです)。

# --------------- Main parameters ---------------------------------------------
DATATYPE 0

# インプットファイルの種類

#indfile なら0、popfile なら1

INDFILE Hoge.indfile   
POPFILE         

#インプットファイルの名前

#インプットファイルの種類によって変える、使っていない方は空欄

#今回はindfile を使った
OUTFILE Hoge.outfile     

#アウトプットファイルの名前
MISCFILE Hoge.miscfile

#サマリーファイルの名前
K 5      #クラスター数
C 67       #個体もしくは集団数
R 10       #run もしくはQ-matrices の数 #inputfileの中を見ればわかる
M 1

# Method to be used (1 = FullSearch, # 2 = Greedy, 3 = LargeKGreedy).

#どのアルゴリズムを使って計算するか

#基本的にはFullSearch を使えってAn overview of STRUCTURE に書いてあった

#manual のセクション2に詳しく書かれている

#読んでもわかりませんでした、ごめんなさい
W 0

 

#母集団数に関係なく各母集団に等しく重み付けを行うなら0、行わないなら1

# indfile をインプットにするときは自動的に0 になる
S 1                                     

# Pairwise matrix similarity statistic 

#使用されるペアワイズ行列類似度統計量 G = 1か2

#マニュアルを読んでもよくわかりませんでした。                                       

#1 = G, 2 = G'.

----------------------- Optional outputs ------------------------

PRINT_PERMUTED_DATA 1 # Print the permuted data (clusters) in INDFILE or POPFILE to  PERMUTED_DATAFILE 0 = don't print, 1 = print into one file, 2 = print into separate files for each run).

#INDFILEまたはPOPFILE内のパーミュレートされたデータ(クラスタ)をPERMUTED_DATAFILEにプリントするかしないか

#0 = 印刷しない

#1 = 1つのファイルに印刷

#2 = 各実行ごとに別々のファイルに印刷
PERMUTED_DATAFILE permuted_data.txt

# The permuted data (clusters) will be printed to this file (if PRINT_PERMUTED_DATA = 2, several files with the extensions "_1" to "_R" will be created).

# このファイルには、パーミューティングされたデータ(クラスタ)が印刷される#PRINT_PERMUTED_DATA = 2の場合、拡張子が"_1 "から"_R "のファイルがいくつか作成される
PRINT_EVERY_PERM 0

# Print every tested permutation of the runs and the corresponding value of H to a file specified by EVERY_PERMFILE

#0 = don't print, 1 = print

#書き込みすると非常に大きなファイルになる場合があるので注意

EVERY_PERMFILE
PRINT_RANDOM_INPUTORDER 0 #必須オプションらしい

PRINT_RANDOM_INPUTORDERFILE 
RANDOM_INPUTORDERFILE               

----------------------- Advanced options ------------------------
OVERRIDE_WARNINGS 1

#このオプションは、ユーザーがプログラムからの重大でない警告を上書きすることを可能にする

#0 は警告を許可し、1 は重大でない警告を発行しない

ORDER_BY_RUN 1

#INDFILEやPOPFILEでのラン数を指定するときは1、指定なしは0

# Permute the clusters of the output files by the specified run.

#これもよくわからなかったのでマニュアルの英文をそのまま載せます

#Permute the clusters according to the cluster order of a specificrun. Set this parameter to a numberrfrom 1 toRto order the clusters by runrfrom theINDFILE or POPFILE; set to 0 to not specify a run. When M=1, setting ORDERBYRUNto 0 will result in the clusters being ordered by run 1, and when M=2 or 3, the clusters willbe ordered by the first run in the first input sequence tested. Setting ORDERBYRUN to anonzero value allows the user to determine how the clusters of the Q-matrices will be orderedafter they have been aligned usingCLUMPP. By reordering the aligned clusters of the runsby the input from one specific run, this option is useful for assessing the consistency of theresults of theGreedyorLargeKGreedyalgorithms across different input sequences of runs.

 

ここまででようやくインプットファイルができます。疲れた。

CLUMPP のソフトウェアファイルが入っているフォルダ(解凍したフォルダの中にあったはず、記憶が不明瞭ですみません)にこのparamfile を移動します。

ここで右クリックでLinux シェルを立ち上げて

./CLUMPP Hoge.paramfile と打ちます。

この後の処理がかなり重くて時間がかかるのでnohup で実行した方がいいです。

今回のわたしのデータでは三日動かず泣く泣く閉じて解析をかけ直しました。

 

で、この記事をここまで書き終えたのが2021年1月13日だったわけですが、1月27日現在なんと、CLUMPPの計算が終了していません。

 

たぶんですが、結局STRUCTURE で計算が収束しなかったものというのはCLUMPP でも収束しないんだな、と思います。

そろそろ修論発表会があるのでいい加減にSTRUCTURE のきれいな図を出したい…。

という訳で、一先ずCLUMPPからdistructという流れは諦めます。 別の方法できれいなバープロットを出してやんよ(別記事に続く…)。

 

【追記】続編です。

nemunemu-nyanko.hatenablog.com