专栏首页生信技能树找个motif嘛,简单

找个motif嘛,简单

上午在朋友圈看到有公众号点名:

好不容易有人搞个公众号活动,就支持一下,毕竟生信技能,都在技能树啊~

原文见:画泡泡图,发CNS

虽然今天整体都在开会ing,但还是抽空上午跑了个程序,晚上出结果,再花30分钟写这个教程,完美!

我在生信菜鸟团发布的自学CHIP-seq分析第八讲就提到过如何寻找motif,motif是比较有特征的短序列,会多次出现的,一般认为它的生物学意义重大,做完CHIP-seq分析之后,一般都会寻找motif 。查找有两种:

  • 一是de novo的,要求的输入文件的fasta序列,一般是根据peak的区域的坐标提取好序列
  • 二是依赖于数据库的搜寻匹配,很多课题组会将现有的ChIP-seq数据进行整合,提供更全面,更准确的motif数据库。

motif的英文定义如下:

motif: recurring pattern. eg, sequence motif, structure motif or network motif

DNA sequence motif: short, recurring patterns in DNA that are presumed to have a biological function.

从上边的定义可以看出,其实motif这个单词就是形容一种反复出现的模式,而序列motif往往是DNA上的反复出现的模式,并被假设拥有生物学功能。而且,经常是一些具有序列特异性的蛋白的结合位点(如,转录因子)或者是涉及到重要生物过程的(如,RNA 起始,RNA 终止, RNA 剪切等等)。

一篇文献列出了2014年以前的近乎所有知名的A survey of motif finding Web tools for detecting binding site motifs in ChIP-Seq data 链接见:https://biologydirect.biomedcentral.com/articles/10.1186/1745-6150-9-4

当然,扯这么多大道理,没用。

直接上实战!


说大事专用分割线~

input是找到的peaks文件,bed格式

上游分析这里略过,我的GitHub里面给了全套流程代码:https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq

反正得到了bed格式的peaks文件,而且一般下载公共数据在GEO里面作者都会顺带上传他们做好的peaks文件,比如:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66581 就可以下载ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE66nnn/GSE66581/suppl/GSE66581_bed_files.tar.gz 文件。

因为是下载作者的数据,所以必须注意;参考基因组版本是:Genome_build: mm9

最常用的是 MEME工具套件 :

http://meme-suite.org/ 输入文件是fasta序列,需要对peaks进行转换,根据bed的基因坐标从基因组里面提取对应的序列咯: http://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

它里面集成了4个寻找motif 的工具,每个工具都是一篇文章,里面有详细描述具体原理,但是整个网页给人的感觉是too busy,让初学者无从下手,这次讲解就略过它咯。

homer软件来寻找motif

这个软件安装当初特别麻烦: https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step8-Homer-findMotif.sh

但是现在有了conda,一句话搞定:conda install -c bioconda homer ,值得提醒的是,如果是在中国大陆,那么需要设置一下:

# https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda config --show

下载软件及数据

conda install -y -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda homer
perl ~/miniconda2/share/homer-4.9.1-5/configureHomer.pl -install mm9 
ls -lh ~/miniconda2/share/homer-4.9.1-5/data/
## 作者使用的是mm9找到的peaks文件
## 数据库下载取决于网速咯
## 下载成功后会多出 ~/miniconda2/share/homer-4.9.1-5/data/genomes/mm9/ 文件夹, 共 4.8G
## 这个文件夹取决于你把homer这个软件安装到了什么地方。


## 或者用下面代码安装:

cd ~/biosoft
mkdir homer &&  cd homer
wget http://homer.salk.edu/homer/configureHomer.pl 
perl configureHomer.pl -install
perl configureHomer.pl -install hg19

## 软件安装之后就下载peaks文件,肯定是bed格式的
mkdir ~/motif
cd ~/motif
wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE66nnn/GSE66581/suppl/GSE66581_bed_files.tar.gz 
tar zxvf GSE66581_bed_files.tar.gz 

下载得到的bed格式的peaks文件如下:

motif/
|-- [3.4M]  1k_es_peaks.bed
|-- [2.5M]  200_es_peaks.bed
|-- [780K]  2cell_early_peaks.bed
|-- [415K]  2cell_early_trx_inhibitor_peaks.bed
|-- [2.3M]  2cell_peaks.bed
|-- [1.7M]  4cell_dba_peaks.bed
|-- [2.2M]  4cell_peaks.bed
|-- [2.1M]  50k_es_peaks.bed
|-- [1.9M]  8cell_peaks.bed
|-- [5.9M]  GSE66581_bed_files.tar.gz
`-- [1.4M]  icm_peaks.bed

homer软件找motif整合了两个方法,包括依赖于数据库的查询,和de novo的推断。

运行homer软件

但是使用起来很简单:http://homer.ucsd.edu/homer/ngs/peakMotifs.html

for id in *_peaks.bed;
do
echo $id
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp 
findMotifsGenome.pl homer_peaks.tmp mm9 ${id%%.*}_motifDir -len 8,10,12
annotatePeaks.pl    homer_peaks.tmp mm9 1>${id%%.*}.peakAnn.xls 2>${id%%.*}.annLog.txt 
done 

不仅仅找了motif,还顺便把peaks注释了一下。

然后为了出图,我还需要解析这个结果,把需要的基因的motif数据抓出来

也是非常简单,只需要对每个peaks注释结果的 knownResults.txt grep即可。

接下来,就看Y叔表现咯。

整个流程我是我的云服务器操作的,其中homer软件+mm9数据库约5G,下载的peaks及注释后的文件不到200M,软件安装及数据库下载,数据下载,注释,找motif的代码都给了,如果你想实践一下,可以考虑购买我的云服务器, 99 元一个月,不仅仅是找motif,更多详情请访问: 古有杨志卖刀,今有jimmy售器

其它资源

还有一些R包可以,直接从BED文件里面记录的基因坐标来找motif,有的需要输入fasta序列,就需要自己根据bed的基因坐标从基因组里面提取对应的序列咯:

rGADEM (motif discovery): http://bioconductor.org/packages/devel/bioc/html/rGADEM.html

MotIV (motif validation): http://bioconductor.org/packages/devel/bioc/html/MotIV.html

http://lgsun.grc.nia.nih.gov/CisFinder/

http://bioinfo.cs.technion.ac.il/drim/

本文分享自微信公众号 - 生信技能树(biotrainee)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-03-29

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 优秀学员的学习方法展示

    诚然,不同环境下成长的大家吸收新知识的习惯和能力千差万别,但总有一些人的经验非常值得借鉴!同样的指点我发出去了31份,能坚持一个月的寥寥无几,甚至能坚持5天的也...

    生信技能树
  • 爬虫那么危险,干嘛不直接基因数据库下载文件呢?

    问了具体后,才知道原来是ncbi上的信息,相当于在ncbi上在gene库中查找,然后爬取目标信息。如下:

    生信技能树
  • Nanostring的表达矩阵分析也是大同小异

    临床样品的特色是:通常是FFPE样本,在保存过程中往往造成RNA的断裂,不论是qPCR还是RNA-seq都难以进行精准的定量,这个时候Nanostring 仪器...

    生信技能树
  • JavaScript第十二弹——ES6(上)

    Hello大家好,最近我们也讲了不少JavaScript的知识了,今天再来点实用的吧,不管是在工作中还是面试中,ES6都是我们会遇到的一个东西,ES6呢,全称是...

    萌兔IT
  • 献给前端er的各种小技巧(纯干货)

    1.Firefox 的查看页面源代码功能,可以一眼发现未闭合的标签、未转义的HTML字符,另一种办法,提交页面代码到 http://validator.w3.o...

    疯狂的技术宅
  • Spring之事务源码理解,Spring4.3.12.RELEASE版本

    1、声明式事务,境搭建环。在pom.xml配置文件中新增依赖的jar包,导入相关依赖,数据源、数据驱动、Spring-jdbc模块。如下所示:

    别先生
  • DFS练习-HDU1010

    题目来源:HDU1010 DFS的基本原则已经差不多了,但是一些技巧仍然比较难想,所以还是加强练习,然后总结一下。

    Enterprise_
  • 【TBase开源版测评】并行处理

    分布式HTAP数据库(TencentDB for TBase)是腾讯自主研发的分布式数据库系统,集高扩展性、高SQL兼容度、完整的分布式事务支持、多级容灾能力以...

    随风飞舞
  • 从学美容、学理发,到艺术、古典乐和诗歌,社会人可以这样学文化

    知晓君
  • C++ OpenCV图像分割之GrabCut分割

    在OpenCV中的图像分割中GrabCut分割算法,该算法可以方便的分割出前景图像,操作简单,而且分割的效果很好。在前我们刚用学了OpenCV中的鼠标回调函数,...

    Vaccae

扫码关注云+社区

领取腾讯云代金券