Rのgggenesパッケージで原核生物の遺伝子マップを描画する方法
- ただし、遺伝子の方向の反映とdummyの併用が下記の方法だとうまくいかない
- エキソン-イントロン構造とかの複雑な遺伝子図の描画にはあまり向かないかも
Rのgggenesは遺伝子マップを簡単に描くことができる便利なパッケージ。
作者のgggenesのページ
使うにはまずはRで
ggplot2
gggenes
ggfittext
パッケージをCRAN経由でインストール。
そして
molecule start end gene strand
genome1 723332 724332 dummy 1
genome1 724333 725097 hogE 1
genome1 725111 726544 hugA -1
genome1 726541 727215 hugU -1
genome2 100 1100 dummy 1
genome2 1101 2104 hogE 1
といった感じで描画したい遺伝子の情報をtxtファイルで記述し保存(下記ではtestdata.txtとして保存したものを読み込んでいる)。
moleculeは、生物種名等、遺伝子図を1ラインで表示したい単位ごとに別名で記述。
strandは遺伝子の方向。1が+、-1が-。
dummyは、遺伝子図の縮尺をゲノム間であわせるためのダミー遺伝子。
描画は、例えば下記のようなコマンドで可能。
library(ggplot2)
library(gggenes)
library(ggfittext)
example_genes <- read.table("testdata.txt", header=T)
dummies <- make_alignment_dummies(
example_genes,
ggplot2::aes(xmin = start, xmax = end, y = molecule, fill = gene, id = gene),
on = "dummy"
)
ggplot2::ggplot(example_genes, ggplot2::aes(xmin = start, xmax = end, y = molecule, fill = gene, label = gene)) +
geom_gene_arrow(arrowhead_height = unit(3, "mm"), arrowhead_width = unit(1, "mm")) +
ggplot2::geom_blank(data = dummies) +
ggplot2::facet_wrap(~ molecule, scales = "free", ncol = 1) +
geom_gene_label(align = "left") +
theme_genes()
残念ながら、dummyがある状態でstrandの情報を参照すると、エラーが出て描画ができないので上記ではstrandの情報を参照していない。dummyが無いならstrandの情報は
ggplot2::ggplot(example_genes, ggplot2::aes(xmin = start, xmax = end, y = molecule, fill = gene, label = gene, forward = strand)) +
で反映可能。
上記問題が少なくとも私の中では解決できていないので、-の遺伝子の反転や遺伝子の色分けは、私はこの後イラレでやっている。
そのために、下記のようにとりあえず遺伝子図をpdf形式で保存するのが基本。
pdf("gggenestest.pdf", width=10,height = 5)
ggplot2::ggplot(example_genes, ggplot2::aes(xmin = start, xmax = end, y = molecule, fill = gene, label = gene)) +
geom_gene_arrow(arrowhead_height = unit(3, "mm"), arrowhead_width = unit(1, "mm")) +
ggplot2::geom_blank(data = dummies) +
ggplot2::facet_wrap(~ molecule, scales = "free", ncol = 1) +
geom_gene_label(align = "left") +
theme_genes()
dev.off()