科研星球

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

之前的我发现基因组不仅存在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结尾的文件,你会得到如下信息。
1.png
第一行表示,序列长度为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形式进行记录,而是将所有序列从头到尾合并之后的结果。
2.png
第一行记录的是碱基总数,染色体数,随机数种子(seed), 后续的每两行为一个记录,第一行记录染色体编号信息,第二行记录染色体的偏移起始和结束位置,以及模糊碱基数。


没有账号?