专栏名称: 生信媛
生信媛,从1人分享,到8人同行。坚持分享生信入门方法与课程,持续记录生信相关的分析pipeline, python和R在生物信息学中的利用。内容涵盖服务器使用、基因组转录组分析以及群体遗传。
目录
相关文章推荐
51好读  ›  专栏  ›  生信媛

BWA源码阅读笔记(二)索引文件amb/ann/pac文件是什么?

生信媛  · 公众号  · 生物  · 2020-03-24 15:54

正文

请到「今天看啥」查看全文


高通量数据比对讲究的就是一个快和准,因此大部分软件都是用C语言实现。BWA是目前基因组序列比对最常用的工具,由于自我感觉已经入门C语言,为了提高自己的水平,因此开始从源码角度学习李恒大神开发的BWA工具。

使用 bwa index 建立索引后,会得到以下后缀的文件, .amb, .ann, .pac, .bwt, .sa. 我们通常也不在乎这些文件是什么,除非你像我一样,想搞懂bwa建立索引的具体过程。这次我们先来介绍.amb, .ann和.pac这三个文件是什么。

先说结论,amb是ambiguous的缩写,也就是模棱两可的意思,也就是除了ATCG/atcg以外的字符. amb和ann用来记录基因组中除了ATCG以外碱基的信息。而pac文件则是碱基信息高度压缩。

接下来的介绍,主要涉及到两个文件,bntseq.h和bntseq.c, 这里面的代码就是用于将fasta转成amb, ann和pac这三个文件。

bntseq . h 的开头,就定义了对应输出文件的三种数据结构, bntann1_t , bntamb1_t bntseq_t

  1. typedef struct {

  2. int64_t offset; //染色体偏移量

  3. int32_t len; //染色体长度

  4. int32_t n_ambs; //多少个模糊碱基

  5. uint32_t gi;

  6. int32_t is_alt;

  7. char *name, *anno; //染色体名名和注释

  8. } bntann1_t ;


  9. typedef struct {

  10. int64_t offset; //偏移量

  11. int32_t len; //长度

  12. char amb; //模式碱基的字符

  13. } bntamb1_t;


  14. typedef struct {

  15. int64_t l_pac;

  16. int32_t n_seqs; //基因组序列数

  17. uint32_t seed; //随机数种子

  18. bntann1_t *anns; // n_seqs elements

  19. int32_t n_holes; //染色体有多少个空缺

  20. bntamb1_t *ambs; // n_holes elements

  21. FILE *fp_pac;

  22. } bntseq_t;

前两种数据结构能够保存成文本,我们以拟南芥基因组为例。ann文件前几行如下,它记录每个染色体信息,以及里面amb数目

  1. 239335500 7 11 # bns->l_pac, bns->n_seqs, bns->seed)

  2. 0 Chr1 (null) # p->gi, p->name, p->anno

  3. 0 30427671 363 # p->offset, p->len, p->n_ambs

  4. 0 Chr2 (null)

  5. 30427671 19698289 58

  6. ...

amb文件前几行如下,它 记录确切的amb位置

  1. 239335500 7 563 # bns->l_pac, bns->n_seqs, bns->n_holes

  2. 12192623 1 Y # p->offset, p->len, p->amb

  3. 13201252 100 N # p->offset, p->len, p->amb

  4. ...

第三种数据结构, bntseq_t ,则是用于打包将前两者的信息进行了整合,同时记录 pac 文件的文件指针地址。

对于pac文件,他通过哈希映射的方法(哈希函数如下),把一个比较大的位置映射到比较小的位置中,保证了最终序列文件大小低于原始文件。

  1. #define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<l)&3)<<1))

  2. #define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)

总结一下,amb,ann,pac这三个文件主要是用于压缩原始的参考基因组序列,将碱基信息分为两部分,ATCG和其他模棱两可的序列。amb记录模棱两可序列的具体为止,ann记录参考基因组中染色体信息,每条染色体的偏移信息和长度信息。我们原本可以去掉染色体名,将所有的序列保存到一个文件中,但是这会导致文件过大,不利于后续的读取操作。因此,李恒采用了一种计算方式,使得只用四分之一的体质保存原来的信息。








请到「今天看啥」查看全文