做基因组学相关的课题估计迟早会在某个分析结果里碰到一组看起来像密码的字母比如 AAGCGT、CACGTC或者带方括号的写法如 [CA]G[TC]A。这就是今天要聊的 motif。它不是什么复杂概念但很多人第一次接触时都会懵这串短序列是干嘛的为什么一堆分析工具都在找它找到之后又能说明什么这篇文章我会把它拆开讲清楚motif 到底指什么在基因调控中扮演什么角色怎么表示、怎么找、怎么验证最后分享一些我在真实项目里踩过的坑。适合刚入门生物信息的同学也适合湿实验背景、想看懂 ChIP-seq 或 ATAC-seq 结果的朋友。保证不讲虚的尽量用大白话和能直接抄的实操步骤。1. 先弄明白motif 到底是个什么东西1.1 用一次 ChIP-seq 分析来理解 motif 的定位设想一个最常见的场景你在研究某个转录因子 XX想知道它在一个特定细胞系里调控哪些基因。你做了 ChIP-seqcall peaks 之后拿到几千个显著的结合区域。这时候你会遇到一个问题——这几千个 peak 分布在基因组各处看起来各有各的序列完全找不到规律。于是你拿这些 peak 的序列去做 motif discovery跑完之后第一个显著结果是一个 8bp 的短序列比如 CACGTG。你立刻发现这个序列反复出现在大量 peak 的中心位置而不是随机分布。这个 CACGTG 就是一个 motif它是转录因子 XX 结合 DNA 时的序列“偏好”也就是结合位点的共同特征。实际结合的时候并不是每个位点都严格等于 CACGTG。有些位点可能是 CACGTG有些是 CATGTG有些是 CACCTG。但整体上它们都带有这个核心特征。motif 就是把这些结合位点放在一起抽象出来的一个共同模式。说得再直白一点转录因子是靠识别 DNA 序列来找到结合位置的motif 就是它识别的那一小段“口令”。1.2 motif 为什么值得单独研究基因调控的“钥匙孔”基因表达调控的一个核心环节是转录因子找到正确的 DNA 片段并贴上去。转录因子在细胞核里搜索 DNA 时不是靠人眼去“看”而是通过蛋白质表面的氨基酸残基与 DNA 碱基之间的氢键、静电作用、形状互补来识别。这种识别依赖的就是序列特征也就是 motif。如果把转录因子比作门禁卡那一堆峰序列里的 motif 就是它要插的钥匙孔。一个基因能不能被某个转录因子调控首先要看它的调控区域里有没有对应的 motif如果 motif 因为突变被破坏转录因子的结合就会受影响基因表达随之改变。这也是为什么大量疾病相关 SNP 落在非编码区的转录因子结合位点里——它们中很多就是直接改变了一个 motif。所以研究 motif 不只是做生信分析时的一个步骤它是连接 DNA 序列与基因调控之间的桥梁。你找到一个新的 motif就等于发现了一个可能的新调控位点你把一个已知 motif 对应到某个疾病位点就等于把突变与调控变化联系起来。1.3 不只是转录因子多种调控信号都叫 motif我在实际教学中发现很多人以为 motif 只存在于转录因子的 DNA 结合位点里。其实不是。广义的 motif 可以出现在任何“蛋白识别核酸序列”的场合转录因子结合 DNA最常见的 TF binding motifRNA 结合蛋白结合 RNA比如剪接因子识别外显子附近的剪接增强子序列这类也有 RBP motif组蛋白修饰酶识别特定组蛋白标记与 DNA 序列无直接关系但对应的“reader”蛋白存在底物偏好有些生信分析也会把它类比成 motifCRISPR 系统中的 PAM 序列Cas9 蛋白结合 DNA 时要求目标区域旁边有特定的 PAM motif比如 NGG复制起点、拓扑异构酶切割位点等都对应特定的序列特征。不过日常口语里如果没有特殊说明大家说“motif 分析”通常默认是转录因子结合基序。你可以把它理解成基因组学里的一类通用概念一段短小的、有生物学意义的、可被识别的序列模式。2. motif 的表示方法从“一串字母”到“一张打分矩阵”2.1 共有序列和正则表达式直观但信息量有限最朴素的表示方法是共有序列consensus sequence就是取每一个位置上出现频率最高的碱基拼起来。比如上面 CACGTG 就是共有序列。它非常直观看一眼就知道大概长什么样但缺点也很明显完全丢失了每个位置上碱基分布的细节。假设一个 motif 有 10 个位置有些位置非常保守比如 95% 的序列都是 G有些位置很松G 只占 40%C 占 30%A 占 20%T 占 10%。共有序列只看最高频那个这两种情况都写成 G但实际的识别强度完全不同。IUPAC 简并码可以表达一定程度的模糊性。比如 R 表示 A 或 GY 表示 C 或 TM 表示 A 或 CS 表示 G 或 C。正则表达式如 A[CGT]{2}A 可以描述更灵活的模式但本质上仍是离散的、不能定量打分。2.2 位置权重矩阵 PWM当前最主流的方式现在数据库和工具里最常用的表示方式是位置权重矩阵Position Weight MatrixPWM也叫位置特异性打分矩阵PSSM。它是一张 4×L 的矩阵L 是 motif 的长度每一行对应 A、C、G、T每一列对应一个位置每个格子是一个得分或概率值。PWM 的好处是把“这个位置 A 出现的概率是 0.7G 是 0.2C/T 分别是 0.07/0.03”这种信息完整保存下来。判断一段序列是否可能是某个转录因子的结合位点时就把该序列的每个碱基按位置找到对应列的分数加起来得到一个总得分。得分越高越可能是真正的结合位点。这样的设计非常合理因为蛋白质和 DNA 的结合本质上就是一个热力学过程允许一定程度的错配错配代价因位置而异。PWM 以量化方式模拟了这个过程比单纯看“序列完全匹配还是不完全匹配”要精准得多。2.3 手算一个 PWM从四条短序列推导打分矩阵这部分我用一个最小例子来演示保证你看完能理解数据库里的矩阵是怎么来的。假设我们克隆了 4 条被某转录因子结合的 6bp 序列seq1: A C G T C A seq2: A C G T C T seq3: T C G T C A seq4: A C G T T A第一步统计每个位置的碱基计数位置1A 3T 1位置2C 4位置3G 4位置4T 4位置5C 3T 1位置6A 3T 1第二步把计数转为频率。比如位置1A 的频率是 0.75T 是 0.25位置2C 是 1.0。第三步构建 PWM。标准做法是取 log2(频率 / 背景频率)背景频率一般取基因组或启动子区域的碱基频率。这里简化假设背景 A/C/G/T 都是 0.25。位置1 的 A 得分就是 log2(0.75/0.25)1.58T 是 log2(0.25/0.25)0。位置2 的 C 是 log2(1.0/0.25)2。最终得到的 PWM 大约这样碱基pos1pos2pos3pos4pos5pos6A1.58-inf-inf-inf-inf1.58C-inf2.00-inf-inf1.58-infG-inf-inf2.00-inf-inf-infT0-inf-inf2.0000-inf 表示该位置从未出现某个碱基实际构建时通常不会用负无穷会加一个小伪计数避免除零比如把 0 的次数替换成 0.01 或 0.05。这个细节后面讲参数时还会提到。注意这里为什么要除以背景频率如果一个位置 A 的概率是 0.75但基因组背景里 A 的概率本身就是 0.75那就说明没有任何偏好A 的出现完全是背景水平只有相对于背景明显升高时才能说明该位置对 A 有选择性。3. motif 发现实操把自己手上的序列变成可解读的结果3.1 输入数据准备这一步不花力气后面全是坑做 motif 分析的第一步永远不是直接调工具而是准备干净且合理的输入序列。如果你从 ChIP-seq 数据出发流程通常是先 call peaks比如用 MACS2拿到 narrowPeak 文件。然后从参考基因组里把这些 peak 区域对应的序列导出来。这里我有一个很关键的建议不要拿整个 peak 长度去跑 motif。ChIP-seq 的峰长度在几百到上千 bp真实转录因子结合位点往往集中在峰中心附近的一个小小窗口里。拿全峰去做等于把大量无关的基因组背景序列混进来会严重稀释 motif 信号。我一般取 peak summit 为中心上下各扩 100bp也就是一个 200bp 的窗口如果测序深度较低扩到 150bp 或 200bp 也够用。提取序列用 bedtools 就可以awk {print $1\t$2-100\t$3100\t$4} peaks.narrowPeak peaks_center200.bed bedtools getfasta -fi hg38.fa -bed peaks_center200.bed -fo peaks_center200.fa如果你是做启动子区域分析就取 TSS 上游 -2000 到下游 200 的区间。注意方向性启动子序列的 fasta 提取建议带上链方向或者为了找 motif 方便统一用参考基因组正链提取后续扫描时也要把两条链都考虑进去。另外还有一个很容易被忽视的步骤过滤。建议先利用基因组注释文件把已知的高重复区域过滤掉或者至少用重复序列掩蔽。基因组里有大量转座子来源的重复序列它们富含各种短模式会让富集结果出现大量与转录因子无关的信号。这一步不做后面会很被动。3.2 de novo discoveryMEME 与 HOMER 两种思路输入序列准备好了就可以做从头发现de novo motif discovery。最常用的两个工具是 MEME 和 HOMER它们的算法思路不太一样。MEME基于期望最大化EM算法核心思想是把每一条输入序列看成一段“背景序列 可能存在的 motif 实例”的混合体然后迭代地优化一个 PWM让它在全部序列中的总似然最大化。使用时需要单独指定几条参数meme peaks_center200.fa -dna -nmotifs 5 -minw 6 -maxw 15 -mod zoops -oc meme_out其中-mod zoops表示假设每条序列最多出现一个 motif 实例-nmotifs 5表示找 5 个候选 motif-minw和-maxw限制 motif 宽度。MEME 适合序列数量不多、但你要精细建模的场景输出里有 motif 的 E-value、位点分布图、motif logo信息量很足。HOMER走的是差异富集路线。它会把你的输入序列和一个背景序列集做比较统计所有 k-mer比如 6~12bp 的短单词在目标集里相对背景集的富集倍数然后从显著富集的 k-mer 出发拼接、扩展成完整 motif。HOMER 我最常用的命令是findMotifsGenome.pl peaks_center200.bed hg38 homer_out -size 200 -len 8,10,12这个命令直接在 bed 文件上操作省去单独提取序列的步骤内部会自动界定背景。HOMER 的一个大优势是它会顺带给出“已知 motif 富集”结果直接把你的输入序列和 JASPAR 等数据库做匹配很多情况下刚跑完就能看到一个熟悉的名字。选哪个我的经验是如果是标准的 ChIP-seq 转录因子分析优先试试 HOMER速度快、结果直观、附带已知注释如果输入的序列是深度较深的结合区间或者你想把一个未知的蛋白结合特征从头刻画出来MEME 往往能给出更精细的 motif。更稳妥的方案是两个都跑交叉印证。3.3 与已知数据库比对我的 motif 到底可能是谁de novo 得到的是一个无名 PWM还需要回答一个问题这个 motif 是不是某个已知转录因子的结合位点比对工具有很多最常用的是 MEME 套件里的 TOMTOM。它把一个 PWN 与数据库里所有已知 motif 做两两比对计算相似性并给出 q-valuetomtom meme_out/meme.html jaspar.meme -o tomtom_out -thresh 0.05JASPAR 是植物和脊椎动物转录因子 motif 的主要公共数据库可以直接下载 meme 格式的矩阵文件。除此之外还有 CIS-BP真核转录因子结合位点数据库、TRANSFAC商业库等。HOMER 在运行时会自动做类似的事情它的输出目录里就有 KnownMotif 相关文件。看比对结果时要注意一点某些不同转录因子家族的 motif 可能长得非常像。比如 KLF 和 SP1 家族都偏好 GC 富集序列很多时候 TOMTOM 会同时匹配上一堆名字相似的因子。这时候你需要结合转录因子的表达量、蛋白结构域、已有文献来判断到底哪个是“真正的”结合蛋白。单一靠 sequence 相似性下结论是不够的。3.4 在基因组上回头扫描FIMO 的典型用法很多时候你已经有了一个 motif是从公共数据库下载的或是自己 de novo 找到的接下来想知道它在某个基因组区域里具体落在哪些位置。这时候用 FIMO 做全区域扫描fimo --text --thresh 1e-4 motif.meme query_sequences.fa fimo_hits.tsvFIMO 会逐条扫描输入序列的正义链和反义链对每个位置计算 PWM 打分然后转成 p-value。得到的结果是一个表格每一行是一个命中包含序列名、起始位置、链方向、打分和 p-value。你可以用这个结果去 overlap 启动子列表、去统计某个基因附近有没有结合位点或者画在一张基因组浏览器图上。这里提醒一个比例PWM 阈值设得严命中数少但精确设得松命中剧增但假阳性也跟着涨。一般做候选位点筛选我会用 1e-4 到 1e-5 之间如果要做全基因组范围的粗扫可能需要结合其他特征比如开放染色质信号、保守性进一步过滤。4. 这几种情况我全踩过帮你提前避开4.1 背景模型出错结果全是噪音我最早做 HOMER 的时候直接用了默认参数结果跑出来最显著的“motif”是一长串 poly-A 和 poly-G。当时我还一愣心想这算什么转录因子后来排查发现问题出在目标序列和背景序列的 GC 含量不一致。基因组不同区域的 GC 含量差异很大启动子区普遍高 GC着丝粒附近低 GC线粒体更夸张。如果你的输入序列是高 GC 的启动子而默认背景是随机基因组区域那么“高 GC”这件事本身就会被当成富集信号。poly-G、poly-C 这类短串会被识别为高度显著但它们只是序列组成偏差的结果而不是真正的结合位点。解决办法是让背景序列与目标序列在长度、GC 含量、染色体分布上尽量匹配。HOMER 有专门的选项可以去匹配 GC 分布MEME 里你也可以提供自定义的 control 序列集。一个很实用的习惯是每次跑完 de novo先瞄一眼输出的 motif logo如果出现一堆单调重复的碱基串第一反应应该是背景模型有问题而不是怀疑数据有异常。4.2 重复序列和低复杂度序列没有处理干净第二个高频坑是输入序列里混入了大量转座子或卫星 DNA。例如 Alu 元件、LINE 元件它们在人类基因组里占据了接近一半的位置。如果你的 peaks 里有一部分落在这些重复区域里这些区域的序列模式高度一致会迅速抢占 motif 富集的排行榜。处理方式是在提取序列前先用基因组上的重复序列注释把相关区域刨掉或者使用 RepeatMasker 对序列做硬掩蔽/软掩蔽。硬掩蔽会把重复碱基直接替换成 N软掩蔽则是改成小写字母很多工具默认会忽略小写字母区域。我自己的流程里只要分析目的是“找结合位点”基本都会做一次 mask。有一个例外要说明如果你研究的转录因子本身参与转座子调控比如某些 KRAB-ZNF 蛋白结合转座子区域的 motif那重复区域恰恰是你要找的。这时候不要盲目删要先想清楚生物学问题。工具是死的问题是活的。4.3 长度、motif 数、显著性阈值该怎么选de novo 工具往往会给你一堆可选参数很多人直接跑默认就完事。实际项目中我建议针对自己的数据特点手动设置。宽度上经典的转录因子 motif 一般是 6~12bp短到 5bp 的也有长到 20bp 以上的也有比如 CTCF 的核心结合区域有较长的接触面。如果你设的窗口太短比如 4bp所有结果都会高度相似因为它们本质上只是随机的短序列组合如果设得太长比如 25bp在常见转录因子分析里找到的往往是两个相邻基序拼在一起形成的组合模式解读起来反而复杂。一般我先用 6,8,10,12 这几种长度跑一遍再挑结果中信号最强的那一个做下游分析。motif 数量上常见的 ChIP-seq 数据一个转录因子通常只会有一个主导 motif。如果设-nmotifs 20后面那些只是逐渐衰减的噪音或者是你这个转录因子的共因子cofactor结合位点不需要太当真。我一般取 5~10 个候选然后按显著性、保守性、已知数据库匹配度做人肉筛选。显著性方面MEME 输出里有 motif E-valueHOMER 输出里有 p-value 和 q-value。有一点必须说清楚q-value 是多重检验校正后的 p-value比原始 p-value 更有参考意义。不要看到一个 p-value10^-8 就觉得稳了如果做了几十万次比较这个值可能经不起校正。另外显著性只是概率判断不能代表生物学意义——一个 p-value 极低的 motif 也可能只是重复序列的假象一定要结合后面的验证步骤。4.4 拿到一堆 motif 后如何筛选出“真”的候选这是很多人最困惑的一步工具输出了一堆候选到底信哪个我总结了一套自己的筛选顺序按性价比排序。第一看位置分布真正的转录因子 motif 会显著富集在 peak 的中心区域。如果一个 motif 平均分布在 200bp 窗口的各个位置它大概率是背景噪音或共因子。第二看数据库匹配用 TOMTOM 算一算它是否与某个已知转录因子显著匹配。如果完全匹配不上任何已知 motif有两种可能一是这是新的、未被注释的因子二是这是噪音。后者概率大得多。第三看对应因子表达你研究的系统里这个转录因子有没有 RNA 或蛋白水平的表达如果它连表达都没有很难相信藏在数据里的这个 motif 有功能。第四做 motif 实例检查手动跑到基因组浏览器里看几个命中位点确认它们确实落在可解释的调控区域比如启动子、增强子而且没有落在奇奇怪怪的重复序列里。这套流程不是严格的统计学检验但它能帮你快速把资源集中在最有可能真实存在的候选上。一个 motif 发现项目里后面那些验证工作的时间成本非常高筛选这步值得多花半小时。5. 确认重要性计算验证与湿实验验证5.1 位置偏好、保守性与随机对照一个计算候选经过了 de novo 和数据库匹配还不能直接写进论文。这件事的本质是你找到了一个在统计上富集的序列模式但还没证明它是不是有功能的调控信号。计算层面最容易做的验证是看位置分布。刚才提到ChIP-seq 数据显示真实结合位点往往富集在 peak 中心。可以用 CentriMo 这样的工具做位置富集分析它会检查 motif 在输入序列中的分布是否在中心有显著峰值。如果你发现 motif 距离 peak 中心的距离呈均匀分布那它在 ChIP-seq 背景下是不可信的。第二个有效手段是保守性分析。功能重要的调控序列通常跨物种保守不会在进化中随机漂变。你可以把 motif 中的每个碱基位点与 phyloP 或 PhastCons 分数做关联看保守位点是否在 motif 核心区域富集。一个只在人类基因组里出现、在黑猩猩里就漂变掉的 motif 位点很难说有多少功能价值。第三个是随机对照把输入序列随机打乱或者从基因组中随机抽取与目标序列长度、GC 含量相同的对照序列重复一遍 motif discovery。如果同样的 motif 在你的随机对照里也能富集出来那就说明它只是反映了序列组成的局部特征不具备特异性。5.2 湿实验端怎么验证一个 motif 的功能如果计算验证都通过了且有条件和意愿做实验验证常用的湿实验手段有三个一是 EMSA凝胶迁移实验。合成一段包含该 motif 的双链 DNA 探针与目标转录因子蛋白一起孵育后跑非变性凝胶。如果蛋白能结合这段 DNA复合物会变重迁移变慢出现一个 shift 条带。再把 motif 中的关键碱基突变掉作为竞争探针如果突变后不再竞争就说明结合确实依赖这个 motif 的核心序列。二是荧光素酶报告基因实验。把包含 motif 的调控序列克隆到荧光素酶报告载体里转染细胞后检测发光强度。再突变 motif做同样的转染。如果突变后荧光活性显著下降说明这段 motif 对基因表达有正向贡献。三是 ChIP 实验的延伸如果你已经有了目标转录因子的 ChIP 抗体可以设计包含 motif 和突变 motif 的等位基因在体内检测结合信号差异。这是最贴近真实调控环境的验证。这些实验不需要全部做根据课题需求选一两个就足够。关键是先想清楚问题你是要证明“这个因子能结合这样的序列”还是“这个序列在某个基因调控中起作用”。两者对应的实验策略完全不同。我自己的一个习惯是无论计算还是湿实验所有结论都要落到一条具体的逻辑链上序列有 motif → 转录因子结合 → 调控区域有活性 → 目标基因表达变化。任何一步缺了都要明确说出来而不是含糊带过。做 motif 分析这几年我最大的体会是它不像机器学习模型那样靠一个准确率就能交差它更像考古——从一堆看似混乱的序列里辨认出规律然后小心地论证这个规律是否真的有生物学意义。下一次你再看到报告里的“CACGTG”或者一张漂亮的 motif logo希望你能想起它背后那个完整的逻辑链条而不只是一串 DNA 字母。如果你的课题里刚好要做转录因子分析建议第一步先把自己的输入序列质量搞好这直接决定后面所有结果的可信度。
