ゲノムシンテニーの可視化

# はじめに

ゲノムアセンブリやパンゲノム系の論文でよく見かけるシンテニーを可視化した図。これを書きたいと思って調べたときの備忘録です。今回は可視化ツールとして**SynVisio**と**GENESPACE**を取り上げました。

# 解析に用いたデータ
練習として、３品種のソルガム（*Sorghum bicolor* (L.) Moench）の遺伝子・タンパク質のアノテーションデータで解析してみた。Phytozome (https://phytozome-next.jgi.doe.gov/) からダウンロードできる以下のデータを用いた。
|Cultivar  |Gene annotation |Protein annotation  |
|---|---|---|
|BTx623  |Sbicolor_454_v3.1.1.gene.gff3.gz  |Sbicolor_454_v3.1.1.protein_primaryTranscriptOnly.fa.gz　 |
|BTx642  |SbicolorBTx642_564_v1.1.gene.gff3.gz  |SbicolorBTx642_564_v1.1.protein_primaryTranscriptOnly.fa.gz  |
|Rio|SbicolorRio_468_v2.1.gene.gff3.gz|SbicolorRio_468_v2.1.protein_primaryTranscriptOnly.fa.gz|

これらのファイルをダウンロードし、`gzip`で解凍する。Protein annotationについては、解凍後それぞれ`BTx623.fa`、`BTx642.fa`、`Rio.fa`として ファイル名を簡略化しておいた。

# 実行環境
* Ubuntu 20.04 LTS
* Python 3.8.10
* R version 4.2.2
* RStudio Desktop 2022.07.2+576

# 【１】 MCScanX と SynVisio によるゲノムシンテニーの可視化
## 0. 必要なツールのインストール　 
### MCScanX
https://github.com/wyp1125/MCScanX#installation

### bedtools
https://bedtools.readthedocs.io/en/latest/content/installation.html

### BLAST+ 
以下のサイトがわかりやすい。
https://heavywatal.github.io/bio/blast.html
## 1. Fastaのマージ　
3品種すべての組み合わせでBLASTを行うので、BLAST用のデータベースには3品種をマージしたfastaを使用する。
```bash
cat BTx623.fa BTx642.fa Rio.fa > 623_642_rio.fa
```
## 2. BLAST
```bash
# データベースの作成
makeblastdb -in 623_642_rio.fa -out 623_642_rio -dbtype prot -parse_seqids
#1品種ずつBLAST
blastall -i 623.fa -d 623_642_rio -p blastp -e 1e-10 -b 5 -v 5 -m 8 -o 623.blast
blastall -i 642.fa -d 623_642_rio -p blastp -e 1e-10 -b 5 -v 5 -m 8 -o 642.blast
blastall -i rio.fa -d 623_642_rio -p blastp -e 1e-10 -b 5 -v 5 -m 8 -o rio.blast
#BLAST結果をマージ
cat 623.blast 642.blast rio.blast > 623_642_rio.blast
```
## 3. gffの作成
MCScanXは以下のような独特な形式のgffを要求する。左から、染色体、遺伝子ID、開始位置、終了位置。
```
os5	Os05t0507200	25150045	25150392
os2	Os02t0774100	33540257	33541838
os5	Os05t0480600	23720056	23723583
os6	Os06t0299300	11200887	11202509
os7	Os07t0587100	24548577	24551667
os5	Os05t0526700	26286838	26287599
```
この形式を目標として、はじめにダウンロードしたアノテーションファイルをもとに整形する。Phytozomeからダウンロードできるgff3の中身は、このようになっている。
```Sbicolor_454_v3.1.1.gene.gff3
##gff-version 3
##annot-version v3.1.1
##species Sorghum bicolor
Chr01	phytozomev12	gene	1951	2616	.	+	.	ID=Sobic.001G000100.v3.2;Name=Sobic.001G000100;ancestorIdentifier=Sobic.001G000100.v2.1
Chr01	phytozomev12	mRNA	1951	2616	.	+	.	ID=Sobic.001G000100.1.v3.2;Name=Sobic.001G000100.1;pacid=37943753;longest=1;ancestorIdentifier=Sobic.001G000100.1.v2.1;Parent=Sobic.001G000100.v3.2
Chr01	phytozomev12	CDS	1951	2454	.	+	0	ID=Sobic.001G000100.1.v3.2.CDS.1;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01	phytozomev12	CDS	2473	2616	.	+	0	ID=Sobic.001G000100.1.v3.2.CDS.2;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01	phytozomev12	gene	11180	14899	.	-	.	ID=Sobic.001G000200.v3.2;Name=Sobic.001G000200;ancestorIdentifier=Sobic.001G000200.v2.1
Chr01	phytozomev12	mRNA	11180	14899	.	-	.	ID=Sobic.001G000200.1.v3.2;Name=Sobic.001G000200.1;pacid=37944019;longest=1;ancestorIdentifier=Sobic.001G000200.1.v2.1;Parent=Sobic.001G000200.v3.2
```
これをbedtoolsでbedに変換。
```bash
sortBed -i Sbicolor_454_v3.1.1.gene.gff3 | gff2bed > ６２３.bed
```
中身を確認すると、
```623.bed
Chr01	1950	2454	Sobic.001G000100.1.v3.2.CDS.1	.	+	phytozomev12	CDS	0	ID=Sobic.001G000100.1.v3.2.CDS.1;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01	1950	2616	Sobic.001G000100.1.v3.2	.	+	phytozomev12	mRNA	.	ID=Sobic.001G000100.1.v3.2;Name=Sobic.001G000100.1;pacid=37943753;longest=1;ancestorIdentifier=Sobic.001G000100.1.v2.1;Parent=Sobic.001G000100.v3.2
Chr01	1950	2616	Sobic.001G000100.v3.2	.	+	phytozomev12	gene	.	ID=Sobic.001G000100.v3.2;Name=Sobic.001G000100;ancestorIdentifier=Sobic.001G000100.v2.1
Chr01	2472	2616	Sobic.001G000100.1.v3.2.CDS.2	.	+	phytozomev12	CDS	0	ID=Sobic.001G000100.1.v3.2.CDS.2;Parent=Sobic.001G000100.1.v3.2;pacid=37943753
Chr01	11179	11510	Sobic.001G000200.1.v3.2.three_prime_UTR.1	.	-	phytozomev12	three_prime_UTR	.	ID=Sobic.001G000200.1.v3.2.three_prime_UTR.1;Parent=Sobic.001G000200.1.v3.2;pacid=37944019
Chr01	11179	14899	Sobic.001G000200.1.v3.2	.	-	phytozomev12	mRNA	.	ID=Sobic.001G000200.1.v3.2;Name=Sobic.001G000200.1;pacid=37944019;longest=1;ancestorIdentifier=Sobic.001G000200.1.v2.1;Parent=Sobic.001G000200.v3.2
```
このようになり、目標に少し近づいた。
これを以下のスクリプトでさらに加工する。
```python:bed2gff.py
cultivar = 'BTx623'

with open(f'./{cultivar}.bed', mode = 'r') as inp:
    with open(f'./{cultivar}.gff', mode = 'w') as out:

        for line in inp:
            line = line.rstrip('\n|\r|\r\n')
            line = line.split('\t')

            if (line[7] == 'mRNA') & line[0].startswith('Chr'):
                Chr = line[0].replace('Chr',f'{cultivar}_0').replace('010', '10')
                gene_id = line[3].split('.')
                protein_id = gene_id[0] + '.' + gene_id[1] + '.' + gene_id[2] + '.p'

                out.write(f'{Chr}\t{protein_id}\t{line[1]}\t{line[2]}\n')

            else:
                pass
```
すると以下の形式になる。
```BTx623.gff 
623_01	Sobic.001G000100.1.p	1950	2616
623_01	Sobic.001G000200.1.p	11179	14899
623_01	Sobic.001G000400.6.p	22390	37967
623_01	Sobic.001G000400.5.p	22390	37972
623_01	Sobic.001G000400.3.p	22390	42443
623_01	Sobic.001G000400.4.p	22489	42443
```
このgffを品種分作成し、最後にマージする。
```bash
cat BTx623.gff BTx642.gff Rio.gff > 623_642_rio.gff
```
## 4. MCScanX の実行
必要なのは`623_642_rio.blast`と`623_642_rio.gff`だが、ファイル名は統一しておく必要がある（コマンドが以下のようになっているため）。
```bash
MacscanX ./623_642_rio
```
これでcollinearityが計算される。
## 5. SynVisio による可視化
SynVisio (https://synvisio.github.io/#/) にアクセス。動画でのチュートリアルも充実している。
![スクリーンショット 2022-12-10 23.31.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/a7d20453-bacf-7ee6-c236-de397252fb6a.png)
上の`UPLOAD OWN DATA TO DASHBOARD`に入り、ファイルを選択して`UPLOAD`。
![スクリーンショット 2022-12-10 23.32.37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/3d341d77-a1b5-d2f6-9269-4955562cc12c.png)
次に、`SYNTENY DASHBOARD`で Multi Lebel Analysisにチェック。今回は３つの品種のゲノムを比較するので、`GO`の隣の`＋`を押して３つのレベルを用意する。それぞれのプルダウンから各レベル（品種）の染色体を選択し、`GO`。
![スクリーンショット 2022-12-10 23.35.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/22118534-97de-f86c-f39e-641d45c3936c.png)
このような図が作成できる。
![synvisio-export-PNG-lbgc0n0q.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/9bdcc11d-84cf-e141-5778-12b53a91ff60.png)

---
# 【２】 GENESPACE によるシンテニーの可視化

## 0. 必要なツールのインストール
### OrthoFinder
https://github.com/davidemms/OrthoFinder#installing-orthofinder-on-linux
### MCScanX 
## 1. データの準備
MCScanX も GENESPACE も、gffのフォーマットをいちいち作り変えないといけないのが難点。我慢して、整形する。ここでは、MCScanX のところでgff3から作成したbedをベースにしている。

```Python
cultivar = 'Rio'

with open(f'./{cultivar}.bed', mode = 'r') as inp:
    with open(f'./{cultivar}.gff', mode = 'w') as out:

        for line in inp:
            line = line.rstrip('\n|\r|\r\n')
            line = line.split('\t')

            if (line[7] == 'mRNA') & line[0].startswith('Chr'):
                Chr = int(line[0].replace('Chr',''))
                gene_id = line[3].split('.')
                protein_id = gene_id[0] + '.' + gene_id[1] + '.' + gene_id[2] + '.p'

                out.write(f'{Chr}\tPhytozome\tgene\t{line[1]}\t{line[2]}\t\t{line[5]}\t\tlocus={protein_id}\n')

            else:
                pass
```
これを実行すると、以下のようなファイルが出来上がる。
```Rio_gene.gff
1	Phytozome	gene	23885	25465		+		locus=SbRio.01G000100.1.p
1	Phytozome	gene	33801	34918		-		locus=SbRio.01G000200.1.p
1	Phytozome	gene	39376	40597		+		locus=SbRio.01G000300.1.p
1	Phytozome	gene	48624	52364		-		locus=SbRio.01G000400.1.p
1	Phytozome	gene	59938	61551		-		locus=SbRio.01G000500.1.p
1	Phytozome	gene	61878	75618		-		locus=SbRio.01G000600.1.p
```
Phytozomeからダウンロードしておいたprotein fastaも用意する（BTx623_protein.fa）。
これらのファイルはgz形式で圧縮しておく。
```bash
gzip BTx623_gene.gff
gzip BTx623_protein.fa
```
ファイル構成は、テストデータと同じく以下のようにしておく。
```
── rawGenomes
   ├── BTx623
   │   └── BTx623
   │       └── annotation
   │           ├── BTx623_gene.gff.gz
   │           └── BTx623_protein.fa.gz
   ├── BTx642
   │   └── BTx642
   │       └── annotation
   │           ├── BTx642_gene.gff.gz
   │           └── BTx642_protein.fa.gz
   └── rio
       └── rio
           └── annotation
               ├── rio_gene.gff.gz
               └── rio_protein.fa.gz
```

## 2. GENESPACE のインストールと実行
ここからは、Rstudioで作業する。実行スクリプトの全体は以下の通りです。
```R: plot_synteny_genespace.R
if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("Biostrings", "rtracklayer"))
devtools::install_github("jtlovell/GENESPACE", upgrade = F)

gids <- c("BTx623","rio","BTx642") 
gpar <- init_genespace(
  genomeIDs = gids, 
  speciesIDs = gids, 
  versionIDs = gids, 
  ploidy = rep(1,3),
  wd = runwd, 
  gffString = "gff", 
  pepString = "protein",
  path2orthofinder = "/mnt/data/software/OrthoFinder/orthofinder", 　
  path2mcscanx = "/mnt/data/software/MCScanX",
  rawGenomeDir = file.path(runwd, "rawGenomes"))

parse_annotations(
  gsParam = gpar, 
  gffEntryType = "gene", 
  gffIdColumn = "locus",
  gffStripText = "locus=", 
  headerEntryIndex = 1,
  headerSep = " ", 
  headerStripText = "locus=")

gpar <- run_orthofinder(gsParam = gpar)
gpar <- synteny(gsParam = gpar)

pdf('synteny.pdf', height = 3, width = 7)
plot_riparianHits(gpar, blackBg=F)
dev.off()

# pngで保存の場合
png('synteny.png',width=6,height=2,units="in",res=1200)
plot_riparianHits(gpar, blackBg=T)
dev.off()
```
以下、補足など。
devtools をインストールする。
```R
if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")
```
必要な外部パッケージをインストールされていないというエラーがでたが、その都度 `sudo apt-get install -y ~` でインストールしていけば大丈夫。

依存パッケージ (Biostrings および rtracklayer) をインストールする。
```R
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("Biostrings", "rtracklayer"))
```
github経由でGENESPACEをインストールする。
```R
devtools::install_github("jtlovell/GENESPACE", upgrade = F)
```
GENESPACE を実行する。
```R
gids <- c("BTx623","rio","BTx642") #出力される図に表示したい順番
gpar <- init_genespace(
  genomeIDs = gids, 
  speciesIDs = gids, 
  versionIDs = gids, 
  ploidy = rep(1,3),
  wd = runwd, 
  gffString = "gff", 
  pepString = "protein",
  path2orthofinder = "/mnt/data/software/OrthoFinder/orthofinder", 　
  path2mcscanx = "/mnt/data/software/MCScanX",
  rawGenomeDir = file.path(runwd, "rawGenomes"))
```
OrthoFinder のパスはバイナリまで。MCScanX のパスはフォルダまで。
次に、gffのパーサを動かす。
```R
parse_annotations(
  gsParam = gpar, 
  gffEntryType = "gene", 
  gffIdColumn = "locus",
  gffStripText = "locus=", 
  headerEntryIndex = 1,
  headerSep = " ", 
  headerStripText = "locus=")
```
Orthofinder で解析後、シンテニー解析を行う。
```R
gpar <- run_orthofinder(gsParam = gpar)
gpar <- synteny(gsParam = gpar)
```
最後に可視化し、pdf で保存する。デフォルトだと背景が黒くなるが、白いほうが何かと使いやすい気がしたので `blackBg=F` を指定している。
```R
pdf('synteny.pdf', height = 3, width = 7)
plot_riparianHits(gpar, blackBg=F)
dev.off()

または、
png('synteny.png',width=6,height=2,units="in",res=1200)
plot_riparianHits(gpar, blackBg=T)
dev.off()
```
これで以下のような図が保存される。
![synteny.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/0ff7e4e4-c189-8c42-ef7d-b28cc1621056.png)
ちなみに背景が黒いとこんな感じになる。マージンは変更できなさそう。
![synteny_bk.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3044655/678750f1-34d8-35ea-69b6-e3edd07d3bf7.png)