【原】基于gff3/gtf文件-批量提取启动子或CDS序列-任何人都可以

写在前面

Emm....最近确实是忙到晕头转向，所以一直没有继续写。今晚似乎还有大半个小时的时间，便尝试写写。这篇推文应是 《-零命令行-生信下游数据分析》的第五贴，主题定为基于gff3/gtf从基因组序列中批量调取序列。
伴随高通量测序的普及，越来多的物种基因组得到测序，为从事响应物种的研究工作者提供了极大的便利，但也带来了不少数据提取和整理的烦恼。
有一句名言说的好，如下：

海量的生物数据，如果无法解读或利用，那么不是生物信息，而是占据硬盘的垃圾。--CJ

而利用这些数据的其中一个常见的热点就是，序列的提取。主要原因是：

1. 有些基因组确实只提供了基因组序列和注释文件

2. 有时候我们需要基因组网站所没有发布的部分，如某些基因的启动子区域

这个时候，就需要自己进行序列提取

正文

这篇推文主要讲三个模式，分别对应三种的需求。当然，任何工具都像一口锅，任何人能介绍的最多是这口锅可能能拿来烧这道菜那道菜，但不代表你不能拿这个锅来煮面条，或者是蒸鱼。
回到主题，本次演示，需要的文件是：

1. 基因组序列，fasta格式

2. 基因组注释信息，gff3或者gtf格式，TBtools会自动识别
   如：

   

提取所有基因的CDS

需要做的事情很简单，

1. 将基因组注释文件拖进去

2. 点击Initialize，此时会弹出gff3文件的预览，关闭即可

3. 选择Feature Tag为CDS

4. 选择Feature ID为 Parent，因为我们最后是使用CDS的parent标签进行CDS串联的，所以选择parent即可。如果是使用gtf文件，这里可能要选择transcript_id

3. 将基因组序列文件拖进去

4. 设置输出文件路径和文件名

5. 点Start

大概十几秒钟后，之所以这么慢，是因为这需要建立一下fasta索引，后续就快了，提取完成，可以直接将cds.fa拖入到TBtools的Fasta stater快速统计序列概况，提取出来35386个cds序列（这肯定是包括了可变剪切）

大致看下序列文件情况

提起所有基因的ATG上游2000bp（所谓的启动子序列）

有时候，我们会将基因的ATG上游2000bp定位基因的启动子序列（当然其实这个是不严谨的，只是没办法）
此时，只需要调整一处即可，设置Up Stream Bases为2000，修改下输出文件名

点击Start之后，几秒钟即可提取结束，因为fasta索引之前已经建立好了
查看下序列情况，并与前面的cds序列进行对比

所以，很多人纠结了很久的基因组已知，如何获取某个基因的启动子序列，就这么简单。。。

只提取基因的CDS上游2000bp

上述我们提取之后，还保留了CDS序列，但是有时候我们并不需要，此时，只需要做一个小的修改即可

这个时候，提取出来的，就只有2000bp的序列。。。

写在最后

其实，CDS只是一个操作，如果基因组注释文件中有UTR，你也可以提取UTR上游2000bp，或者1000bp，或者100bp...
当然，TBtools其实还提供了很多参数，比如你只想要没有内含子的基因的CDS序列

功能有太多，写是写不完了
索性就不写了。
对TBtools感兴趣的朋友，欢迎加入QQ群

或者加我微信，带上微信红包（避免广告），拉你进微信群