2020 11/2 誤字修正
先日紹介したneedleallやvsearchによるall versus allの配列比較のテキスト出力をもとに、ヒートマップで視覚化する。ここではggplot2パッケージを使う。
EMBOSS needleallによるall versus allの配列比較
１、配列の準備
all versus allの配列比較に使う塩基配列は、ある程度長くて全長がよく似ていれば何でも良いが、ここでは大腸菌と赤痢菌の16S rRNA遺伝子を比較してみる。（このブログでは何でも良いが）一般論としてよく分からない株を比較しても議論できないので、まずLPSNでEscherichia coli とShigella dysenteriae のタイプストレインを調べ、その遺伝子を使う。
Escherichia coli (Migula 1895) 、1895年に"Bacillus coli"として再分類（wiki））は
ATCC 11775; CCUG 24; CCUG 29300; CIP 54.8; DSM 30083; JCM 1649; LMG 2092; NBRC 102203; NCCB 54008; NCTC 9001、がタイプストレインとして登録されている。
Shigella dysenteriae (Shiga 1897)は
ATCC 13313; CIP 57.28; DSM 4781; NCTC 4837、がタイプストレインとして登録されている。
NCBIのbrowse by organism機能（紹介）を使って、これらのタイプストレインの利用可能なゲノムデータを探した。Escherichia coliで探すといきなりATCC 11775の完全長ゲノムが見つかった（accession: NZ_CP033092のバージョン２であるNZ_CP033092.2）。これをダウンロードした（*2）。
（K-12は病原性がなくて弱い大腸菌）
次にShigella dysenteriaeだが、ATCC 13313の完全長ゲノムが登録されている（NZ_CP026774.1）。これをダウンロードする。
ここでは2つのゲノムファイルだけ取ってくれば良いので問題ないが、タイプストレインのゲノムを自動で大量に取ってくるのは意外にも難しい（*1）。
ダウンロードしたゲノムからbarrnap（紹介）を使ってSSU rRNA（以後16S rRNA）を取り出す（16S rRNA配列も複数登録されており、アノテーションファイル中にも当然存在すると思われるが、ここでは演習としてゲノムから取り出す）。トラブルの元なのでファイル名に空白は入れないこと。
#E.coli_ATCC11775
barrnap --outseq E.coli_ATCC11775_rRNA.fna --kingdom bac --threads 4 E.coli_ATCC11775.fna > E.coli_ATCC11775_rRNAs.gff
#S.dysenteriae_ATCC13313
barrnap --outseq S.dysenteriae_rRNAs.fa --kingdom bac --threads 4 S.dysenteriae_ATCC13313.fna > S.dysenteriae_rRNAs.gff
- --kingdom Kingdom: euk mito bac arc (default 'bac')
- --lencutoff [n.n] Proportional length threshold to label as partial (default '0.8')
- --reject [n.n] Proportional length threshold to reject prediction (default '0.25')
- --evalue [n.n] Similarity e-value cut-off (default '1e-06')
barrnapの出力から16S rRNAを取り出す。S.dysenteriae_ATCC13313ゲノムからは16S rRNAが７コピー見つかった。
（barrnapで撮り出されたrRNAは向きが統一されており、向きを考えずそのまま配列比較に使える）
E.coli_ATCC11775_rRNAsゲノムからも16S rRNAは7コピー見つかった。
16S rRNAの配列だけを取り出し、16S.faとして１つのファイルに保存した。長さを確認する。
$ s 16S.fa
file format type num_seqs sum_len min_len avg_len max_len
16S.fa FASTA DNA 14 21,532 1,538 1,538 1,538
全て1,538 bp
ヘッダを確認する。
$ grep -n ">" 16S.fa
$ grep -n ">" 16S.fa
1:>16S_rRNA::NZ_CP026774.1:2536626-2538164(+)
3:>16S_rRNA::NZ_CP026774.1:2636479-2638017(+)
5:>16S_rRNA::NZ_CP026774.1:2768885-2770423(+)
7:>16S_rRNA::NZ_CP026774.1:2810847-2812385(+)
9:>16S_rRNA::NZ_CP026774.1:3061594-3063132(+)
11:>16S_rRNA::NZ_CP026774.1:3646780-3648318(+)
13:>16S_rRNA::NZ_CP026774.1:1662010-1663548(-)
15:>16S_rRNA::NZ_CP033092.2:458561-460099(+)
17:>16S_rRNA::NZ_CP033092.2:1285125-1286663(+)
19:>16S_rRNA::NZ_CP033092.2:3780640-3782178(-)
21:>16S_rRNA::NZ_CP033092.2:4551515-4553053(-)
23:>16S_rRNA::NZ_CP033092.2:4591684-4593222(-)
25:>16S_rRNA::NZ_CP033092.2:4844587-4846125(-)
27:>16S_rRNA::NZ_CP033092.2:4726193-4727731(-)
grepで検索するとき、grep -n > 16S.fa としてはならない。
２、16S配列の総当たり比較
ここではvsearchを使う。
vsearch --allpairs_global 16S.fa --id 0.9 --uc vsearch_output
$ head vsearch_output
vsearch_outputの4,9,10列目を取り出す。不要な行（先頭がアスタリスク）があるので、パイプしてgrepで除く（grep -vでマッチは除外）。
cut -f 9,10,4 vsearch_output | grep -v ^* - > seq_ide
$ head -n 5 seq_ide
99.2 16S_rRNA::NZ_CP026774.1:1662010-1663548(-) 16S_rRNA::NZ_CP033092.2:4726193-4727731(-)
99.1 16S_rRNA::NZ_CP026774.1:1662010-1663548(-) 16S_rRNA::NZ_CP033092.2:458561-460099(+)
99.1 16S_rRNA::NZ_CP026774.1:1662010-1663548(-) 16S_rRNA::NZ_CP033092.2:1285125-1286663(+)
99.1 16S_rRNA::NZ_CP026774.1:1662010-1663548(-) 16S_rRNA::NZ_CP033092.2:3780640-3782178(-)
99.1 16S_rRNA::NZ_CP026774.1:1662010-1663548(-) 16S_rRNA::NZ_CP033092.2:4551515-4553053(-)
遺伝子名に冗長な情報が多く見えづらい。
見づらいのでヘッダの不要な情報を消す。アクセッションIDも株名に置換する。
sed -e "s/16S_rRNA:://g" \
-e s/$'NZ_CP033092.2:'/S.dysenteriae_/g \
-e s/$'NZ_CP026774.1:'/E.coli_/g \
seq_ide > seq_rename
$ head -n 5 seq_rename
99.2 E.coli_1662010-1663548(-) S.dysenteriae_4726193-4727731(-)
99.1 E.coli_1662010-1663548(-) S.dysenteriae_458561-460099(+)
99.1 E.coli_1662010-1663548(-) S.dysenteriae_1285125-1286663(+)
99.1 E.coli_1662010-1663548(-) S.dysenteriae_3780640-3782178(-)
99.1 E.coli_1662010-1663548(-) S.dysenteriae_4551515-4553053(-)
seq_renameを開き、先頭行にidentity<tab>genome1<tab>genome2をつけて上書き保存した。
３、結果の視覚化
Rのggplot2を使う。
Rのコンソールに切り替え、ファイルを読み込む。
#R
#データの読み込み（カレントにないならフルパスで指定）
input <- read.table("seq_rename", header=T, sep="\t")
> head(input)
identity genome1 genome2
1 99.2 E.coli_1662010-1663548(-) S.dysenteriae_4726193-4727731(-)
2 99.1 E.coli_1662010-1663548(-) S.dysenteriae_458561-460099(+)
3 99.1 E.coli_1662010-1663548(-) S.dysenteriae_1285125-1286663(+)
4 99.1 E.coli_1662010-1663548(-) S.dysenteriae_3780640-3782178(-)
5 99.1 E.coli_1662010-1663548(-) S.dysenteriae_4551515-4553053(-)
6 99.1 E.coli_1662010-1663548(-) S.dysenteriae_4591684-4593222(-)
視覚化する。
#ggplot2
library("ggplot2")
#default settings
ggplot(data = input, aes(x=genome1, y=genome2, fill=identity)) + geom_tile()
#カスタム。99-100%、色は青。内部に数値。figのraioを1
ggplot(data = input, aes(x=genome1, y=genome2, fill=identity)) + geom_tile(color = "white")+
scale_fill_gradient2(low = "#ffffff", high = "#0000cd", mid = "#ffffff",
midpoint = 99.5, limit = c(99,100), space = "Lab",
name="16S rRNA\nsequence identity (%)") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, vjust = 1,
size = 9, hjust = 1))+
coord_fixed() + theme(aspect.ratio = 1) +
geom_text(aes(label = identity), size = 3)
出力された。一番低い組み合わせでも99.1%の配列同一性があった。
一部欠損があるが、これはvsearchが計算していない配列の組み合わせがあるため。EMBOSS neddleallではこうゆうことはなくて、14x14比較なら196行出力（& 最後に4行のコメント行）される。また、色が二層性で綺麗に分かれてしまったのは、シゲラ のゲノムとE.coliのゲノムの16S rRNAはいずれのコピーも互いにペアワイズ比較すると99.1程度の配列同一性しかなく、同じゲノム内のコピーはほぼ100%合致したため。比較するデータがヒートマップで視覚化するには不向きだったということになる（= 総当たりの表を示せば議論できる。あえてヒートマップで視覚化する必要はない）。
2020 10/31 追記
もう少し距離がある配列間の相同性を視覚化してみる。脊索動物のピルビン酸デヒドロゲナーゼ複合体サブユニット遺伝子５つを使う。
#1 配列のall vs all比較
needleall PDHA1proteins.faa PDHA1proteins.faa out -gapopen 10.0 -gapextend 0.5 -aformat pair
#2 scirptを書いてidentitywの数値を取り出す("# Identity"で始まるの行の()内がIdentity、"# 1:"が比較元の名前、"# 2:"が比較先の名前
perl extracct_best_hit.pl -i out -o out
2の出力
| genome1 | genome2 | identity |
| NP_000275.1[Homo | NP_000275.1[Homo | 100 |
| XP_010641961.1[Fukomys | NP_000275.1[Homo | 99.2 |
| XP_006911688.1[Pteropus | NP_000275.1[Homo | 97.4 |
| XP_021570775.1[Carlito | NP_000275.1[Homo | 95.8 |
| PNI15832.1 | NP_000275.1[Homo | 93.8 |
| NP_000275.1[Homo | XP_010641961.1[Fukomys | 99.2 |
| XP_010641961.1[Fukomys | XP_010641961.1[Fukomys | 100 |
| XP_006911688.1[Pteropus | XP_010641961.1[Fukomys | 98.2 |
| XP_021570775.1[Carlito | XP_010641961.1[Fukomys | 95.5 |
| PNI15832.1 | XP_010641961.1[Fukomys | 93 |
| NP_000275.1[Homo | XP_006911688.1[Pteropus | 97.4 |
| XP_010641961.1[Fukomys | XP_006911688.1[Pteropus | 98.2 |
| XP_006911688.1[Pteropus | XP_006911688.1[Pteropus | 100 |
| XP_021570775.1[Carlito | XP_006911688.1[Pteropus | 93.8 |
| PNI15832.1 | XP_006911688.1[Pteropus | 91.3 |
| NP_000275.1[Homo | XP_021570775.1[Carlito | 95.8 |
| XP_010641961.1[Fukomys | XP_021570775.1[Carlito | 95.5 |
| XP_006911688.1[Pteropus | XP_021570775.1[Carlito | 93.8 |
| XP_021570775.1[Carlito | XP_021570775.1[Carlito | 100 |
| PNI15832.1 | XP_021570775.1[Carlito | 90 |
| NP_000275.1[Homo | PNI15832.1 | 93.8 |
| XP_010641961.1[Fukomys | PNI15832.1 | 93 |
| XP_006911688.1[Pteropus | PNI15832.1 | 91.3 |
| XP_021570775.1[Carlito | PNI15832.1 | 90 |
| PNI15832.1 | PNI15832.1 | 100 |
Rのコンソールに移動する。上の表をコピペ読み込み。
#mac
input <- read.table(pipe("pbpaste"), header = TRUE)
#windows
input <- read.table("clipboard"), header = TRUE)
#またはファイルinput.tsvに保存し、データフレームに読み込む。
input <- read.table("input.tsv", header=T, sep="\t")
#視覚化する
#ggplot2
library("ggplot2")
#default settings
ggplot(data = input, aes(x=genome1, y=genome2, fill=identity)) + geom_tile()
#custom1 - 90-100%、青のグラデーション。内部にidentityの数値、figの比率は1
ggplot(data = input, aes(x=genome1, y=genome2, fill=identity)) + geom_tile(color = "white")+
scale_fill_gradient2(low = "#ffffff", high = "#0000cd", mid = "#ffffff",
midpoint = 93, limit = c(90,100), space = "Lab",
name="PDHA\nsequence identity (%)") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, vjust = 1,
size = 9, hjust = 1))+
coord_fixed() + theme(aspect.ratio = 1) +
geom_text(aes(label = identity), size = 5)
#default settings
#custom1
クラスタリングしたり本格的に計算してから視覚化する場合、データフレームをas.matrixで行列に変換してから進める。
参考
大腸菌研究の歴史
https://www.sbj.or.jp/wp-content/uploads/file/sbj/9010/9010_yomoyama-1.pdf
関連
全配列間の多重整列を実行し、その結果から距離行列ファイルを出す。
*1
最近出たGTDBの論文（DOI: https://doi.org/10.1038/s41587-020-0501-8）では、LPSNやGTDB、bacdiveとの識別子の互換性が乏しい点についての苦労話がある。現状、タイプストレインのゲノムだけ公共データベースから大量にfetchするのは骨が折れる。
*2
タイプストレインは有用だが、完全長ゲノム情報が利用できない株がたくさん存在する。
乖離している理由だが、それには第一に、いわゆる細菌学は百年よりずっと昔から続く歴史ある研究分野であるのに対し、細菌ゲノムをサンガーシークエンシングやHTSのシークエンサーで読んでゲノムベースで取り組む研究手法は、登場から数十年程度の歴史しかない、相対的にずっと新しい研究分野であることを考慮する必要がある。
