R語言-制作motif的PWM

一、PWM與PFM的介紹

motif指的是轉(zhuǎn)錄因子偏好結(jié)合的DNA序列模式或RNA結(jié)合蛋白偏好結(jié)合的序列模式，一般使用PWM來表示motif。

制作PWM的過程如下：

  1. 首先計(jì)算所有序列每個(gè)位置的堿基頻數(shù)，可以得到PFM(Position Frequency Matrix)，這個(gè)可以作為ggseqlogo包的輸入來繪制motif圖。
  2. 通過PFM可以進(jìn)一步計(jì)算得到PPM(Position Probability Matrix)，表示位置頻率矩陣，是計(jì)算了每個(gè)位置四種堿基的出現(xiàn)頻率后得到的矩陣，這個(gè)PPM矩陣可以作為TOMTOM的輸入去與其他motif數(shù)據(jù)庫比對(duì)。
  3. 將PPM中的頻率除以堿基背景頻率然后取對(duì)數(shù)(log2)以后就得到了PWM(Position Weight Matrix)矩陣。

制作PWM的一個(gè)示例數(shù)據(jù)可參考我之前寫的博客。

下面介紹一下如何使用R語言來計(jì)算motif的PFM和PPM，代碼參考了ggseqlogo的源代碼并進(jìn)行了修改與注釋。

二、根據(jù)堿基序列手工制作PFM與PPM

要實(shí)現(xiàn)的目標(biāo)：輸入含有n條相同長度的DNA或RNA序列，返回可以表征該motif的PFM或PPM。

1. 主要實(shí)現(xiàn)函數(shù)

letterMatrix <- function(input){
  # Ensure kmers are the same length characters(ggseqlogo)
  # 首先要確保輸入的堿基序列的長度都是一致的
  seq.len = sapply(input, nchar) # 計(jì)算每條序列的堿基數(shù)目
  num_pos = seq.len[1] # 第一條序列的堿基數(shù)目
  if(! all(seq.len == num_pos)) { # 所有序列的堿基數(shù)目必須一致，不一致則報(bào)錯(cuò)
    stop('Sequences in alignment must have identical lengths')
  }

  # Construct matrix of letters(ggseqlogo)
  # 接下來構(gòu)建一個(gè)矩陣，每個(gè)元素是一個(gè)堿基
  split = unlist( sapply(input, function(seq){strsplit(seq, '')}) ) # strsplit可以將字符串切割成單個(gè)字符

  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第二個(gè)參數(shù)為2，表示對(duì)列進(jìn)行操作
    # Get frequencies (ggseqlogo)
    t = table(pos.data) # 計(jì)算該位置四種堿基的頻數(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")
# 計(jì)算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

# 計(jì)算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. 實(shí)現(xiàn)效果

ggseqlogo包可以接受一個(gè)存儲(chǔ)堿基序列的字符串向量或者PFM來繪制motif logo，下面就依次使用這兩種方法對(duì)5條DNA和RNA序列繪制motif logo，以驗(yàn)證手工計(jì)算出的PFM與該包計(jì)算的結(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ù)來計(jì)算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é)果是一樣的。

該代碼也可以計(jì)算出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ù)計(jì)算PFM，在對(duì)RNA進(jìn)行計(jì)算時(shí)需要指定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

上面實(shí)現(xiàn)了根據(jù)堿基序列制作PFM與PPM，但如果已經(jīng)有了PFM呢？有了PFM以后就可以依次計(jì)算PPM和PWM了，下面寫一下PFM->PPM->PWM的代碼。

  1. PFM->PPM

代碼如下，

pfm2ppm <- function(pfm) {
  ppm <- apply(pfm, 2, function(col) {col / sum(col)} ) # 對(duì)
  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

# 兩種方法計(jì)算出來的PPM是一致的。
  1. PPM->PWM

由PPM到PWM需要先將PPM除以每種堿基的背景頻率(一般情況下四種堿基的背景頻率均為0.25)，然后取對(duì)數(shù)(log2)就得到了PWM，代碼如下：

ppm2pwm <- function(ppm) {
  pwm <- log2(ppm / 0.25)
  pwm[is.infinite(pwm)] <- 0 # 頻率為0的堿基取對(duì)數(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
  1. 直接讀入PFM矩陣并計(jì)算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)
# 計(jì)算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)化

  1. 優(yōu)化代碼，使其可自動(dòng)判斷輸入類型是DNA還是RNA，即無需指定seq_type參數(shù)。
  2. 博客文章需要改一下數(shù)字，順便將代碼放上去。
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

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

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