一褪子、PWM與PFM的介紹
motif指的是轉(zhuǎn)錄因子偏好結(jié)合的DNA序列模式或RNA結(jié)合蛋白偏好結(jié)合的序列模式辙售，一般使用PWM來表示motif。
制作PWM的過程如下：
- 首先計算所有序列每個位置的堿基頻數(shù)飞涂，可以得到PFM(Position Frequency Matrix)旦部，這個可以作為ggseqlogo包的輸入來繪制motif圖。
- 通過PFM可以進(jìn)一步計算得到PPM(Position Probability Matrix)较店，表示位置頻率矩陣士八，是計算了每個位置四種堿基的出現(xiàn)頻率后得到的矩陣，這個PPM矩陣可以作為TOMTOM的輸入去與其他motif數(shù)據(jù)庫比對梁呈。
- 將PPM中的頻率除以堿基背景頻率然后取對數(shù)(log2)以后就得到了PWM(Position Weight Matrix)矩陣婚度。
制作PWM的一個示例數(shù)據(jù)可參考我之前寫的博客。
下面介紹一下如何使用R語言來計算motif的PFM和PPM官卡，代碼參考了ggseqlogo的源代碼并進(jìn)行了修改與注釋蝗茁。
二、根據(jù)堿基序列手工制作PFM與PPM
要實現(xiàn)的目標(biāo)：輸入含有n條相同長度的DNA或RNA序列寻咒，返回可以表征該motif的PFM或PPM哮翘。
1. 主要實現(xiàn)函數(shù)
letterMatrix <- function(input){
# Ensure kmers are the same length characters(ggseqlogo)
# 首先要確保輸入的堿基序列的長度都是一致的
seq.len = sapply(input, nchar) # 計算每條序列的堿基數(shù)目
num_pos = seq.len[1] # 第一條序列的堿基數(shù)目
if(! all(seq.len == num_pos)) { # 所有序列的堿基數(shù)目必須一致，不一致則報錯
stop('Sequences in alignment must have identical lengths')
}
# Construct matrix of letters(ggseqlogo)
# 接下來構(gòu)建一個矩陣毛秘，每個元素是一個堿基
split = unlist( sapply(input, function(seq){strsplit(seq, '')}) ) # strsplit可以將字符串切割成單個字符
t( matrix(split, seq.len, length(split)/num_pos) )
}
make_ppm <- function(seqs, ppm=TRUE, seq_type="dna") {
# seqs: A vector of strings, each string is a DNA or RNA sequence
# ppm: Whether to return PPM, default is PPM, else return PFM
# seq_type: Sequence type, can be "dna" of "rna"
letter_mat = letterMatrix(seqs) # 構(gòu)建堿基矩陣饭寺，每一行是一條序列阻课，每一列是堿基位置
# Get namespace(ggseqlogo)
if(seq_type == "dna") {
namespace = c("A", "T", "G", "C")
} else if (seq_type == "rna" ) {
namespace = c("A", "U", "G", "C")
} else {
stop('Wrong seq_type! Must be one of "dna" and "rna".')
}
# Construct PWM(ggseqlogo)
pfm_mat = apply(letter_mat, 2, function(pos.data){ # apply第二個參數(shù)為2，表示對列進(jìn)行操作
# Get frequencies (ggseqlogo)
t = table(pos.data) # 計算該位置四種堿基的頻數(shù)
# Match to aa(ggseqlogo)
ind = match(namespace, names(t)) #
# Create column(ggseqlogo)
col = t[ind] #
col[is.na(col)] = 0
names(col) = namespace
if(ppm) { # 若返回PPM艰匙，則將堿基頻數(shù)除以該列堿基總數(shù)
col = col / sum(col)
}
col # 函數(shù)返回值col
})
num_pos = nchar(seqs[1])
colnames(pfm_mat) = 1:num_pos
pfm_mat
}
用法如下：
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
# 計算PPM(TOMTOM的輸入格式)
ppm <- make_ppm(seqs, ppm=TRUE)
ppm
## 1 2 3 4 5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
# 計算PFM(ggseqlogo的輸入格式)
pfm <- make_ppm(seqs, ppm=FALSE) # ppm=FALSE則輸出PFM
pfm
## 1 2 3 4 5
## A 2 1 1 5 3
## T 0 3 4 0 0
## G 0 1 0 0 2
## C 3 0 0 0 0
2. 實現(xiàn)效果
ggseqlogo包可以接受一個存儲堿基序列的字符串向量或者PFM來繪制motif logo限煞，下面就依次使用這兩種方法對5條DNA和RNA序列繪制motif logo，以驗證手工計算出的PFM與該包計算的結(jié)果一致员凝。
2.1 制作DNA的motif logo
首先用ggseqlogo畫一下這5條序列的motif logo署驻，如下圖所示
# install.packages('ggseqlogo')
library('ggseqlogo')
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
ggseqlogo(seqs)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
image
接下來使用上面自定義的函數(shù)來計算PFM然后繪制motif logo：
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
pfm <- make_ppm(seqs, ppm=FALSE)
pfm
## 1 2 3 4 5
## A 2 1 1 5 3
## T 0 3 4 0 0
## G 0 1 0 0 2
## C 3 0 0 0 0
ggseqlogo(pfm)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
image
可以看到畫出來的圖與直接使用ggseqlogo的結(jié)果是一樣的。
該代碼也可以計算出PPM(TOMTOM需要的格式)绊序，如下所示：
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
pfm <- make_ppm(seqs)
pfm
## 1 2 3 4 5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
2.2 制作RNA的motif logo
首先用ggseqlogo畫motif logo如下：
# install.packages('ggseqlogo')
library('ggseqlogo')
seqs <- c("CGUAA", "AUUAG", "CUAAG", "AUUAA", "CAUAA")
ggseqlogo(seqs)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
image
使用make_ppm函數(shù)計算PFM硕舆，在對RNA進(jìn)行計算時需要指定seq_type="rna"，代碼如下：
seqs <- c("CGUAA", "AUUAG", "CUAAG", "AUUAA", "CAUAA")
pfm <- make_ppm(seqs, ppm=FALSE, seq_type='rna')
pfm
## 1 2 3 4 5
## A 2 1 1 5 3
## U 0 3 4 0 0
## G 0 1 0 0 2
## C 3 0 0 0 0
ggseqlogo(pfm)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
image
上述logo圖與ggseqlogo畫出來的結(jié)果是一樣的骤公。
三抚官、PFM->PPM->PWM
上面實現(xiàn)了根據(jù)堿基序列制作PFM與PPM，但如果已經(jīng)有了PFM呢阶捆？有了PFM以后就可以依次計算PPM和PWM了凌节，下面寫一下PFM->PPM->PWM的代碼。
- PFM->PPM
代碼如下洒试，
pfm2ppm <- function(pfm) {
ppm <- apply(pfm, 2, function(col) {col / sum(col)} ) # 對
return(ppm)
}
# 示例
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
pfm <- make_ppm(seqs, ppm=FALSE)
pfm
## 1 2 3 4 5
## A 2 1 1 5 3
## T 0 3 4 0 0
## G 0 1 0 0 2
## C 3 0 0 0 0
ppm_ori <- make_ppm(seqs)
ppm_ori
## 1 2 3 4 5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
ppm <- pfm2ppm(pfm)
ppm
## 1 2 3 4 5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
ppm_ori == ppm
## 1 2 3 4 5
## A TRUE TRUE TRUE TRUE TRUE
## T TRUE TRUE TRUE TRUE TRUE
## G TRUE TRUE TRUE TRUE TRUE
## C TRUE TRUE TRUE TRUE TRUE
# 兩種方法計算出來的PPM是一致的倍奢。
- PPM->PWM
由PPM到PWM需要先將PPM除以每種堿基的背景頻率(一般情況下四種堿基的背景頻率均為0.25)，然后取對數(shù)(log2)就得到了PWM垒棋，代碼如下：
ppm2pwm <- function(ppm) {
pwm <- log2(ppm / 0.25)
pwm[is.infinite(pwm)] <- 0 # 頻率為0的堿基取對數(shù)以后是負(fù)無窮卒煞，需要將其置為0
return(pwm)
}
seqs <- c("CGTAA", "ATTAG", "CTAAG", "ATTAA", "CATAA")
ppm <- make_ppm(seqs)
ppm
## 1 2 3 4 5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
pwm <- ppm2pwm(ppm)
pwm
## 1 2 3 4 5
## A 0.6780719 -0.3219281 -0.3219281 2 1.2630344
## T 0.0000000 1.2630344 1.6780719 0 0.0000000
## G 0.0000000 -0.3219281 0.0000000 0 0.6780719
## C 1.2630344 0.0000000 0.0000000 0 0.0000000
- 直接讀入PFM矩陣并計算PPM，然后使用ggseqlogo畫圖
seqs_text <- ' 1 2 3 4 5
A 2 1 1 5 3
T 0 3 4 0 0
G 0 1 0 0 2
C 3 0 0 0 0'
seqs_pfm <- read.table(text=seqs_text)
# 計算PPM
pfm2ppm(seqs_pfm)
## X1 X2 X3 X4 X5
## A 0.4 0.2 0.2 1 0.6
## T 0.0 0.6 0.8 0 0.0
## G 0.0 0.2 0.0 0 0.4
## C 0.6 0.0 0.0 0 0.0
# 直接繪圖
library(ggseqlogo)
ggseqlogo(as.matrix(seqs_pfm))
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
image
四叼架、待優(yōu)化
- 優(yōu)化代碼畔裕，使其可自動判斷輸入類型是DNA還是RNA，即無需指定seq_type參數(shù)乖订。
- 博客文章需要改一下數(shù)字扮饶，順便將代碼放上去。
