まずは、ドメイン構造を書きたいタンパク質のアミノ酸配列をとってきて、fasta形式で保存します。ここではパラスペックルの形成を制御する代表的なRNA結合タンパク質、Fus、Nono、Pspc1、Rbm14、Sfpq、Tardbpを集めてみました。このファイル、一応paraspeckles_fasta.txtでこのページの下のやる気のないリンク、探すのに一苦労のリンク先に置いときます。

次に、このファイルをHMMERに投げて、アミノ酸の何番目から何番目がどのドメインか、という情報をとってきます。HMMERのインストールは毎度おなじみターミナルからのhomebrewで簡単にできるようです。

brew install hmmer

次に、Kazmaxneoさんのページにあるように、タンパク質ドメインのデータベースをとってきて、インデックスをつけときます。ここまでの操作はターミナルで行います。今回はncRNA_reviewというフォルダを作業フォルダにしてやっています。

cd /Users/Nakagawa_imac/GoogleDrive/ncRNA_review
wget ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam32.0/Pfam-A.hmm.gz
gzip -dv Pfam-A.hmm.gz
hmmpress Pfam-A.hmm

あとは、同じフォルダに入れたparaspeckles_fasta.txtをhmmerに投げるだけ。アウトプットファイルはRで読み込みやすいように--domtbloutのオプションで指定しています。他はCPUをどれだけ使うかとか、なんちゃらで、意外と高速に結果を出してくれます。

hmmscan --domtblout paraspeckles_domain.txt --cpu 4 -E 1e-10 Pfam-A.hmm paraspeckles_fasta.txt

さて。これで、ドメインの情報が入ったファイルparaspeckles_domain.txtが作られました。ここから先はRの作業です。 ファイルの読み込みはdata.tableのfreadで行ってます。これがやはり一番便利ですね。ただ、#で始まるコメント行が付いていてなんかややこしそう。コメント行を除くには、、、そう。困った時は検索検索！「R fread comment」あたりで検索するとstack overflowの何やらそれらしいページがヒットしてきます。いつぞやトロント大のモリスさんにベンチ屋なんですがなんちゃってレベルで良いのでインフォを勉強したかったらどうしたらいいんですか？とおそるおそる聞いたら、「ググれ」と即答だったのですが、この手の業界の人って、本当にネットとの相性が良いと言いますか、親切なページがたくさんあって本当に助かります。こいつを読んでいると、どうやらgrepと正規表現でチャチャっとできるみたいです。コメント行を除くとヘッダーがなくなってしまうのでヘッダーなしオプションを指定して、あとはなんかRに怒られたのでfill = TRUEをつけています。ここで取り込んだデータテーブル、たくさん列があるのですが、とりあえずドメイン構造を書くためにはドメインの名前、queryのタンパク質の名前、タンパク質の長さ、ドメインの開始と終わりの情報だけあれば良いので、それらの5列を取り出して名前を付け直しています。

library(data.table)
library(ggplot2)
dt <- fread("grep -v '^#' paraspeckles_domain.txt", header = FALSE, fill = TRUE)
dt <- dt[, c(1, 4, 6, 20, 21)]
colnames(dt) <- c("domain_name", "query_name", "length", "query_from", "query_to")

ここまでデータが整形できれば後は我らが強い味方、ggplotに投げるだけです。geom_segmentでタンパク質の全長を書いておいて、geom_rectでドメインの構造を書いて、それぞれのタンパク質の図をfacetで並べて、後はガチャガチャいろいろある軸やらタイトルのラベルをひたすら消しまくってるだけです。最後は見た目が良くなるサイズでpngファイルに書き出してます。

g <- ggplot(data = dt) +
  geom_segment(aes(x = 0, xend = length, y = 5, yend = 5), color="gray", size = 4) +
  geom_rect(aes(xmin = dt$query_from, xmax = dt$query_to, ymin = 0, ymax = 10, fill = domain_name)) +
  facet_grid(vars(query_name)) +
  theme(strip.text.y = element_text(angle = 0)) +
  #erase miscellaneous elements
  xlab("") +
  ylab("") +
  theme(axis.text.y = element_blank()) +
  theme(axis.text.x = element_blank()) +
  theme(panel.grid.major = element_blank()) +
  scale_y_continuous(breaks=NULL) +
  scale_x_continuous(breaks=NULL, expand = c(0, 0)) +
  theme(strip.text.y = element_text(angle = 0))
g
ggsave(filename = "domain.png", g, width = 5, height = 2)

うんうん。なんだかいい感じ。イラレのトレースの時代にはもう戻れません。正確だし、綺麗だし、何より手軽。情報提供の各サイトを作ってくださった方々に感謝です！！