geom_logoによる配列ロゴの表示

目的

アライメントした塩基配列を得た時に、配列のどの部位にどんな塩基が保存されているかをgeom_logoを使ってr上で視覚的に表示させる。今回は、アライメントされたfasta形式のファイルを使う。下は結果の一例、位置4-6においてCAT配列が保存されているのが分かる。こういうのを配列ロゴ(Sequence logo)と呼ぶらしいのだが何でlogo?。

準備

必要なパッケージをインストール。今回は、fasta形式のファイルを読み込むread.fastaも使うのでseqinrもインストール。

install.package("ggplot2");install.package("ggseqlogo");install.package("seqinr")

Alignmentされたfasta形式のファイルを適当に用意

>seq1
-GTCGA
>seq2
GGACG-

サンプルプログラム

require(ggplot2)
require(ggseqlogo)
require(seqinr)

setwd("C:/Users/hogehoge/") #適当な位置を指定
data <- read.fasta("test_align.fas", seqtype = "DNA", as.string = TRUE)#自分で用意したfastaファイル名を入力
data_m <- list(toupper(c(data[1:length(data)])))

#plot
ggplot() + 

##graphic parameter
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +

geom_logo(data_m, seq_type = "dna", method = "probability")

表示方法

geom_logoは表示形式を、bits形式とprobability形式の二つから選べる。

bits形式

bits形式では、配列の偏りを表示の高さで示してくれるため、どの位置で配列が偏っているか、つまりは保存されているかが一目で分かる。各位置での高さは、エントロピーの最大値から各位置でのエントロピーを引いた値となる。エントロピーの最大値は、DNAやRNA配列などの塩基4種類であれば$log_2 (4)$、タンパク質などでアミノ酸20種類の場合は$log_2 (20)$となる。各位置でのエントロピー$H_i$は、$H_i=-\sum_{i=a,t,g,c}P_ilog_2 Pi$で表される。$P_i (0\leqq P_i\leqq1)$は、各位置での各塩基の出現頻度。

probability形式

probability形式は、各位置での各塩基の出現確立に応じた表示をしてくれる。なので、bit形式で表示させると文字が潰れてしまうような位置10などにおいてどんな塩基が使われているのかも分かる。