堕络的小屋

  • 首页
  • 小工具
    • 百度
    • 武器库
    • 堕络的图床
    • 子域名爆破
    • 音乐搜索器
    • 手绘相片制作
    • 微博图片找博主
    • 社会主义核心价值观编码转换
  • 值得一看
    • 黄色
    • 天天优惠
    • 剑灵小助手
  • 系统
    • 高清壁纸
    • 全网优惠券
    • 付费音乐解锁
    • 自动签到框架
    • 我们的足迹系统
    • 网易云音乐签到打卡
    • 全自动网页生成系统
    • 自动采集活动线报
堕络哥哥
一个专业打杂的程序猿
  1. 首页
  2. Linux
  3. 正文

C++使用htslib库读入和写出bam文件

2017年11月16日

有时候我们需要使用C++处理bam文件,比如取出read1或者read2等符合特定条件的序列,根据cigar值对序列指定位置的碱基进行统计或者对序列进行处理并输出等,这时我们可以使用htslib库。htslib可以用来处理SAM, BAM,CRAM 和VCF文件,是samtools、bcftools的核心库。

#include <stdio.h>
#include <stdlib.h>
#include <htslib/sam.h>

using namespace std; 

#define bam_is_read1(b) (((b)->core.flag&BAM_FREAD1) != 0)

uint8_t Base[16] = {0,65,67,0,71,0,0,0,84,0,0,0,0,0,0,78};

int main(int argc, char **argv)
{
    bam_hdr_t *header;
    bam1_t *aln = bam_init1();

    samFile *in = sam_open(argv[1], "r");
    htsFile *outR1 = hts_open(argv[2], "wb");
    header = sam_hdr_read(in);
    if (sam_hdr_write(outR1, header) < 0) {
        fprintf(stderr, "Error writing output.\n");
        exit(-1);
    }
    uint8_t *seq;
    int32_t lseq;
    uint32_t *cigar;
    char* qname;
    while (sam_read1(in, header, aln) >= 0) {
        if (bam_is_read1(aln)){
            sam_write1(outR1, header, aln);
        }
        else {
            seq = bam_get_seq(aln);
            lseq = aln->core.l_qseq;
            qname  = bam_get_qname(aln);
            printf("%s\n",qname);
            cigar = bam_get_cigar(aln);
            for(int i=0; i < aln->core.n_cigar;++i){
                int icigar = cigar[i];
                printf("%d%d\n",bam_cigar_op(icigar),bam_cigar_oplen(icigar));
            }
            for(int i=0; i < lseq;++i){
               printf("%c", Base[bam_seqi(seq, i)]);
            }
            printf("\n");

        }
    }
    sam_close(in);
    sam_close(outR1);
}

cigar值存储形式

32位int,通过bam_get_cigar获得地址,aln->core.n_cigar返回cigar operation的个数

  • 低 4位存储cigar operation;通过函数bam_cigar_op()获得operation
  • 高28位存储cigar值的长度;通过函数,bam_cigar_oplen获得

seq存储形式

8位int,4位存储一个碱基,1,2,4,8,15分别代表A、C、G、T、N,高四位存储坐标数较小的碱基,可通过bam_seqi(seq,i)遍历。

参考资料

htslib sam.h文件:https://github.com/samtools/htslib/blob/develop/htslib/sam.h
htslib sam文件格式说明:https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf

标签: Linux
最后更新:2017年11月16日

chenxing

'

点赞
< 上一篇
下一篇 >

文章评论

razz evil exclaim smile redface biggrin eek confused idea lol mad twisted rolleyes wink cool arrow neutral cry mrgreen drooling persevering
取消回复

COPYRIGHT © 2024 堕络的小屋. ALL RIGHTS RESERVED.

Theme Kratos Made By Seaton Jiang