R包ggseqlogo 繪製seq logo圖

簡介

在生物信息分析中，經常會做序列分析圖（sequence logo），這裡的序列指的是核苷酸(DNA/RNA鏈中)或胺基酸(在蛋白質序列中)。sequence logo圖是用來可視化一段序列某個位點的保守性，據根提供的序列組展示位點信息。常用於描述序列特徵，如DNA中的蛋白質結合位點或蛋白質中的功能單元。

實現以上可視化過程的工具有很多，本文介紹一個使用起來非常簡單，不拖泥帶水的R包ggseqlogo，只要你根據此包要求的數據格式上傳一堆DNA序列或者胺基酸序列，再根據現成的命令流程就能畫出logo圖。

安裝到作圖的代碼如下：

安裝

安裝方式有兩種

#直接從CRAN中安裝install.packages("ggseqlogo")#從GitHub中安裝devtools::install.github("omarwagih/ggseqlogo")

數據加載

ggseqlogo提供了測試數據ggseqlogo_sample。

#加載包library(ggplot2)library(ggseqlogo)#加載數據data(ggseqlogo_sample)

ggseqlogo_sample數據集是一個列表，裡面包含了三個數據集：

seqs_dna:12種轉錄因子的結合位點序列

pfms_dna:四種轉錄因子的位置頻率矩陣

seqs_aa:一組激動酶底物磷酸化位點序列

#seqs_dnahead(seqs_dna)[1]

## $MA0001.1##  [1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG"##  [6] "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG"

#pfms_dnahead(pfms_dna)[1]

## $MA0018.2##   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]## A    0    0   11    0    1    0    2    8## C    1    1    0    9    0    3    7    0## G    1   10    0    2   10    0    1    1## T    9    0    0    0    0    8    1    2

#seqs_aahead(seqs_aa)[1]

## $AKT1##   [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG"##   [4] "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"

外部數據讀入

也可以用自己的數據集，支持兩種格式，序列和矩陣。

# 長度為7的motif。每一行為一條序列，長度相同，每一列的鹼基組成代表對應位置的鹼基偏好性。fasta = "ACGTATGATGTATGACGTATGACATATGACGTACG"fasta_input <- read.table(fasta, header=F, row.names=NULL)fasta_input <- as.vector(fasta_input$V1)# 長度為5的motif矩陣示例，每一列代表一個位置，及鹼基在該位置的出現次數。共4行，每一行代表一種鹼基matrix <- "Base    1    2    3    4    5A    10    2    0    8    1C    1    12    1    2    3G    4    0    9    1    1T    0    0    0    1    9"matrix_input <- read.table(matrix, header=T, row.names=1)matrix_input <- as.matrix(matrix_input)

可視化

ggplot()+geom_logo(seqs_dna$MA0001.1)+theme_logo()

ggseqlogo提供了一個直接繪圖的函數ggseqlogo(),這是一個包裝函數。下面命令結果同上面的。

ggseqlogo(seqs_dna$MA0001.1)

輸入格式

ggseqlogo支持以下幾種類型數據輸入：

下面是使用數據中的位置頻率矩陣生成的seqlogo

ggseqlogo(pfms_dna$MA0018.2)

方法

ggseqlogo通過method選項支持兩種序列標誌生成方法：bits和probability。

p1 <- ggseqlogo(seqs_dna$MA0001.1, method="bits")p2 <- ggseqlogo(seqs_dna$MA0001.1, method="prob")gridExtra::grid.arrange(p1,p2)

序列類型

ggseqlogo支持胺基酸、DNA和RNA序列類型，默認情況下ggseqlogo會自動識別數據提供的序列類型，也可以通過seq_type選項直接指定序列類型。

ggseqlogo(seqs_aa$AKT1, seq_type="aa")

自定義字母

通過namespace選項來定義自己想要的字母類型

#用數字來代替鹼基seqs_numeric <- chartr("ATGC", "1234", seqs_dna$MA0001.1)ggseqlogo(seqs_numeric, method="prob", namespace=1:4)

配色

ggseqlogo可以使用col_scheme參數來設置配色方案，具體可參考?list_col_schemes

ggseqlogo(seqs_dna$MA0001.1, col_scheme="base_pairing")

自定義配色

ggseqlogo提供函數make_col_scheme來自定義離散或者連續配色方案

離散配色

csl <- make_col_scheme(chars = c("A","T", "C", "G"), groups = c("gr1","gr1", "gr2","gr2"), cols = c("purple","purple","blue","blue"))ggseqlogo(seqs_dna$MA0001.1,col_scheme=csl)

連續配色

cs2 <- make_col_scheme(chars = c("A", "T", "C", "G"), values = 1:4)ggseqlogo(seqs_dna$MA0001.1, col_scheme=cs2)

同時繪製多個序列標誌

ggseqlogo(seqs_dna, ncol = 4)

上述命令實際上等同於

ggplot()+geom_logo(seqs_dna)+theme_logo()+  facet_wrap(~seq_group,ncol = 4,scales = "free_x")

自定義高度

通過創建矩陣可以生成每個標誌的高度，還可以有負值高度

set.seed(1234)custom_mat <- matrix(rnorm(20), nrow = 4, dimnames = list(c("A","T","C", "G")))ggseqlogo(custom_mat,method="custom",seq_type="dna")+  ylab("my custom height")

字體

可以通過font參數來設置字體，具體可參考?list_fonts

fonts <- list_fonts(F)p_list <- lapply(fonts, function(f){  ggseqlogo(seqs_dna$MA0001.1,font=f)+ggtitle(f)})do.call(gridExtra::grid.arrange,c(p_list, ncol=4))

注釋

注釋的話跟ggplot2是一樣的

ggplot()+  annotate("rect", xmin = 0.5, xmax = 3.5, ymin = -0.05, ymax = 1.9, alpha=0.1, col="black", fill="yellow")+  geom_logo(seqs_dna$MA0001.1, stack_width = 0.9)+  annotate("segment", x=4, xend = 8, y=1.2, yend = 1.2, size=2)+  annotate("text", x=6, y=1.3, label="Text annotation")+  theme_logo()

圖形組合

將ggseqlogo生成的圖形與ggplot2生成的圖形組合在一起。

p1 <- ggseqlogo(seqs_dna$MA0008.1)+theme(axis.text.x = element_blank())aln <- data.frame(  letter=strsplit("AGATAAGATGATAAAAAGATAAGA", "")[[1]],  species=rep(c("a","b","c"), each=8),  x=rep(1:8,3))aln$mut <- "no"aln$mut[c(2,15,20,23)]="yes"p2 <- ggplot(aln, aes(x, species)) +  geom_text(aes(label=letter, color=mut, size=mut)) +  scale_x_continuous(breaks=1:10, expand = c(0.105, 0)) + xlab('') +  scale_color_manual(values=c('black', 'red')) +  scale_size_manual(values=c(5, 6)) +  theme_logo() +  theme(legend.position = 'none', axis.text.x = element_blank()) bp_data <- data.frame(  x=1:8,  conservation=sample(1:100, 8))p3 <- ggplot(bp_data, aes(x, conservation))+  geom_bar(stat = "identity", fill="grey")+  theme_logo()+  scale_x_continuous(breaks = 1:10, expand = c(0.105, 0))+  xlab("")suppressMessages(require(cowplot))plot_grid(p1,p2,p3,ncol = 1, align = "v")

嚴濤：浙江大學，作物遺傳育種在讀博士。愛鼓搗各種可視化軟體，愛折騰。點擊閱讀原文跳轉其博客。

精品回顧

畫圖三字經 生信視頻 生信系列教程 

心得體會 癌症資料庫 Linux Python 

高通量分析 在線畫圖 測序歷史 超級增強子

培訓視頻 PPT EXCEL 文章寫作 

色彩搭配 圖形排版 互作網絡

後臺回復「生信寶典福利第一波」獲取教程合集

（錯誤矯正基金：如果您在閱讀過程中發現文字或命令錯誤，請留言或加小編微信指出，獲取紅包或累積獎勵。希望大家多監督，反饋。適用於所有原創文章。）