ChIP-Seq数据挖掘系列-2: Motif 分析(2) - HOMER Motif 分析基本步骤

Posted 2023-04-15

参考技术A

在基因组调控元件分析中，HOMER 可以用于发现新的motif。HOMER 通过比较两个序列集，再使用ZOOPS scoring (zero or one occurrence per sequence)和超几何检验进行富集分析。HOMER主要被用于 ChIP-Seq 和 promoter 分析，但是核酸序列motif寻找问题都可以尝试使用HOMER。

HOMER预测Motif 需要的两个序列集

HOMER 分析基本步骤：
1. 预处理
1.1 提取序列 (findMotifs.pl/findMotifsGenome.pl)
提供的数据是基因组位置信息，就需要提取对应的DNA信息；提供基因号时，需要选择启动子区域。

1.2 背景选择 (findMotifs.pl/findMotifsGenome.pl)
未指定背景序列时，HOMER 会自动选择。
对基因组某些区域进行分析时，从基因组随机选择GC含量一致的序列作为背景序列。
对启动子进行分析时，除用来分析外的所有启动子将被作为背景。
自定义背景使用参数"-bg <file>"。

1.3 GC 标准化 (findMotifs.pl/findMotifsGenome.pl)
目标序列和背景序列会基于GC含量按5%作为bin 查看GC含量的分布。背景序列会得到权值，从而使得其GC含量分布与目标序列一致。
ChIP-Seq 实验得到序列GC含量。

1.4 自动标准化 (New with v3.0, homer2/findMotifs.pl/findMotifsGenome.pl)
需要分析的序列除了GC含量会带来误差，其他的生物学现象，外显子中密码子偏好性或测序实验中偏好性都会影响分析。对于足够强的偏差，HOMER 会自动追踪目标序列和背景中显著差异的特征序列，并通过调整背景序列的权重来平衡输入数据和背景中短寡聚核酸序列不平衡。短寡聚核酸序列长度可以通过参数"-nlen <#>"指定。

2. 重头预测Motifs (homer2)
默认情况下，HOMER 调用homer2 进行motif 分析；通过参数"-homer1" 可以指定老版本工具。

2.1 将输入序列解析为寡聚核苷酸序列
将输入序列按照motif 长度期望值解析为寡聚核苷酸序列，以及创建Oligo 数据表。Oligo 数据表中记录着每条oligo 在目标序列和背景中被发现的次数。

2.2 Oligo 自动标准化 (可选)
2.3 全局搜索阶段
Oligo 表格信息构建好之后，HOMER 对富集的Oligo 进行全局搜索。如果一个Motif是富集的，那么属于这个Motif的Oligo 也应该会富集。首先，HOMER 会搜索可能富集的Oligo 。HOMER 允许错配 ，使用参数"-mis <#>" 调节允许的错配数目。

2.3.1 Motif 富集分析
Motif 富集分析使用超几何分布和二项式分布。一般情况下，序列较多或者背景序列远远多于目标序列，二项式分布计算比较快，因此findMotifsGenome.pl默认使用二项式分布；当自定义背景序列时，这时序列较少，使用超几何检验比较好("-h")。findMotifs.pl用于启动子分析，并且默认使用超几何检验。

2.4 矩阵优化
2.5 Mask and Repeat
当最优oligo被优化成motif后，motif 对应的序列从要分析的数据中移除，接下来再分析最优的.....直到 25(默认值，"-S <#>")个motifs 被发现。

3. 计算已知Motifs是否富集 (homer2)
3.1 导入Motif库
为了搜索输入数据中已知Motifs ，HOMER 可以输入已知Motifs 数据，可以时HOMER 默认的 ("data/knownTFs/known.motifs")，也可以是自己构建("-mknown <file>") 。

3.2 筛选每一个Motif
对于每个motif，HOMER 计算丰度（包含motif的序列/background sequences）， ZOOPS (zero or one occurence per sequence)计数以及使用超几何检验或二项式计算显著性。

4. Motif 分析结果
4.1 Motif Files (homer2, findMotifs.pl, findMotifsGenome.pl)
" .motif"包含motifs的信息
" .motif"文件格式：

一个motif 的信息分为一块。motif 信息首行是motif 各种统计信息；其他行对应各个A/C/G/T的占比。
motif 信息首行解析：

4.2 重头预测的 motif (findMotifs.pl/findMotifsGenome.pl/compareMotifs.pl)
首先会对motif进行去冗余，将每个motif 的概率矩阵转换为向量，求motif之间的Pearson 相关性。
html 结果：

4.2 已知 motif 的富集情况

参考：
Homer

ChIP-Seq 数据挖掘系列文章目录：
ChIP-Seq数据挖掘系列-1:Motif 分析(1)-HOMER 安装
ChIP-Seq数据挖掘系列-2: Motif 分析(2) - HOMER Motif 分析基本步骤
ChIP-Seq数据挖掘系列-3: Motif 分析(3) - 利用ChIP-Seq结果在基因组区域中寻找富集的Motifs
ChIP-Seq数据挖掘系列-4: liftOver - 基因组坐标在不同基因组注释版本间转换
ChIP-Seq数据挖掘系列-5.1: ngs.plot 可视化ChIP-Seq 数据
ChIP-Seq数据挖掘系列-5.2: ngs.plot 画图工具ngs.plot.r 和 replot.r 参数详解

r 导入荷马ChIP-Seq Motif数据

# Daniel Cook 2014
# Danielecook.com
#
# Use this function to import ChIP Seq Data generated by Homer. This data is generated using homers findMotifsGenome.pl
# command with the '-find <motif file>' argument. Generate output
# looks like this:
#
# 1. Peak/Region ID
# 2. Chromosome
# 3. Start
# 4. End
# 5. Strand of Peaks
# 6-18: annotation information
# 19. CpG%
# 20. GC%
# 21. Motif Instances <distance from center of region>(<sequence>,<strand>,<conservation>)
#
# Getting the information from column 21 for plotting, for example, frequency can be tricky. This funciton aims to help that
# by converting the data format from wide to long. Currently, annotation columns are ignored but you can modify the function to
# fix.

library(stringr)
library(splitstackshape)
library(reshape2)
library(dplyr)
suppressMessages(library(data.table))

import_peak_file <- function(filename, peak_type) {
  df <- fread(filename)
  setnames(df, names(df)[1],"Peak")
  setnames(df, names(df)[9:length(names(df))] , str_extract(names(df)[9:length(names(df))], "[A-Za-z0-9]+"))
  motifs <- names(df)[9:length(names(df))]
  # Reshape latter columns

  # Remove extra columns if present
  df$Annotation <- NULL
  df$"Focus Ratio/Region Size" <- NULL
  df$Detailed <- NULL
  df$Distance <- NULL
  df$Nearest <- NULL
  df$Entrez <- NULL
  df$Nearest <- NULL
  df$Nearest <- NULL
  df$Nearest <- NULL
  df$Gene <- NULL
  df$Gene <- NULL
  df$Gene <- NULL
  df$Gene <- NULL
  
  # Melt Again
  df <- melt(df, id.vars = 1:8)
  df$value <- str_replace_all(df$value, "),", "|")
  df <- filter(df, value != "")
  df <- concat.split(as.data.frame(df), split.col = c("value"), sep="|")
  df$value <- NULL
  
  df <- rename(df, c("variable"="motif_name"))
  df <- melt(df, id.vars = 1:9)
  df$value_1 <- NULL
  df$variable <- NULL
  df <- filter(df, value != "")
  
  df$MOTIF <- str_extract(df$value,"[ATCG]+")
  df$POS <- as.integer(str_extract(df$value,"[-0-9]+"))
  df$STRAND <- str_extract(df$value,"(,[+|-],)")
  df$STRAND <- str_replace_all(df$STRAND,",","")
  df$value <- NULL
  df <- filter(df, !is.na(POS))
  df$peak_type <- peak_type
  df
}