motif分析——从实战到原理(Homer篇)

一、不求甚解系列

软件下载及配置

  1. conda安装: conda install -c bioconda homer

  2. 使用configureHomer.pl完成HOMER软件的配置

# 下载配置文件
wget http://homer.salk.edu/homer/configureHomer.pl 
​
# 使用配置文件进行软件配置
perl configureHomer.pl -install
​
# 下载 hg19 人的参考基因组
# 将hg19替换为mm10可下载mm10 小鼠的参考基因组
# 因为课题组人和小鼠的研究居多,这里只以这两者举例
perl configureHomer.pl -install hg19 

Motif 分析及peak注释

此处以 MACS2 软件call peak得到的bed文件举例:

HOMER软件接受两种文件格式,这里不展开介绍文件的具体格式,细节可看使用详解系列

findMotifsGenome.pl H3K4Me3.bed hg19 H3K4Me3_motif -len 8,10,12
# peak文件:ChIP-Seq_H3K4Me3_1_homer.bed
# 基因组版本:hg19
# 输出路径:ChIP-Seq_H3K4Me3_1_motifDir
# motif长度:-len 8,10,12 
# motif的软件默认长度为8,10,12
​
annotatePeaks.pl H3K4Me3.bed hg19 1> ChIP-Seq_H3K4Me3_1_peakAnn.xls 2>H3K4Me3_annLog.txt </pre>

结果展示

软件运行后的结果全部放在我们制定的输出文件路径下,将该目录下所有文件下载到本地,打开homerResults.html 文件,即可得到如下结果展示:

随便找的一个例子,不一定是H3K4me3的motif

二、使用详解

HOMER 一个常用的motif分析软件。它通过比较两个序列集,并使用ZOOPS scoring和超几何分布(或者负二项分布)进行motif的富集分析。它主要用于ChIP-seq和promoter分析,但也可以用于核酸序列的motif分析问题。

HOMER软件可以进行多种类型的motif分析,如 promoter motif analysis ,基因组位置motif分析(ChIP-seq分析中的motif分析),利用自定义的fasta文件进行motif分析,RNA序列的motif分析(分析CLIP-seq数据中的RNA binding elements)

HOMER进行motif分析时,需要两个数据集:

  1. 感兴趣的目标序列,如ChIP-seq分析中的peak文件

  2. 背景序列,如ChIP-seq分析中的物种全基因组序列

HOMER包的介绍与安装

HOMER包分为四种:

  1. SOFTWARE: 软件包,里面包括所有的代码以及一些常见的数据文件,如motif矩阵

  2. ORGANISMS:物种特异性数据包,包括一些指定物种的accession conversion数据,gene description数据,GO 分析数据。大部分数据都是以NCBI的gene数据库为基础的

  3. PROMOTERS:启动子数据包,对promoter进行motif分析时所需要的相关数据和promoter序列信息。大部分数据都是基于RefSeq 转录组。

  4. Genomes:基因组数据包,包括一些基因组序列和注释信息

HOMER软件的安装可看不求甚解部分。HOMER软件刚安装之后,是不包含任何序列信息的。我们可以使用 configureHOMER.pl 脚本下载所需数据。 我们可以使用 perl configureHOMER.pl -list 命令得到已有的HOMER包。

有些软件包的名称后面会跟有 -o, -p, -g等参数,这是HOMER为了区分同一物种不同类型数据所添加的。如 human-p。其中 -o 表示organism,-p 表示promoter,-g 表示genome。 这里给出几个软件包安装与删除的例子:

# install 表示安装
perl configureHomer.pl -install mouse # 安装小鼠的promoter数据集
perl configureHomer.pl -install hg19 # 安装人的hg19版本的基因组数据集
​
# remove 表示删除
perl configureHomer.pl -remove mm10 # 删除小鼠mm10版本的基因组数据

此外,需要注意的是,当我们下载基因组数据时,同时会自动下载物种包的数据(比如,下载hg19数据时,会自动下载human数据)。每次我们下载某个物种的基因组数据或者promoter数据时,它都会自动替我们检查是否下载了对应的物种数据。

motif分析

正如前面所说,HOMER可以进行多种类型的motfi分析。这里我们以最常见的ChIP-seq流程中的peaks motfi分析举例。该分析需要使用 findMotifsGeome.pl 脚本。

该脚本的输入文件支持两种类型的文件,一种是常见的bed文件格式(MACS2软件call peak得到的bed文件即可直接使用),另一种是HOMER软件指定使用的peak文件格式。 下面详细讲一下这两种文件格式:

  1. bed文件格式:findMotifsGeome.pl 脚本用到的信息为bed文件的前六列,分别是 chr, start, end, peak ID, score(HOMER用不到该信息,但是必须有), strand

  2. peak文件格式:该文件格式需要五列信息(使用Tab分隔),分别是 peak ID , chr , start , end ,strand

当我们准备分析一个peak/bed文件的motif时,需要运行以下命令:

findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size # [options]
# 比如: findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput -size 200 -mask
  • -mask : 该参数告诉motif分析程序,在得到一个可能的motif之后,在后续的motif分析中是否排除该motif的影响。有点类似于抽样调查中的无放回抽样。

  • -size: 指定用于motif分析的片段长度,默认为200; -size given 设置片段大小为目标序列长度。越大,motif分析所需要的计算资源越多。

下面介绍一些常用的其他参数:

  • -p : 设置线程数-S : 所寻找的motif数目,默认为25。25已经不算少了,如果自定义数目,推荐设置更少而不是更多。

  • -mis : 所允许的错配数,默认为2

  • -cpg : 在对背景序列和目标序列进行normalization时,使用CpG%矫正而非GC%。

  • norevopp : 只在peak所在的链搜索motif,不进行反义链搜索。HOMER默认会在两条链上都进行搜索。

  • -rna : 搜索RNA上的motif(如使用CLIP-seq数据进行分析的时候)

  • -bg : 指定自定义的背景序列。HOMER默认随机选取基因组序列作为背景。

  • -h: findMotifsGenome.pl脚本默认使用二项分布进行motif富集分析,这在背景序列多于目标序列时是十分有用的。但是有时我们使用 -bg 参数自定义背景序列时,其数目可能会小于目标序列,此时使用 -h 参数选择超几何分布会更加适合。

  • -len : motif的长度设置,默认8,10,12,越大越消耗计算资源。

  • -N: 用于motif分析时所需的序列数目,通常当我们设置 -len过大,内存不够时,会选择减小 -N参数或者 -size 参数。

三、原理

前面我们说到,HOMER软件可以进行多种类型的motif分析,如 promoter motif analysis ,基因组位置motif分析(ChIP-seq分析中的motif分析),利用自定义的fasta文件进行motif分析,RNA序列的motif分析(分析CLIP-seq数据中的RNA binding elements)。 但是,不管我们如何调用HOMER,对于发现motif的基本步骤都是一致的,下面主要介绍一下各步骤的原理:

原理理解即可,HOMER已经将各个步骤封装起来,不用一步一步分别进行。

预处理

1. 序列提取

  • 若给出的基因组位置信息,则提取出来的是对应的基因组序列;

  • 若给出的是基因accession number,给需要选择适当的promoter区域

2. 背景选择

如果没有自定义背景序列信息(即: -bg 参数),HOMER将自动为用户选取。如果使用的是基因组位置,将从整个基因组序列抓取GC含量一致的序列作为背景序列。若进行的是promoter分析,则将所有的promoter作为背景。

3. GC含量矫正

将目标序列和背景序列对GC含量进行分bin(5%区间),背景序列通过调节权重得到和目标序列同样的GC含量分布。这可以使得HOMER在分析来自CpG岛时的序列时,不会简单的找到那些GC富集的motif。下图是一个ChIP-seq分析中GC含量分bin的例子


4. 自动标准化

除了GC含量之外,目标序列还会存在其他的情况影响分析结果。如外显子上的密码子偏好等生物学现象,由A富集序列优先排列引起的实验偏差。当这些偏差显著时,将会影响HOMER的motif分析。HOMER通过给背景序列分配适当的权重,减少寡核苷酸序列频率(如AA)在背景与目标的差异。

重头预测motif

默认情况下HOMER调用homer2版本,若想使用旧版本,在命令行中添加 -homer1 参数即可

5. 解析输入序列

根据设置的motif长度,将输入序列解析为寡核苷酸,并生成一个寡核苷酸频度表。该表只记录出现的寡核苷酸和其出现的次数,可以增加motif搜索时的效率,但是也会破坏寡核苷酸与其原始序列的关系。

6. 寡核苷酸自动均一化

HOMER除了在第四步中对全局序列进行了自动均一化矫正,还在寡核苷酸频度表中进行了矫正。该方法的主要思想是将那些较短的寡核苷酸与较长的寡核苷酸相平衡,但是该方法由于存在许多与motif同长的寡核苷酸,需要分配数量众多的权重。对于那些极端的序列偏好,该方法不会移除它们。

7. 全局搜索

在构建好寡核苷酸频度表后,HOMER进行全局motif搜索。其基本原理就是若某个motif富集,则其存在的寡核苷酸也同样富集。为了增加灵敏度,HOMER允许搜索motif时的错配,同时跳过多个错配的寡核苷酸节省计算资源。

8. 矩阵优化

9. Mask and Repeat

当最优oligo被优化成motif后,motif 对应的序列从要分析的数据中移除,接下来再分析最优的…..直到设定个数的motif被发现。

搜索已知的motif(homer2版本)

10. 加载moitf数据库

HOMER从以前的数据中加载一个已知的motif列表从而搜索已知的motif还可以通过在命令行中指定的motif( -mknown)或通过编辑homer包中的文件( data/knownTFs/known.motifs )来添加motifs。

11. 扫描所有的motif

HOMER通过扫描所有的序列,确定每个motif的富集度,并计算其显著性。

motif输出的结果分析

12. Motif文件

HOMER的输出有许多后缀名为 motif 的文件,其内容如下:

每个motif信息是一块,均以>开头。>所在行的信息以tab分隔。 motif首行信息解释:

  • > + 一致性序列:如图上的>ASTTCCTCTT

  • Motif名称:如图上的1-ASTTCCTCTT

  • 比值检测概率的log值:8.059752

  • P-value得log值:-23791.535714

  • 占位符:如上图得0,不具有任何信息

  • 逗号分隔得富集信息,如:T:17311.0(44.36%),B:2181.5(5.80%),P:1e-10317

    • T表示带有该motif的目标序列在总的目标序列(target)中得百分比

    • B表示带有该motif的背景序列在总的背景序列(background)中得百分比

    • P表示最终富集的p-value

  • 逗号分隔的motif统计信息,如:Tpos:100.7,Tstd:32.6,Bpos:100.1,Bstd:64.6,StrandBias:0.0,Multiplicity:1.13

    1. Tpos:motif在目标序列中的平均位置

    2. Tstd:motif在目标序列中位置的标准差

    3. Bpos:motif在背景序列中的平均位置

    4. Bstd:motif在背景序列中位置的标准差

    5. StrandBias: 链偏好性,在正义链上的motif数与反义链motif数的比值的log

    6. Multiplicity:具有一个或多个结合位点的序列中每个序列的平均出现次数

13 denovo预测的motif

首先对产生的motif进行去冗余,然后将motif转换为向量值从而计算皮尔逊相关系数。HOMER产生的motif结果使用网页展示,该文件名称为 homerResults.html ,而目录homerResults/ 存放着所有该网页依赖的文件和图片。生成的motif首先去冗余。以下是一个简单的输出示例:

shi

输出结果按照p-value排序,最后一列是一个链接到motif文件的超链接,可以从这个文件中找到包含此motif的其他序列。在Best Match/Details 列中,HOMER将会展示与denovo motif最匹配的已知motif(不可全信,最匹配不一定是最优)。该列还包含有一个名为 More information 的超链接,点击后会出现以下页面:

该页面包含motif的一些基本信息,如链接到motfi文件的超链接,且可以生成motif logo的pdf格式。

接着是与已知motif的match。这部分主要看看denovo motif和已知的motif的相似性,需要检查一下是否合理,如下:


如果点击 Simiar motifs 超链接,将会展示在denovo motif搜寻过程中与该motif相似的其他motifs(这些motif的enrichment value相对较低)。通常我们会点击查看一下是否有一些motifs被不正确的分类,比如几个motif具有几个相同的碱基。具体界面如下:

14. 已知motif的输出

该输出也是以网页形式展现,页面结果根据enrichement进行排序,具体如下:


©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 157,198评论 4 359
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,663评论 1 290
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 106,985评论 0 237
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,673评论 0 202
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,994评论 3 285
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,399评论 1 211
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,717评论 2 310
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,407评论 0 194
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,112评论 1 239
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,371评论 2 241
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,891评论 1 256
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,255评论 2 250
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,881评论 3 233
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,010评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,764评论 0 192
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,412评论 2 269
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,299评论 2 260

推荐阅读更多精彩内容