motif富集分析（关于mRNA上的m6A富集分析文献解读）

本文目录

• 关于mRNA上的m6A富集分析文献解读
• ATAC-seq专题---生信分析流程
• ChIP-Seq数据挖掘系列-6: 怎么选择HOMMER结果中的motif
• 文献精读 | 基因组6mA修饰相关文章解读及部分分析重现
• ChIP-Seq数据挖掘系列-5.1: ngs.plot 可视化ChIP-Seq 数据
• ChIP-Seq数据挖掘系列-2: Motif 分析(2) - HOMER Motif 分析基本步骤
• meme suite —— Motif分析百宝箱（二）
• 外泌体多组学19-外泌体miRNA分选机制motif研究
• 一篇2区的基因家族分析文章带你飞
• ChIP-Seq数据挖掘系列-5.2: ngs.plot 画图工具ngs.plot.r 和 replot.r 参数详解

关于mRNA上的m6A富集分析文献解读

我今天给大家带来一篇2012年cell杂志发表的 关于 m6A在mRNA上富集的综合分析 文献解读和文献中研究方法的介绍。

Comprehensive Analysis of mRNA Methylation Reveals Enrichment in 3’ UTRs and Near Stop Codons

作者首先通过dot blot的方法验证m6A抗体是否能与RNA上的m6A产生免疫反应，先分别将m6A修饰的寡聚核苷酸和非修饰的寡聚核苷酸点在尼龙膜上，然后用m6A抗体进行免疫反应，再将m6A修饰的寡聚核苷酸点在膜上，用于免疫反应的m6A抗体分别用不同浓度的修饰RNA和非修饰的RNA竞争结合后与膜上的m6A修饰的RNA进行免疫反应，说明该m6A抗体只和m6A修饰的RNA产生免疫反应。

图注：A是不同浓度的M6A修饰和未修饰的RNA与m6A抗体的免疫反应，B是用M6A修饰的RNA和未修饰的RNA与膜上的M6A修饰的RNA竞争性结合m6A抗体，表明m6A抗体与m6A修饰的RNA结合而与未修饰的RNA不结合。

随后作者通过检测不同组织中m6A说明神经组织的RNA甲基化比其他组织高，并且在整个大脑发育过程中m6A含量丰富。作者用寡聚核苷酸（dT）进行RNApull-down钓取细胞Mrna后用m6A抗体进行免疫检测说明mRNA的甲基化，然后通过寡聚核苷酸（dT）杂交获取的mRNA用RNaseH（一种降解polyA尾的RNA水解酶）处理后再进行免疫检测说明Mrna的m6a修饰不是在poly(A)尾上。

图注：A是小鼠不同组织RNA的m6A检测，B是小鼠生长不同阶段的大脑组织RNA的m6A检测，C是用寡聚核苷酸Dt序列钓取细胞的mRNA然后检测m6A，D是钓取的Mrna用RnaseH去除Mrna的poly(A)尾后检测m6A。

为了检测RNA上的m6A，作者通过MeRIP先将细胞总RNA片段化为大约100bp左右，用m6A抗体钓取细胞含m6A的RNA片段，然后进行高通量测序分析m6A的位置，表明大部分的m6A在mRNA的3’端。随后作者用其中的Ldlr mRNA的反向互补序列作为探针钓取细胞的Ldlr mRNA进行m6A的免疫反应验证。

图注：A为MeRIP-Seq的数据分析中mRNA甲基化的信息，B为Ldlr mRNA的RNA-RNApull-down后进行m6A免疫检测验证Ldlr mRNA甲基化。

MeRIP-Seq数据分析表明m6A集中在在motif GAC上，并且在U-rich的motif上几乎不含m6A。

先将m6A修饰的RNA固定在尼龙膜上，然后用m6A抗体进行免疫反应，分别用m6A修饰的RNA和未修饰的RNA作为m6A抗体的竞争结合探针，在加入m6A修饰的RNA时，随着加入的RNA量增加，抗体与膜上的RNA反应降低；而加入非修饰的RNA时，抗体与膜上的RNA反应不受影响

MeRIP-Seq实验方法：先将提取的细胞总RNA进行片段化处理，然后与m6A抗体一起孵育，再用磁珠吸附m6A抗体与m6A-RNA的复合体，洗脱的RNA片段进行高通量测序。

作者用Ldlr mRNA的互补序列探针进行pull-down实验钓取细胞的Ldlr Mrna，然后用m6A抗体进行WB检测pull-down产物，表明Ldlr Mrna上含有m6A富集。CTL探针的pull-down产物作为对照。

ATAC-seq专题---生信分析流程

ATAC-seq信息分析流程主要分为以下几个部分：数据质控、序列比对、峰检测、motif分析、峰注释、富集分析，下面将对各部分内容进行展开讲解。 下机数据经过过滤去除接头含量过高或低质量的reads，得到clean reads用于后续分析。常见的trim软件有Trimmomatic、Skewer、fastp等。fastp是一款比较新的软件，使用时可以用--adapter_sequence/--adapter_sequence_r2参数传入接头序列，也可以不填这两个参数，软件会自动识别接头并进行剪切。如： fastp \ --in1 A1_1.fq.gz \ # read1原始fq文件 --out1 A1_clean_1.fq.gz \ # read1过滤后输出的fq文件 --in2 A1_2.fq.gz  \ # read2原始fq文件 --out2 A1_clean_2.fq.gz \ # read2过滤后输出的fq文件 --cut_tail  \ #从3’端向5’端滑窗，如果窗口内碱基的平均质量值小于设定阈值，则剪切 --cut_tail_window_size=1 \ #窗口大小 --cut_tail_mean_quality=30 \ #cut_tail参数对应的平均质量阈值 --average_qual=30 \ #如果一条read的碱基平均质量值小于该值即会被舍弃 --length_required=20  \ #经过剪切后的reads长度如果小于该值会被舍弃 ***隐藏网址*** 如： fastqc A1_1.fq.gz fastqc A1_2.fq.gz ***隐藏网址*** 经过trim得到的reads可以使用BWA、bowtie2等软件进行比对。首先需要确定参考基因组fa文件，对fa文件建立索引。不同的软件有各自建立索引的命令，BWA软件可以参考如下方式建立索引： bwa index genome.fa 建立好索引后即可开始比对，ATAC-seq推荐使用mem算法，输出文件经samtools排序输出bam： bwa mem genome.fa  A1_clean_1.fq.gz A1_clean_2.fq.gz | samtools sort -O bam -T A1 》 A1.bam 值得注意的是，在实验过程中质体并不能完全去除，因此会有部分reads比对到质体序列上，需要去除比对到质体上的序列，去除质体序列可以通过samtools提取，具体方法如下：首先将不含质体的染色体名称写到一个chrlist文件中，一条染色体的名称写成一行，然后执行如下命令即可得到去除质体的bam samtools view -b A1.bam $chrlist 》 A1.del_MT_PT.bam 用于后续分析的reads需要时唯一比对且去重复的，bwa比对结果可以通过MAPQ值来提取唯一比对reads，可以用picard、sambamba等软件去除dup，最终得到唯一比对且去重复的bam文件。 比对后得到的bam文件可以转化为bigWig（bw）格式，通过可视化软件进行展示。deeptools软件可以实现bw格式转化和可视化展示。首先需要在linux环境中安装deeptools软件，可以用以下命令实现bam向bw格式的转换： bamCoverage -b A1.bam -o A1.bw 此外，可以使用deeptools软件展示reads在特定区域的分布，如： computeMatrix reference-point   \ # reference-pioint表示计算一个参照点附近的reads分布，与之相对的是scale-regions，计算一个区域附近的reads分布 --referencePoint TSS   \#以输入的bed文件的起始位置作为参照点 -S  A1.bw \ #可以是一个或多个bw文件 -R  gene.bed \ #基因组位置文件 -b 3000   \ #计算边界为参考点上游3000bp -a 3000   \ #计算边界为参考点下游3000bp，与-b合起来就是绘制参考点上下游3000bp以内的reads分布 -o  A1.matrix.mat.gz \ #输出作图数据名称 #图形绘制 plotHeatmap \ -m  new_A1.matrix.mat.gz \ #上一步生成的作图数据 -out A1.pdf \ # 输出图片名称 绘图结果展示： MACS2能够检测DNA片断的富集区域，是ATAC-seq数据call peak的主流软件。峰检出的原理如下：首先将所有的reads都向3’方向延伸插入片段长度，然后将基因组进行滑窗，计算该窗口的dynamic λ，λ的计算公式为：λlocal = λBG（λBG是指背景区域上的reads数目），然后利用泊松分布模型的公式计算该窗口的显著性P值，最后对每一个窗口的显著性P值进行FDR校正。默认校正后的P值（即qvalue）小于或者等于0.05的区域为peak区域。需要现在linux环境中安装macs2软件，然后执行以下命令： macs2 callpeak \ -t A1.uni.dedup.bam \ #bam文件 -n A1 \ # 输出文件前缀名 --shift -100 \ #extsize的一半乘以-1 --extsize 200 \ #一般是核小体大小 --call-summits #检测峰顶信息 注：以上参数参考文献（Jie Wang，et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications） ATAC分析得到的peak是染色质上的开放区域，这些染色质开放区域常常预示着转录因子的结合，因此对peak区域进行motif分析很有意义。常见的motif分析软件有homer和MEME。以homer软件为例，首先在linux环境中安装homer，然后用以下命令进行motif分析： findMotifsGenome.pl \ A1_peaks.bed \ #用于进行motif分析的bed文件 genome.fa  \ #参考基因组fa文件 A1  \ #输出文件前缀 -size  given \ #使用给定的bed区域位置进行分析，如果填-size -100,50则是用给定bed中间位置的上游100bp到下游50bp的区域进行分析 ***隐藏网址*** 根据motif与已知转录因子的富集情况可以绘制气泡图，从而可以看到样本与已知转录因子的富集显著性。 差异peak代表着比较组合染色质开放性有差异的位点，ChIP-seq和ATAC-seq都可以用DiffBind进行差异分析。DiffBind通过可以通过bam文件和peak的bed文件计算出peak区域标准化的readcount，可以选择edgeR、DESeq2等模型进行差异分析。 在科研分析中我们往往需要将peak区域与基因联系起来，也就是通过对peak进行注释找到peak相关基因。常见的peak注释软件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker为例，需要在R中安装ChIPseeker包和GenomicFeatures包，然后就可以进行分析了。 library(ChIPseeker) library(GenomicFeatures) txdb《- makeTxDbFromGFF(‘gene.gtf’)#生成txdb对象，如果研究物种没有已知的TxDb,可以用GenomicFeatures中的函数生成 peakfile 《-readPeakFile(‘A1_peaks.narrowPeak’)#导入需要注释的peak文件 peakAnno 《- annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb) # 用peak文件和txdb进行peak注释，这里可以通过tssRegion定义TSS区域的区间 对于peak注释的结果，也可以进行可视化展示，如： p 《- plotAnnoPie(peakAnno) 通过注释得到的peak相关基因可以使用goseq、topGO等R包进行GO富集分析，用kobas进行kegg富集分析，也可以使用DAVID在线工具来完成富集分析。可以通过挑选感兴趣的GO term或pathway进一步筛选候选基因。

ChIP-Seq数据挖掘系列-6: 怎么选择HOMMER结果中的motif

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

文献精读 | 基因组6mA修饰相关文章解读及部分分析重现

Zhu S , Beaulaurier J , Deikus G , et al. Mapping and characterizing N6-methyladenine in eukaryotic genomes using single molecule real-time sequencing. Genome Research, 2018, 28(7):gr.231068.117.

Xiao C L , Zhu S , He M , et al. N6-Methyladenine DNA Modification in the Human Genome. Molecular cell, 2018, 71(2).

简单来说，这篇文章从“华夏1号”的PacBio三代测序数据出发， 前6个results 介绍了PacBio测序寻找6mA修饰的精准度，并对基于二代测序的6mA-IP-seq，基于质谱的LC-MS/MS，以及基于qPCR的6mA-IP-qPCR和6mA-RE-qPCR等方法进行验证和补充说明，结论是基于二代测序的6mA-IP-seq能够得到的6mA修饰位点在SMRT的数据中都有，但是位点要显著少于SMRT测序数据能找到的6mA位点。 Result-7和Result-8 分别介绍了基于猜想寻找人类系统中6mA修饰的甲基化和去甲基化酶，分别是N6MT1和ALKBH1。这里主要是进行了一系列细胞水平的过表达或者干扰/敲除实验，以及体外甲基化实验。首次鉴定了人类基因组6mA修饰的甲基化酶——N6AMT1。最后两个 结果9和10 主要是粗略研究了一下一小部分gastric cancer 和liver cancer病人中的6mA修饰情况以及N6AMT1/ALKBH1的表达情况，得出 癌组织比癌旁组织的6mA修饰少 ， 值得注意的一点是，这篇文章当中，对癌和癌旁组织中的6mA进行检测用的方法是免疫组织化学（IHC），而不是测序的方法 。并且发现6mA修饰水平和预后有关。通过对肿瘤细胞系做N6AMT1和ALKBH1过表达、沉默实验发现6mA修饰水平的高低会影响肿瘤细胞的增殖、迁移、侵袭等表型。

个人觉得这篇文章值得参考的地方是利用PacBio的测序数据去提取修饰信息的Bioinformatics分析部分。

Cell. 2018 Nov 15;175(5):1228-1243.e20. doi: 10.1016/j.cell.2018.10.006. Epub 2018 Nov 1 N6-methyladenine DNA Modification in Glioblastoma.

我个人感觉 这篇文章思路比较迷 （也许应该说比较丰富，所以一遍读不明白），发表于第二篇文章之后，出现很多有意思的，与第二篇文献的结论完全相反的结果。在看了Author Contribution之后我的感觉更加强烈——即这篇文章应该是两个竞争课题组后来整合成果共同发表为一篇文章的套路。（只是猜测，不一定对哈）由于相对来说这篇文章比肖传乐的文章更复杂，所以在附上完整思路整理图的同时，我对单个result进行拆分梳理如下。前三个结果主要谈论的是N6-mA修饰；后三个结果主要谈论的是ALKBH1这个基因。

首先是 背景 ，值得一说的是为什么作者选择了神经胶质母细胞瘤（Glioblastoma）来研究DNA的6mA修饰。Glioblastoma是一种发病率很高，并且很具侵略性的原发性脑肿瘤， 这个癌种被证明有广泛的表观遗传修饰失调 ，例如DNA 5mC的高甲基化、染色质重塑酶基因如EZH2和BMI1的表观修饰变化。然后作者们就选择了用来自病人的胶质母细胞瘤细胞对DNA 6mA修饰及其功能进行研究。

这部分主要用dot blot, MS和免疫细胞荧光(ICC)三种技术对Glioblastoma病人衍生的干细胞GSCs和两个primary tissues进行研究(primary tissues只做了dot blot)，然后发现 在GSCc和primary tissues中的6mA修饰都比human astrocytes (人类形形胶质细胞，作为normal control)中要高。这个结论与第二篇肖传乐等人的文章中得出的‘normal组织中的6mA修饰水平更高’刚好相反 。

这个部分主要是对两个GSCs细胞系和一个primary组织进行了N6-mA DIP-seq测序。分析测序结果发现(1) 每个样本能检测到7,282-17,263个N6-mA peak, 这一结果与2018-molecular cell上发表的肖传乐等人的文章中用的6mA-IP-seq在HuaXia1基因组中检测到的(21,129)差不多？ (2) 每条染色体上都检测到了N6mA，其中7号染色体、19号染色体和21号染色体上最多； 而关于分布的染色体的特征，肖传乐等人的文章似乎表示6mA在常染色体上的分布是比较均一的？ (PS: 在文章中进行了Genome Background比较得出的这个结论，由于我是生信菜鸟，我不太懂，希望有专业人士可以指导一下~)(3) GO富集分析发现m6A修饰的基因主要富集在神经发生和神经发育通路；(4) motif分析发现6mA的常见修饰寡核苷酸是GGAAT，肖传乐等人的文章中得到的motif是 。

作者提出 “N6mA是一种抑制性的表观修饰” 。怎么说呢，因为在Result-2中鉴定到的富含6mA修饰的DNA motif- GGAAT-是人类大部分微卫星重复序列的结构，而人类微卫星重复序列主要在异染色质区域。所以Result-3就开始探索N6mA与异染色质之间的关系。首先对Result-2部分的三个样本进行了H3K9me3，H3K27me3和H3K4me3的ChIP-seq，发现80%的N6-mA的峰能和异染色质marks，即H3K9me3, H3K27me3重合。之后又对387这个GSC做了WGBS (全基因组甲基化测序)，发现有轻微差异。（小声：连显著性都没标是不是因为没有显著差异，23333）

这个结果部分首先对ALKBH1这个基因进行体外和体内实验，均证明能影响6mA的量； 最重要的是，这篇文章也做了N6AMT1这个基因，通过CRISPR构建N6AMT1敲除的细胞系的6mA修饰水平没有发生变化。体外的实验也没有变化。这里有几个地方我个人认为比较不明确：首先，做CRISPR-Cas9敲除的GSCs细胞系未告知明确是哪个；然后做dot blot的DNA用量也没有标明。因此不能和肖传乐等人的文章进行更精确的比较。 _________________________________________________________________________ 还有不解的是：后面他们先是做了CRISPR敲除ALKBH1基因的细胞系，然后去做了RNAseq，得到表达变化了的差异基因；然后他们用shRNA干扰ALKBH1表达的细胞做N6-mA DIP-seq去找ALKBH1被knockdown之后N6mA修饰发生改变的位点，然后比较ALKBH1被敲除后表达下调的基因的N6-mA的修饰的变化情况。为什么不直接用ALKBH1-KO的细胞系做N6-mA DIP-seq呢？？？关键是他文写的还是’after ALKBH1 deleption’？？？感觉略混乱啊。。 ________________________________________________________________________ 继续‘震惊’的是，他们又做了ATAC-seq，发现ATAC-seq的峰刚好和N6-mA DIP-seq的峰呈负相关。。。这思路真(kan)清(bu)奇(dong)。。。难道前面不是已经证明了N6-mA喜欢在异染色质区域吗？异染质区域难道不是开放性弱的区域吗？？那么不是理论上就是在ATAC-seq不能拉到的区域吗？？ ________________________________________________________________________ 然后，为了搞明白ALKBH1到底是怎么去掉N6-mA修饰的，他们做了ALKBH1 pull-down实验。发现ALKBH1喜欢跑到带有N6-mA修饰的寡聚核苷酸那里去。。为了解释这种现象，他们又做了ALKBH1的ChIP-seq。。。然后发现ALKBH1 ChIP-seq的峰和N6mA的峰重叠。。然后得到的结论是：ALKBH1是一个转录因子，它跑到N6-mA富集的地方去并且移除N6-mA修饰的抑制作用。。思路绕吗？？可能CNS级文章就是这样吧。。

第五个结果比较简单，但是到了第五个结果这里，ALKBH1, Histone modification (eg, H3K9me3), N6-mA modification和Gene expression这4者之间的关系其实相对来说似乎就明朗了。 (1) ALKBH1升高会减少N6-mA的修饰； (2) N6-mA修饰减少会去掉基因表达的抑制性修饰，所以理论上来说正常组织中带有较多N6-mA修饰的基因在ALKBH1表达升高之后它们的表达也会升高，即ALKBH1的水平与基因表达呈正相关； (3) N6-mA修饰与H3K9me3的修饰peaks高度重叠 ； 所以在第五个结果里，通过对ALKBH1 knockdown的细胞系做H3K9me3 ChIP-seq，弥补上了ALKBH1与H3K9me3的关系。本文的研究结果是： 检测到ALKBH1 knockdown的细胞系基因组中有更多H3K9me3 peaks，因此认为ALKBH1可以调控组蛋白修饰的形成。

这篇文章内容很丰富，思路很丰满，太累了。后面主要是设计了相对完整的对照实验做RNAseq等实验证明N6-mA与缺氧调控通路的关系，以及去做了个ALKBH1敲除后的肿瘤细胞表型实验等。

Wu T P , Wang T , Seetin M G , et al. DNA methylation on N6-adenine in mammalian embryonic stem cells. Nature, 2016.

Yao B , Cheng Y , Wang Z , et al. DNA N6-methyladenine is dynamically regulated in the mouse brain following environmental stress. Nature Communications, 2017, 8(1):1122.

这篇文章的第二个结果，通过6mA-IP-seq测了3对stress/control鼠脑PFC区的6mA变化，得到主要变化——不论gain 6mA还是loss-of 6mA都是发生在intergenic基因间区。当关注基因内的6mA变化时，发生主要是发生在introl区。 而2018-Moleculer Cell-肖传乐等人的文章的结果是6mA在基因内intragenic的修饰区域主要是在exon区。 在第四个结果里，通过对RNAseq的数据进行分析，并且结合6mA-IP-seq数据，把所有基因分为以下四类：① gain-of-6mA/upregulated; ② gain-of-6mA/downregulated; ③ loss-of-6mA/upregulated; ④ loss-of-6mA/downregulated。GO分析之后发现只有第③类loss-of-6mA/upregulated的基因能够富集到与神经发育及功能相关的通路上。

Liu J , Zhu Y , Luo G Z , et al. Abundant DNA 6mA methylation during early embryogenesis of zebrafish and pig. Nature Communications, 2016, 7:13052.

这篇文章比较短，只有7页pdf，接下来按照文章给出的result来理一下整个项目的思路。 首先第一个结果部分 (1) DNA 6mA modification during early embryogenesis of zebrafish ，作者主要利用UHPLC-QQQ-MS/MS超高分辨率-三重四级杆-串联质谱，以斑马鱼为模型，非常细致地研究了斑马鱼的精细胞、卵子、早期胚胎发育的各个细胞时期（受精卵、2细胞、4细胞、8细胞、16细胞、32细胞、64细胞、128细胞、256细胞、512细胞时期）以及成年斑马鱼的不同组织器官（包括脑、眼睛、心脏、卵巢、精囊、肌肉和肠）进行了6mA丰度研究，平均每一个实验组需要约10,000个细胞提取的基因组DNA。结果表明精细胞中6mA/A的比例约在0.003%，卵母细胞中约为0.015%；在32-64细胞时期6mA/A的比例达到最高，约为0.1%，之后慢慢下降，直至回归原始基因组水平。 同时成年斑马鱼的不同组织器官中6mA/A的丰度也是符合一般水平，无明显变化。

ChIP-Seq数据挖掘系列-5.1: ngs.plot 可视化ChIP-Seq 数据

基因组基础包只是包含 genebody, CGI, exon ；对于 hg19 和 mm9 ， ngs.plot 准备了额外的Enhancers 和 DHSs注释；

ngs.plot.r 参数设置，参考文章 ： ngs.plot 画图工具ngs.plot.r 和 replot.r 参数详解

ngs.plot.r 参数设置，参考文章 ：

H3K36me3配置文件 ： config.hesc.k36.txt

H3K27me3配置文件 ： config.hesc.k27.txt

做热图时，ngs.plot 有不同的基因或区域排序算法，也有层次聚类和k-means两种聚类方法。为例避免高测序深度样本带来的偏差; 在聚类时,将值转化为排名。

配置文件 : config.k4k27.inp.txt

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

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

meme suite —— Motif分析百宝箱（二）

Motif Discovery中还包括MEME-ChIP，可对ChIP-seq或CLIP-seq数据的DNA序列进行一系列的Motif分析。该方法整合了： 1、MEME&STREME功能，用以预测Motif（de novo Motif）； 2、CentriMo用以寻找输入序列中间区段内已知的Motif，适合ChIP-Seq数据检出峰所在序列上富集的已知 Motif； 3、通过工具Tomtom功能比较已知的motif进行相似性分析，并对重要的motif进行分组； 4、Spamo与CentriMo功能相似，也是对Motif进行富集； 5、FIMO功能旨在输出Motif在基因的位置信息。 MEME-ChIP具体参数可以使用默认值，或按需求更改。示意图如下： 我们以网站Sample Output中MEME-ChIP example结果为例进行说明。

外泌体多组学19-外泌体miRNA分选机制motif研究

外泌体和其他小的胞外囊泡(sev)提供了一种独特的细胞间通信模式，在这种模式中，从一个细胞产生和释放的microrna(miRNAs)被远处的细胞吸收，在那里它们可以引起基因表达的变化。然而，mirna被分选到外泌体/sev或保留在细胞中的机制仍在很大程度上未知。 在这里，我们证明了miRNA具有分选序列，决定其在sEV（EXOmotifs）或细胞保留(CELLmotifs)中的分泌，并且不同的细胞类型，包括白色和棕色脂肪细胞、内皮、肝脏和肌肉，优先使用特异性分选序列，从而定义该细胞类型的sEVmiRNA谱。 为了研究miRNA细胞保留和外泌体/sEV释放的具体特征，我们从5个小鼠细胞系中分离出了外泌体/sEV，这些细胞系代表了参与代谢调节的重要组织： 不同类型的细胞释放不同数量的sev，3T3-L1脂肪细胞的产生和释放率最高，C2C12肌管的释放率最低；这与观察到的脂肪组织是体内循环外泌体/sEVmirna的主要贡献者相一致 在所有病例中，sev的直径为50-200nm。 富集经典的外泌体标记物ALIX、TSG101、CD9和CD63，而细胞标记物GM130和CANX基本没有 这些sev的RNA含量与释放的囊泡总数相平行 使用基于实时荧光定量PCR(qPCR)的阵列评估分泌的sev及其来源的细胞的miRNA组成(图1a)。与预期的一样，细胞miRNA的主成分分析(PCA)显示，尽管棕色和白色的脂肪细胞和内皮细胞聚集在一起，而肝细胞和肌细胞更为不同，但每种细胞类型都有不同的miRNA谱。 当PCA同时包括细胞和sEVmiRNA时，sEVmiRNA彼此之间以及与其来源细胞中的miRNA有很大的不同(图1b)，突出了miRNA分泌的特殊性 然后，我们比较了每种细胞或sEV类型中每个miRNA与其他类型的水平。在细胞体中评估的664个miRNAs中，有210个miRNAs（32%）在一种细胞类型中的表达明显高于其他四种细胞类型 同样，与其他细胞类型的囊泡相比，大约三分之一的sEVmiRNAs（218/660）在一种细胞类型的囊泡中富集，可用于预测来自不同器官的混合sev中的组织来源，如在血液中。 一些sev特异性的miRNA反映了miRNA的细胞类型特异性，例如，miR-133a/b，它在C2C12肌管中特异性表达。 然而，对于每种细胞类型，73-92%的被认为是针对特定细胞类型的sEVmiRNAs也在其他细胞类型的细胞体中以类似或甚至更高的水平表达，提示了一种细胞类型特异性的分类机制。同样地，在两种不同细胞类型的sev中同样丰富的miRNA在两种分泌细胞类型中可能有相当不同的表达水平。因此，确定sEVmiRNA的组织起源并不像知道miRNA在哪个组织中高表达那么简单（上图e）。 通过比较细胞体中每个miRNA与sev的相对水平可以揭示miRNA分泌和保留的独特模式。 与细胞体相比，一些mirna在所有五种细胞类型的sev中均富集，而其他的则在只有一种或两种细胞类型的sev中表现出选择性富集，还有一些发现在sev中很少或根本没有发现，尽管它们存在于细胞体中。 将sEV中miRNA的相对丰度除以其在细胞体中的相对丰度(sEV富集)，表明sEV和细胞中miRNA的分类存在大量差异，一些miRNAs在sEV中明显富集，而另一些在细胞体中明显富集（保留）. 为了确定miRNA分选的潜在机制，我们分析了显示sEV或细胞体富集的miRNA的miRNA序列和结构。5p和3pmiRNAs没有普遍富集（补充表6）。然而，与细胞中保留的序列相比，sEV富集水平更高的miRNA序列具有更高的G+C含量和更低的吉布斯自由能(ΔG)(扩展数据图2c，d). 识别潜在的差异分选miRNAs序列，我们分析在sEV或细胞体中富集的miRNAs的序列，将它们与未显示出优先sEV分选或细胞保留的mirna序列进行比较。

一篇2区的基因家族分析文章带你飞

BZR（BRASSINAZOLE-RESISTANT）家族基因是编码参与油菜素内酯信号转导的植物特异性转录因子，在植物生长中起着至关重要的作用。 今天我就给大家带来一篇甜菜中BZR基因家族分析的文章。文章于2019年5月9日在线发表在BMC Plant Biology（影响因子3.93，中科院分区二区）。 具体分析内容如下： 一、甜菜中 BvBZR 基因的鉴定 通过鉴定，共鉴定出6个BvBZR基因： Bv5_cuzi 、 Bv_epwr 、 Bv1_fxre 、 Bv6_nyuw 、 Bv1_qnjn 、 Bv_yfzt 。 二、Motif分析和系统发育分析 为了阐明BZR家族的进化关系，作者基于来自甜菜、拟南芥、水稻和大白菜的41个BZR家族成员的氨基酸序列构建了系统发育树，并进行了motif分析。 三、 BvBZR 基因染色体分布和基因结构分析 作者将鉴定的6个 BvBZR 基因定位到了甜菜基因组的5条染色体上，并对基因结构进行了分析。 四、 BvBZR 基因的顺势作用原件分析 五、不同甜菜品种根茎的生长特征规律统计 作者统计了包括主根的生长曲线（根重）、主根的生长速度、主根的含糖量以及主根含糖量的增加速率4个指标。 六、与甜菜生长特征相关的基因表达模式和相关性分析 七、 BvBZR 基因在E型和Z型根、茎、叶组织中表达模式分析 八、 BvBZR 基因对植物激素响应的基因表达模式分析 为了研究 BvBZR 基因的表达水平是否受外源植物激素的调节，作者对甜菜根喷洒了IAA、ABA、MeJA、GA3共4种植物激素，并检测了 BvBZR 基因的表达水平。 九、 BvBZR 基因的亚细胞定位 作者首先使用Wolf PSORT软件对 BvBZR 基因进行亚细胞定位预测，并采用实验手段对预测结果进行了验证。 总结 到此为止，这篇基因家族类文章的所有分析就完成了，在内容上还是比较常规的，只是在实验方面补充了因的亚细胞定位实验和一些生长指标，并没有复杂的实验操作和分析内容，值得大多数研究者借鉴！ ***隐藏网址*** 更多生物信息课程： 1. 文章越来越难发？是你没发现新思路，基因家族分析发2-4分文章简单快速，学习链接： 基因家族分析实操课程 、 基因家族文献思路解读 2. 转录组数据理解不深入？图表看不懂？点击链接学习深入解读数据结果文件，学习链接： 转录组（有参）结果解读 ； 转录组（无参）结果解读 3. 转录组数据深入挖掘技能-WGCNA，提升你的文章档次，学习链接： WGCNA-加权基因共表达网络分析 4. 转录组数据怎么挖掘？学习链接： 转录组标准分析后的数据挖掘 、 转录组文献解读 5. 微生物16S/ITS/18S分析原理及结果解读 、 OTU网络图绘制 、 cytoscape与网络图绘制课程 6. 生物信息入门到精通必修基础课，学习链接： linux系统使用 、 perl入门到精通 、 perl语言高级 、 R语言画图 7. 医学相关数据挖掘课程，不用做实验也能发文章，学习链接： TCGA-差异基因分析 、 GEO芯片数据挖掘 、 GSEA富集分析课程 、 TCGA临床数据生存分析 、 TCGA-转录因子分析 、 TCGA-ceRNA调控网络分析 8.其他课程链接： 二代测序转录组数据自主分析 、 NCBI数据上传 、 二代测序数据解读 。

ChIP-Seq数据挖掘系列-5.2: ngs.plot 画图工具ngs.plot.r 和 replot.r 参数详解

ProgramArguments101

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 参数详解 ChIP-Seq数据挖掘系列-6: 怎么选择HOMMER结果中的motif