はじめに
ゲノムアセンブリやパンゲノム系の論文でよく見かけるシンテニーを可視化した図。これを書きたいと思って調べたときの備忘録です。今回は可視化ツールとして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
bedtools
BLAST+
以下のサイトがわかりやすい。
https://heavywatal.github.io/bio/blast.html
1. Fastaのマージ
3品種すべての組み合わせでBLASTを行うので、BLAST用のデータベースには3品種をマージしたfastaを使用する。
cat BTx623.fa BTx642.fa Rio.fa > 623_642_rio.fa
2. BLAST
# データベースの作成
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に変換。
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
このようになり、目標に少し近づいた。
これを以下のスクリプトでさらに加工する。
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を品種分作成し、最後にマージする。
cat BTx623.gff BTx642.gff Rio.gff > 623_642_rio.gff
4. MCScanX の実行
必要なのは623_642_rio.blastと623_642_rio.gffだが、ファイル名は統一しておく必要がある（コマンドが以下のようになっているため）。
MacscanX ./623_642_rio
これでcollinearityが計算される。
5. SynVisio による可視化
SynVisio (https://synvisio.github.io/#/) にアクセス。動画でのチュートリアルも充実している。
上のUPLOAD OWN DATA TO DASHBOARDに入り、ファイルを選択してUPLOAD。
次に、SYNTENY DASHBOARDで Multi Lebel Analysisにチェック。今回は３つの品種のゲノムを比較するので、GOの隣の＋を押して３つのレベルを用意する。それぞれのプルダウンから各レベル（品種）の染色体を選択し、GO。
このような図が作成できる。
【２】 GENESPACE によるシンテニーの可視化
0. 必要なツールのインストール
OrthoFinder
MCScanX
1. データの準備
MCScanX も GENESPACE も、gffのフォーマットをいちいち作り変えないといけないのが難点。我慢して、整形する。ここでは、MCScanX のところでgff3から作成したbedをベースにしている。
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形式で圧縮しておく。
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で作業する。実行スクリプトの全体は以下の通りです。
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 をインストールする。
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
必要な外部パッケージをインストールされていないというエラーがでたが、その都度 sudo apt-get install -y ~ でインストールしていけば大丈夫。
依存パッケージ (Biostrings および rtracklayer) をインストールする。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("Biostrings", "rtracklayer"))
github経由でGENESPACEをインストールする。
devtools::install_github("jtlovell/GENESPACE", upgrade = F)
GENESPACE を実行する。
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のパーサを動かす。
parse_annotations(
gsParam = gpar,
gffEntryType = "gene",
gffIdColumn = "locus",
gffStripText = "locus=",
headerEntryIndex = 1,
headerSep = " ",
headerStripText = "locus=")
Orthofinder で解析後、シンテニー解析を行う。
gpar <- run_orthofinder(gsParam = gpar)
gpar <- synteny(gsParam = gpar)
最後に可視化し、pdf で保存する。デフォルトだと背景が黒くなるが、白いほうが何かと使いやすい気がしたので blackBg=F を指定している。
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()
これで以下のような図が保存される。
ちなみに背景が黒いとこんな感じになる。マージンは変更できなさそう。
