#目的
アライメントした塩基配列を得た時に、配列のどの部位にどんな塩基が保存されているかを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などにおいてどんな塩基が使われているのかも分かる。
