Analysis: MCScanX

このページはゲノム間のシンテニー解析に有用なソフトウェアMCScanXの使い方（細菌ゲノム特化）の備忘録です．

必要PCスペック: Blastが快適にできる程度のLaptopで十分可能．

＜注意＞下記の内容は間違いを含む可能性がありますので，参考程度に自己責任でお願いします．

目次

MCScanXとは？

進化や種分化に伴う遺伝子の獲得/喪失に関する洞察を得るために同種または近縁種間で染色体ゲノムにおける遺伝子の順序（シンテニー）を比較することは有用である．

2012年に公表されたソフトウェアMCScanX（Wang et al., 2012; https://github.com/wyp1125/MCScanX）では、複数のゲノム（または染色体）間で相同性のある遺伝子をBLASTP（all-against-all）で検出した結果をもとに「collinear block」としてシンテニーを整列する．オンラインツールを用いて，gffファイルおよびMCScanXが出力するcollinearityファイルをアップロードすることでシンテニーの可視化ができる．

Wang et al. (2012)のこれまでの被引用数は5270件（2025年2月15日時点），また，発表から12年経過した2024年に解説論文（Wang et al., 2024）が出版されるなど，かなり息の長いソフトウェアといえる．インストール，実行，および可視化の流れを解説した日本語サイト（下記）は複数存在する．しかし，入力データファイルの用意についてこと細かく解説しているサイトはあまりない．そこで，データ整形に重点を置きながら使い方を解説してみようと思う．

Wang et al. (2012)  https://doi.org/10.1093/nar/gkr1293

MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity, Nucleic Acids Research, 40, e49,

Wang et al. (2024)   https://doi.org/10.1038/s41596-024-00968-2

Detection of colinear blocks and synteny and evolutionary analyses based on utilization of MCScanX. Nature Protocols, 19, 2206–2229.

可視化に使用するオンラインツール

SynVisio https://synvisio.github.io/#/

AccuSyn（環状の図のみ） https://accusyn.usask.ca/

本ページをまとめる上で参考にさせていただいたウェブサイト

https://kazumaxneo.hatenablog.com/entry/2020/06/20/144410

https://kazumaxneo.hatenablog.com/entry/2023/02/15/100040

https://qiita.com/Shumpei_Hashimoto/items/ea1f00b32d040d11cd7a

0．そもそもシンテニー解析をする必要がある（細菌）ゲノムなのか？

事例１: 単一ゲノムのCircos plotを作りたい場合

最近はあまり見なくなってきた？ような気がするCircos plotで細菌の環状ゲノムを図示するような場合にはMCScanXは不適である．

例えば，GenoVi（https://github.com/robotoD/GenoVi）などを使えば，そのまま論文に使えるクオリティの図が出来上がる．

MCScanXの結果を利用して同等の図を作るためには，GC含量，GC skew，CDSの位置・向き等のTrackファイルを別途用意しておく必要がある（番外編を参照）．

GenoViの使い方の参照　https://kazumaxneo.hatenablog.com/entry/2023/04/10/020916

事例２: なんとなくアタリを付けたいだけのとき

MCScanXを使うには，後述の通り，そこまで大変ではないが，入力ファイルの整形が必要である．そのため，結果がいまいちな場合は徒労となる．

もし不安な場合は，複数ゲノムのシンテニーおよびDot plotを出力してくれるDigAlign（https://www.genome.jp/digalign/）を使ってアタリをつけると良いかもしれない．

参考: TogoTVによる解説動画 「DiGAlign を使ってゲノムの変化を可視化する」　https://togotv.dbcls.jp/en/20241022.html

１．入力データファイルの用意

必要なファイルは，各ゲノムのタンパク質配列（.faa）および簡略化されたgff3ファイル（.gff）

gff3のフォーマット（特に遺伝子IDは.faaのヘッダーと必ず一致させなければいけない）　下の完成形gff3の画像も参照

1列目 2列目 3列目 4列目

ゲノム配列名 遺伝子ID 開始位置 終了位置

※細菌ゲノムで利用する場合の事例を紹介するが，ほかの生物でも流れは同じはず．

※３つ以上の細菌ゲノムも比較可能だが，あまりに多くしすぎると見えづらい図になってしまう

（ただし，複数のゲノムを含む解析をした場合でも可視化の時点で使用するゲノムは選択可能．その際に数多く選びすぎなければよい）．

A: NCBIに登録されているゲノムのアノテーション結果を使用する場合

A１．比較したいゲノムをNCBIのGenomeデータベースで検索

細菌ゲノムの場合，基本的にはコンプリートゲノムでないとシンテニーの比較は難しい．コンプリートでない場合は，最長のコンティグに絞るなど工夫はできると思う．

A２．タンパク質配列の取得

Downloadのアイコンから.faaをダウンロード，解凍する．

このときgff3ファイルも取得可能だが，MCScanXの入力ファイルフォーマットにするための整形が大変なので，別の方法でMCScanX用のgff3ファイルを作る．

※下記のサイトを参考にDownloadから取得したgff3ファイルでもbedtools等を使って整形することも可能．

https://qiita.com/Shumpei_Hashimoto/items/ea1f00b32d040d11cd7a

A3．タンパク質配列のヘッダー名の整理

.faaを正規表現（grep）で置換可能なテキストエディタで開く．※ファイルをいちいち開かなくてもコマンドでできるはず．．．

grep [.1 .+] → [.1] （.1スペース以降の文字列を.1に置換する）．遺伝子IDのみが残っているはずなので，その状態で上書き保存．

＜置換前＞　赤字が不要

>WP_005504750.1 MULTISPECIES: 30S ribosomal protein bS22 [Bacteria]

MGSVIKKRRKRMSKKKHRKLLRKTRHQRRNKK

>WP_029067950.1 MULTISPECIES: type B 50S ribosomal protein L36 [Micrococcales]

MKVRNSLKSLKNQPGSQVVRRRGRVFVINKKNPRLKARQG

＜置換後＞

>WP_005504750.1

MGSVIKKRRKRMSKKKHRKLLRKTRHQRRNKK

>WP_029067950.1

MKVRNSLKSLKNQPGSQVVRRRGRVFVINKKNPRLKARQG

＜追記１＞ヘッダー名の整理

友人のツボカビ博士がヘッダー整理はseqkitでできるよ！と教えてくれました．

seqkit seq -i replace_test.faa > replace_test_seqkit.faa

-i, --only-id               print IDs instead of full headers（1つ目のspace以降を削除）

Chatgptに聞いてみたら教えてくれた方法

 sed -E 's/(>WP_[0-9]+\.[0-9]+) .*/\1/' replace_test.faa > replace_test_output.faa

コマンドの解説：

sed -E：拡張正規表現を有効化

s/(>WP_[0-9]+\.[0-9]+) .*/\1/：>WP_[0-9]+\.[0-9]+ → WP_ で始まり、数字と . を含む識別子をキャプチャ

.* → その後のスペース以降の文字列を削除

\1 → キャプチャした識別子のみを残す

A4．アノテーションデータの取得

View annotated genesをクリックし，すべての項目にチェックを入れてからDownloadをクリックし，Download Tableで.tsvファイルを取得する．

ちなみにRefSeqの場合（GCF_），A２でダウンロードした.faaのヘッダーは，Proteins列にある名前（WP_）と一致する．

A5．アノテーションデータの整形

A4で取得した.tsvファイルをExcelで開く．必要なのは赤四角の列のみ

※こちらの整形もちゃんとスクリプトを作れば一括処理できるはず．．．

＜Excel上の手順＞

・Excelで，データ→フィルターでProtein accession（K列）をソート

・protein-codingではないLocus（tRNA, rRNA, pseudogene）が含まれているので，それらの行を削除する

・Locus tagで再度ソート（たぶん必要ないけど揃えたほうが気分が良い）

・Accession（A列）を適宜置換する

今のままではNZ_CP080203.1が最終的な名前として使用される．染色体とプラスミドが両方あるゲノムの場合は区別可能な名前に変えておかないと可視化した際にアクセッション番号の文字列をみて使用するゲノムを選ぶことになるので大変面倒である．自分で識別可能な簡単な名前に変える．例: NZ_CP080203.1 -> Mspo_Chr1またはChr1など

※SynVisioでは名前は割と自由に設定できる一方，AccuSynでは短い名前で最後は数字が付いている必要があるなど制約が多い．

AccuSymを使いたい場合はあらかじめChr1, Chr2等シンプルな名前にしておいたほうが無難．

・K列をA列とB列の間に挿入

・ヘッダー行を削除

・セルの全選択→メモ帳に貼り付け→拡張子.gffで保存

B: 自分のゲノムを使用する場合（DFASTでアノテーションした場合）

B1: DFAST（https://dfast.ddbj.nig.ac.jp/）のアノテーション結果を取得

あらかじめSeq-names/Gene-IDsのところからSequence NameおよびLocus tag prefix（BioSample登録可能な文字列にしておくと尚良し）を分かりやすくしてSaveしておく．

結果のページからprotein.faaおよびannotation.gffを取得する（annotation.zipをダウンロードでもよい）．

B2．タンパク質配列のヘッダー名の整理（DFAST）

.faaを正規表現（grep）で置換可能なテキストエディタで開く．※ファイルをいちいち開かなくてもコマンドでできるはず．．．

grep [ .+] → [入力なし] （スペース以降の文字列を削除）．Locus名のみが残っているはずなので，その状態で上書き保存．

＜置換前＞　赤字が不要

>MLACM10_00010 hypothetical protein

MDKPVNNFWGDMSNQLDETQQSWQLVLDAMSQDSRVTPPLHGFLNLVIPKGVMAG

>MLACM10_00020 DNA polymerase III subunit beta

MRFSVNRDVFSEAVSFAVKLLPQRTTLPILSGVLIETNDTGLTLSSFDYEVSAQTRIVAEIEE

＜置換後＞

>MLACM10_00010

MDKPVNNFWGDMSNQLDETQQSWQLVLDAMSQDSRVTPPLHGFLNLVIPKGVMAG

>MLACM10_00020

MRFSVNRDVFSEAVSFAVKLLPQRTTLPILSGVLIETNDTGLTLSSFDYEVSAQTRIVAEIEE

B3．アノテーションデータの整形

B1で取得した.tsvファイルをExcelで開く．必要なのは赤四角の列のみ

※こちらの整形もちゃんとスクリプトを作れば一括処理できるはず．．．

＜Excel上の手順＞

・I列を選択，データ→区切り位置セミコロン指定

・I列に「ID=MLACM10_00010」のように残るはず．「ID=」の部分を置換で削除

・Excelで，データ→フィルターでH列をソート

・protein-codingではないLocus（tRNA, rRNA, pseudogene）が含まれているので，それらの行を削除する（０と．のうち，．のほうの行）

その際，末尾の##FASTA以下にゲノム配列が載っているのでそれらも行ごとすべて削除

・Locus tag（I列）で再度ソート（たぶん必要ないけど揃えたほうが気分が良い）

・Accession（A列）を適宜置換する（理由は上記）

・BおよびC列を削除，A列および新しいB列の間にLocus tagの列（元のI列）を挿入

・ヘッダー行を削除，A-D以外の不要な列を削除

・セルの全選択→メモ帳に貼り付け→拡張子.gffで保存

完成形のgff3が完成！

＜追記2＞gffreadを利用した入力ファイルの作成

今回もツボカビ博士からの情報提供（多謝）．

gffreadの使い方　https://kazumaxneo.hatenablog.com/entry/2019/08/26/073000

https://ccb.jhu.edu/software/stringtie/gff.shtml#gffread

1.  gffreadを使ってgff3ファイルとgenome.fastaファイルからprotain.fastaを作る

gffread XXX.gff3 -g genome.fasta -y gffread_protein.faa

-x    write a fasta file with spliced CDS for each GFF transcript

 -y    write a protein fasta file with the translation of CDS for each record

2. gffreadを使ってgff3ファイルをbedファイルに変換

gffread XXX.gff3 -g genome.fasta --bed -o YYY.bed

 --bed output records in BED format instead of default GFF3

3. awkコマンドを使ってbedファイルからscaffold名(列1), 遺伝子名(列4), start position(列2), end position(列3)の順番，かつタブ区切りで出力

awk -v 'OFS=\t' '{print $1,$4,$2,$3}' YYY.bed > YYY_edit.bed

-v 'OFS=\t'　OFS（Output Field Separator, 出力のフィールド区切り文字）をタブ区切り (\t) に設定

'{print $1,$4,$2,$3}'　print コマンドを使って、YYY.bed の各行から特定の列を取り出して並べ替え

コマンド実行後 のYYY_edit.bed

chr1    geneA   1000    2000

chr2    geneB   1500    2500

chr3    geneC   3000    4000

２．ファイルの結合

用意したゲノムの.faaファイルおよび.gff3ファイルは結合しておく．その際，catではなくawk 1を使用する．

catで連結した場合は改行の処理をミスるのでこの後の処理でエラーが出る．

catコマンドで連結した場合の問題点　https://t-arae.blog/posts/2023/2023-08-31-concatenate-fastq-files-into-one/

コマンド例

awk 1 Mspo_protein.faa Mlac_protein.faa > MCscanX_test.faa

awk 1 Mspo.gff Mlac.gff > MCscanX_test.gff

３．MCScanXの導入

下記のサイトを参考にして，gitからファイルを複製，コンパイルする．

https://kazumaxneo.hatenablog.com/entry/2020/06/20/144410

Username@PCname: /mnt/c/XXX/$ git clone https://github.com/wyp1125/MCScanX.git

Username@PCname: /mnt/c/XXX/$ cd MCScanX/

Username@PCname: /mnt/c/XXX/MCScanX$ make -j 8

テストを実行してうまくいけばok

Username@PCname: /mnt/c/XXX/MCScanX$ ./MCScanX data/os_sb 

※私は使用頻度が高くないので，PATHを通さず，MCScanXディレクトリ直下にディレクトリを作り，そこにファイルを移動して実行している．

４．BLASTP

・BLASTPを使える環境を用意する．

Windowsの場合はwsl2の設定

miniforgeのインストールを先にしておく

condaでインストールする場合の例

(base) Username@PCname: /mnt/c/XXX/MCScanX$  conda create -n MCscanX

(base) Username@PCname: /mnt/c/XXX/MCScanX$  conda activate MCscanX

(MCscanX) Username@PCname: /mnt/c/XXX/MCScanX$  conda install -c conda-forge -c bioconda blast

※業務でcondaを使う場合は有償化されたdefaultリポジトリを使わないように要注意

【Python】Anaconda、MinicondaではなくMiniforgeを使おう

https://qiita.com/Torahugu/items/0efe216da26ee65db547

次のサイトを参考にする　https://kazumaxneo.hatenablog.com/entry/2020/06/20/144410

・makeblastdb

連結した.faaを指定する．-outの名前は任意．諸々のファイルがカレントディレクトリにできる．

(MCscanX) Username@PCname: /mnt/c/XXX/MCScanX$ makeblastdb -in MCscanX_test.faa -dbtype prot -out MCscanX_test

・BLASTP（ゲノムの数が多いとちょっと時間がかかる）-outは任意の名前，拡張子は.blast

(MCscanX) Username@PCname: blastp -query MCscanX_test.faa（連結ファイル） -db MCscanX_test -outfmt 6 -out MCscanX_test.blast -num_threads [threads] -evalue 1e-10

５．MCScanXの実行

MCScanX実行ファイルのある場所にディレクトリを作り，その中に.blastと連結した.gffを移動させておく

MCScanXの実行（e-valueはデフォルトで1e-05，好みで変更する）

※今回はMCScanX_test1というディレクトリの中に.gffおよび.blastのファイルが入っている（MCScanX_test1/MCScanX_test1　前半部分）．

出力ファイルの名前はMCScanX_test1になる（MCScanX_test1/MCScanX_test1　後半部分）．

(MCscanX) Username@PCname: /mnt/c/XXX/MCScanX$ ./MCScanX MCScanX_test1/MCScanX_test1

Reading BLAST file and pre-processing

Generating BLAST list

25029 matches imported (28854 discarded)

6 pairwise comparisons

201 alignments generated

Pairwise collinear blocks written to MCScanX_test1/MCScanX_test1.collinearity [0.835 seconds elapsed]

Tandem pairs written to MCScanX_test1/MCScanX_test1.tandem

Writing multiple syntenic blocks to HTML files

Chr1.html

Chr2.html

Plas1.html

Done! [1.033 seconds elapsed]

こんな感じの出力が得られれば成功．

MCScanX_test1.collinearityというファイルが出力されているはず．

.htmlをブラウザで開くとBLASTで相同性のあった遺伝子IDが分かる．

６．可視化: SynVisio

SynVisioにアクセスする　https://synvisio.github.io/#/

Upload OWN Data to Dashboardから連結した.gffおよび.collinearityをアップロードする→Dashboardに移動

Filter Panelというところから比較するゲノムを選択→GO

色々いじっているうちに使い方が分かるはず．解説動画もあり→　https://www.youtube.com/watch?v=C4fTi9bVHEY

※SynVisioではBLASTPで1つも自他の遺伝子に相同性がなかったゲノムは比較の選択肢に出てこない．

右下のDownload ImagesからSynteny plotまたはdot plotがSVGでダウンロード可能．

７．可視化: AccuSyn

AccuSynにアクセスする　https://synvisio.github.io/#/

Get startedから連結した.gffおよび.collinearityをアップロードする→Dashboardに移動

Connectionsというところから比較するゲノムを選択

こちらも色々いじっているうちに使い方が分かるはず．

※AccuSynでは相同性がなかったゲノムも表示可能．下図の円の上部にある短い配列はプラスミド，Chr1およびChr2の両方と相同性のある遺伝子がないことが分かる．

右下にある↓アイコンからSVGでダウンロード可能．

番外編　Trackを追加する

※試験的な段階です．もっと良いやり方があるといいのですが．．．

・GC contentを追加する（Heatmapで指定）

ゲノム配列（.fna）から5000塩基ずつスライドしながらGC contentを計算する（seqkit）．

cat Mspo.fna | seqkit sliding -s 5000 -W 5000 -g | seqkit fx2tab -n -g >> Mspo_gc_content.txt

Excel上でファイルを区切り位置などで整形して上書き保存

＜整形前＞

1列目 2列目

NZ_CP080203.1_sliding:1-5000 65.66

NZ_CP080203.1_sliding:5001-10000 64.56

NZ_CP080203.1_sliding:10001-15000 67.18

NZ_CP080203.1_sliding:15001-20000 67.1

＜整形後＞

1列目 2列目 3列目 4列目

Chr1 1 5000 65.66

Chr1 5001 10000 64.56

Chr1 10001 15000 67.18

Chr1 15001 20000 67.1

・GC skewを追加する（Line chartで指定）

ゲノム配列（.fna）から5000塩基ずつスライドしながらGC skewを計算する．

cat Mspo.fna | seqkit sliding -s 5000 -W 5000 -g | seqkit fx2tab -n -G >> Mspo_gc_skew.txt

Excel上でファイルを整形して上書き保存

＜整形前＞

1列目 2列目

NZ_CP080203.1_sliding:1-5000 8.07

NZ_CP080203.1_sliding:5001-10000 -3.78

NZ_CP080203.1_sliding:10001-15000 3.72

NZ_CP080203.1_sliding:15001-20000 2.53

＜整形後＞

1列目 2列目 3列目 4列目

Chr1 1 5000 8.07

Chr1 5001 10000 -3.78

Chr1 10001 15000 3.72

Chr1 15001 20000 2.53

※Trackファイルは1列目でゲノム名を識別する（.gffの名前と一致させる必要あり）．

複数のゲノムのTrackを追加する場合は，上記の計算・整形したファイルを連結することで単一のtxtファイルにする．

※SeqkitでFlagを間違えないように注意

https://bioinf.shenwei.me/seqkit/usage/#fx2tab-tab2fx

seqkit fx2tab [flags]

Flags:

-g, --gc print GC content

-G, --gc-skew print GC-Skew

・CDS，tRNA，rRNAの追加および遺伝子の向きによる色分け（Heatmapで指定）

※色分けはヒートマップでわざと差をつけることで識別できるようにできる（絶対的な色の指定のやり方がわからない）

※AccuSynの場合は，色分けができないのでplusとminusでtxtファイルを分けたほうが良い

・A4で取得した.tsvファイルをExcelで開く．

A列，B列，C列，E列，K列が今回は重要．

・Accessionの名前をgffファイルに合わせて変更

・I列でソート

CDS→protein-coding

tRNA，rRNA→そのまま書いてある

それぞれ別のファイルとして保存して整形していく

・CDSの場合: Orientationの列のうち，minus→0，plus→10000に置換（数字は逆でもよい），ヘッダー以外をメモ帳に張り付けてtxtで保存．

1列目 2列目 3列目 4列目

Chr1 861 2372 0

Chr1 2411 3160 0

Chr1 3357 3929 0

Chr1 3982 4659 0

Chr1 4835 5626 10000

Chr1 5659 7140 10000

Chr1 7192 7794 10000

Chr1 7831 9201 10000

Chr1 9210 9989 10000

Chr1 9982 11112 10000

※任意の遺伝子IDのみを抜き出してtrackファイルにしてアップロードするとそこだけ目立つ

transposonやファージなどを表示するとよいかも？

・tRNA，rRNAの場合: Gene Typeの列のうち，tRNA→0，rRNA→10000に置換（数字は逆でもよい），ヘッダー以外をメモ帳に張り付けてtxtで保存．

1列目 2列目 3列目 4列目

Chr1 2018637 2020156 10000

Chr1 2100756 2100872 10000

Chr1 2101027 2104148 10000

Chr1 2104564 2106083 10000

Chr1 156393 156468 0

Chr1 158816 158891 0

Chr1 190275 190348 0

Chr1 190579 190652 0

可視化: Track 

Track File（txtファイル）を指定して一緒にアップロードする（SynVisio）

シンテニープロットを表示した後，上の方にあるSelect TrackでTrackのプロットの種類・カラーパレットを1つずつ設定していく

下記の図の例

Track1　 Heatmap Plasma

Track2 Line Chart Blue

Track3 　 Heatmap Plasma

Track4　 Heatmap Viridis

AccuSynはシンテニープロットを表示後にAdditional tracksからTrackを追加していく

Track追加の操作感はSynVisioに似るが，RのviridisパッケージにあたるカラーパレットがAccuSynでは使えず，０と10000の間で色分けができない

（Track3のように無色と色ありになってしまう）

Show all chromosomesの✓を外し，１つの染色体のみを選択することで単一ゲノムのCircos plot風の表示も可能．

Version1 (20250216) 　作成

Version1.1 (20250217)　seqkit seqによるfastaヘッダーの修正を追記．誤字修正．　

Version2 (20250223)　gffreadについてを追記．