R語言繪圖包09--序列分析圖的繪制ggseqlogo

R語言繪圖包系列：

序列分析圖（sequence logo）的序列指的是核苷酸(DNA/RNA鏈中)或氨基酸(在蛋白質(zhì)序列中)。sequence logo圖是用來可視化一段序列某個位點的保守性，根據(jù)提供的序列組展示位點信息。常用于描述序列特征，如DNA中的蛋白質(zhì)結(jié)合位點或蛋白質(zhì)中的功能單元。

ggseqlogo包提供了一種簡便的序列分析圖繪制方法。其主函數(shù)為ggseqlogo()，該函數(shù)可以直接接受已經(jīng)比對好的DNA或者氨基酸序列，但是有一點需要注意：比對好的序列中可以有g(shù)ap，但是所有序列的長度必須相同，否則會報錯。

1. 安裝
#直接從CRAN中安裝
install.packages("ggseqlogo")
#從GitHub中安裝
devtools::install.github("omarwagih/ggseqlogo")
2. 加載數(shù)據(jù)集
  • 演示數(shù)據(jù)集
#加載包
library(ggplot2)
library(ggseqlogo)
#加載數(shù)據(jù)
data(ggseqlogo_sample)

ggseqlogo_sample數(shù)據(jù)集包含了三個列表：
pfms_dna：四種轉(zhuǎn)錄因子的位置頻率矩陣
seqs_aa：一組激動酶底物磷酸化位點序列
seqs_dna：12種轉(zhuǎn)錄因子的結(jié)合位點序列

#pfms_dna
head(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_aa
head(seqs_aa[[1]])
# [1] "VVGARRSSWRVVSSI" "GPRSRSRSRDRRRKE" "LLCLRRSSLKAYGNG" "TERPRPNTFIIRCLQ" "LSRERVFSEDRARFY" "PSTSRRFSPPSSSLQ"
#seqs_dna
seqs_dna$MA0001.1
#  [1] "CCATATATAG" "CCATATATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATAAATAG" "CCATATATGG" "CCATATATGG" "CCAAATATAG" "CCAAATATAG" "CCAAATATAG" "CCATAAATGG"
# [14] "CCATAAATGG" "CCATAAATGG" "CCATAAATGG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAAAATAG" "CCAAATATGG" "CCAAAAATGG" "CCAAAAATGG" "CCATTTATAG" "CCATTTATAG"
# [27] "CCATTTATAG" "CCATTAATAG" "CCATTAATAG" "CCATTAATAG" "CCAATTATAG" "CCATTAATGG" "CCAATAATAG" "CCAATTATGG" "CCAAAATTAG" "CCAAAATTAG" "CCAAAATTAG" "CCACATATAG" "CCATTTTTAG"
# [40] "CCATTTTTAG" "CCATTTTTAG" "CCATAAGTGG" "CTATATATAG" "CCATACATGG" "CTATAAATAG" "CTATAAATAG" "CCTAATATAG" "CCTAATATAG" "CCTAATATAG" "CCAAAACTAG" "CTAAATATAG" "CTAAATATAG"
# [53] "CCATAAACAG" "CTAAAAATAG" "CTAAAAATAG" "CCTAAAATGG" "CCTAAAATGG" "CCAGATTTAG" "CCTTTTATAG" "CCTTTTATAG" "CCAAAAAAAG" "CCAAATAAGG" "CCAGAATTAG" "CCATTTCTGG" "CTAAATTTAG"
# [66] "CCCGATATAG" "CCAAAAATAC" "CCATATTGGG" "CCGTTTTTAG" "CAATTTATAG" "CTAAAAGTAG" "CTAAAAGTAG" "CCATTTTAAG" "CCATTTTAAG" "CCATTTATGT" "CTGTAAATAG" "CCATTTTGGG" "CAATAATTAG"
# [79] "CTCAAAATAG" "CCAATTTAGG" "CCAATTTAGG" "CCTAATACGG" "CTATATTAAG" "CTATAAATGT" "CCTTTTATGT" "CCCTTTTTCG" "CTAAAAGAAG" "CTATATTTGT" "CTATAATTGT" "CTACTTAAAG" "CCGAAAAATG"
# [92] "CCTAAAGTGT" "CCTTTTCTAC" "CCCTCACAAG" "TAATTTCTCG" "GCACACAGAC" "TTAGGATTCG"

可以看到seqs_dna中的MA0001.1對象其實是一個含有97個堿基序列的向量，而且，這97個堿基序列其實都是獨立的字符串元素。因此，如果大家利用R隨機生成一串字符數(shù)相同的字符串向量，也是可以用來繪制seqlogo的。

  • 也可以導(dǎo)入自己的數(shù)據(jù)
    支持兩種格式，序列和矩陣。
# 長度為7的motif。每一行為一條序列，長度相同，每一列的堿基組成代表對應(yīng)位置的堿基偏好性。
fasta = "ACGTATG
ATGTATG
ACGTATG
ACATATG
ACGTACG"

fasta_input <- read.table(fasta, header=F, row.names=NULL)
fasta_input <- as.vector(fasta_input$V1)

# 長度為5的motif矩陣示例，每一列代表一個位置，及堿基在該位置的出現(xiàn)次數(shù)。共4行，每一行代表一種堿基
matrix <- "Base    1    2    3    4    5
A    10    2    0    8    1
C    1    12    1    2    3
G    4    0    9    1    1
T    0    0    0    1    9
"
matrix_input <- read.table(matrix, header=T, row.names=1)
matrix_input <- as.matrix(matrix_input)
3. 可視化
  • 使用數(shù)據(jù)中的位置頻率矩陣生成seqlogo
    ggseqlogo通過method選項支持兩種序列標(biāo)志生成方法：bits和probability。
p1 <- ggseqlogo(seqs_dna$MA0001.1, method="bits")
p2 <- ggseqlogo(seqs_dna$MA0001.1, method="prob")
gridExtra::grid.arrange(p1,p2)
##或者直接繪圖（不使用ggseqlogo），結(jié)果和p1是一樣的
ggplot()+geom_logo(seqs_dna$MA0001.1)+theme_logo()
在上面的seqlogo中，橫軸是每一個堿基的位點，縱軸默認(rèn)情形下是Bits (也就是在不設(shè)置method時，默認(rèn)出的是上面那張圖)。關(guān)于bits和probability的區(qū)別，在此不做贅述，大家可以去維基百科查詢?？傊涀∫稽c，不管是什么指標(biāo)，都是用來衡量某個位點某個堿基出現(xiàn)的可能性的大小。顯然，堿基的字符越高（寬度是相同的），其出現(xiàn)的可能性越大。
  • 使用矩陣?yán)L圖
df <- data.frame(seq=sample(c('A','T','C','G'),1000,replace = T),type=rep(c('gt1','gt2','gt3','gt4'),each=250))
mat <- table(df$seq,df$type)
mat   
#     gt1 gt2 gt3 gt4
#   A  74  60  61  64
#   C  55  59  58  50
#   G  61  61  59  61
#   T  60  70  72  75
ggseqlogo(mat)
  • 序列類型
    ggseqlogo支持氨基酸、DNA和RNA序列類型，默認(rèn)情況下ggseqlogo會自動識別數(shù)據(jù)提供的序列類型，也可以通過seq_type選項直接指定序列類型。
ggseqlogo(seqs_aa$AKT1, seq_type="aa")
  • 自定義字母
    通過namespace選項來定義自己想要的字母類型
#用數(shù)字來代替堿基
seqs_numeric <- chartr("ATGC", "1234", seqs_dna$MA0001.1)
ggseqlogo(seqs_numeric, method="prob", namespace=1:4)
  • 配色
    ggseqlogo可以使用col_scheme參數(shù)來設(shè)置配色方案，具體可參考?list_col_schemes
ggseqlogo(seqs_dna$MA0001.1, col_scheme="base_pairing")

自定義配色：ggseqlogo提供函數(shù)make_col_scheme來自定義離散或者連續(xù)配色方案

# 離散配色
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)

# 連續(xù)配色
cs2 <- make_col_scheme(chars = c("A", "T", "C", "G"), values = 1:4)
ggseqlogo(seqs_dna$MA0001.1, col_scheme=cs2)
  • 同時繪制多個序列標(biāo)志
ggseqlogo(seqs_dna, ncol = 4)
##直接繪圖（不使用ggseqlogo），結(jié)果是一樣的
ggplot()+geom_logo(seqs_dna)+theme_logo()+  
facet_wrap(~seq_group,ncol = 4,scales = "free_x")
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容