如何找到基因组中的模糊碱基呢

之前的我发现基因组不仅存在ATCG和表示Gap的N,还有一些模糊序列,用来表示两种或两种以上的碱基。
对于这些序列,我有以下几个问题:
  1. 这些非ATCG的序列在基因组的哪些位置?
  2. 这些非ATCG的序列长度分别是多少?
  3. 基因组上存在多少个gap序列?
解决这个问题通常有以下两个思路:
  1. 通过检索,找到能够回答以上问题的工具
  2. 自己编写脚本,写一段代码进行分析。
而这里,我会用一个大家都想不到的工具来解答这些问题,这个工具就是我们经常用于二代序列回帖的BWA。
通常而言,我们使用bwa index建立索引,建立完索引之后,我们就直接用索引进行比对,而不会在乎索引文件。bwa index建立完索引之后,会得到5个文件,分别是amb, ann, bwt, pac,sa. 而amb就是我们所需的文件,amb是ambiguous的缩写,也就是模糊之意。
当我们对拟南芥的基因组建立索引之后,然后用less打开以amb结尾的文件,你会得到如下信息。
第一行表示,序列长度为119667750, 一共有7条,共563个模糊序列。
从第二行开始,分别记录着模糊序列的位置,它的长度,以及它的碱基信息。
通过Linux的基本文本操作,我们可以知道拟南芥一共有165个gap,占185,738 bp
cat Athaliana.amb | grep 'N$' | wc -l#165cat Athaliana.amb | grep 'N$' | datamash -t ' ' sum 2# 185738
这些碱基的具体位置,需要用ann文件进行解读,因为amb并非以chr:start-end形式进行记录,而是将所有序列从头到尾合并之后的结果。
第一行记录的是碱基总数,染色体数,随机数种子(seed), 后续的每两行为一个记录,第一行记录染色体编号信息,第二行记录染色体的偏移起始和结束位置,以及模糊碱基数。
生物信息学

WGCNA分析是如何找出基因模块的?

2020-6-25 20:25:02

生物信息学

GEO的数据注释文件没有基因名怎么办?

2020-6-25 20:31:07

声明 本网部分文章源于互联网,转载出于传递更多信息和学习之目的,并不意味着赞同其观点。
如转载稿涉及版权等问题,请立即联系管理员,我们会予以修改、删除相关文章,请留言反馈
When your legal rights are being violated, please send an email to: sci666net@qq.com
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索