PhyloNetを使って系統推定/visualize(RAD-seq)

**PhyloNet**というのは系統樹の推定・作成ソフトです。このソフトはIQ-TREEなどで推定された複数の遺伝子ごとの系統樹をもとに種間のネットワークを推定することができます。が、これを使おうとして公式サイトのTutorialを見ると、インストールの次にあるのがPhyloNetの起動コマンドとinputファイルの中身が貼られており、inputファイルの構造について説明があるのみで面食らいます。食らいました。


複数の系統樹を考慮して系統ネットワークを作成する、というところがミソ、というかそれさえ理解しておけばどうにかなると思いますが、せっかくなのでここではRAD-seqのシーケンスデータからipyrad, IQ-TREEを用いてPhyloNetのinputファイル(NEXUS形式)を作成するまでの流れを備忘録的に書いておこうと思います。
説明はipyradで処理したデータがあることを前提として進めますが、stacksなどで処理したデータでも可能と思われます。

なお、ここで記載しているコードの多くはChatGPTに投げたり修正したりで構成していますので、具体的によくわかっていないところも存在しますというか、うp主はプログラミング初心者なので大目に見てください。
##### 間違いあったら教えて下さい！

### 環境

```
macOS 14.2.1
ipyrad  0.9.95
IQ-TREE2 2.3.4
PhyloNet 3.8.2
```

## インストール
[公式サイト](https://phylogenomics.rice.edu/html/phylonet.html)にいって書いてある通りやるだけです。
これを書いているときの最新バージョンは3.8.2。
とりあえず最新のものをDLしておけばよいと思います。`.Jar`ファイルをてきとうなディレクトリに配置すればおっけー。

## Inputファイル作成の流れ
ipyradを用いてアライメントを行った場合、step7ですべてのローカスの配列が1つにまとめられたlociファイルが出力されます。これを1ローカスごとのFASTAファイルに分割します。
この分割したFASTAファイルごとにIQ-TREEで系統の推定を行い、その結果を`.treefile`として出力します。この`.treefile`群に記述されている系統樹を一つのNEXUSファイルにまとめ、Tutorialにあるようなinputファイルを作ります。
ではやっていきましょう。

:::note 
やり口は上記の方法以外にも色々ありそうなような、そんなことはないような…
:::

### ipyrad出力のlociファイルから切り出し
ipyradで処理されたデータがあることを前提として進めます。ipyradは1~7のステップがありますが、このうちのstep7で出力される`.loci`ファイルを使用します。`.loci`ファイルには各ローカスの配列情報が入っているため、この配列情報をローカスごとに切り出し、FASTAファイルとして保存します。

以下のコードを`任意の名前.py`で保存してターミナル上にて`python3 任意の名前.py`などと打って実行します。すると、`locus_fasta`という名前のディレクトリができ、その中にlocusごとに分割されたFASTAファイルが生成されます。
```ruby:ターミナル
from pathlib import Path
input_file = "input_data.loci"  # 実際のパスに合わせる
output_dir = Path("locus_fasta")
output_dir.mkdir(exist_ok=True)
with open(input_file) as f:
    locus = []
    i = 1
    for line in f:
        if line.startswith('//'):
            if locus:
                with open(output_dir / f"locus_{i}.fasta", "w") as out:
                    for entry in locus:
                        name, seq = entry
                        out.write(f">{name}\n{seq}\n")
                i += 1
                locus = []
        elif line.strip():
            parts = line.strip().split()
            if len(parts) == 2:
                locus.append((parts[0], parts[1]))

```

### IQ-TREEで各FASTAファイルの系統推定

`cd`で`locus_fasta`に移動し、IQ-TREEで系統の推定をします。
```ruby:ターミナル
cd <パス名/locus_fasta>

for file in locus_*.fasta; do
    iqtree2 -s "$file" -m MFP -nt AUTO -pre "${file%.fasta}"
done

#IQ-TREEの細かい設定は調べてください。すまん。
#-m MFP はmodel finder。いい感じのモデルを勝手に選んでくれってやつです。
```
すると、`.bjonj` `ckp.gz` `.iqtree` `.log` `.mldist` `model.gz` `.treefile` `.uniqueseq.phy`などのファイルがoutputされますが、PhyloNetのinputファイルであるNEXUS形式のファイルを作成するにあたって使用するのは`.treefile`のみです。この時点で`.treefile`には樹長のデータが入っており、これがあるとPhyloNetは文句を言って解析をしてくれないので削除する必要があります。

以下のコードを`任意の名前.py`で保存してターミナル上で`python3 任意の名前.py`的な感じに実行すると樹長が削除された`.treefile`が入った`tree_nobranch`といったディレクトリが作業ディレクトリ内に作製されます。
```ruby:ターミナル
import re
from pathlib import Path
input_dir = Path("<絶対パスで指定した方がよさげ/locus_fasta>")  # treefile のあるフォルダに合わせて変更
output_dir = Path("tree_nobranch")
output_dir.mkdir(exist_ok=True)
for treefile in input_dir.glob("*.treefile"):
    with open(treefile) as f:
        tree = f.read().strip()
    # 枝長情報 :0.123 のような部分を削除
    tree_no_branch = re.sub(r":-?\d+\.?\d*(e[-+]?\d+)?", "", tree)
    # 出力ファイル名
    out_file = output_dir / treefile.name
    with open(out_file, "w") as f:
        f.write(tree_no_branch + "\n")
```


### PhyoNetのinputファイル作成
Inputファイルの作製はIQ-TREEから出力された`.treefile`内の系統樹をPhyloNet公式サイトにあるようなNEXUSファイルの形にまとめることとなります。先ほど樹長を削除した`.treefile`の中身をまとめるわけです。

以下のコードを`任意の名前.py`で保存し、`python3 任意の名前.py`で実行するという流れを行うと、`phylonet_input.nexus`というinputファイルができます。
```ruby:ターミナル
cd <パス名/tree_nobranch>

from pathlib import Path
# .treefile があるフォルダ
tree_dir = Path("パス名/tree_nobranch")  
# 出力ファイル
output_file = Path("phylonet_input.nexus")
# すべての .treefile を取得
treefiles = sorted(tree_dir.glob("*.treefile"))
with open(output_file, "w") as out:
    out.write("#NEXUS\n")
    out.write("BEGIN TREES;\n")
    for i, tf in enumerate(treefiles, 1):
        with open(tf) as f:
            newick = f.read().strip()
            # 空でなければ書き出し（空行やエラー回避）
            if newick:
                out.write(f"  Tree tree_{i} = {newick}\n")
    out.write("END;\n\n")
    out.write("BEGIN PHYLONET;\n")
    out.write("  InferNetwork_ML (tree_1")
    for i in range(2, len(treefiles)+1):
        out.write(f", tree_{i}")
    out.write(") 1 -pl 6 ;\n")
    out.write("END;\n")
print(f": {len(treefiles)} 本の系統樹を含む NEXUS ファイルを生成しました: {output_file}")
```
細かい中身は違えど、[TutorialのInferNetwork_ML](https://phylogenomics.rice.edu/html/commands/InferNetwork_ML.html)にもあるような中身のファイルになるはずです。
↓こんな感じのやつ

```
#NEXUS

BEGIN TREES;

Tree geneTree1 = ((C,((B,D),A)),E);
Tree geneTree2 = (B,(D,(C,(A,E))));
Tree geneTree3 = (D,(B,((C,E),A)));
Tree geneTree4 = (D,((B,E),(C,A)));

END;


BEGIN PHYLONET;

InferNetwork_ML (geneTree1,geneTree2,geneTree3,geneTree4) 1 -bl 6;

END;
```
`BEGIN PHYLONET;`ブロックでコマンドを打ちます。
手法を指定（ここでは`InferNetwork_ML`）して、推定に使用する系統樹のリストを書きます。
` 1 -bl;`この1はreticulation数（交雑・遺伝子流動イベントの回数）です。想定されるreticulation数を書いて解析すればよいと思います。
PhyloNetはなにも指定しなかった場合コア数1で計算をするので、`-bl 6`など指定してコア数を指定するとよいです。

:::note warn
`InferNetwork_ML (geneTree1,geneTree2,geneTree3,geneTree4) 1 -bl;`ここの`(geneTree1,geneTree2,geneTree3,geneTree4)`部分にスペースを入れるとエラーを吐くので注意。
:::

### PhyloNetの実行


めでたくinputファイルができたところでPhyloNetの実行です。
特にむずかしいことはなく、`PhylonNet.jarのパス` `inputファイルのパス` `outputファイルのパス`といった順序で並べるだけといった雰囲気。
以下が実行のコードです。
```ruby:ターミナル
java -jar パス名/PhyloNetvx_y_z.jar パス名/phylonet_input.nexus > パス名/phylonet_output.txt
```
みたいな感じです。
ここは[Tutorial](https://phylogenomics.rice.edu/html/tutorials.html)を見てもわかると思います。実行すると指定したパスにoutputファイルができるはず。

### visualize(図示)
無事にoutputファイルができたら、もう一仕事必要です。
outputファイルは以下のような内容と思われますが、このうちの`inferred Network:`のブロックにある系統樹を切り取って1行のみの`.txt`ファイルを作ります。
```
InferNetwork_ML_(tree_1, tree_2, tree_3) 1 -pl 6
Results after run #1
-123.456:((((A,B),C),D),E);
Running Time (min): 1.2345
============================

Results after run #2
-123.456:((((A,B),C),D),E);
Running Time (min): 1.2345
============================

Inferred Network #1:
((A:0.2,((B:0.1)#H1:0.3):0.3,(#H1:0.2,C:0.3):0.2):0.4,D:0.5,E:0.6); 
Total log probability: -12.3456789Running time: 123.456789
```
この中の

```
((A:0.2,((B:0.1)#H1:0.3):0.3,(#H1:0.2,C:0.3):0.2):0.4,D:0.5,E:0.6); 
```
ここだけコピーして新たな`.txt`ファイルを作ります。
#H1の部分が交雑(hybrid)の表記のようです。

図示にあたっては**Dendroscope**というソフトを使う方法と **icytree** ブラウザ上で動かせるツールを使う方法があります。Dendroscopeを使う場合は[公式サイト](https://uni-tuebingen.de/fakultaeten/mathematisch-naturwissenschaftliche-fakultaet/fachbereiche/informatik/lehrstuehle/algorithms-in-bioinformatics/software/dendroscope/)に飛んでインストールしてください。

:::note
Dendroscopeはインストールが必要ですが、icytreeはオンライン上で使えるらしいです。
:::

なお、ここで系統樹内の`#H`の後に`<NUMBER>::<NUMBER>`といったような2重コロンが確認されるのですが、どうもこれがあるとDendroscopeが読み込んでくれません。公式Tutorialの **”Visualizing a Phylogenetic Network”** を見ると、`::<NUMBER>`の部分をしゅどうで消去し、1つ目の樹長のみを残すといいよ、というようなことが書いてあるのでそうしましょう。

:::note
この数値は "inheritance probabilities(遺伝確率？)"らしいです。
手動で消去するか、オプション`-di`を使用してくれ、といったことが書かれています。
icytreeの場合は認識してくれたり認識してくれなかったりするようです。
なお、樹長の方は双方問題ないとのこと。
:::


あとは系統樹のみ抜き出した`.txt`ファイルをDendroscopeで開くとFigが出ます。
細かい表示の設定はまあ、てきとうに……。

#### 2025.06.02 追記
PhyloNetの解析はざっくり触っている感じかなり遅いです。
ところで、公式のサイトを見るとPhyloNetPyなるものがあります。どうもこれはPhyloNetと比べてかなり速いようなのですが、現状はサイト内にあるTutorialに存在する解析しかできないらしい。つまり、ML法はできなさそう＆導入がかなりやっかいそうなので私自身は触れていません。
ちゃんと調べてみたら使えるよ！とｋいい方法など知っている方がいましたら教えてください……。


## 引用文献・サイト
・PhyloNetのTutorial：https://phylogenomics.rice.edu/html/phylonetOverview.html
・Dendroscopeのマニュアル（リンク先PDF）：https://software-ab.cs.uni-tuebingen.de/download/dendroscope/manual.pdf