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

在基因组调控元件分析中，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 参数详解

HOMER 是一套用于Motif查找和二代数据分析的工具。hommer结果中一般包含已知motif富集情况，并且也会对用户提供的序列进行重头预测motif。

很多同学在拿到这个结果后都是一脸懵，尽管Homer motif Results是有排序的，排名第一的可能并不是用户的期望，后面的结果也是可以选择的，但是如何评价及选择这些结果？

这段序列比对到了已知motif(YY1)，但是我们可以发现比对上的位置并不是YY1 motif (CAAGATGGC)的中心区域。因此YY1 motif 并不是完整的motif在用户的数据中富集，因此这个结果是不太可信的。

在许多情况下，HOMER 结果有很显著的p值，但是motifs 却不是好的。

因此，在选择motif的时候，用户需要注意以下的原则：

低复杂度的motif序列的核苷酸倾向于都是同一种核苷酸，从而导致GC含量异常。

当目标序列和背景库中序列之间存在系统性偏差时会导致这样的结果。通它们的GC含量非常高。在这种情况下，您可以在motif分析命令中添加参数“-gc”，从而使运算按总GC含量而不是CpG岛含量进行标准化。

其它情况，当分析多样的基因组序列时，这是很难在背景中控制的。例如，将一个启动子序列比对到某些物种随机基因组背景序列，结果会对嘌呤或嘧啶的偏好。HOMER非常敏感，所以如果序列的组成有偏差，HOMER 很可能会发现。新版本中的Autonormalization可以尽量减小这个问题的发生。

motifs 有时候会出现一些序列模式的重复

这种motifs 一般会有数个差不多序列的motifs。除非有充分的理由相信这些可能是真实的，否则背景序列可能有问题。如果你的目标序列在外显子和其他类型的序列上高度富集，就会出现这种情况；并且如果" -gc "参数也不能改善结果，用户就需要考虑自己正在分析序列的类型以及怎么去匹配他们。

这种发生在motif看起来很靠谱，但是在序列中出现的百分比缺失很低的。例如，寡核苷酸和重复序列在用户序列中出现从而导致极高的显著性。统计上使显著的，但是事实上却并不是。一些调节基因的启动子序列会发生这样的事情。原则上，motif 在不到5%的靶序列中存在的话，这个motif不太可信。

排名第一的motif：

后续的motif：

这种情况不一定是坏的，但是也需要考虑；常常发生于ChIP-Seq数据中，免疫沉淀的蛋白高表达以及与大量的结合位点紧密结合。这些motifs 可能结合PU.1，但是亲和性不高。处理这种情况的方法是重复motif 分析的过程，但是丢掉 the top motif（排名靠前的），添加参数 -mask <motif file> 就可以在motif分析过程中忽略这些top motif 。

How to Judge the Quality of the Motifs Found

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 参数详解

2017年发表在Nature Methods杂志上的 SCENIC 算法，利用单细胞RNA-seq数据，同时进行基因调控网络重建和细胞状态鉴定，应用于肿瘤和小鼠大脑单细胞图谱数据，提出并证明了 顺式调控网络分析能够用于指导转录因子和细胞状态的鉴定 。SCENIC通过使用生物学驱动的features自动清除肿瘤样本特异性等批次效应。

我们要深刻了解SCENIC并应用在我们的项目上就要对 基因调控网络(GRN, gene regulatory network) 的背景有一些认识。细胞的转录状态来自潜在的基因调控网络，GRN由数量有限的转录因子（TFs）和辅因子相互调节及调节下游靶基因构成。

单细胞转录组分析的最新进展为高分辨率识别转录状态和状态之间的过渡提供了令人兴奋的机会。例如，在分化过程中针对单细胞RNA-seq进行优化的统计和生物信息学方法带来了新的生物学见解。但是， 稳定细胞状态下的特异性和强健的（robust）GRN是否能够被确定还尚不清楚 。

鉴于在单细胞水平来做这件事情，这可能确实具有挑战性。由于转录爆发和其他来源的基因表达的随机变化，基因表达可能会与TF输入的动力学部分断开（不相关）。目前已经有一些从单细胞RNA-seq数据推断共表达网络的方法发表，但是，这些方法未使用调控序列（regulatory sequence）分析来预测TF与靶基因之间的相互作用。

作者认为，将顺式调控序列与单细胞基因表达连接可以克服缺失（基因表达丰度检测不到）和技术多样性，从而优化细胞状态的发现和表征。鉴于此，作者开发了单细胞调控网络推断和聚类（SCENIC）来绘制GRN，通过评估每个细胞中GRN的活性，来识别稳定的细胞状态。

SCENIC workflow 包含3个主要步骤：

使用RcisTarget是SCENIC不同于大多共表达算法的重要区别 。由于GENIE3模块仅基于共表达，因此结果可能包含许多误报和间接target，为了鉴定推断的直接结合的靶标基因，使用RcisTarget对每个共表达模块进行顺式调控基序（motif）分析。 仅保留具有正确基因上游调节子且显着富集TF motif的模块，并对它们进行修剪以除去缺乏基序支持的间接靶标，这些处理后的模块才称为regulon。

作为SCENIC的一部分，作者开发了AUCell算法来对每个细胞中的每个regulon的活性进行评分。

对于一个给定的regulon，通过比较所有细胞间的AUCell打分值，我们可以识别哪些细胞具有更显著高的regulon活性。通过卡阈值得到的二元活性矩阵使矩阵维数减少（可理解为只有 0|1，on|off），对于下游分析很有用。 例如，基于regulon二元活性矩阵的聚类，可以根据某个调控子网络（regulon）的活性来识别细胞群类型和细胞状态。由于regulon是整体评分的，而不是使用单个基因的表达，因此这种方法对于个别基因的dropouts很有效。

作者对SCENIC进行了 准确性 和 鲁棒性 评估。

作者拿成年小鼠大脑的已知细胞类型的scRNA-seq数据集来测试SCENIC，来评估其准确性。

分析了1,046个初始共表达模块中的151个Regulons，这些Regulons包含对应TF的显着富集的motif（占初始TFs的7％）。每个细胞的Regulon活性评分确定了预期的细胞类型（ 上图d，e ）以及每种细胞类型的潜在主调控因子的列表（例如下图中的小胶质细胞网络），细胞聚类的结果比一些专用的单细胞聚类方法( SC3 )要准确。

作者用以下方式分析了小鼠大脑单细胞RNA-seq图谱数据：

SCENIC鉴定出了仅由少量细胞代表的细胞类型（例如，来自小胶质细胞，星形胶质细胞或中间神经元的2至6个细胞；下图）。此外，预测的TFs与细胞类型的关联也与先前已知的一致，并且此准确性优于标准分析流程。

为了验证鉴定小鼠interneurons的Dlx1/2调控网络，作者分析了人脑的sNuc-Seq(Single nuclei RNA-Seq)数据集。

在人类脑sNuc-Seq数据集上，SCENIC也鉴定出由DLX1/2强烈驱动的interneurons细胞群，该群具有与小鼠相同的motif，并且识别出一组保守的靶标，包括DLX1本身。

接下来，作者将这种跨物种分析扩展到其他细胞类型。 与基于归一化表达的标准聚类（产生强大的物种驱动聚类）不同，SCENIC分析有效地按细胞类型对细胞进行了分组（下图）。 这表明网络活性的评分是可靠的（robust），可以用来克服批处理或测序方法等技术影响。

由于肿瘤特异性突变和复杂的基因组畸变，癌细胞状态的鉴定比正常细胞状态更具挑战性。一般的标准聚类会通过表达矩阵将细胞按其肿瘤起源分组（各个样品聚成一类），但SCENIC的结果揭示了不同的图景。

以下是作者将SCENIC用在 少突胶质细胞瘤 （来自 6个肿瘤的4,043个细胞 ）和 黑素瘤 （来自 14个样本的1,252个细胞 ）的scRNA-seq数据集上的结果。

对于少突胶质细胞瘤（oligodendroglioma），在所以肿瘤细胞中鉴定出三种癌细胞状态（上图c-e），每种状态均由预期的TF驱动，包括：

此外，作者将扩散映射（Diffusion Maps）应用于二元SCENIC矩阵（上上图）重建了从stem-like到oligodendrocyte-like 和 astrocyte-like分支的分化轨迹。值得注意的是，与正常oligodendrocyte分化相比，此路径代表不同的“轨迹”。

在黑素瘤（melanoma）数据上观察到了类似的肿瘤效应校正，其中SCENIC识别了跨肿瘤的细胞群。包括一群与少突胶质细胞瘤中类似的TF驱动的周期细胞（例如，E2F1/2/8 和 MYBL2）

黑色素瘤细胞大致分为两组：

SCENIC在MITFlow state的细胞群下鉴定到了两个新的TFs：

NFATC2是JNK/MAPK途径中的转录阻遏物，参与黑色素瘤去分化和免疫逃逸；

NFIB与毛囊和黑素细胞干细胞的干细胞行为有关，它在小细胞肺癌转移过程中起着重要作用。

为了进一步探索NFATC2和NFIB在MITFlow状态下的潜在作用，作者对25个具有不同肿瘤进展的黑色素瘤标本进行了免疫组织化学分析。

作者发现NFIB和NFATC2在前哨淋巴结（sentinel lymph nodes）中表达最高，这与ZEB1表达共定位，这与ZEB1表达共定位，这表明这些markers的表达与最早的转移事件之间存在关系。

在A375黑色素瘤细胞系，NFATC2和NFIB的表达很高，当使用siRNA敲除NFATC2时，作者发现NFATC2调节子中的基因被显着上调。这与先前证实的NFATC2作为阻遏物是一致的。此外参与细胞粘附、细胞外基质和几个先前发表的代表黑色素瘤浸润状态调控的基因也被上调。这表明NFATC2确实可能在疾病进展中起重要作用。作为黑色素瘤regulons的第二次验证，作者使用ChIP-seq数据确定了MITF和STAT的预测靶标。

随着单细胞数据集大小的增加，作者建议两种补充方法来扩展网络推断：

SCENIC是一种普遍适用的分析scRNA-seq数据的方法，利用TF和顺式调控序列来指导细胞状态的发现。文章的结果表明，GRNs是确定细胞状态的可靠方法，并且scRNA-seq数据非常适合跟踪基因调控过程，在基因调控过程中特定组合的TFs驱动细胞特异性的转录组。

作者开发的pyscenic使用教程请参考： https://pyscenic.readthedocs.io/en/latest/ 或者查看 pyscenic github

最近发表的 单细胞转录组测序绘制小鼠内皮细胞精细图谱 分析研究转录因子在不同ECs中的调控作用，发现了转录因子调控网络呈现了组织特异性。

---