做转录因子和motif分析这几年我最大的感受是网上资料多但能一口气把转录因子是什么—motif是什么—怎么拿到motif—拿到之后怎么用这条线讲透的文章太少。多数人要么背了一堆数据库名要么在软件参数里打转最后拿着一个sequence logo却不知道怎么解释生物学意义。这篇我按自己的实操经验整理尽量不绕弯子把转录因子、结合基序、以及分析中常踩的坑一次性说清楚适合刚入门的生信同学也适合做实验想自己跑一轮motif分析的湿实验研究者。1. 转录因子和motif先搞清楚这两个概念到底什么关系1.1 转录因子是做蛋白-DNA互作的开关转录因子Transcription Factor, TF本质是能够识别特定DNA序列并调控基因转录的蛋白质。一个基因要不要表达、表达量多少很大程度上取决于它的启动子和增强子上结合了哪些转录因子。可以把它想象成一把带识别能力的钥匙DNA上被识别的那段序列就是锁孔钥匙插对了锁孔转录机器才被招募过来基因才开闸放水。这个锁孔就是转录因子结合位点Transcription Factor Binding Site, TFBS而描述这类结合位点共性规律的短序列模式就是我们常说的motif。所以motif不是某个具体的基因组位置而是一种序列偏好的抽象表达——它回答的问题是这个转录因子更愿意结合在什么样的DNA序列上。一个典型的转录因子结合位点长度通常只有6-20个碱基远远短于一个基因的编码区。也正因为短它在基因组里天然就会大量随机出现所以转录因子要发挥作用还必须配合染色质可及性、共调控因子、DNA甲基化等机制避免乱开枪。这个背景知识在做motif富集分析时特别重要后面我会专门讲为什么不能只靠一个motif就下结论。1.2 转录因子结构域与motif的对应关系转录因子在结构上是模块化的。绝大多数转录因子至少包含一个DNA结合域DNA-Binding Domain, DBD和一个调控功能域如激活域或抑制域。其中DNA结合域负责读DNA序列而激活域/抑制域负责喊话给转录机器。我们讨论的motif本质上对应的是DNA结合域识别的序列模式而不是整个蛋白序列的保守基序——这是很多人混淆的点。按DNA结合域的结构特征转录因子可以分成好几个大家族锌指蛋白Zinc Finger最常见的一类C2H2型锌指蛋白是哺乳动物基因组里数量最多的转录因子家族识别的motif通常比较长而且往往相邻多个锌指协同识别。bZIP碱性亮氨酸拉链以二聚体形式结合DNA识别ACGT核心类序列典型如AP-1家族。bHLH碱性螺旋-环-螺旋识别E-box类序列CANNTG在发育和神经分化中极其重要。核受体家族Nuclear Receptor识别激素响应元件通常以反向重复inverted repeat形式出现。Homeodomain家族识别TAAT核心基序在胚胎发育中扮演关键角色。这个分类之所以要了解是因为motif分析到了后期总要做家族归类比如你找到一个E-box富集大概率要和bHLH家族挂钩找到一个ACGT核心往往往bZIP方向看。家族先验知识能帮你快速缩小候选转录因子范围不至于在几百个motif里大海捞针。1.3 动辄数千个转录因子背后的逻辑很多初学者会困惑人类基因组里到底有多少个转录因子为什么不同数据库数出来的数量不一样根据Human TFDB等数据库的统计人类有大约1600-1900个转录因子相关基因占了全部蛋白编码基因的7%-8%左右。但真正常见、研究比较透的也就几百个。剩下很多是组织特异性表达、发育阶段特异性表达或者只在特定刺激下才出现。这也是做motif分析时一个非常重要的思维转变数据库里那个motif能从几十个物种里比对上不等于在你的实验体系里那个转录因子真的在发挥作用。转录因子要在某个细胞里干活它首先得表达其次染色质得打开再次它得有协同因子帮忙。所以看motif结果时一定要回头结合RNA-seq或蛋白表达数据做交叉验证。我在帮朋友看数据时见过太多次富集出的motif是某个发育关键因子但样本是成年肝组织那个因子根本不表达——这种结果就只能当计算线索不能当机制结论。2. motif的表达方式与核心概念从consensus sequence到PWM2.1 千万不要只记consensus sequence最直观的motif表达方式是一致序列consensus sequence比如GATA1的motif可以写成WGATAA其中W代表A或TR代表A或GY代表C或TN代表任意碱基。这种记法直观但丢掉了大量的定量信息。它只告诉你每个位置最常见的碱基却不告诉你这个位置上其他碱基允许到什么程度。举个例子如果某个位置A出现60%、C出现40%A是consensus里的碱基但这个位置其实有接近一半的序列是C。你要是拿一个C的位点去比对用consensus的方式就会漏掉。所以在正式分析里我们几乎不用consensus做计算它是给人看概览的不是给软件做判断的。真正做计算用的是位置权重矩阵Position Weight Matrix, PWM也叫PSSMPosition-Specific Scoring Matrix。PW M的每一列对应motif一个位置每一行对应A、C、G、T四个碱基数值是该位置上出现某碱基的对数几率比。简单说它把每个位置的碱基偏好变成了一个可加减的分数。扫描某条序列时把窗口内每个碱基对应位置的分数加起来得到该窗口与motif的匹配得分。顺带一提从PWM可以计算每个位置的信息量information content单位是bit。信息量越高的位置说明该位置碱基越固定对结合特异性贡献越大。这正是sequence logo图里每个字母高度的由来。2.2 怎么读sequence logo做过motif分析的人都会遇到sequence logo图——就是那种每个位置堆叠A、C、G、T字母的图。很多人只扫一眼哪个字母最大其实logo里的信息密度远超consensus。在读logo时重点看两个维度总高度该位置的保守程度越高表示越保守结合特异性贡献越大。通常是0-2bit2bit对应完全保守只有一种碱基。字母顺序同一位置多种碱基按频率从上到下排列主碱基在最上面。如果某个位置两个碱基高度接近说明这个位置兼容两类碱基这类位置往往是决定转录因子结合特异性的可塑性位置。在做项目汇报时我习惯把logo图、PWM信息和对应转录因子家族一起放上这样审稿人或合作者一眼就能看清楚这个motif的特异性结构比单独贴一段consensus有说服力得多。2.3 motif的方向、位置与二聚化一个经常被忽视的motif特性是方向性问题。许多转录因子以二聚体形式结合DNA这使得它们的结合位点呈现**回文palindrome或直接重复direct repeat**结构。举个例子核受体识别的往往是AGGTCA的某种重复排列两个半位点之间隔着数个碱基AP-1识别的是TGACTCA类回文。方向为什么重要因为在扫描基因组找结合位点时软件的搜索策略里有一个是否考虑反义链的选项。如果你研究的转录因子识别的是非回文motif你必须在正负链上都扫描如果是回文motif正链扫到就等同于负链扫到。HOMER、FIMO这类工具有默认处理但你要理解这个逻辑尤其是做自定义脚本时漏掉反向互补会导致结果直接少一半。另外motif在启动子区的位置也有学问。经典观点认为转录因子结合位点大多落在转录起始位点TSS上游几百碱基内但随着增强子研究深入我们知道很多功能性结合位点位于远端调控区域甚至基因内部。所以做motif扫描时区域选择直接影响结论。拿到一个motif信号先问一句我扫的是启动子、增强子还是全基因组开放区域这三类区域的enrichment解读逻辑差别很大。3. motif从哪来实验方法与计算方法两条路线3.1 实验定义motif的主流方法现在数据库里那些motif源头大致可以分两类一类来自低通量的体外实验比如SELEXSystematic Evolution of Ligands by Exponential Enrichment另一类来自高通量测序实验主要是ChIP-seq。SELEX的原理简述把转录因子蛋白和随机DNA文库混合洗掉不结合的DNA留下结合的那部分PCR扩增后再进入下一轮筛选。多轮富集后测序这些DNA片段比对出富集的序列模式。它的优点是体外定义、干净利落直接反映蛋白本身的结合偏好缺点是脱离了染色质环境可能高估或低估了体内的某些结合事件。ChIP-seq则是体内实验用抗体把和目标转录因子结合的DNA片段拽下来然后测序。把测到的peaks取中心区域序列再做de novo motif分析得到的就是该转录因子在细胞内真实结合位置的序列偏好。它的优点是生理相关性强缺点是受抗体质量、peak calling参数、染色质状态影响很大。实际使用中数据库里同一个转录因子的motif可能同时来自SELEX和ChIP-seq两种来源两者可能长得不完全一样。SELEX出来的motif通常更锐利信息量高因为体外没有其他因子干扰ChIP-seq出来的motif往往包含一些附属的次级结合位点因为pull down下来的片段可能还带着与目标蛋白形成复合物的其他因子的足迹。所以你看JASPAR里同一个因子有时有多个条目就是这个道理。3.2 计算上怎么de novo找motif当手里有一批peaks或一组共调控基因的启动子序列时想从头发现富集的motif有很多工具可选MEME、HOMER、DREME、Weeder、STREME等。它们的算法逻辑不完全一样但核心思路都是看哪些短序列模式在目标序列集里出现的频率显著高于背景模型预期。以MEME为例它用的是EM算法期望最大化假设每条序列里可能含0个或1个motif实例也可以配置其他假设迭代优化PWM的各个参数直到收敛。HOMER则更偏快速扫描统计检验适合对大量peak序列做超几何富集分析。STREME是MEME Suite里较新的工具专为找短motif设计速度比MEME快不少。选工具没有唯一标准。我的一般策略是先用STREME跑一个快速初筛再用MEME精细建模最后用HOMER交叉验证一遍。如果两个工具都能富集出同一个motif这个结果的信心度就高很多。别指望一个工具一步到位motif发现这件事本身就带有统计推断的不确定性多工具互证是最实用主义的方法。3.3 常用motif数据库怎么选做motif比对或注释时数据库选择直接决定你查出来的因子准不准。最常用的几个JASPAR非冗余、开放获取的转录因子motif数据库覆盖脊椎动物、植物、昆虫、真菌等多个分类。质量高格式规范是绝大多数分析的首选。日常用JASPAR 2024的非冗余脊椎动物库基本够用。Cis-BP覆盖面极广包含大量物种的转录因子DNA结合域预测motif很多JASPAR里查不到的冷门物种/冷门因子能在Cis-BP里找到。HOCOMOCO专注于人和小鼠整合了大量ChIP-seq和SELEX数据注释质量较高适合做人和小鼠项目时参考。TRANSFAC老牌商业库内容全但收费个人用户一般用不起。我的建议是主力用JASPAR冷门物种查Cis-BP人和小鼠兼看HOCOMOCO。另外注意版本号别用旧版数据库的motif去注释新数据不同版本同一个因子的motif矩阵可能存在差异审稿人有时会揪这个。4. 实操用MEME Suite和JASPAR完成一轮完整的motif分析4.1 准备输入序列格式和区域选择无论用什么工具第一步都是准备FASTA格式的输入序列。以分析一批ChIP-seq peaks为例你需要把peak的基因组坐标转换成序列。这一步通常用bedtools getfasta完成。假设你有一个peaks.bed参考基因组是hg38bedtools getfasta -fi hg38.fa -bed peaks.bed -fo peaks.fa这里有一个细节值得注意peaks通常有几百碱基宽而motif只有6-20bp如果直接把整个peak塞进去背景噪声会比较大。更常用的做法是把peak先收缩到以peak summit为中心、左右各50-100bp的窗口再做motif分析。很多软件如HOMER会自动做类似处理但如果自己控制流程这一步能显著提升信噪比。窗口大小没有标准但我个人经验是**100bpsummit±50和200bpsummit±100**是比较好的起点。太短会丢掉部分真实结合位点太长会稀释信号。另外如果做的是启动子motif分析比如一组差异表达基因的启动子建议取TSS上游-1000到100这个范围。取序列时注意正负链基因的TSS方向最简单的方法是只取上游-1000到TSS这一个方向然后用软件的-revcomp选项同时搜负链。4.2 de novo找motifMEME参数怎么设有了FASTA文件后就可以用MEME找motif了。命令长这样meme peaks_100bp.fa -dna -mod zoops -nmotifs 5 -minw 6 -maxw 20 -revcomp -oc meme_out几个关键参数的含义-dna输入是DNA序列。-mod zoops每条序列里每个motif出现0次或1次。这是最常用的假设。如果你怀疑某些序列里同一个motif会出现多次可以考虑-anrany number of repetitions如果确信每条序列必含一个用-oops。默认其实是zoops我一般保持默认。-nmotifs 5让程序输出排名前5的motif。第一次跑可以多设几个比如10个跑完后看富集排名再筛选。-minw 6 -maxw 20motif宽度范围。转录因子motif大部分落在这个区间如果研究的是特别短的基序比如4-5bp的某些辅助因子可以把minw调低但短motif假阳性会显著上升。-revcomp同时搜正负链。运行完成后meme_out目录里会生成meme.html、meme.txt等结果文件。HTML文件可以在浏览器里直观查看每个motif的logo、E-value、在每个输入序列中的位置分布。MEME结果的解读关键看两个指标E-value这个motif在随机序列中出现的期望次数。越接近0越好一般小于1e-3是比较可信的信号。Motif的宽度和位点覆盖motif宽度太宽比如15bp但明显是两段拼起来的有可能是两个相邻motif融合了后续建议拆开看。4.3 用FIMO扫描新序列找到具体的结合位点MEME找到的是这条motif长什么样但要回答在我的目标序列里具体是哪些位置有结合位点需要做motif扫描。FIMO是MEME Suite里的扫描工具它拿一个motif或一整个motif集合逐一扫描输入序列的每个位置给出每个位置的匹配得分和显著性。用法示例fimo --oc fimo_out JASPAR2024_CORE_vertebrates_non-redundant_pfms_jaspar.meme target_sequences.fa这里第一个文件是motif数据库可以从JASPAR网站下载MEME格式的PFM文件第二个文件是你想扫描的序列。FIMO的输出是fimo.tsv每一行代表一个扫描到的结合位点列包括motif_id、sequence_name、start、stop、strand、score、p-value、q-value。实际分析时我会做一个过滤步骤awk $8 0.05 fimo_out/fimo.tsv fimo_significant.tsv把q-value小于0.05的位点留下。q-value是多重检验校正后的p-value只盯着p-value看会在扫描了几万个位点时积累一堆假阳性。FIMO扫描还有一个派生功能值得做拿到显著位点在基因组上的分布后可以统计它跟TSS的距离分布、在不同注释区间的富集程度甚至可以结合ATAC-seq数据看这些位点是否落在开放染色质区域。这一步能把纯序列层面的预测结合提升到潜在功能性结合解释力强很多。4.4 用TomTom把de novo motif跟已知库比对MEME找出来的motif是个无名氏得给它验明正身。把de novo的motif文件拿去和JASPAR里的已知motif比对看它最像哪个转录因子的motif这件事TomTom可以完成。用法tomtom -oc tomtom_out meme_out/meme.txt JASPAR2024_CORE_vertebrates_non-redundant_pfms_jaspar.meme输出里最重要的是tomtom.tsv里面每个待查询motif会给出最匹配的数据库条目、匹配得分和q-value。当q-value小于0.05时基本可以认为两者的位置权重矩阵高度相似。但这里有个经验提醒TomTom匹配的是矩阵相似性不是功能等价性。同一个家族里的不同成员比如FOXO1和FOXO3矩阵几乎长得一样TomTom在它们之间往往给出一堆并行的低q-value结果这是正常的不代表结论错误。在报告里建议写成该motif与FOXO家族如FOXO1/FOXO3的已知结合基序高度一致而不是斩钉截铁说结合FOXO1。如果TomTom里没有任何已知库匹配上先不要急着宣布发现了新motif。先做反向检查去JASPAR里按转录因子家族搜一下这个核心序列比如你发现的motif里有个E-box核心CAGCTG看是不是属于语法变换或方向变化导致的未匹配。另外考虑可能目标物种的ortholog数据库里确实没有这个因子的条目这时可以去Cis-BP里碰碰运气。5. 实操中我踩过的坑motif分析的典型问题排查5.1 结果太多或太少怎么排查motif分析最常让人抓狂的两类情况一类是de novo跑出来一大串motif全是低复杂度的AT富集序列另一类是跑出来几乎什么都没有E-value差到没法看。先说结果太多、全是AT富集的情况。这通常不是转录因子的真实偏好而是序列组成偏倚或背景模型不匹配的产物。比如你用GC含量很低的基因组区域做输入MEME默认的背景模型可能会把这些低复杂度区域当作显著motif。解决办法是给MEME提供自定义背景模型或者先屏蔽简单重复序列用RepeatMasker处理再跑。也可以在后续过滤时直接丢掉那些consensus里NNN比例过高、或者信息量主要集中在一两个位置的低复杂度motif。再说结果太少的问题。最常见的原因是输入序列数量太少。MEME的统计检验很依赖样本量你塞个20条序列进去想得到显著motif确实比较难。另一种可能是窗口取得太长信号被稀释。把100bp扩到300bp时信号密度变化不大但背景序列量大幅增加显著性自然下降。还有一种是motif太保守、太短——比如4bp的poly(A)结合位点这种短motif在基因组里到处都是背景模型里本来就有很高背景频率很难被检验为显著。实操里我通常这样应对如果输入序列少于50条先不做de novo直接跳到用已知motif做富集分析如果序列多但结果空先调整窗口大小再考虑换工具比如STREME和MEME对短motif的灵敏度不同。5.2 数据库版本和格式转换的坑JASPAR、Cis-BP、HOCOMOCO提供的motif文件格式不同有的是PFMraw counts矩阵有的是PWM有的是transfac格式。MEME Suite的扫描工具通常需要MEME motif格式直接拿PFM文件喂给FIMO不一定会报错但可能被错误解析。如果你需要转换格式可以用MEME Suite自带的脚本chen2meme JASPAR2024_CORE_vertebrates_non-redundant_pfms_jaspar.txt jaspar.meme这个命令把JASPAR的PFM文件简单序列格式转成MEME格式。另一种常见情况是拿motif位置文件BED格式的TFBS去算motif在启动子区上的分布密度这里需要留意BED坐标是0-based还是1-basedbedtools默认0-based错过了这个细节结果差一个碱基对一些短motif来说一个碱基偏移足够反转结论。还有一点JASPAR的CORE和non-redundant版本区别。non-redundant集合里同一个转录因子家族只保留一个代表性motif适合做全基因组扫描避免多个相似motif重复洗数据CORE里则包含更多变体条目适合做精细比对。选择哪个取决于你的目的不要无脑用CORE。5.3 生物学解释才是最后一步最后必须强调一句motif富集是关联不是因果。整个计算流程——de novo发现、已知库比对、位点扫描、富集分析——得出的结论都是这些序列模式在数据里显著富集至于这个富集有没有生物学功能必须结合其他证据链。我通常要求自己在下结论前至少回答三个问题这个转录因子在我的样本里表达吗如果不表达motif富集有什么替代解释比如它可能是另一个共存因子的影子位点这些结合位点是否落在转录起始位点或增强子附近的可及染色质区域没有开放染色质支持的结合位点很可能是存在但不工作的。是否有多组学证据互相印证比如转录因子表达量、motif富集、目标基因表达变化三者之间有没有一致性的逻辑链。很多人在最后一步翻车就是跳过了表达量验证直接说富集到FOXO motif所以FOXO通路激活。这种claim做生信的自己心里有数审稿人迟早也会问。5.4 别忘了考虑家族内共享motif的问题基因重复在进化中非常普遍同一家族的转录因子往往识别几乎一样的DNA序列。比如FOXO家族的FOXO1/FOXO3/FOXO4它们的DBD区高度保守motif矩阵几乎相同。你在数据分析里永远无法单靠motif区分是哪个FOXO成员在结合。这不是软件的局限而是DNA结合层面的事实同一个家族的不同成员通过蛋白-蛋白互作、翻译后修饰、组织表达模式来实现调控特异性而不是靠DNA识别序列的差异。所以在报告里写结果时我的习惯是写到家族层面为止先报motif显著匹配FOXO家族FOXO1/FOXO3的结合基序再结合表达数据和文献说其中FOXO3在样本中高表达提示其可能为主要调控者。这样写既准确又经得起追问。5.5 几个实用技巧速查跑MEME之前用fasta-shuffle-letters生成对照序列看你的motif是否只在目标序列里富集而对照里没有这是排除背景偏倚的快速方法。大规模扫描全基因组时建议把JASPAR的non-redundant库拆成几个大类如家族分类分别跑FIMO输出结果更便于下游层析。可视化motif时R的ggseqlogo包比WebLogo更灵活可以方便地把多个motif排列对比做图也更符合发表要求。记录分析过程时一定要保留motif的版本标识和来源数据库版本。同一个JASPAR条目在不同版本里矩阵可能微调这部分信息在方法学描述里必写。做生物信息学尤其是转录因子结合这种信号弱、背景料杂的分析心态上要有点做侦探的意思数据给你线索但不能只凭一条线索定罪。把motif这个线索拿到了关键是找到表达、染色质、功能实验这几条旁证一起钉进去结论才立得住。这也是我从一开始坚持motif只是起点不是终点的原因。