大范围查找序列模式所在位置，用TBtools啊

写在写在前面的前面

中午一个师弟说，手上有一个需求，需要确定某个模式在荔枝所有基因的启动子区域内是否存在。事实上，这很明显就是要看某个转录因子是否是调控哪些基因，我一看模式就知道，必定是ERF，因为是典型的GCC-box序列。

一提到这个，我的想法是：

  1. 复杂的操作，那么你最好用MEME，MAST，所以可能需要稍微写个小教本，但是事实上，你也可以用TBtools，我早就打包了Windows, Linux, MacOS的执行器。所以。。。有那么复杂吗？

  2. 简单的操作，事实上大部分人需要的是Perfect Match，当然你也可以自己写正则。我在半年前就开放了TBtools的一个序列模式定位功能（附加在下方）。

如果是我直接在命令行下操作，那么我一两分钟就搞定，如下

我在写这个命令的时候，考虑了：

  1. 基因组存储序列是会存在换行的，往往60个字符一行，模式匹配的时候需要考虑这个

  2. 使用-0777一次性读取整个文件进行检索，因为这个正则几乎不会出现在header上，没必要搞得太复杂。

  3. 基因组序列中可能会有大小写（当然，手上这个文件，必然是小写，因为使用TBtools一次性提取的所有基因的启动子区域，我当然知道只有小写）

最后，我们得到上述的一行perl命令。

问题解决，但是还是存在另外两个问题：

  1. 如果用户不是需要数目，而是需要每个序列，甚至是匹配位置？

  2. 如果基因组很大，比如16G的小麦基因组，那么上述这个命令至少主要32Gb的内存，运行慢且不说...

  3. 无论基因组大小如何，模式匹配本身就比直接indexOf()慢很多，这是基本的编程基础知识

  4. ....

上述这些问题，在TBtools里面都不是问题，即是你用的是1个T的基因组，只要你愿意花一点点时间等待，那么TBtools都能帮你找出所有匹配的位置，而且还很快，比上述那行perl快；而且给出的信息还很全面，包括所有匹配位置；当然还支持正则等等等等

输出结果还很优美

同样是736个....

最后，还是补充一句：

无论你是做实验的，还是做数据分析的，TBtools都值得你推荐给身边的人。因为你是前者，那么TBtools减少你的烦恼，给你补充一些基础数据分析能力；如果你是后者，那么TBtools可以减少其他成员占用你的工作时间，你应该去做一些有深度的分析，而不是帮其他人做一些任何人都能简单搞定的事情，比如做个热图，有意思吗？

当然了，你在高校是如此；你在公司亦是如此，减少你客户的烦恼；降低公司技术运维的工作量。

CJchen，公众号：生信札记7000+人教我写TBtools - 他值得任何人推荐给身边的朋友

附加半年前TBtools的推文如下：

写在前面

TBtools开发至今三年有余，让开发工作的持续开展，事实上主要来源于用户朋友的信任与支持。总的来说，我对多数用户存在某种意义上的感谢，毕竟使用一个未正式发表的工具，并在漂亮的发表工作引用，是对我个人能力在某种程度的认可。
与开发早期极少数负面评价相同，同样存在一部分非常支持TBtools的开发工作。部分朋友，甚至将TBtools应用于学校课堂教学。这是我从未敢想的事，毕竟TBtools事实上只是我们课题组（前面硕士与现在博士）分析需求的部分体现。
我是一个潮汕人，潮汕人（或者可能是我家人给我的感觉）有一个传统，即知恩图报，更或者说，滴水之恩，涌泉相报。

开发这一功能的起因

在非常支持TBtools开发的用户朋友中，其常常与我提起，

能否开发出一个支持正则表达式检索的工具？用于查找结构域或者特定的模式，如微卫星....

总的来说，我觉得这个功能似乎我并不一定用得到，所以我或者会建议其去找另外的工具，比如MEME suite的MAST，或者自己编程实现。但是最后，似乎他还是没有按照我说的去操作。
昨晚我在等待数据下载的过程中，再次看到其提起。我想这已经是其超过三次的提起，或许确实可以花时间试试，毕竟辜负一个一直支持你的人的期待，似乎是一种不讲义气的操作。我大体想了下，这个事情似乎还不是我以前想的那么简单。

  1. MAST是不允许gaps的，其最大的作用是允许misM

  2. 自己变成时，那么多行的序列的回车需要处理；而如果一次将整条序列读入内存，那么遇到chromosome时，可能会有内存不足。

  3. 序列模式的Overlap的情况需要处理

  4. 不定长度的正则要求，可能确实不好找到工具来处理

总的来说，似乎确实找不到合适的工具可以很好的支持这些需求。所以，我写了一个。

功能的使用

首先是打开TBtools，找到最新的功能，Fasta Sequence Pattern Locator

需要设置的是三个参数：

对于输入的序列文件，事实上，需要保证的只是Fasta格式，而与其序列长度与类型（全基因组序列或者蛋白序列等）没有关系。
对于序列模式，那么需要用户对正则表达式有一定的了解，比如挖掘微卫星(AT){6,}表示ATATATATATAT....；或者从mRNA中预测ORF，ATG(?:\w{3})+T(?:AG|AA|GA)，随后对每个序列取最长的一个最长ORF即是。当然，也可以使蛋白的某些序列模式...
如果确实不了解正则表达式，其实还是比较简单，大概可能知道

ATCG 对应的就是 四个碱基，ATCG
[ATCG] 对应的是 一个碱基，A或T或C或G
[^AT] 对应的是 一个碱基，但不会是A或者T
(AT) 对应的是 两个碱基，把AT定位一个单元
(AT){6} 对应的是 2x6，一共 12 个碱基，也就是AT重复正好6次，如ATATATATATAT
(AT){,6} 对应的是(AT)重复不多于6次，如AT, ATAT, ATATAT, ATATATAT, ATATATATAT, ATATATATATAT
(AT){,6} 对应的是(AT)重复不低于6次, ATATATATATAT, ATATATATATATATATATATATAT.....
等等....

Overlap参数，大体对应的就是模式之间是否允许Overlap, 比如ATG\w{111}TGA，那么在这个模式捕获出来的序列中间是否有可能出现新的ATG?
Max Sequence Pattern(kb)，查找出来的片段大小，不能长于这个长度，如2即2kb。这个参数事实上，也对模式查找存在一定的影响。对于Overlap模式的，没有影响。对于不Overlap的，可能会丢失一小部分模式。

当所有参数设置好之后，设置输出文件，注意补齐文件名。随后点击Start即可。

输出结果

输出的结果，主要分为四列

分别是：

  1. 序列ID

  2. 起始坐标

  3. 终止坐标

  4. 匹配到模式的序列

其实是可以支持多个模式一次查找的，但我确实没有感觉到这个需求的大小，所以暂时也懒得支持。

赞 (0)

相关推荐